DEB-IBM for Abatus cordatus

; the indivduals are given their own agentsets depending on lifestage, and sex for adults:
breed [embryos embryo]  ;; there shouldn't be any: if an embryo appears, there is an issue in the individual energetic dynamics.
breed [juveniles juvenile]
breed [females female]
breed [males male]

;_-_-_-_-_-_-_-_-_
; ------------------------------------------------------------------------------------------------------------------------------------------
globals[  ; global parameters: are accessible for patches and turtles, and the same for all (one change affects all)

;.....................................CALENDARS - TIMERS (from the start of the simulation)........................................................
;repro_time        ;  7     8     9    10     11     0     1     2     3     4     5     6,     7     8 .......based on reproduction timings
;GSI_time          ;  5     6     7     8     9     10    11    12     1     2     3     4,     5     6 .......for accumulation or decrease of GSI
;birth_time        ;  1     0     0     0     0      8     7     6     5     4     3     2,     1     0 .......only runs if triggered by reproduction, for the DEBbirth of individuals
time               ;  1     2     3     4     5      6     7     8     9    10    11    12,    13    14 .......no reseting: time since the start of the run
;;;                ; Oct   Nov   Dec   Jan   Feb   Mar   Apr   May   Jun    Jul   Aug   Sep,   Oct   Nov.......corresponding Gregorian calendar for AH
month_time         ; 10     11    12    01    02    03    04    05    06    07    08    09,    10    11 .......month numeration for the Gregorian calendar

year ; year count from the beginning of the run

;;; Global parameters and their dimensions:

; ...............parameters linked to  temperature and resources:
  f        ; - , scaled functional response
  TC       ; temperature correction factor used to adapt the standard DEB parameters to species-specific site temperature

; ...............population:
  pop_init          ; #, initial population (in a one-patch simulation configuration)
  prop_fem          ; proportion of females on site
  pop_density       ; indiv/m^2, current population density on site
  car_cap           ; carrying capacity for the site
  density_adults    ; #/m^2, density of adults on the patch
  density_juveniles ; #/m^2, density of juveniles on the patch
  surface           ; number of 1m^2 patches, useful in case the model is used with a spatial component ;;;RMV?
  food_compet       ; quantification of the scale of competition

; ...............mortality accounts:
  olddeath          ; count of the number of death due to ageing  (individual-specific)
  bgdeath           ; count of the number of death due to background death (fixed monthly percentage of random death on the population)
  bgdeath_adu       ; same as bgdeath but only counting the adults
  bgdeath_juv       ; same as bgdeath but only counting the juveniles
  starvedeath       ; count of the number of death due to starvation  (individual-specific)
  tempdeath         ; count of the number of death due to temperature (percentage of random death on the population depending on the SST)

; ............... for environmental input parameters:
  ;;for making the in-program LISTS (order according to the position of the column in the file):
  ;Temperature:
  timeMark    ; month number on the first column of the file (0 to 71)
  avgTemp     ; temperature value ;; average over a month for 2012-2018 from PROTEKER data
  avgTemp_fut ; future temperature list (new list created from present list)

  Tbefore    ; for checking mortality: temperature of the month preceding the current tick

  ;Resources:
  infoMonth  ; month number on the first column of the file (1 to 12)
  inputf     ; f value  ;; from Delille1989
  inputf_fut ; future resources list (new list created from present list)


;;; Variables used only in the setup:

; ............... for calculation of the cost of an egg (specific variables used for the simulation of embryo development in calc-embryo):
  e_scaled_embryo
  ref_e_scaled
  U_E_embryo
  S_C_embryo
  U_H_embryo
  L_embryo
  dU_E_embryo
  dU_H_embryo
  dL_embryo
  lower_bound
  upper_bound
  simulation_embryo   ; to count the number of loops

  U_E^start     ; t L^2, initial embryo reserves calculated from the bisection method
  L_0           ; cm, initial structural volume of embryo
  estim_embryo  ; t L^2, initial embryo reserve in the simulation

; ............... for calculation of initial parameters for each age 1 to 5yo (simulation of individual development in calc-init):
  simul          ; to count the number of loops
  timestep_simul ; a different timestep is used in the setup development simulation
  f_simul        ; constant f value for the simulation

  L_simul        ; L, structural length used in the simulation

  U_E_simul
  U_H_simul
  U_H^p_simul
  U_R_simul
  h_rate_simul
  q_acceleration_simul
  dh_rate_simul
  dq_acceleration_simul
  e_scaled_simul
  S_C_simul
  S_A_simul
  dU_E_simul
  dU_H_simul
  dU_R_simul
  dL_simul
  age_months_simul


  U_E_1yo
  U_H_1yo
  U_R_1yo
  L_1yo
  h_rate_1yo
  q_acceleration_1yo

  U_E_2yo
  U_H_2yo
  U_R_2yo
  L_2yo
  h_rate_2yo
  q_acceleration_2yo

  U_E_3yo
  U_H_3yo
  U_R_3yo
  L_3yo
  h_rate_3yo
  q_acceleration_3yo

  U_E_3yo_juv
  U_H_3yo_juv
  U_R_3yo_juv
  L_3yo_juv
  h_rate_3yo_juv
  q_acceleration_3yo_juv

  U_E_4yo
  U_H_4yo
  U_R_4yo
  L_4yo
  h_rate_4yo
  q_acceleration_4yo

  U_E_5yo
  U_H_5yo
  U_R_5yo
  L_5yo
  h_rate_5yo
  q_acceleration_5yo

;;; other global variables (used in GO):
  U_E_0yo   ; Reserve at birth, calculated from the bisection method (embryo simulation)
  L_b       ; (L) structural length at birth

; ............... for plotting and data exploitation: mean values for the state variable for all turtles
  avg_UR_fem     ; lists the mean value of females U_R for each tick
  min_UR_fem     ; lists the min value of females U_R for each tick
  max_UR_fem     ; lists the max value of females U_R for each tick
  avg_dUR_fem    ; lists the mean value of females dU_R for each tick
  min_dUR_fem
  max_dUR_fem
  avg_dUE_fem
  min_dUE_fem
  max_dUE_fem
  avg_UE_fem
  min_UE_fem
  max_UE_fem
  avg_Lphy
  min_Lphy
  max_Lphy
  avg_UE_adu     ; lists the mean value of adults U_E for each tick
  min_UE_adu
  max_UE_adu
  avg_dUE_adu
  min_dUE_adu
  max_dUE_adu
  avg_Lphy_adu
  min_Lphy_adu
  max_Lphy_adu
  avg_UE_juv     ; lists the mean value of juveniles U_E for each tick
  min_UE_juv
  max_UE_juv
  avg_dUE_juv
  min_dUE_juv
  max_dUE_juv
  avg_UE         ; lists the mean value of inidividuals (juveniles & all adults included) U_E for each tick
  min_UE
  max_UE
  avg_dUE
  min_dUE
  max_dUE
  avg_dL
  min_dL
  max_dL
]

patches-own[   ; parameters for the environment (patch-specific):

; ............... DEB parameters for calculation of TC:
  T_ref    ; reference temperature in Kelvin
  T_A      ; Arrhenius temperature
  T        ; actual temperature (in Celsius)
  f-scaled ; functional response f from the input file
]

turtles-own[ ;variables that are turtle-only and turtle-specific: each individual has to change their own value for the variable, as changes only applies to the turtle that is asked to modify it

;;;;.....................................CALENDARS - TIMERS (from the start of the simulation)........................................................
repro_time        ;  7     8     9    10     11     0     1     2     3     4     5     6,     7     8 .......based on reproduction timings
GSI_time          ;  5     6     7     8     9     10    11    12     1     2     3     4,     5     6 .......for accumulation or decrease of GSI
;birth_time       ;  1     0     0     0     0      8     7     6     5     4     3     2,     1     0 .......only runs if triggered by reproduction, for the DEBbirth of individuals
;time             ;  1     2     3     4     5      6     7     8     9    10    11    12,    13    14 .......no reseting: time since the start of the run
;;;               ; Oct   Nov   Dec   Jan   Feb   Mar   Apr   May   Jun    Jul   Aug   Sep,   Oct   Nov.......corresponding Gregorian calendar for AH
;month_time       ; 10     11    12    01    02    03    04    05    06    07    08    09,    10    11 .......month numeration for the Gregorian calendar


; ............... State variables (dimensions):
  L           ; (L) structural length
  dL          ; change of structural length in time
  U_H         ; (t L^2) scaled maturity
  dU_H        ; change of scaled maturity in time
  U_E         ; (t L^2) scaled reserves
  dU_E        ; change of scaled reserves in time
  e_scaled    ; (-) scaled reserves per unit of structure
  U_R         ; (t L^2) scaled energy in reproduction buffer
  dU_R        ; change of energy in reproduction buffer
  Lphy        ; (L), physical length (L/del_M)

  age_months   ; time in month from the birth of the indiv
  age_years    ; time in years from the birht of the indiv


; ............... for ageing:
  q_acceleration  ; (t^-2) ageing acceleration
  dq_acceleration ; change of ageing acceleration in time
  h_rate          ; (t-1) specific death probability rate
  dh_rate         ; change of hazard rate in time
  h_a             ; (t^-2) Weibull ageing acceleration

; ............... Fluxes:
  S_A         ; (L^2) assimilation flux (scaled)
  S_C         ; (L^2) mobilisation flux (scaled)

; ............... Standard DEB parameters (dimensions):
  g           ; (-) energy investment ratio
  v_rate      ; (L/t) energy conductance
  kap         ; (-) fraction of mobilised reserve allocated to soma
  kap_R       ; (-) reproduction efficiency
  k_M_rate    ; (1/t) somatic mainitenance rate coefficient
  k_J_rate    ; (1/t) maturity mainitenance rate coefficient
  U_H^b       ; (t L^2) scaled maturity at birth
  U_H^p       ; (t L^2) scaled maturity at puberty

  scatter-multiplier  ; parameter used to put a random variation in individual parameters

; ............... specific flags for dead or starving state.
  deceased       ; flag to indicate the individual died during the tick (0/ON)
  deceased_old   ; flag to indicate the individual died of senescence (0/ON)
  deceased_bg    ; flag to indicate the individual died due to population bacckground mortality rate (0/ON)

  starvation     ; flag to indicate the individual is starving (ON) or not (0)
]

females-own [
; ............... for the reproduction and laying of offsprings:
  reproduction   ; flag switch to allow reproduction (OFF/ON)
  GSI            ; Gonado-somatic index, needed to switch reproduction on/off
  GSI_init       ; minimum GSI (outside of reproduction)
  dGSI           ; delta GSI
  GSI_max        ; GSI at the beginning of reproduction
  U_R_max        ; U_R at the beginning of reproduction
  eggs           ; number of eggs produced
  Ri             ; reproductive output (for the number of juveniles born)

;.....................................CALENDARS - TIMERS (from the start of the simulation)........................................................
;repro_time        ;  7     8     9    10     11     0     1     2     3     4     5     6,     7     8 .......based on reproduction timings, backbone timer
;GSI_time          ;  5     6     7     8     9     10    11    12     1     2     3     4,     5     6 .......for accumulation or decrease of GSI
birth_time         ;  1     0     0     0     0      8     7     6     5     4     3     2,     1     0 .......only runs if triggered by reproduction, for the DEBbirth of individuals
;time              ;  1     2     3     4     5      6     7     8     9    10    11    12,    13    14 .......no reseting: time since the start of the run
;;;                ; Oct   Nov   Dec   Jan   Feb   Mar   Apr   May   Jun    Jul   Aug   Sep,   Oct   Nov.......corresponding Gregorian calendar for AH
;month_time        ; 10     11    12    01    02    03    04    05    06    07    08    09,    10    11 .......month numeration for the Gregorian calendar

]


; -----------------------------------------------SETUP-----------------------------------------------------------------------------
;_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-INITIALIZATION OF THE MODEL_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

to setup
  clear-all

  ;setup the timers and calendars:
  set time 0
  set month_time 10
  set year 0


  ask patches [
    set T_ref 293.15   ; reference temperature for which TC = 1
    set T_A  9000      ; Arrhenius temperature in tolerance range
    input-temperature  ; get the temperature data from external file, depending on the site chosen
    input-resources    ; procedure to get f data from external files, depending on the site chosen

    if projection != "present" [  ; and modify it if future scenarios are chosen for the projection
      future-temp
      future-rsces
    ]

    setup-temperature  ; to get the environmental varb for the month of the setup
    setup-resources

    set pop_init 120            ; initial population for the site
    set car_cap 200           ; the carrying capacity of the site.
    set prop_fem 0.5   ; proportion of females in the pop, from Poulin1996
 ]

  if add-my-pet? [convert-parameters]  ; to use the model with another species with DEB parameters taken from the DEBtool database add-my-pet

  ; calcluation of initial individuals parameters:
  calc-embryo
  calc-init

  ; setup the pop:
    ask patches [
      sprout pop_init [          ; each patch creates a number of individuals equal to pop_init
        init-birth               ; and each of these individuals initialize their parameters.
      ]
    ]
    init-pop                     ; the population is divided into age classes with age-specific parameters


  ; setup of calendars for reproduction cycles:
  ask turtles [
      set GSI_time 5
      set repro_time 7
    ]

  ask females [
    set GSI 0.03
    set birth_time 0
  ]

reset-ticks
end 

; ----------------------------------------------------------GO--------------------------------------------------------------------------------
;_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-PROCEEDINGS OF THE MODEL_-__-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
;_-_-_-Each round of the model realizes the following procedures in this order and then starts again from the top:

to go
  reset-countdeath      ; reset the death counts. Disabble to obtain a cumulative death plot.

  update-year_timer     ; timers for the month of the year are updated at the beginning for the input of environmental data

;;;Environmental settings:
  ask patches [
    update-resources    ; take the data for f and temperature for this month (calendar is the month_time calendar)
    update-temperature
    ;set f-scaled 1     ; uncomment this line and the next to run the model with DEB standard and constant values for f and T
    ;set T 20

    if competition? [calc-compet]   ; if the competition is switched on in the interface, calculate the effect of competition on f
    set f (f-scaled + food_compet)  ; set the functional response f based on resources availability (external file) and competition
    if f > 1 [set f 1]  ; f is 1 at maximum
    if f < 0 [set f 0] ; f is 0 at minimum
  ]

;;;Individual settings:
  ask turtles [
    if deceased = "on" [die]  ; the individuals that died during the previous timestep are removed
      ;(they are not removed in the previous timestep so that they appear in the plots to monitor their death, this varb was made to monitor ageing effect)

    convert-TC        ; calculate the temperature correction and convert the parameters affected by temperature
    calc-dU_E         ; calculate the changes in state variables that happened over the month
    calc-dU_H
    calc-dU_R
    calc-dL
    calc-ageing       ; calculate the changes in ageing parameters
  ]

  update-individual        ; update state variables of individuals by applying the calculated changes

  ask turtles [update-repro_timers]   ; the timers for reproduction are updated

;;;Females settings (for reproduction):
  ask females [    ; mature females
    update-birth-timer
    if repro_time = 0 and reproduction = "on" [         ; if it is the timing for reproduction, and if the reproduction has been allowed by the GSI the month before,
      set GSI_max GSI           ; set the current GSI as the GSI from which calculations will be made
      set eggs 27
    ]

    calc-GSI                    ; in any case they calculate their GSI each month

    ifelse GSI >= 0.07          ; if the GSI is enough,
     [set reproduction "on"       ; allow reproduction
        if repro_time = 11 [
          set U_R_max U_R         ; keep the last U_R before reproduction period in mind, to force the gradual decrease during reproduction
          set birth_time 9        ; make it known that reproduction is happening
        ]
    ]

     [set reproduction "off"]     ; if the GSI is not higher than the required threshold, the possibility to reproduce is cancelled

    if (birth_time >= 6) and (birth_time <= 8) [  ; if it's one of the three reproduction month,
      reproduce ]                                  ; females produce embryos (by adapting the reprobuffer)

      ; 5 months after conception, the juveniles are DEBborn
    if (birth_time <= 3) and (birth_time >= 1) [
      release-offsprings ]                        ; let the offspring out and give it initial parameters to start dvpmt
  ]

;;;Population data:
  update-pop                     ; update the population by calculating the density
  update-time                    ; update the calendars and timers that haven't been updated yet

  tick

  if count turtles = 0 [         ; stop the model if all individuals died
    user-message ("What have you done ?! The population went extinct! (O_O') ")
    stop]

  if time = run_time * 12 [      ; stop the model after the number years indicated in the interface
    stop]
end 


; ---------------------------------------DESCRIPTION OF THE PROCEDURES---------------------------------------------------------------------------------------------------
; ---------------------------------------Procedures for SETUP: ---------------------------------------------------------------------------------------------------

; ===============
; ====================================Input of environmental data from external file:
; ===============

;_-_-_-_-_-_-_-_-_for temperature:

to input-temperature     ; creating lists for each parameter
  if Sites = "Anse du Halage" [file-open "temp_time_monthavg_Halage.txt"]   ; open the text file in the folder

  ;; the following is for when other sites are chosen: in that case, careful that there needs to have the appropriate resources files for the site in the folder:
  ;/!\for now choosing sites other than Anse du Halage, Port Couvreux or Ile Haute is not possible, the files are not available/!\/!\/!\/!\/!\/!\/!\
  if Sites = "Ile Suhm" [file-open "temp_time_monthavg_Suhm.txt"]
  if Sites = "Port aux Français" [file-open "temp_time_monthavg_PaF.txt"]
  if Sites = "Ile Haute" [file-open "temp_time_monthavg_Haute.txt"]
  if Sites = "Fjord des Portes Noires" [file-open "temp_time_monthavg_PortesNoires.txt"]
  if Sites = "Port Couvreux" [file-open "temp_time_monthavg_Couvreux.txt"]

  set timeMark []    ; create empty list
  set avgTemp []

  while [file-at-end? = false] [                    ; concatenate the empty list and a list containing the next value in the file (the file is read from left to right, and the values taken one by one to add to the list)
    set timeMark sentence timeMark (list file-read)   ; first list with the values for the month
    set avgTemp sentence avgTemp (list file-read)     ; second list (parallel to the first) with the values for temp
  ]
  file-close       ; created a duo of parallel lists with the timestamp and the corresponding temperature
end 

;_-_for the future temperature data:

to future-temp
  ifelse future != "food only" [     ; for projections including the effect of temp changes:
    if projection = "RCP 2.6" [                      ; depending on the chosen scenario in the interface
      set avgTemp_fut map [i -> i + 1.1 ] avgTemp ]    ; create a new list from the current lists by adding the expected changes

    if projection = "RCP 8.5" [
      set avgTemp_fut map [i -> i + 1.7 ] avgTemp ]
  ]
  [                                  ; for projections isolating the effect of food changes:
    if projection = "RCP 2.6" [
      set avgTemp_fut avgTemp ]

    if projection = "RCP 8.5" [
      set avgTemp_fut avgTemp ]
  ]
end 

;_-_then get the values for the month used for setup:

to setup-temperature    ; see 'to-report select-rsces' procedure for explanation
  ifelse projection = "present"
   [set T select-temp timeMark avgTemp 0]     ; gives to T the value for temperature corresponding to the line 0 in the temp list.
   [set T select-temp timeMark avgTemp_fut 0]
end 

to-report select-temp [list1 list2 timing]   ; see 'to-report select-rsces' procedure for explanation:
  let m timing
  if ((timing / 72) >= 1) [              ; (if the number of ticks is higher than the number of values in the temperature list,
    set m (timing - (72 * floor (timing / 72))) ; get the corresponding value counting from the start of the list.)
  ]
  let Q (position m list1)
  report item Q list2
end 


;_-_-_-_-_-_-_-_-_for resources:

to input-resources
  if Sites = "Anse du Halage" [file-open "inputRSces_M_f_Delille.txt"]

  ;; the following is for when other sites are chosen: in that case, careful that there needs to have the appropriate resources files for the site in the folder:
  ;/!\for now the files are not available for other sites: choosing Port Couvreux of Ile Haute will use the file for Anse du Halage, while other site choices will not work at all/!\/!\/!\/!\/!\/!\/!\
  if Sites = "Ile Suhm" [file-open "inputRSces_M_f_Suhm.txt"]
  if Sites = "Port aux Français" [file-open "inputRSces_M_f_PaF.txt"]
  if Sites = "Ile Haute"  [file-open "inputRSces_M_f_Delille.txt"] ; [file-open "inputRSces_M_f_Haute.txt"] if the site-sepcific data is available, exchange the two commands
  if Sites = "Fjord des Portes Noires" [file-open "inputRSces_M_f_PortesNoires.txt"]
  if Sites = "Port Couvreux"  [file-open "inputRSces_M_f_Delille.txt"] ; [file-open "inputRSces_M_f_Couvreux.txt"] if the site-sepcific data is available, exchange the two commands

  set infoMonth []
  set inputf []
  while [file-at-end? = false] [
    set infoMonth sentence infoMonth (list file-read)
    set inputf sentence inputf (list file-read)
  ]
  file-close    ; created a duo of parallel lists with (month) and (f)
end 

;_-_for future resources data:

to future-rsces
  ifelse future != "temp only"          ; for projections including the effect of food changes:
  [
    if projection = "RCP 2.6" [                       ; depending on the chosen scenario in the interface
      set inputf_fut map [i -> i - i * 0.1 ] inputf ] ; create new lists from the current lists by adding the expected changes

    if projection = "RCP 8.5" [
      set inputf_fut map [i -> i - i * 0.2] inputf ]
  ]
  [                                ; for projections isolating the effect of temperature changes:
     if projection = "RCP 2.6" [
      set inputf_fut inputf ]

    if projection = "RCP 8.5" [
      set inputf_fut inputf ]
  ]
end 

;_-_then get the values for the month used for setup:

to setup-resources   ; see 'to report select-rsces' procedure for explanation
  ifelse projection = "present"
  [set f-scaled select-rsces infoMonth inputf month_time ]    ; get the f-scaled from the f list created in input-resources using the month list as an index.
  [set f-scaled select-rsces infoMonth inputf_fut month_time]
end 

;_-_What 'select-rsces' does:

to-report select-rsces [list1 list2 timeval]   ;list1 (infoMonth), list2 (inputf) and timeval (month_time) are indicators of where to take the informations
  let month timeval   ; (month is a temporary variable existing only within the procedure)

  let R (position month list1) ; R temporarily marks the position of the value corresponding to month in the 'list1' (eg. it finds the value of month_time in the list infoMonth and gives to R the line number of this value in the list)
  report item R list2          ; gives the item with index R in list2 = gives the corresponding f value indexed at this month
end 


; ===============
; ====================================Parameters conversion procedures:
; ===============

;_-_-_-_-_-_-_-_-_Conversion of parameters from add-my-pet parameters (for adaptation to other species):

to convert-parameters
  let p_am p_m * zoom / kap_init
  set U_H^b_init E_H^b / p_am
  set U_H^p_init E_H^p / p_am
  set k_M_rate_init p_m / E_G
  set g_init (E_G * v_rate_init / p_am) / kap_init
end 


; ===============
; =============================================Initialisation of the population:
; ===============

;_-_-_-_-_-_-_-_- Bisection method for the calculation of energy of an egg:

to calc-embryo
  set L_0 0.00001              ; (embryos start with almost no structure)
  set lower_bound 0
  set upper_bound 1000        ; the upper bound is given an arbitrary first value,

  set simulation_embryo 0
  loop [
    set simulation_embryo simulation_embryo + 1     ; (to keep track of the number of loops)
    set estim_embryo 0.5 * (lower_bound + upper_bound)
    set L_embryo L_0
    set U_E_embryo estim_embryo
    set U_H_embryo 0
    set e_scaled_embryo v_rate_init * (U_E_embryo / L_embryo ^ 3)
    set ref_e_scaled 1

    while [U_H_embryo < U_H^b_init and e_scaled_embryo > ref_e_scaled]  ; simulation of embryo development until DEBbirth:
    [ set e_scaled_embryo v_rate_init * (U_E_embryo / L_embryo ^ 3)
      set S_C_embryo L_embryo ^ 2 * (((g_init * e_scaled_embryo) / (g_init + e_scaled_embryo)) * (1 + ((L_embryo * k_M_rate_init) / v_rate_init)))
      set dU_E_embryo (-1 * S_C_embryo)                ; (embryo doesn't feed and only mobilizes reserves)
      set dU_H_embryo ((1 - kap_init) * S_C_embryo - k_J_rate_init * U_H_embryo)
      set dL_embryo ((1 / 3) * (((v_rate_init / (g_init * L_embryo ^ 2)) * S_C_embryo) - k_M_rate_init * L_embryo))

      set U_E_embryo U_E_embryo + dU_E_embryo
      set U_H_embryo U_H_embryo + dU_H_embryo
      set L_embryo L_embryo + dL_embryo
      set e_scaled_embryo v_rate_init * (U_E_embryo / L_embryo ^ 3)
    ]
   ; when the embryo reaches DEBbirth (first condition in the 'while' loop),
    if (e_scaled_embryo < ref_e_scaled + 0.05) and (e_scaled_embryo > ref_e_scaled - 0.05) and (U_H_embryo >= U_H^b_init) [
        ; if e is around the reference, the initial reserve is set and the loop stopped:
      set U_E^start estim_embryo  ; the cost of an egg is set.
      set U_E_0yo U_E_embryo      ; the U_E is set for the juveniles initialization.
      stop]

        ; if the reserve is not around 1,
    ifelse U_H_embryo > U_H^b_init
     [set upper_bound estim_embryo]   ; the upper or lower bound is modified accordingly for a new loop.
     [set lower_bound estim_embryo]

    if simulation_embryo > 200            ; the simulation is stopped if it an initial reserve is not found
    [user-message ("embryo calc did not converge")
      stop]
  ]
end 


;_-_-_-_-_-_-_-_-_ Calculation of individual parameters at each age:

to calc-init
  set L_b 0.276 * del_M    ; (phy length from Schatt1985)
  set f_simul 1
  set timestep_simul 2
  set simul 0

  set L_simul L_b
  set U_E_simul U_E_0yo    ; get the initial reserve from the bisection method
  set U_H_simul U_H^b_init
  set U_H^p_simul U_H^p_init
  set U_R_simul 0
  set h_rate_simul 0
  set q_acceleration_simul 0
  set age_months_simul 0

loop [        ; Simulation of an individual development from birth till 5yo:
    set simul simul + 1    ; (to keep track of the number of loops)

    set e_scaled_simul (v_rate_init * (U_E_simul / L_simul ^ 3))
    set S_C_simul L_simul ^ 2 * (g_init * e_scaled_simul / (g_init + e_scaled_simul)) * (1 + ((L_simul * k_M_rate_init) / v_rate_init))
    set S_A_simul f_simul * L_simul ^ 2
    set dU_E_simul (S_A_simul - S_C_simul)
    ifelse U_H_simul < U_H^p_init
     [set dU_H_simul ((1 - kap_init ) * S_C_simul - k_J_rate_init * U_H_simul)
      set dU_R_simul 0 ]
     [set dU_H_simul 0
      set dU_R_simul (( 1 - kap_init) * S_C_simul - k_J_rate_init * U_H^p_init) ]
    set dL_simul ((1 / 3) * (((v_rate_init /( g_init * L_simul ^ 2 )) * S_C_simul) - k_M_rate_init * L_simul))
    set dq_acceleration_simul ((q_acceleration_simul * (L_simul ^ 3 / (v_rate_init / ( g_init * k_M_rate_init)) ^ 3) * sG + h_a_init) * e_scaled_simul * ((v_rate_init / L_simul) - ((3 / L_simul) * dL_simul )) - ((3 / L_simul ) * dL_simul ) * q_acceleration_simul)
    set dh_rate_simul (q_acceleration_simul - ((3 / L_simul) * dL_simul ) * h_rate_simul)

    set U_E_simul U_E_simul + (dU_E_simul * 30.5 / timestep_simul )
    set U_H_simul U_H_simul + (dU_H_simul * 30.5 / timestep_simul )
    set U_R_simul U_R_simul + (dU_R_simul * 30.5 / timestep_simul )
    set L_simul L_simul + (dL_simul * 30.5 / timestep_simul)
    set q_acceleration_simul (q_acceleration_simul + (dq_acceleration_simul * 30.5 / timestep_simul ))
    set h_rate_simul (h_rate_simul + (dh_rate_simul * 30.5 / timestep_simul ))

    set age_months_simul age_months_simul + 1            ; update the age of individuals

; get the individual parameters for each age as the loop goes:
    if age_months_simul = 12 [
      ifelse turtles with [U_H_simul < U_H^p_init] = nobody
       [user-message ("no juveniles at 1 yo, calc-init failed")]
       [set U_E_1yo U_E_simul
        set U_H_1yo U_H_simul
        set U_R_1yo U_R_simul
        set L_1yo L_simul
        set h_rate_1yo h_rate_simul
        set q_acceleration_1yo q_acceleration_simul
      ]
    ]

    if age_months_simul = 24 [
      ifelse turtles with [U_H_simul < U_H^p_init] = nobody
       [user-message ("no juveniles at 2 yo, calc-init failed")]
       [set U_E_2yo U_E_simul
        set U_H_2yo U_H_simul
        set U_R_2yo U_R_simul
        set L_2yo L_simul
        set h_rate_2yo h_rate_simul
        set q_acceleration_2yo q_acceleration_simul
      ]
    ]

    if age_months_simul = 36 [
      ifelse turtles with [U_H_simul < U_H^p_init] = nobody
       [user-message ("no juveniles at 3 yo, calc-init failed")]
       [set U_E_3yo_juv U_E_simul
        set U_H_3yo_juv U_H_simul
        set U_R_3yo_juv U_R_simul
        set L_3yo_juv L_simul
        set h_rate_3yo_juv h_rate_simul
        set q_acceleration_3yo_juv q_acceleration_simul
      ]
    ]

    if age_months_simul = 36 [
      ifelse turtles with [U_H_simul >= U_H^p_init] = nobody
       [user-message ("no adults at 3 yo, calc-init failed")]
       [set U_E_3yo U_E_simul
        set U_H_3yo U_H_simul
        set U_R_3yo U_R_simul
        set L_3yo L_simul
        set h_rate_3yo h_rate_simul
        set q_acceleration_3yo q_acceleration_simul
      ]
    ]

    if age_months_simul = 48 [
      ifelse turtles with [U_H_simul >= U_H^p_init] = nobody
       [user-message ("no adults at 4 yo, calc-init failed")]
       [set U_E_4yo U_E_simul
        set U_H_4yo U_H_simul
        set U_R_4yo U_R_simul
        set L_4yo L_simul
        set h_rate_4yo h_rate_simul
        set q_acceleration_4yo q_acceleration_simul
      ]
    ]

    if age_months_simul = 60 [
      ifelse turtles with [U_H_simul >= U_H^p_init] = nobody
       [user-message ("no adults at 5 yo, calc-init failed")]
       [set U_E_5yo U_E_simul
        set U_H_5yo U_H_simul
        set U_R_5yo U_R_simul
        set L_5yo L_simul
        set h_rate_5yo h_rate_simul
        set q_acceleration_5yo q_acceleration_simul
      ]
    ]

    ifelse  age_months_simul >= 60
     [stop]
     [if simul > 70 [print age_months_simul user-message "calc-init did not end" stop]]
  ]
end 


;_-_-_-_-_-_-_-_- Give initial parameters for newborn individuals:

to init-birth  ; juveniles are initialized as they just reached DEBbirth (=start feeding)
  set v_rate v_rate_init
  set kap kap_init
  set kap_R kap_R_init
  set k_M_rate k_M_rate_init
  set k_J_rate k_J_rate_init

;input a slight variation in parameters between individuals:
  set scatter-multiplier e ^ (random-normal 0 cv)  ; scatter-multiplier is the exp of a number taken randomly on a normal distribution of mean 0 and sd cv(from interface)
  set g g_init / scatter-multiplier
  set U_H^b U_H^b_init / scatter-multiplier
  set U_H^p U_H^p_init / scatter-multiplier
  set U_E U_E_0yo / scatter-multiplier      ; get the initial reserve from the bisection method
  set U_H U_H^b
  set U_R 0
  set L L_b / scatter-multiplier
  set h_rate 0
  set q_acceleration 0

  set age_months 0
  set age_years 0
  set breed juveniles
end 


;_-_-_-_-_-_-_-_- Give initial parameters for all ages:

to init-pop

  let nbtot count turtles     ; at first, all the individuals that were created by 'sprout' are juveniles with 0yo parameters (given by init-birth procedure)

  ask n-of (pop_init / 6) turtles [                      ; divide the initial population into 6 age-groups in a cascade-like manner and give them their parameters.
   ; 5 year-olds initialisation:
    set U_E U_E_5yo / scatter-multiplier     ; parameters come from the individual development simulation calc-init,
    set U_H U_H_5yo / scatter-multiplier     ; including a slight variation in parameters from individual to individual.
    set U_R U_R_5yo / scatter-multiplier
    set L L_5yo / scatter-multiplier
    set h_rate h_rate_5yo / scatter-multiplier
    set q_acceleration q_acceleration_5yo / scatter-multiplier
    set age_months 60
    set age_years 5
  ]

  let Five_yo turtles with [age_months = 60]
  let nb_5yo (count Five_yo)                                       ; for age-groups above puberty,
  ask n-of floor (prop_fem * nb_5yo) Five_yo [set breed females]   ; set the sex-ratio of the population.
  ask Five_yo [if breed != females [set breed males] ]

   ; 4 year-olds initialisation:
  let nb4 count turtles with [age_months < 60]              ; remaining individuals which have not yet been given initial parameters
  ask n-of (pop_init / 6) turtles with [age_months < 60] [
    set U_E U_E_4yo / scatter-multiplier
    set U_H U_H_4yo / scatter-multiplier
    set U_R U_R_4yo / scatter-multiplier
    set L L_4yo / scatter-multiplier
    set h_rate h_rate_4yo / scatter-multiplier
    set q_acceleration q_acceleration_4yo / scatter-multiplier
    set age_months 48
    set age_years 4
 ]

  let Four_yo turtles with [age_months = 48]
  let nb_4yo (count Four_yo)
  ask n-of floor (prop_fem * nb_4yo) Four_yo [set breed females]
  ask Four_yo [if breed != females [set breed males] ]

   ; 3 year-olds initialisation:
  let nb3 count turtles with [age_months < 48]              ; remaining individuals which have not yet been given initial parameters
  ask n-of (pop_init / 6)turtles with [age_months < 48] [
    set U_E U_E_3yo / scatter-multiplier
    set U_H U_H^p
    set U_R U_R_3yo / scatter-multiplier
    set L L_3yo / scatter-multiplier
    set h_rate h_rate_3yo / scatter-multiplier
    set q_acceleration q_acceleration_3yo / scatter-multiplier

    set age_months 36
    set age_years 3
  ]

  let Three_yo turtles with [age_months = 36]
  let nb_3yo (count Three_yo)
  ask n-of floor (prop_fem * nb_3yo) Three_yo [set breed females]
  ask Three_yo [if breed != females [set breed males] ]

   ; 2 year-olds initialisation:
  let nb2 count turtles with [age_months < 36]              ; remaining individuals which have not yet been given initial parameters
  ask n-of (pop_init / 6) turtles with [age_months < 36] [
    set U_E U_E_2yo / scatter-multiplier
    set U_H U_H_2yo / scatter-multiplier
    set U_R U_R_2yo / scatter-multiplier
    set L L_2yo / scatter-multiplier
    set h_rate h_rate_2yo / scatter-multiplier
    set q_acceleration q_acceleration_2yo / scatter-multiplier

    set age_months 24
    set age_years 2
   ]

   ; 1 year-olds initialisation:
  let nb1 count turtles with [age_months < 24]              ; remaining individuals which have not yet been given initial parameters
  ask n-of (pop_init / 6) turtles with [age_months < 24] [
    set U_E U_E_1yo / scatter-multiplier
    set U_H U_H_1yo / scatter-multiplier
    set U_R U_R_1yo / scatter-multiplier
    set L L_1yo / scatter-multiplier
    set h_rate h_rate_1yo / scatter-multiplier
    set q_acceleration q_acceleration_1yo / scatter-multiplier

    set age_months 12
    set age_years 1
   ]

   ; 0 year-olds initialisation:
   ; the remaining individuals keep their 0yo parameters given to them by the init-birth procedure.
end 




; ---------------------------------------Procedures for GO: ---------------------------------------------------------------------------------------------------

;_-_-_-_-_-_-_-_-_Resetting the death counts:

to reset-countdeath
  set starvedeath 0
  set tempdeath 0
  set olddeath 0
  set bgdeath_adu 0
  set bgdeath_juv 0
end 


; ===============
; =============================================Setting the environmental parameters up-to-date:
; ===============

;_-_-_-_-_-_-_-_-_First updating of calendar:

to update-year_timer
  ifelse ticks = 0
   [set month_time month_time]
   [ifelse month_time != 12
    [set month_time month_time + 1]
    [set month_time 1]]            ; resetting the month_time calendar after 12 months
end 


;_-_-_-_-_-_-_-_-_To get the environmental data for the month:

to update-resources   ; see 'to-report select-resources' procedure for explanation
  ifelse projection = "present"      ; depending on the selected projection, take a different set of lists
   [set f-scaled select-rsces infoMonth inputf month_time]
   [set f-scaled select-rsces infoMonth inputf_fut month_time]
end 

to update-temperature
  ifelse projection = "present"       ; depending on the selected projection, take a different set of lists
   [set T select-temp timeMark avgTemp ticks]  ; get the T at the corresponding time of the year, from the temp list created in input-resources using the month list as an index.
   [set T select-temp timeMark avgTemp_fut ticks]
end 


;_-_-_-_-_-_-_-_-_For the effect of competition on the f:

to calc-compet

  calc-density    ; calculate the current month's population density

  ifelse pop_density < 1.9 * car_cap   ; this because if = 2, div by 0 error and if > 2, then the formula gives the untrue result (bigger pop less compet)
   [set food_compet ((1 - f-scaled) * (1 - (pop_density / (2 * car_cap - pop_density))))]
   [set food_compet ((1 - f-scaled) * (1 - (pop_density / (car_cap / 10))))]

;;explanation of the procedure:
 ;(1 - f_scaled) : compet is only effective if the food availability is less than the max.
  ; in which case, proportionnally to how much the f is lessened compared to max, the size of the current pop has an influence on how important the compet is:
   ;  if the pop is below car_cap, then the food is more available for the individuals present, but if the pop is over cc, the availability of food is lessened.
  ; Therefore, if pop > carrying capacity, food_compet < 0 <=> f decreases
   ;           if pop = cc, food_compet = 0 <=> f constant
   ;           if pop < cc, food_compet > 0 <=> f increases

 ;Meaning, the compet is actually calculated depending on how far from the carrying capacity the pop density is, and how far from max f it is.
 ; the pop that isn't there times the food not available  --> the less food there is and the bigger the pop, the higher the compet is.
end 


; ===============
; =============================================Changing individuals states and parameters:
; ===============

;_-_-_-_-_-_-_-_-_Calcualtion of the temperature correction factor:

to convert-TC
  set TC e ^ ((T_A / T_ref) - (T_A / (T + 273.15))) ; standard DEB parameters are for 20°C while A. Cordatus lives in an average SST of 5°C,
  set k_M_rate k_M_rate_init * TC                     ; so TC corrects parameters that are time dependent
  set v_rate v_rate_init * TC
  set k_J_rate k_J_rate_init * TC
  set h_a h_a_init * TC ^ 2
end 


;_-_-_-_-_-_-_-_-_Calculation of the changes in state variables:

;_-_change in energy contained in the reserve:

to calc-dU_E
  if U_H < U_H^b [set f 0]                   ; individuals don't feed before birth
  set e_scaled (v_rate * (U_E / L ^ 3))      ; calculation of scaled reserve
  set S_C L ^ 2 * (g * e_scaled / (g + e_scaled)) * (1 + ((L * k_M_rate) / v_rate))  ; caculation of mobilisation rate (rate at which energy is leaving the reserve)
  set S_A f * L ^ 2                          ; calculation of assimilation rate (rate at which energy is entering the reserve)
  set dU_E (S_A - S_C)
end 

;_-_change in energy contained in the maturity compartment:

to calc-dU_H
  ifelse U_H < U_H^p                        ; individuals supply their maturity compartment only until puberty
   [set dU_H ((1 - kap ) * S_C - k_J_rate * U_H)]
   [set dU_H 0]                              ; after puberty, quantity of energy stored in the maturity compartment stays the same.
end 

;_-_change in energy contained in the reproduction buffer:

to calc-dU_R
  if U_H >= U_H^p [                          ; for mature individuals only.
    ifelse (breed = males) or (birth_time < 6)   ; check if it is a male or a female outside of reproduction time,
     [set dU_R (( 1 - kap) * S_C - k_J_rate * U_H^p)]   ; in which case the change in energy stored in repro buffer is calculated as usual.
     [set dU_R 0]                            ; for females during reproduction, repro buffer calculated in the 'reproduce' procedure.
  ]
  if U_H < U_H^p [set dU_R 0]                ; for non mature individuals, no change in the repro buffer
end 

;_-_change in structural length:

to calc-dL
  set dL ((1 / 3) * (((v_rate /( g * L ^ 2 )) * S_C) - k_M_rate * L))
  if e_scaled < (L / (v_rate / (g * k_M_rate))) [   ; when scaled reserve density is insufficient to maintain the structural length,
    set starvation "on"   ;(starvation flag)
    starve                                      ; the individual is starving (see below for procedure).
  ]
end 

;_-_-_-_-Starvation strategy:

to starve
  set dL 0                  ; allocation of energy only for somatic maintenance
  ifelse U_H < U_H^p        ; all other flux are stopped
   [set dU_H 0]
   [set dU_R 0]
                                                    ; knowing dL = 0, the mobilisation flux is recalculated as S_C = [p_M]L^3 / {p_Am}
  set S_C (kap * (k_M_rate * g * L ^ 3) / v_rate)   ; with [p_M] = [E_G]k_M and [E_G] = gkap{p_Am}/v
  set dU_E S_A - S_C                                ; then dU_E is recalculated and will be applied in the next procedure
  set e_scaled v_rate * U_E / L ^ 3                 ; this new dU_E will affect e_scaled in the next tick

  if e_scaled <= 0 [                             ; when all energy is redirected to soma maintenance, if the scaled reserve still cannot recover
    set starvedeath starvedeath + 1 ; (for mortality accounts)
    set deceased "on"                            ; then the individual dies (flag which will get the individual removed in the next tick).
  ]
end 

;_-_calculation of the damage due to ageing:

to calc-ageing
    set dq_acceleration (q_acceleration * (L ^ 3 / (v_rate / ( g * k_M_rate)) ^ 3) * sG + h_a) * e_scaled * ((v_rate / L) - ((3 / L) * dL)) - ((3 / L ) * dL) * q_acceleration
    set dh_rate (q_acceleration - ((3 / L) * dL) * h_rate)
end 


;_-_-_-_-_-_-_-_-_Apply calculated changes to the individuals:

to update-individual
    ask turtles [
    set U_E U_E + (dU_E * 30.5)   ; changes are based on daily calculations upscaled for a month
    set U_H U_H + (dU_H * 30.5)
    set U_R U_R + (dU_R * 30.5)
    set L L + (dL * 30.5)
    if U_H > U_H^b [               ; damage due to ageing only starts from DEBbirth
      set q_acceleration (q_acceleration + dq_acceleration * 30.5)
      set h_rate (h_rate + dh_rate * 30.5)
    ]

    if deceased != "on"[
      let x random 2
      let y random 1
      if x = y [  ; with the current input of x = random 2 and y = random 1 (i.e. y = 0), this has 50% chance of being true
        if random-float 1 <  ( 1 - (1 - h_rate) ^ 30.5) [    ; damage induces death in a random manner, with a higher probability as it accumulates: here the individual checks that probability
          set olddeath olddeath + 1
          set deceased_old "on" set deceased "on"
        ]
      ]
    ]
  ]

   ;; setting the lifestage label of each individual:
  if any? juveniles with [U_H >= U_H^p] [                      ; juveniles who just reached puberty this round
    let newadults juveniles with [U_H >= U_H^p]                ; are called newadults temporarily
    let nb count newadults                                     ; count how many new adults there are
    ask n-of (prop_fem * nb) newadults [set breed females]     ; then they turn into female or male 50/50 (male priority: eg. in case of one individual, it becomes male)
    ask newadults [if breed != females [set breed males] ]]

  if any? embryos with [U_H >= U_H^b and U_H < U_H^p] [
    ask embryos with [U_H >= U_H^b and U_H < U_H^p] [set breed juveniles]  ; embryos turn into juveniles if they reach DEBbirth
  ]

  ask turtles with [U_H < U_H^b] [set breed embryos]            ; turtles that aren't DEBborn yet are called embryos (but shouldn't appear in the model: if they do, there is an issue in the ODE)

  ask turtles [                     ; update the age of individuals
    set age_months age_months + 1
    set age_years age_years + (1 / 12)
  ]
end 


; ===============
; =============================================Procedures involved in reproduction:
; ===============

;_-_-_-_-_-_-_-_-_Updating of calendars/timers that are used to time reproduction events:

;_-_reproduction and GSI calendars based on reproductive cycle:

to update-repro_timers
  ifelse ticks = 0
   [set repro_time repro_time]  ; at the start, keep the timing of the setup for the first round
   [ifelse repro_time != 11     ; for the rest of the run, the repro_time evolves
    [set repro_time repro_time + 1]
    [set repro_time 0 ]         ; reset when the end of the cycle is reached
  ]

  ifelse ticks = 0
   [set GSI_time GSI_time]      ; at the start, keep the timing of the setup for the first round
   [ifelse GSI_time != 12       ; for the rest of the run, the repro_time evolves
    [set GSI_time GSI_time + 1]
    [set GSI_time 1]            ; reset when the end of the cycle is reached
  ]
end 

;_-_birth and reproduction period timer:

to update-birth-timer               ; birth time, only for females.
  ifelse birth_time = 0             ; only update birth_time if it was triggered by reproduction (if not triggered, it should be 0)
    [set birth_time 0]
    [set birth_time birth_time - 1] ; this timer is a countdown
end 


;_-_-_-_-_-_-_-_-_Calculations of the gonado-somatic index for reproduction:

to calc-GSI
  ifelse GSI_time < 10                      ; if it is outside of the 3 month of egg-laying, calculate the GSI as usual:
   [set f f-scaled
    set GSI (TC * GSI_time * 30.5 * ((k_M_rate_init * g / f ^ 3) / (f + kap * g * y_VE)) * ((1 - kap) * f ^ 3 - k_J_rate_init * U_H^p / zoom ^ 2 / s_M ^ 3))]

  ; if it is in the 3 months of egg-laying, force the gradual decrease of GSI:
   [set GSI_init (GSI_max - GSI_max * 0.52)  ; Magniez1983 females invest 52% into reproduction
    set dGSI ((GSI_max - GSI_init) / 3)      ; change in GSI over a month
    set GSI (GSI - dGSI)]
end 


;_-_-_-_-_-_-_-_-_Reproduction event (conception of eggs):

to reproduce
  set U_R (U_R - ((U_R_max * 0.52) / 3))     ; adapt the reproduction buffer every month of reproduction according to data in Magniez1983
end 


;_-_-_-_-_-_-_-_-_How offsprings are born into the model:

to release-offsprings
  set Ri (0.65 * eggs)             ; 65% of embryos survive (Poulin 1996) until birth
  hatch-juveniles (Ri / 3) [       ; hatch a third of the calculated number of offsprings (1/3 per month over 3 months)
    init-birth                     ; and give them their initial parameters
  ]
end 


; ===============
; =============================================Updates:
; ===============

;_-_-_-_-_-_-_-_-_Update the state of the population:

to update-pop
  background-mortality  ; apply background mortality rates on the population
  ask patches [check-temperature]     ; check mortality induced by critical temperature.
  calc-density   ; calculations for the pop structure
  calc-pop-varb  ; calculations for all individuals variables in the pop
end 

;_-_-_-_-_-_-_-_-_Apply background mortality:

to background-mortality
  let juvzombies count juveniles                           ; temporary variable with the number of juveniles in the model
  ask n-of (round ((0.41 / 12) * juvzombies)) juveniles [  ; 41% of juveniles die each year
    if deceased != "on" [                                  ; check that the individual hasn't already died from another cause
      set bgdeath_juv bgdeath_juv + 1                      ; add to the account of deaths
      set deceased "on"                                    ; indicate the individual is dead
    ]
  ]

  let aduzombies count turtles with [U_H >= U_H^p]         ; temporary variable with the number of adults in the model
  ask n-of (round ((0.24 / 12) * aduzombies)) turtles with [U_H >= U_H^p] [  ; 24% of adults die each year
    if deceased != "on" [                                  ; check that the individual hasn't already died from another cause
      set bgdeath_adu bgdeath_adu + 1                      ; add to the account of deaths
      set deceased_bg "on"                                 ; indicate the individual is dead of this specific cause
      set deceased "on"                                    ; indicate the individual is dead
    ]
  ]
end 

;_-_-_-_-_-_-_-_-_Check mortality by temperature:

to check-temperature
   ; to get the temperature from the month before on the list:
  ifelse projection = "present"                         ; for projections under current conditions
   [ifelse ticks = 0
     [set Tbefore select-temp-before timeMark avgTemp 1]     ; for setup
     [set Tbefore select-temp-before timeMark avgTemp ticks] ; fot the run
  ]
                                                        ; for projections under one of the other scenarios
   [ifelse ticks = 0
     [set Tbefore select-temp-before timeMark avgTemp_fut 1]     ; for setup
     [set Tbefore select-temp-before timeMark avgTemp_fut ticks] ; for the run
  ]

  let zombies count turtles                  ; temporary variable with the number of individuals in the model

;Choice of population with:
 ;;; a lower proportion of individuals sensitive to high temperatures (can withstand for 1 but not 2 months):
  if sensitivity = "resistant" [

  ifelse T >= 12
    [if Tbefore >= 12 [ask turtles [set tempdeath tempdeath + 1 die]]]                              ; if avg temperature is >=12 for two months, 100% die
    [ifelse  T >= 11
     [if Tbefore >= 11 [ask n-of (0.30 * zombies) turtles [set tempdeath tempdeath + 1 die]]]       ; if avg temperature is >=11 for two months, 30% die
     [ifelse T >= 9.5
      [if Tbefore >= 9.5 [ask n-of (0.20 * zombies) turtles [set tempdeath tempdeath + 1 die]]]     ; if avg temperature is >=11 for two months, 20% die
      [if T >= 8
        [if Tbefore >= 8 [ask n-of (0.10 * zombies) turtles [set tempdeath tempdeath + 1 die]]]]    ; if avg temperature is >=11 for two months, 10% die
      ]
    ]
  ]

 ;;; individuals not able to withstand high temperatures for 1 month:
  if sensitivity = "intermediate" [

    ifelse T >= 12
      [ask turtles [set tempdeath tempdeath + 1 die]]                               ; if avg temperature is >=12 for one month, 100% die
      [ifelse T >= 11
        [ask n-of (0.30 * zombies) turtles [set tempdeath tempdeath + 1 die]]         ; if avg temperature is >=12 for one month, 30% die
        [ifelse T >= 9.5
           [ask n-of (0.20 * zombies) turtles [set tempdeath tempdeath + 1 die]]        ; if avg temperature is >=12 for one month, 20% die
           [if T >= 8 [ask n-of (0.10 * zombies) turtles [set tempdeath tempdeath + 1 die]]]]]          ; if avg temperature is >=12 for one month, 10% die
  ]

 ;;; a higher proportion of individuals sensitive to high temperatures (can withstand for 1 but not 2 months):
  if sensitivity = "vulnerable" [

  ifelse T >= 12
    [if Tbefore >= 12 [ask turtles [set tempdeath tempdeath + 1 die]]]                           ; if avg temperature is >=12 for two months, 100% die
    [ifelse  T >= 11
     [if Tbefore >= 11 [ask n-of (0.45 * zombies) turtles [set tempdeath tempdeath + 1 die]]]    ; if avg temperature is >=11 for two months, 45% die
     [ifelse T >= 9.5
      [if Tbefore >= 9.5 [ask n-of (0.35 * zombies) turtles [set tempdeath tempdeath + 1 die]]]  ; if avg temperature is >=10 for two months, 35% die
      [if T >= 8
        [if Tbefore >= 8 [ask n-of (0.25 * zombies) turtles [set tempdeath tempdeath + 1 die]]]] ; if avg temperature is >=9 for two months, 25% die
      ]
    ]
  ]
end 

;_-_reporting procedure to get the temperature from the month before on the list:

to-report select-temp-before [list1 list2 timing]   ; get the avg temp for the month before, that is at position m - 1
  let m timing                          ; create a temporary variable m with the value of the tick
  let P 0                               ; create a temporary variable P
  if ((timing / 72) >= 1) [    ; if the end of the file has been reached,
    set m (timing mod 72)]     ; start counting from the beginning of the file to obtain the temp:

  ifelse m = 0                       ; when it's the first line of the file,
   [set P (position 71 list1)]       ; get the temp from the last line.
   [set P ((position m list1) - 1) ] ; otherwise simply get it from the line before

  report item P list2          ; give the corresponding temp
end 


;_-_-_-_-_-_-_-_-_Get the density data for the month:

to calc-density
  let indiv count turtles with [U_H >= U_H^b]    ; count the number of individuals that are DEBborn
  let space count patches                        ; count the number of 1m2 patches (here, only 1 in the current implementation)
  set pop_density (indiv / space)                ; calculate density of total population
  let adu ((count females) + (count males))
  set density_adults (adu / space)               ; calculate density of adults only
  let juve count juveniles
  set density_juveniles (juve / space)           ; calculate density of juveniles only
end 


;_-_-_-_-_-_-_-_-_Get the state variables data for the month:

to calc-pop-varb
  ifelse any? females [
    set avg_UR_fem mean [U_R] of females              ; reproduction buffer of females
    set min_UR_fem min [U_R] of females
    set max_UR_fem max [U_R] of females

    set avg_dUR_fem mean [dU_R * 30.5] of females
    set min_dUR_fem min [dU_R * 30.5 ] of females
    set max_dUR_fem max [dU_R * 30.5 ] of females

    set avg_UE_fem mean [U_E] of females              ; energy reserve of females
    set min_UE_fem min [U_E] of females
    set max_UE_fem max [U_E] of females

    set avg_dUE_fem mean [dU_E * 30.5] of females
    set min_dUE_fem min [dU_E * 30.5] of females
    set max_dUE_fem max [dU_E * 30.5] of females
  ]
  [    ; if there is no female, set a mock value that can be removed when processing the outputs
    set avg_UR_fem 500
    set min_UR_fem 500
    set max_UR_fem 500

    set avg_dUR_fem 50
    set min_dUR_fem 50
    set max_dUR_fem 50

    set avg_UE_fem 1500
    set min_UE_fem 1500
    set max_UE_fem 1500

    set avg_dUE_fem 100
    set min_dUE_fem 100
    set max_dUE_fem 100
  ]

  ifelse any? turtles with [U_H >= U_H^p] [
    set avg_UE_adu mean [U_E] of turtles with [U_H >= U_H^p]             ; energy reserve of adults only
    set min_UE_adu min [U_E] of turtles with [U_H >= U_H^p]
    set max_UE_adu max [U_E] of turtles with [U_H >= U_H^p]

    set avg_dUE_adu mean [dU_E * 30.5] of turtles with [U_H >= U_H^p]
    set min_dUE_adu min [dU_E * 30.5] of turtles with [U_H >= U_H^p]
    set max_dUE_adu max [dU_E * 30.5] of turtles with [U_H >= U_H^p]


    ask turtles [set Lphy L / del_M]             ; physical length of all individuals
    set avg_Lphy_adu mean [Lphy] of turtles with [U_H >= U_H^p]
    set min_Lphy_adu min [Lphy] of turtles with [U_H >= U_H^p]
    set max_Lphy_adu max [Lphy] of turtles with [U_H >= U_H^p]
  ]
  [     ; if there is no adult, set a mock value that can be removed when processing the outputs afterwards
    set avg_UE_adu 1500
    set min_UE_adu 1500
    set max_UE_adu 1500

    set avg_dUE_adu 100
    set min_dUE_adu 100
    set max_dUE_adu 100

    set avg_Lphy_adu 10
    set min_Lphy_adu 10
    set max_Lphy_adu 10
  ]


  ifelse any? turtles with [U_H >= U_H^b and U_H < U_H^p] [
    set avg_UE_juv mean [U_E] of turtles with [U_H >= U_H^b and U_H < U_H^p]          ; energy reserve of juveniles only
    set min_UE_juv min [U_E] of turtles with [U_H >= U_H^b and U_H < U_H^p]
    set max_UE_juv max [U_E] of turtles with [U_H >= U_H^b and U_H < U_H^p]

    set avg_dUE_juv mean [dU_E * 30.5] of turtles with [U_H >= U_H^b and U_H < U_H^p]
    set min_dUE_juv min [dU_E * 30.5] of turtles with [U_H >= U_H^b and U_H < U_H^p]
    set max_dUE_juv max [dU_E * 30.5] of turtles with [U_H >= U_H^b and U_H < U_H^p]
  ]
  [     ; if there is no juvenile, set a mock value that can be removed when processing the outputs afterwards
    set avg_UE_juv 1500
    set min_UE_juv 1500
    set max_UE_juv 1500

    set avg_dUE_juv 100
    set min_dUE_juv 100
    set max_dUE_juv 100
  ]

  ifelse any? turtles [
    set avg_UE mean [U_E] of turtles             ; energy reserve of all individuals
    set min_UE min [U_E] of turtles
    set max_UE max [U_E] of turtles

    set avg_dUE mean [dU_E * 30.5] of turtles
    set min_dUE min [dU_E * 30.5] of turtles
    set max_dUE max [dU_E * 30.5] of turtles

    set avg_dL mean [dL * 30.5] of turtles
    set min_dL min [dL * 30.5] of turtles
    set max_dL max [dL * 30.5] of turtles

    ask turtles [set Lphy L / del_M]             ; physical length of all individuals
    set avg_Lphy mean [Lphy] of turtles
    set min_Lphy min [Lphy] of turtles
    set max_Lphy max [Lphy] of turtles
  ]
  [     ; if there is no individual, set a mock value that can be removed when processing the outputs afterwards
    set avg_UE 1500
    set min_UE 1500
    set max_UE 1500

    set avg_dUE 100
    set min_dUE 100
    set max_dUE 100

    set avg_dL 1
    set min_dL 1
    set max_dL 1

    set avg_Lphy 10
    set min_Lphy 10
    set max_Lphy 10
  ]
end 


;_-_-_-_-_-_-_-_-_Updating of the other timings:

to update-time
  set time time + 1
  set year year + (1 / 12)
end 





  ;;;;;;;;;;;END;;;;;;;;;;;;;;