﻿  !*********************************************************************************************************************
  !  AUTHOR : Tom Langbehn | tom.langbehnb@uib.no
  !  PROJECT: Limits to northern range shifts in marine fish
  !  MODEL  : Dynamic optimization model of vertical migration and fitness trade-offs in mesopelagic fish
  !  TYPE   : Main program
  !  PURPOSE: Optimize vertical positioning
  !  VERSION: June 2020   | v2.0
  !*********************************************************************************************************************
  PROGRAM main_program_vassild

  use GLOBAL_PARAMETERS ! The global prarameters module contains all global parameters and  declarations
  use IFPORT            ! The ifport module is needed to return elapsed time
  use HELPERS           ! The helpers module contains several utility functions, ::proust()
  use VISRANGE          ! The visual range calculation module itself
  use REAL_PRECISION    ! The precision module is used in visual range calculation module
  use SET_ENV           ! The set env module sets the temperature and light environment for a given latitude
  use BINARY_IO         ! The binary_IO module allows to write and read ninary .bin files
  use BIOENERGETICS     ! The bioenergetics module provdes a suit of funtions related to fish bioenergetics
  use FORAGING          ! The foraging interactions module, contains functions for clearance rate
  use SET_PREYFIELD     ! The set preyfield module sets the preyfield for a given latitude, based on generic assumptions

  IMPLICIT NONE

  !declarations
  character(len=1024):: filename
  character(2) :: percent = ']'// char(13)
  integer n
  real(SRP), dimension(1:no_depth,1:E_max, 1:G_max, 1:L_max, 1:int(t_max), 1:no_days):: strategy
  real(SRP), dimension(1:no_depth, 1:E_max, 1:G_max, 1:L_max) :: V , Vnext
  real(SRP), dimension(1:no_depth):: visr_meso_m, clear_meso_m3, encr_meso_ind, ingr_meso_ind, ingr_meso_J
  real(SRP), dimension(1:no_depth):: visr_pred_m, clear_pred_m3, encr_pred_ind, predr_indts
  real(SRP), dimension(1:no_depth, 1:no_ts, 1:8):: output_bwd
  real(SRP), dimension(1:30, 0:no_ts, 1:L_max):: output_fwd
  real(SRP), dimension(1:E_max):: survival, reserves, surplus, dE
  integer  , dimension(1:E_max):: intE
  real(SRP), dimension(1:no_depth, 1:E_max):: fitnessatdepth
  real(SRP), dimension(1:no_days,1:no_lat,1:L_max):: cohort_fitness

  real(SRP), dimension(1:no_ts, 1:2):: lcz

  real(SRP) gut, gut_new, visr_tmp, GER, SMR, RMR, F, U, S, Mt, digestion_max, digestion, netenergy
  real(SRP) newG, dG, newE, dE_fwd, reserves_fwd, surplus_fwd, survival_fwd
  real(SRP) depth, excess, temp_u, temp_l
  real(SRP) surplus_cum, survival_cum
  integer G, intG, currentz, newz, E, z_opt(1), v_opt(1),intE_fwd

  real(SRP) daily_survival(1:no_days), daily_surplus (1:no_days), expct_survival(1:no_days), expct_surplus(1:no_days)
  real(SRP) survival_useday, surplus_useday, survival_tmp, surplus_tmp
  integer cohortday,useday, lat_index
  real(SRP) t1,t2

  character(len=8):: systime

  !> Initalize and warning messages
  !! elapsed time start
  run_start = dclock()

  !> show some number son the screen
  call time(systime)
  print '(A4,A21,A9)','>>', 'Model run started at', systime
  print*, ' '

  print'(A4,A17)', '>', 'Model resolution'
  print'(A4,A)', '','----------------------------------------------------------------------------'
  print'(A24,f5.2,A,f6.2,A,f5.2,A,A,I4,A)', ' | Warming : from    ', warm_min         ,' to  '    , warm_max         ,' by  '   , warm_inc           ,' degrees'    ,' \ ' , no_warm             , ' scenarios    |'
  print'(A24,f5.2,A,f6.2,A,f5.2,A,A,I4,A)', ' | Latitude: from    ', lat_min          ,' to  '    , lat_max          ,' by  '   , lat_inc ,' degrees',' \ ' , no_lat              , ' latitudes    |'
  print'(A24,f5.2,A,f7.2,A,f7.2,A,A,I4,A)', ' | Depth   : from    ', dble(z_min)      ,' to '     , dble(z_max)      ,' by'     , dble(z_inc)        ,' meters '    ,' \ ' , no_depth            , ' depth bins   |'
  print'(A24,I2,A,I3,A,I2,A,A,I4,A)'      , ' | Time    : from day', day_start        ,'    to  ' , day_end          ,'    by  ', int(ts_dur*60)     ,'    minutes' ,' \ ' , no_ts               , ' timesteps    |'
  print'(A4,A)', '','----------------------------------------------------------------------------'
  print*, ' '
  print*, ' '

  !> if not existent create new sub directory for different types of data output
  call create_directory(sub_dir_bin)
  call create_directory(sub_dir_txt )

  !> initate files to collect cumulative survival and surplus
  open (510, file= trim(sub_dir_txt) // '\output_fitness_by_lat__' // trim(parametrization) // '.txt')  ! output to write test files
  write(510,'(7A)') 'Temp_inc ', 'Latitude ', 'Length_m ', 'Survival_cum ', 'Surplus_JJ ', 'expct_surplus_JJ ', 'parametrization'

  !open (511, file= trim(sub_dir_txt) // '\output_fitness_by_lat_' // trim(parametrization) // '_in_J.txt')  ! output to write test files
  !write(511,'(6A)') 'Temp_inc ', 'Latitude ', 'Length_m ', 'Survival_cum ', 'Surplus_J ', 'expct_surplus_J'

  !open (512, file= trim(sub_dir_txt) // '\antipredationwindow_by_lat_' // trim(parametrization) // '.txt')
  !write(512,'(7A)') 'Lat ', 'ts ',  'depth_m ', 'predR ', 'cleaR ', 'enctR ', 'eb'

  !> MAIN BODY                                                                                [UNITS]
  ## Set the agents
  write(*, '(A4,A69)',advance='no')'>>', 'Setting the agents..................................................'
  !!> set CALANUS properties
  ! calanus width                                                                             [m]
  calanus%width_m(1) = calanus%shape(1) * calanus%length_m(1)
  ! calanus image area                                                                        [m-2]
  calanus%area_m2(1) = calanus%width_m(1) * calanus%length_m(1) * 0.75_SRP
  ! calanus energy content                                                                    [J ind-1]
  calanus%energy_J_ind(1) = calanus%weight_gWW(1) * calanus%energy_J_gWW(1)

  !> set MESOPELAGIC FISH properties
  !> loop over length classes
  do L = L_min, L_max
    ! calculate weight from length [m] for a given species ("BenGla" = Benthosema glaciale)   [g]
    mesopf%weight_g(L) = f_weight(mesopf%length_m(L), "BenGla")
    ! maximum gut capacity for a given length class                                           [J calanus-1]
    mesopf%gut_max(L) =  (mesopf%weight_g(L) * mesopf%gut_capacity(L)) * calanus%energy_J_gWW(1)
    ! maximum reserve capacity for a given length class                                       [J]
    mesopf%reserves_max(L) =  (mesopf%weight_g(L) * mesopf%res_capacity(L)) * mesopf%energy_J_gWW(L)
    ! swimming speed                                                                          [m s-1]
    mesopf%swimvel_ms(L) = mesopf%swimvel_BLs(L)*mesopf%length_m(L)
    ! eye sensitivity                                                                         [m]
    mesopf%eyesens = ((mesopf%length_m(L) * mesopf%maxvisr_BL(L))**2.)/(calanus%contrast(1)*calanus%area_m2(1))
    ! mesopelagic fish width                                                                  [m]
    mesopf%width_m   = mesopf%shape(L)*mesopf%length_m(L)
    ! mesopelagic fish image area                                                             [m-2]
    mesopf%area_m2 = mesopf%width_m(L)*mesopf%length_m(L)*0.75
    ! energy content per indidividual                                                         [J ind-1]
    mesopf%energy_J_ind(L) = mesopf%weight_g(L) * mesopf%energy_J_gWW(L)
  enddo

  write(*, '(A7)', advance='yes')'[100%]'
  print*,''

  !!> set PISCIVOROUS FISH properties
  ! swimming speed [m s-1]
  pred%swimvel_ms(1) = pred%swimvel_BLs(1)*pred%length_m(1)
  ! eye sensitivity (calculated base don the largest mesoeplagic fish size class)
  pred%eyesens = ((pred%length_m(1) * pred%maxvisr_BL(1))**2.)/(mesopf%contrast(2)*mesopf%area_m2(2))

  ## Set the environment
  !> set the environment for all days and latitudes

  ### INDEXING OF LATITUDES
  call set_environment()
  !! the environment is now available as multi-dimensional arrays of the form "Lat x Timestep (x Depth)"
  !! set_environment() returns: "sun_height_2d", "light_surf_2d", "light_amb_3d", Temp: "temp_3d"

  ## set the prey field
  !> set the prey field for all days and latitudes
  call set_zoop_preyfield()
  !>set_zoop_preyfield() returns zoop_3d(50:85, 1:no_days, 1:1001), prey is calculated for every meter between 0-1000m depth

  !> fix the environment for debugging
  IFTRUE_fixenv: if (fix_env == .TRUE.) then
    temp_3d(:,:,:)      =  7._SRP
    sun_height_2d(:,:)  = 10._SRP
    light_surf_2d(:,:)  = 10._SRP
    light_amb_3d(:,:,:) = 10._SRP
    zoop_3d(:,:,:)      = 10._SRP
  endif IFTRUE_fixenv


  ENV_ONLY: if (run_env_only == .FALSE.) then
    print '(A4,A24)','>>', 'Find optimal behaviours'

    ## Warming
    !> loop over warming scenarios
    WARMING_LOOP: do itemp = 1, no_warm
      warming = all_warm_degC(itemp)


      ## Latitude
      !> loop over latitudes
      LAT_LOOP: do ilat = 1, no_lat
        !> latitude as real number
        latitude_real   =  all_latitudes(ilat)
        !> relativ eposition of the chosen latitude between 50-85°N
        lat_index = int(latitude_real-50._SRP)/1 + 1

        print*,''
        write(*, '(A4,A9,F4.1,A2,I3.3,A1,I3.3,A1,A4,A9,F5.1,A2,I3.3,A1,I3.3,A1)',advance='yes')'>' ,'Warming:' , warming , '[', itemp, '/', no_warm,']',' -- ','Latitude:', latitude_real, '[', ilat, '/' , no_lat ,']'
        print*,''

        !> calculate depth of the "light comfort zone" for every timestep
        ! upper boundary
        lcz(:,1) = f_isodepth(1.6E-5_SRP, light_amb_3d(ilat,:,1), DiffAtt_H2O_upr, DiffAtt_H2O_lwr)
        ! lower boundary
        lcz(:,2) = f_isodepth(2.5E-7_SRP, light_amb_3d(ilat,:,1), DiffAtt_H2O_upr, DiffAtt_H2O_lwr)

        !> set the terminal fitness criterion
        ### maybe 0?
        V(:,:,:,:) = 1.

        ## -- BACKWARDS ------
        IFTRUE_bwd: if (run_bwd == .TRUE.) then

          ## Year
          !> loop over years to avoid terminal effects
          CONV_LOOP_bwd: do conv_bwd = 1, conv_max_bwd
            !> tic
            t1 = dclock()
            ! the dynamic optimization works backwards in time
            ! therefore, we reset the timestep counter to the terminal step for every new year
            ts = no_ts + 1

            ## Day
            !> loop over days
            DAY_LOOP_bwd: do d = no_days, 1, -1
              day = all_days(d)

              !> calculate warming in near surface waters
              temp_inc(:) = (abs(temp_z(:)-maxval(temp_z(:)))**9 / maxval(abs(temp_z(:)-maxval(temp_z(:)))**9))* warming
              !> add warming to the temperature environemnt
              temp_3d(ilat,:, d) =  temp_3d(ilat,:, d) + temp_inc(:)

              ## Timestep
              !> loop over timesteps (t_max * no_days)
              TIMESTEP_LOOP_bwd: do t = t_max, 1, -1
                ! convert t to hours [0-23]
                hour = (1 + (t - 1)*ts_dur) - 1
                ! move one step back in time
                ts = ts - 1
                ! set Vnext
                Vnext(:,:,:,:) = V(:,:,:,:)

                ## Length
                !> loop over length class
                LENGTH_LOOP_bwd: do L = L_min, L_max

                  ## Depth
                  !> loop over depth
                  DEPTH_LOOP_bwd: do z = 1, no_depth
                    depth_m = dble(z_min + (z - 1)*z_inc)

                    !> MESOPELAGIC FISH FORAGING
                    !> calc mesopelagic fish visual range at given depth                          [m]
                    call srgetr(visr_tmp, BeamAtt_H2O, calanus%contrast(1), calanus%area_m2(1), mesopf%eyesens(L), mesopf%Ke(L), light_amb_3d(ilat, ts, z))
                    !> make sure prey is always detected when close                               [m]
                    visr_meso_m(z) = max(visr_tmp, mesopf%minvisr_m(L))
                    !> calc mesopelagic fish clearance rate                                       [m^3 s-1]
                    clear_meso_m3(z) = f_clear(mesopf%scndarea(L), visr_meso_m(z), mesopf%swimvel_ms(L))
                    !> calc encounter rate of mesopelagic fish and their prey                     [ind s-1]
                    encr_meso_ind(z) = f_encr(clear_meso_m3(z), zoop_3d(lat_index, day, int(depth_m)+1), mesopf%handlingtime(L))
                    !> convert envounters from [ind s-1] to                                       [ind timestep-1]
                    encr_meso_ind(z) = encr_meso_ind(z) * 3600._SRP * ts_dur
                    !> calculate encounters where the prey is successfully ingested               [ind timestep-1]
                    ingr_meso_ind(z) = encr_meso_ind(z) * mesopf%pc(L)

                    !> total ingested energy per timestep                                         [J timestep-1]
                    ingr_meso_J(z) = ingr_meso_ind(z) * calanus%energy_J_ind(1)

                    !> EPIPELAGIC PISCIVOROUS FISH FORAGING (mesopleagic fish mortality) at a given depth
                    !> calc predatory fish visual range                                           [m]
                    call srgetr(visr_tmp, BeamAtt_H2O, mesopf%contrast(L), mesopf%area_m2(L), pred%eyesens(1), pred%Ke(1), light_amb_3d(ilat, ts, z))
                    !> make sure prey is always detected when close                               [m]
                    visr_pred_m(z) = max(visr_tmp, pred%minvisr_m(1))
                    !> calc predatory fish clearance rate                                         [m^3 s-1]
                    clear_pred_m3(z) = f_clear(pred%scndarea(1), visr_pred_m(z), pred%swimvel_ms(1))
                    !> calc encounter rate of predatory fish and mesopelagic fish                 [ind s-1]
                    encr_pred_ind(z) = clear_pred_m3(z) * pred%no_ind_m3(1)
                    !> convert envounters from [ind s-1] to                                       [ind timestep-1]
                    encr_pred_ind(z) = encr_pred_ind(z) * 3600._SRP * ts_dur
                    !> calculate predation rate of mesopelagic fish                               [ind timestep-1]
                    predr_indts(z) = encr_pred_ind(z) * pred%pc(1)

                    !> write anitpredation windwo to txt.
                    !if(latitude_real == 55. .AND. day == 172 .AND. L == 2 .OR. latitude_real == 80. .AND. day == 202 .AND. L == 2) then
                    !  write(512,'(F5.2,I5, F7.1,3F15.5, ES10.2)') latitude_real , ts , depth_m, predr_indts(z), clear_meso_m3(z), encr_meso_ind(z), light_amb_3d(ilat, ts, z)
                    !endif

                    !> save output from backwards iteration
                    if(save_bwd == .TRUE. .AND. L == 2) then
                      !!(ts,z,i )                 (ts,z,L)
                      output_bwd(z, ts, 1) =  visr_pred_m (z)
                      output_bwd(z, ts, 2) =  predr_indts (z)
                      output_bwd(z, ts, 3) =  visr_meso_m  (z)
                      output_bwd(z, ts, 4) =  ingr_meso_ind(z)
                      output_bwd(z, ts, 5) =  zoop_3d(lat_index, day, int(depth_m)+1)
                      output_bwd(z, ts, 6) =  light_amb_3d(ilat, ts, z)
                      output_bwd(z, ts, 7) =  temp_3d(ilat, z, d)
                      output_bwd(z, ts, 8) =  light_surf_2d(ilat, ts)
                    end if

                  enddo DEPTH_LOOP_bwd

                  ## Gut
                  GUT_LOOP_bwd: do G = 1, G_max
                    ! gut content                                                               [J]
                    gut = ((real(G) - 1.)/(real(G_max) - 1.)) * mesopf%gut_max(L)

                    ## Current depth
                    CURRENTDEPTH_LOOP_bwd: do currentz = 1, no_depth

                      !!STRATEGIES
                      ! intentionally unrealistic number to fill
                      fitnessatdepth(:,:) = -1.d20

                      ## New depth
                      !> new depth constrained to be somewhat near current depth but between z_min+1 and z_max+1
                      NEWDEPTH_LOOP_bwd: do newz = max(1, currentz - maxvertdist_zh), min(currentz + maxvertdist_zh, no_depth)

                        !! BIOENERGETICS
                        !> gastric evacuation rate                                              [fraction of gut h-1]
                        GER = f_GER(temp_3d(ilat,newz,d), "BenGla")
                        !> convert from [fraction of gut h-1] to                                [fraction of gut timestep-1]
                        GER = GER*ts_dur
                        !> maximum digestion during any timestep, how much can be used          [J]
                        digestion_max = mesopf%gut_max(L) * (1. - exp(-GER))
                        !> digestion during any timestep, how much can be used                  [J]
                        digestion = gut * (1. - exp(-GER))
                        !> fish might encounter more food than the gut can handle               [J]
                        gut_new = max(0._SRP, min(gut + ingr_meso_J(newz) - digestion, mesopf%gut_max(L)))
                        !> as in Hewett & Johnson 1988
                        !> calculate F : EGESTION                                               [J ind-1 timestep-1]
                        F = mesopf%Fa(L) * digestion
                        !> calculate U : EXCRETION                                              [J ind-1 timestep-1]
                        U = mesopf%Ua(L) * (digestion - F)
                        !> calculate standart metabolic rate                                    [J ind-1 h-1]
                        SMR = f_SMR(mesopf%weight_g(L), temp_3d(ilat,newz,d), "BenGla")
                        !> convert from [J ind-1 h-1] to                                        [J ind-1 timestep-1]
                        SMR = SMR*ts_dur
                        !> scale SMR with an activity coeff. to get routine metabolic rate      [J ind-1 timestep-1]
                        RMR = SMR*activity

                        !> calculate energy accounted for by specific dynamic action            [J ind-1 ts-1]
                        S = mesopf%Sa(L)*(digestion - F)
                        !> calculate netenergy                                                  [J ind-1 ts-1]
                        netenergy = digestion - (RMR + S) - (F + U)
                        !> calculate mortality and survival                                     [surv. probability timestep-1]
                        Mt = predr_indts(newz)
                        survival(:) = exp(-Mt)

                        !> record fitness at every depth, basis for making optimal decision
                        !> interpolate new gut state for inbetween gut states
                        newG = (gut_new/mesopf%gut_max(L)) * (dble(G_max) - 1.) + 1._SRP
                        intG = max(1, min(idint(newG), G_max - 1))
                        dG   = max(0._SRP, min(newG - dble(intG), 1._SRP))

                        ## RESERVES
                        ENERGY_RESERVES_LOOP: do E = 1, E_max
                          !> level of reserves                                                  [J]
                          reserves(E) = (real(E) - 1._SRP)/(E_max - 1._SRP) * mesopf%reserves_max(L)
                          !> level or reserves after feeding                                    [J]
                          reserves(E) = reserves(E) + netenergy
                          !> aquired energy not needed to fill the reserves is surplus energy   [J]
                          surplus (E) = max(0._SRP, reserves(E) - mesopf%reserves_max(L) )
                          !> if energy reserves are depleted fish will starve                   [J]
                          if(reserves(E) < 0._SRP) survival(E) = 0._SRP
                          !> interpolate new reserve state for inbetween reserve states         [J]
                          newE    = (reserves(E)/mesopf%reserves_max(L)) * (dble(E_max) - 1._SRP) + 1._SRP
                          intE(E) = max(1, min(idint(newE), E_max - 1))
                          dE  (E) = max(0._SRP, min(newE - dble(intE(E)), 1._SRP))

                        enddo ENERGY_RESERVES_LOOP

                        !> dynamic programming equation, fitness at depth for all E's
                        fitnessatdepth(newz,:) = survival(:)*(surplus(:) + &
                          (1. - dE(:)) * (1. - dG) * Vnext(newz, intE(:)    , intG, L) +  (1. - dE(:)) * dG * Vnext(newz, intE(:)    , intG + 1, L) + &
                          dE(:)  * (1. - dG) * Vnext(newz, intE(:) + 1, intG, L) +        dE(:)  * dG * Vnext(newz, intE(:) + 1, intG + 1, L))

                      enddo NEWDEPTH_LOOP_bwd

                      do E = 1, E_max
                        !> find optimal depth where fitness is maximized
                        z_opt(:) = maxloc(fitnessatdepth(:, E))
                        !z_opt(:) = 30
                        !> write optimal depth to strategy matrix
                        strategy(currentz, E, G, L, t, d) = z_min + (z_opt(1) - 1)*z_inc
                        !> write fitness at optimal depth to fintess matrix
                        V(currentz, E, G, L) = fitnessatdepth(z_opt(1), E)
                      enddo


                    enddo CURRENTDEPTH_LOOP_bwd
                  enddo GUT_LOOP_bwd
                enddo LENGTH_LOOP_bwd
              enddo TIMESTEP_LOOP_bwd
            enddo DAY_LOOP_bwd


            !> toc
            write(*,'(A6, A18, A2,I3.3,A1,I3.3,A1)', advance='no')' | ','Backward iteration', '[', conv_bwd, '/' , conv_max_bwd ,']'
            t2 = dclock()
            write(*,'(A4,A,A34)', advance='yes')' -- ',f_hhmmss(idint(t2-t1)), '|'

          enddo CONV_LOOP_bwd
        endif IFTRUE_bwd

        ## -- FORWARD ------
        !> tic
        t1 = dclock()
        IFTRUE_fwd: if (run_fwd == .TRUE.) then
          !> INITIATE FORWARD
          !> array collects output
          output_fwd(:,:,:)  = 0._SRP

          ## LENGTH
          LENGTH_LOOP_fwd: do L = L_min, L_max
            !> initial gut fulness    , arbitrary number                     [J]
            gut = 1.0 * mesopf%gut_max(L)
            !> initial reserve fulness, arbitrary number                     [J]
            reserves_fwd = 1.0 * mesopf%reserves_max(L)
            !> initial depth bin      , arbitrary number                     [z]
            z = no_depth/2

            ## Year
            !> run multiple years to allow convergence
            CONV_LOOP_fwd: do conv_fwd = 1, conv_max_fwd
              !> we simulate forward in time, starting at timestep 0
              ts = 0
              !> cumulative surplus is 0 at the start and builds up over time
              surplus_cum = 0._SRP
              !> cumulative survival is 1 at the start and declines over time
              survival_cum = 1._SRP

              ## DAY
              DAY_LOOP_fwd: do d = 1, no_days
                day = all_days(d)
                !> at the beginning of each day, reset arrays to collect sum of hourly value for given day
                survival_tmp = 1.0
                surplus_tmp  = 0.0

                ## TIMESTEP
                TIMESTEP_LOOP_fwd: do t = 1, t_max
                  ! convert ts to hours [0-23]
                  hour = 1._SRP + (t - 1._SRP)*ts_dur
                  !> move one step forward in time
                  ts = ts + 1

                  !> FIRST, look up OPTIMAL DEPTH
                  ! calc gut state
                  newG = (gut/mesopf%gut_max(L))*(dble(G_max) - 1.) + 1.
                  intG = max(1,  min(idint(newG), G_max - 1))
                  dG   = max(0._SRP, min(newG - dble(intG), 1._SRP))

                  ! calc reserve state
                  newE = (reserves_fwd/mesopf%reserves_max(L))*(dble(E_max) - 1._SRP) + 1._SRP
                  intE_fwd = max(1, min(idint(newE), E_max - 1))
                  dE_fwd   = max(0._SRP, min(newE - dble(intE_fwd), 1._SRP))

                  !> look up optimal depth strategy in matrix
                  depth = (1. - dE_fwd)*(1. - dG)*strategy(z, intE_fwd    , intG , L, t, d) + (1. - dE_fwd)*dG*strategy(z, intE_fwd    , intG + 1, L, t, d) + &
                    dE_fwd *(1. - dG)*strategy(z, intE_fwd + 1, intG , L, t, d) +       dE_fwd *dG*strategy(z, intE_fwd + 1, intG + 1, L, t, d)

                  !> we place fish in nearest z bin.
                  z = idnint((depth-dble(z_min))/dble(z_inc)) + 1

                  !> convert to depth
                  depth = dble(z_min + (z - 1)*z_inc)

                  !> SECOND, calculate MESOPELAGIC FISH FORAGING
                  !> calc mesopelagic fish visual range at given depth                          [m]
                  call srgetr(visr_tmp, BeamAtt_H2O, calanus%contrast(1), calanus%area_m2(1), mesopf%eyesens(L), mesopf%Ke(L), light_amb_3d(ilat, ts, z))

                  !> make sure prey is always detected when close                               [m]
                  visr_meso_m(z) = max(visr_tmp, mesopf%minvisr_m(L))

                  !> calc mesopelagic fish clearance rate                                       [m^3 s-1]
                  clear_meso_m3(z) = f_clear(mesopf%scndarea(L), visr_meso_m(z), mesopf%swimvel_ms(L))

                  !> calc encounter rate of mesopelagic fish and their prey                     [ind s-1]
                  encr_meso_ind(z) = f_encr(clear_meso_m3(z), zoop_3d(lat_index, day, int(depth)+1), mesopf%handlingtime(L))

                  !> convert envounters from [ind s-1] to                                       [ind timestep-1]
                  encr_meso_ind(z) = encr_meso_ind(z) * 3600._SRP * ts_dur

                  !> calculate encounters where the prey is successfully ingested               [ind timestep-1]
                  ingr_meso_ind(z) = encr_meso_ind(z) * mesopf%pc(L)

                  !> total ingested energy per timestep                                         [J timestep-1]
                  ingr_meso_J(z) = ingr_meso_ind(z) * calanus%energy_J_ind(1)

                  !> THIRD, calculate  EPIPELAGIC PISCIVOROUS FISH FORAGING (mesopleagic fish mortality) at a given depth
                  !> calc predatory fish visual range                                           [m]
                  call srgetr(visr_tmp, BeamAtt_H2O, mesopf%contrast(L), mesopf%area_m2(L), pred%eyesens(1), pred%Ke(1), light_amb_3d(ilat, ts, z))

                  !> make sure prey is always detected when close                               [m]
                  visr_pred_m(z) = max(visr_tmp, pred%minvisr_m(1))

                  !> calc predatory fish clearance rate                                         [m^3 s-1]
                  clear_pred_m3(z) = f_clear(pred%scndarea(1), visr_pred_m(z), pred%swimvel_ms(1))

                  !> calc encounter rate of predatory fish and mesopelagic fish                 [ind s-1]
                  encr_pred_ind(z) = clear_pred_m3(z) * pred%no_ind_m3(1)

                  !> convert envounters from [ind s-1] to                                       [ind timestep-1]
                  encr_pred_ind(z) = encr_pred_ind(z) * 3600._SRP * ts_dur

                  !> calculate predation rate of mesopelagic fish                               [ind timestep-1]
                  predr_indts(z) = encr_pred_ind(z) * pred%pc(1)

                  !! BIOENERGETICS
                  !> gastric evacuation rate                                                    [fraction of gut h-1]
                  GER = f_GER(temp_3d(ilat,z,d), "BenGla")

                  !> convert from [fraction of gut h-1] to                                      [fraction of gut timestep-1]
                  GER = GER*ts_dur

                  !> maximum digestion during any timestep, how much can be used                [J]
                  digestion_max = mesopf%gut_max(L) * (1._SRP - exp(-GER))

                  !> digestion during any timestep, how much can be used                        [J]
                  digestion = gut * (1._SRP - exp(-GER))

                  !> fish might encounter more food than the gut can handle                     [J]
                  gut_new = max(0._SRP, min(gut - digestion + ingr_meso_J(z), mesopf%gut_max(L)))

                  !> how much more does it encounter than it can process                        [fraction]
                  excess = (gut - digestion + ingr_meso_J(z) - mesopf%gut_max(L))/mesopf%gut_max(L)

                  !> update gut
                  gut = gut_new

                  !> as in Hewett & Johnson 1988
                  !> calculate F : EGESTION                                                     [J ind-1 timestep-1]
                  F = mesopf%Fa(L) * digestion

                  !> calculate U : EXCRETION                                                    [J ind-1 timestep-1]
                  U = mesopf%Ua(L) * (digestion - F)

                  !> calculate standart metabolic rate                                          [J ind-1 h-1]
                  SMR = f_SMR(mesopf%weight_g(L), temp_3d(ilat,z,d), "BenGla")

                  !> convert from [J ind-1 h-1] to                                              [J ind-1 timestep-1]
                  SMR = SMR*ts_dur

                  !> scale SMR with an activity coeff. to get routine metabolic rate            [J ind-1 timestep-1]
                  RMR = SMR*activity

                  !> calculate energy accounted for by specific dynamic action                  [J ind-1 ts-1]
                  S = mesopf%Sa(L)*(digestion - F)

                  !> calculate netenergy                                                        [J ind-1 ts-1]
                  netenergy = digestion - (RMR + S) - (F + U)
                  reserves_fwd = reserves_fwd + netenergy

                  !> to avoid numerical errors where reserves turn slightly negative in forward part, despite being positive in optimization
                  if ((reserves_fwd < 0._SRP) .and. (reserves_fwd > -0.00001_SRP)) reserves_fwd = 1.E-10_SRP

                  !> all energy not needed to fill reserves is allocated to surplus
                  surplus_fwd  = max(0._SRP, reserves_fwd - mesopf%reserves_max(L))                                                                             !(J)
                  !> cumulative surplus at every timestep
                  surplus_cum = surplus_cum + surplus_fwd
                  !> summed for every day
                  surplus_tmp = surplus_tmp + surplus_fwd


                  !> calculate mortality and survival                                           [surv. probability timestep-1]
                  Mt = predr_indts(z)
                  survival_fwd = exp(-Mt)
                  !> cumulative survival at every timestep
                  survival_cum = survival_cum * survival_fwd
                  !> summed for every day
                  survival_tmp = survival_tmp * survival_fwd

                  !> reserve size is constraint by an upper limit reserves_max
                  reserves_fwd = min(reserves_fwd, mesopf%reserves_max(L))
                  !> if energy reserves are depleted fish will starve
                  if(reserves_fwd < 0._SRP) survival_fwd = 0._SRP


                  ! migration behaviour across lats and time
                  if(L == L_max) then
                    mig_behav_2d_real(ilat, ts) = depth
                  end if

                  !> write output from forward to array
                  output_fwd( 1,ts,L) = depth                                                                                                      !(m)

                  output_fwd( 2,ts,L) = light_surf_2d(ilat, ts)                    ! surface irradiance                                       !(W m-2 s-1)
                  output_fwd( 3,ts,L) = light_amb_3d(ilat, ts, z)                  ! light intensity at depth                                      !(W m-2 s-1)
                  output_fwd( 4,ts,L) = light_amb_3d(ilat, ts, z)/light_surf_2d(ilat, ts) ! light intensity at depth                               !(% of surface irradiance)
                  output_fwd( 5,ts,L) = sun_height_2d(ilat, ts)

                  output_fwd( 6,ts,L) = visr_meso_m(z)*100.                        ! visual range (from m to cm)                                   !(cm)
                  output_fwd( 7,ts,L) = clear_meso_m3(z)                           ! clearance rate                                                !(m-3 ts-1)

                  output_fwd( 8,ts,L) = ingr_meso_ind(z)                           ! ingestion rate (successful capture events)                    !(prey items ts-1)
                  output_fwd( 9,ts,L) = (ingr_meso_ind(z)*calanus%weight_gWW(1))/mesopf%weight_g(L)! ingestion rate                                !(g g-1 ts-1)
                  output_fwd(10,ts,L) = ingr_meso_J  (z)                           ! ingestion rate                                                !(J ts-1)

                  output_fwd(11,ts,L) = digestion                                  ! digestion rate                                                !(J ind-1 ts-1)
                  output_fwd(12,ts,L) = digestion_max                              ! max. digestion rate                                           !(J ind-1 ts-1)
                  output_fwd(13,ts,L) = excess                                     !                                                               !(frac.)

                  output_fwd(14,ts,L) = RMR                                        ! active metabolic rate                                         !(J ind-1 ts-1)
                  output_fwd(15,ts,L) = surplus_fwd                                ! surplus rate                                                  !(J ind-1 ts-1)
                  output_fwd(16,ts,L) = surplus_cum                                ! cumulative surplus sum over season                            !(J ind-1 yr-1)
                  output_fwd(17,ts,L) = surplus_cum/mesopf%energy_J_ind(L)         ! cumulative surplus                                            !(J J-1 yr-1)

                  output_fwd(18,ts,L) = predr_indts(z)                             !                                                               !(encounters ts-1)
                  output_fwd(19,ts,L) = survival_fwd                               !                                                               !(survival ts-1)
                  output_fwd(20,ts,L) = survival_cum                               ! cumulative survival sum over season                           !(survival prob.)

                  output_fwd(29,ts,L) = survival_cum*surplus_fwd                   ! annual expected surplus => SUM(surv(t0->t)) * surplus(t)      !(J)

                  output_fwd(21,ts,L) = gut                                        !                                                               !(J)
                  output_fwd(22,ts,L) = gut/mesopf%gut_max(L)                      !                                                               !(fraction of maximal gut)

                  output_fwd(23,ts,L) = RMR + S                                    ! metabolic losses (RMR + S)                                    !(J)
                  output_fwd(24,ts,L) = F + U                                      ! waste losses (F + U)                                          !(J)

                  output_fwd(24,ts,L) = F + U                                      ! waste losses (F + U)                                          !(J)

                  output_fwd(25,ts,L) = reserves_fwd/mesopf%reserves_max(L)        ! reserve fulness (E - state)
                  output_fwd(26,ts,L) = netenergy                                  ! netenergy = digestion - (RMR + S) - (F + U)                   !(J ind-1 ts-1)

                  output_fwd(27,ts,L) = temp_3d(ilat,z,d)
                  output_fwd(28,ts,L) = zoop_3d(lat_index, day, int(depth)+1)

                enddo TIMESTEP_LOOP_fwd

                daily_survival(day) = survival_tmp
                daily_surplus (day) = surplus_tmp

              enddo DAY_LOOP_fwd
            enddo CONV_LOOP_fwd

            !> toc
            write(*,'(A6, A18, A2,I3.3,A1,I3.3,A1)', advance='no')' | ','Forward simulation', '[', conv_bwd, '/' , conv_max_bwd ,']'
            t2 = dclock()
            write(*,'(A4,A,A34)', advance='yes')' -- ',f_hhmmss(idint(t2-t1)), '|'

            !> calculate cohort fitness
            ! write expected surplus counting from every day of the season forward
            COHORT_LOOP: do cohortday = day_start, day_end
              survival_useday = 1.0
              surplus_useday  = 0.0
              do day = cohortday, cohortday+364
                if(day > 365) then
                  useday = day-365
                else
                  useday = day
                end if
                survival_useday =                  survival_useday*daily_survival(useday)
                surplus_useday  = surplus_useday + survival_useday*daily_surplus (useday)
              enddo
              expct_survival(cohortday) = survival_useday
              expct_surplus (cohortday) = surplus_useday
            enddo COHORT_LOOP

            cohort_fitness(:,ilat,L) = expct_surplus(:)/mesopf%energy_J_ind(L)

            !> write cumulative surplus and survival to txt file
            write(510,'(F5.2, 3F8.3, 2F10.3, A20)') warming , latitude_real , mesopf%length_m(L) , output_fwd(20,no_ts,L) , output_fwd(17,no_ts,L) , (sum(expct_surplus(:))/no_days)/mesopf%energy_J_ind(L), ' '// parametrization
            !write(511,'(F5.2, 3F8.3, 2F10.3)') warming , latitude_real , mesopf%length_m(L) , output_fwd(20,no_ts,L) , output_fwd(16,no_ts,L) ,  sum(expct_surplus(:))/no_days


          enddo LENGTH_LOOP_fwd
        end if IFTRUE_fwd

        !> write output to file
        if(save_bwd == .TRUE.) then
          write (filename, '(A,F5.2,A,F3.1,A)'), '\output_bwd__', latitude_real, 'N_warm', warming, 'C.bin'
          call array_to_binary(output_bwd(:,:,:), trim(sub_dir_bin)//trim(filename))
        endif

        if(save_fwd == .TRUE.) then
          write (filename, '(A,F5.2,A,F3.1,A)'), '\output_fwd__', latitude_real, 'N_warm', warming, 'C.bin'
          call array_to_binary(output_fwd(:,:,:), trim(sub_dir_bin)//trim(filename))
        endif

        ! lcz across lats and time
        lcz_3d_real(ilat, :, :) = lcz(:,:)

      enddo LAT_LOOP
    enddo WARMING_LOOP
  end if ENV_ONLY

  !> write to binary
  !! solar elevation
  call array_to_binary(sun_height_2d, trim(sub_dir_bin) // "\sun_height_2d.bin")
  !! surface light
  call array_to_binary(light_surf_2d, trim(sub_dir_bin) // "\light_surf_2d.bin")
  !! light at depth
  call array_to_binary(light_amb_3d , trim(sub_dir_bin) // "\amblight_3d.bin")
  !! temp at depth
  call array_to_binary(temp_3d      , trim(sub_dir_bin) // "\temp_3d.bin")
  !! zooplankton preyfield
  call array_to_binary(zoop_3d      , trim(sub_dir_bin) // "\zoop_3d.bin")
  !! migration behaviour across lats and time
  call array_to_binary(mig_behav_2d_real, trim(sub_dir_bin) // "\mig_behav_2d.bin")
  !! lcz across lats and time
  call array_to_binary(lcz_3d_real, trim(sub_dir_bin) // "\lcz_3d.bin")

  !! TERMINAL SEQUENCE
  
  !> close open connectiond
  close(510)
  !close(511)
  !close(512)
  
  !! elapsed time finish
  run_finish = dclock()

  write(*, '(A)')''
  write(*, '(A4,A39)')'>>', 'Model finished, hit ENTER to terminate'
  write(*, '(A4,A10,A9)') '>', 'Run time:',f_hhmmss(idint(run_finish-run_start))
  print'(A4,A)', '','----------------------------------------------------------------------------'

  write(*,*) achar(7)
  read *

  END PROGRAM main_program_vassild



