m_fish.f90 (16601)

Go to the documentation of this file.

!> @file m_fish.f90
!! The FishMet Model: Definition of fish model.
 
! ------------------------------------------------------------------------------
! Copyright notice:
!   Copyright (c) 2015-2024 Sergey Budaev & Ivar Rønnestad
! FishMet model code is distributed according to GNU AGPL v3, see LICENSE.md
! ------------------------------------------------------------------------------
 
!> This module defines the fish and its components.
!! The model is discrete, time is based on time steps, in s.
module the_fish
 
  use realcalc
  use commondata
  use base_random, only : rand_uni => rand_r4
  use base_utils, only : cspline_scalar
  use environ
  use def_salmon
 
  implicit none
 
  !> The size of the history that saves the fish data back in time.
  integer, parameter, public :: history_size = 20
 
  !> This *type* defines history arrays that are saved for the FISH class.
  type, public :: history_array_fish
    !> Total number of food items ingested
    integer, dimension(HISTORY_SIZE)   :: total_n_food_ingested
    !> Total mass of food ingested
    real(srp), dimension(HISTORY_SIZE) :: total_mass_food_ingested
    !> Total mass of excrements evacuated
    real(srp), dimension(HISTORY_SIZE) :: total_mass_evacuated
    !> Instantaneous absorption history, midgut
    real(srp), dimension(HISTORY_SIZE) :: food_mass_absorb_curr
    !> Cumulativa absorption history, midgut
    real(srp), dimension(HISTORY_SIZE) :: food_mass_absorb_cum
    !> Current energy balance of the fish, history
    real(srp), dimension(HISTORY_SIZE) :: energy_balance_curr
    !> Current body mass, with latest history dynamics.
    real(srp), dimension(HISTORY_SIZE) :: body_mass_current
  end type
 
  !> This defines the root type for the FISH *class*.
  !! @dot
  !!   digraph{
  !!   edge [arrowhead = "empty"]
  !!   node [shape = "record", fontname=Arial]
  !!   FISH_SKEL [ label = "{FISH_SKEL|| new() }" ]
  !!   STOMACH [ label = "{STOMACH|| + st_step() }" ]
  !!   MIDGUT [ label = "{MIDGUT|| + mg_step() }" ]
  !!   FISH [ label = "{FISH|| + decide_eat() }" ]
  !!   FISH -> MIDGUT -> STOMACH -> FISH_SKEL
  !!   }
  !! @enddot
  type, public :: fish_skel
    !> ID number of the fish
    integer :: idnum
    !> A limited history of past indicators (of size ::history_array_fish)
    type(history_array_fish) :: history
 
    ! Energetics parameters (maybe ti a separate class?)
    !> Body mass of the fish, g
    real(srp) :: body_mass0
 
    contains
      !> Initialise an empty fish base. See the fish::fish_init_new_zero().
      procedure, public :: new => fish_init_new_zero
      !> Calculate SMR, see the_fish::fish_skel::smr() and the_fish::smr().
      procedure, public :: smr => fish_smr_temp_bodymass_evan1990_scalar
      !> Calculate whole body energy of the fish,
      !! see the_fish::fish_whole_body_energy_trout()
      procedure, public :: body_energy => fish_whole_body_energy_trout
      !> Get the current body mass of the fish, note that this function is
      !! a shortcut to obtain the last body mass from the history array
      !! (the_fish::fish_skel::history::body_mass_current).
      !! See the_fish::fish_body_mass_current_get().
      procedure, public :: mass => fish_body_mass_current_get
  end type fish_skel
 
  !> Calculate SMR, with generic interface to non-object oriented procedure
  !> @note Using both the_fish::fish_skel::smr() and the_fish::smr() allows
  !!        to call the same `smr` function with both the_fish::fish class
  !!        objects (e.g. `fish_agent%smr(18.0_SRP)`) and raw data
  !!        (e.g. `smr(6560.0_SRP, 18.0_SRP)`).
  interface smr
    module procedure fish_smr_temp_bodymass_evan1990_scalar
    module procedure fish_smr_temp_bodymass_evan1990_value
  end interface
 
  !> Calculate stress pattern based on cubic spline interpolation over
  !! the grid defined by commondata::global_stress_factor_hour and
  !! commondata::global_stress_fact_suppress.
  interface stress_factor
    module procedure stress_suppress_factor_s
    module procedure stress_suppress_factor_v
  end interface stress_factor
 
  ! This type defines a food items that is within the fish stomach
  type, public, extends(food_item) :: food_item_eaten
    !> The time that the food item has spent in the fish stomach
    integer :: time_in_process
    !> Ingestion delay in s (time steps) between ingestion of the food item
    !! by the fish end the start of its ingestion. Water uptake occurs during
    !! this period.
    integer   :: ingestion_delay
 
    contains
      !> Initialise an **empty** food item when it got eaten by the fish
      !! See the_fish::food_item_eaten_set_zero_init_state().
      procedure, public :: zero => food_item_eaten_set_zero_init_state
      !> Initialise a food item when it got eaten by the fish
      !! See the_fish::food_item_eaten_eat_init_at_ingest().
      procedure, public :: init_ingest => food_item_eaten_eat_init_at_ingest
      !> Update the time a food item spent in the fish gastrointestinal tract.
      !! See the_fish::food_item_eaten_update_age_in_stomach().
      procedure, public :: time_update => food_item_eaten_update_age_in_stomach
      !> Check if the food item in the stomach shrunk to zero mass and
      !! disappears.
      !! See the_fish::food_item_eaten_check_destruct().
      procedure, public :: check_disappear => food_item_eaten_check_destruct
  end type food_item_eaten
 
  type, public, extends(food_item_eaten) :: food_item_eaten_mg
    !> The starting mass of the food item at the point it is ingested.
    real(srp) :: mass_0
    !> The cumulative mass processed (absorbed) in midgut.
    real(srp) :: absorped
    !> Duration of the time the food item was in the midgut.
    integer :: time_in_mg
    !> Duration of the time since full absorption occurred.
    integer :: time_absorbed
    contains
      procedure, public :: zero => food_item_eaten_mg_set_zero_init_state
 
      procedure, public :: is_absorbed_full => food_item_eaten_is_absorbed_full
 
      procedure, public :: check_disappear => food_item_eaten_mg_check_destruct
 
      procedure, public :: is_evacuate =>                                     &
                                  food_item_eaten_mg_exceed_max_is_to_destruct
 
  end type food_item_eaten_mg
 
  !> Defines the stomach of the fish
  type, public, extends(fish_skel) :: stomach
    !> Total mass capacity of fish stomach.
    real(srp) :: mass_stomach
    !> Number of food items within the stomach of the fish.
    !! @note Note that commondata::max_food_items_index parameter that
    !!       define the size of the array keeping all food items should be
    !!       big enough to fit all food items.
    integer   :: n_food_items_stomach
    !> An array of ::food_item_eaten objects representing food items in the
    !! stomach
    type(food_item_eaten),dimension(MAX_FOOD_ITEMS_INDEX) :: food_items_stomach
 
    contains
      !> Initialise an empty fish stomach object.
      !! See the_fish::fish_stomach_init().
      procedure, public :: st_init => fish_stomach_init
      !> Update one time step (s) in the fish stomach life cycle.
      !! See the_fish::fish_stomach_update_step().
      procedure, public :: st_step => fish_stomach_update_step
      !> Calculate total mass of food in the fish stomach.
      !! the_fish::fish_stomach_total_mass_food().
      procedure, public :: st_food_mass => fish_stomach_total_mass_food
      !> Calculate the fish appetite factor based on relative stomach fullness.
      !! the_fish::fish_appetite_factor_stomach().
      procedure, public :: appetite_stomach => fish_appetite_factor_stomach
      !> Object-oriented frontend for ::stomach_emptying_time(): Calculate the
      !! stomach emptying time for a fish with specific body mass at a specific
      !! temperature. Calculation is based on spline interpolation of
      !! experimental data.
      !! See the_fish::fish_stomach_emptying_time().
      procedure, public :: emptying => fish_stomach_emptying_time
 
  end type stomach
 
  !> Defines the midgut of the fish
  type, public, extends(stomach) :: midgut
    !> Total mass capacity of the fish midgut
    real(srp) :: mass_midgut
    !> Number of food items within the midgut of the fish.
    !! @note Note that commondata::max_food_items_index parameter that
    !!       define the size of the array keeping all food items should be
    !!       big enough to fit all food items.
    integer :: n_food_items_midgut
    !> An array of ::food_item_eaten objects representing food items in the
    !! midgut
    type(food_item_eaten_mg),dimension(MAX_FOOD_ITEMS_INDEX)::food_items_midgut
 
    contains
      !> Initialise an empty fish stomach object.
      !! See the_fish::fish_midgut_init().
      procedure, public :: mg_init => fish_midgut_init
      !> Update the time the food item spent in the **midgut** of the fish.
      !! See the_fish::fish_midgut_update_from_stomach().
      procedure, public :: mg_step => fish_midgut_update_from_stomach
      !> Calculate the cumulative history of absorption in midgut
      !! See the_fish::fish_midgut_history_update()
      procedure, public :: history_update => fish_midgut_history_update
      !> Calculate the total mass of food in the fish midgut.
      !! See the_fish::fish_midgut_total_mass_food().
      procedure, public :: mg_food_mass => fish_midgut_total_mass_food
      !> Shift food items in the midgut after a new food item is ingested and
      !! therefore food items are shifted in stomach adding the new item.
      !! See the_fish::fish_midgut_shift_items_after_new_ingest().
      procedure, public :: add_midgut=>fish_midgut_shift_items_after_new_ingest
      !> Calculate the fish appetite factor based on relative midgut fullness.
      !! See the_fish::fish_appetite_factor_midgut().
      procedure, public :: appetite_midgut => fish_appetite_factor_midgut
      !> Calculate instantaneous absorption rate.
      !! See the_fish::fish_midgut_absorption_rate().
      procedure, public :: absorption_rate => fish_midgut_absorption_rate
      !> Appetite factor multiplier based on the total absorption rate.
      !! See the_fish::fish_midgut_absorption_factor().
      procedure, public :: absorption_factor => fish_midgut_absorption_factor
      !> Determine the index of the latest absorption item in the history
      !! stack array. Note that in most cases the function should return
      !! commondata::history_size as the history stack is filled with previous
      !! values. See the_fish::midgut_absorp_history_last_past_idx()
      procedure, public :: absorp_last_i => midgut_absorp_history_last_past_idx
      !> Calculate the total mass of food evacuated from the midgut
      !!see the_fish::fish_midgut_total_mass_evacuated_step()
      procedure, public :: evacuated_step=>fish_midgut_total_mass_evacuated_step
      !> Calculate the number of the total mass of food evacuated from the
      !! midgut after full absorption and evacuation delay (see
      !! commondata:global_maximum_duration_midgut_min) that occurs at
      !! one time step.
      !! See the_fish::fish_midgut_n_evacuated_step()
      procedure, public :: evacuated_n => fish_midgut_n_evacuated_step
 
  end type midgut
 
  !> Defines an umbrella class for the complete fish organism
  !! An outline of the ::fish class inheritance is presented below
  !! @dot
  !!   digraph{
  !!     edge [arrowhead = "empty"]
  !!     node [shape = "record", fontname=Arial]
  !!     FISH_SKEL [ label = "{FISH_SKEL|| + new() }" ]
  !!     STOMACH [ label = "{STOMACH|| + st_step() }" ]
  !!     MIDGUT [ label = "{MIDGUT|| + mg_step() }" ]
  !!     FISH [ label = "{FISH|| + decide_eat() }" ]
  !!     FISH -> MIDGUT -> STOMACH -> FISH_SKEL
  !!   }
  !! @enddot
  type, public, extends(midgut) :: fish
    contains
      !> Calculate the baseline locomotor activity of the fish, it is fixed
      !! to the values of commondata::global_baseline_activity_day and
      !! commondata::global_baseline_activity_night for day and night
      !! respectively
      !! See the_fish::fish_activity_baseline_component().
      procedure, public :: activity_baseline => fish_activity_baseline_component
 
      procedure, public :: activity_appetite => fish_activity_appetite_factor
 
      !> Calculate AMR, see the_fish::fish_amr_locomotion_energy_cost()
      procedure, public :: amr => fish_amr_locomotion_energy_cost
      !> Calculate and update the current energy balance of the fish.
      !! See the_fish::fish_energy_balance_update()
      procedure, public :: energy_update => fish_energy_balance_update
      !> Calculate energy energy expenditures and O2 uptake
      !! See the_fish::fish_energy_consumption_step() and
      !! the_fish::fish_oxygen_consumption_step())
      procedure, public :: uptake_energy => fish_energy_consumption_step
      !> Calculate osygen consumption (l) of the fish depending on
      !! temperature (O2 = O2_SMR + O2_AMR + O2_SDA)
      !! See the_fish::fish_oxygen_consumption_step()
      procedure, public :: uptake_o2 => fish_oxygen_consumption_step
      !> Calculate the scaling factor for SDA (specific dynamic action) that
      !! depends on the absorption rate. See the_fish::fish_energy_sda_factor()
      procedure, public :: sda_fact => fish_energy_sda_factor
      !> Calculate bronchial and urine energy
      !! see the_fish::fish_energy_branchial_step()
      procedure, public :: ue_ze => fish_energy_branchial_step
      !> Calculate appetite factor that depends on the energy balance
      !! See the_fish::fish_appetite_factor_energy_balance()
      procedure, public :: appetite_energy=>fish_appetite_factor_energy_balance
      !> Calculate the overall appetite of the fish, which defines the
      !! probability to ingest any food item.
      !! See the_fish::fish_appetite_probability_ingestion().
      procedure, public :: appetite => fish_appetite_probability_ingestion
      !> Making decision to eat a single food item that is available.
      !> See the_fish::fish_decide_eat_food_item_ptr().
      !! @warning This method must be applied to the top-level object
      !!          hierarchy.
      procedure, public :: decide_eat =>  fish_decide_eat_food_item_ptr
      !> The fish does eat a single food item.
      !> See the_fish::fish_ingest_food_item().
      procedure, public :: do_ingest => fish_ingest_food_item
      !> Grow the fish body mass: increment mass based on energy
      !! See the_fish::fish_grow_increment_energy()
      procedure, public :: do_grow => fish_grow_increment_energy
 
  end type fish
 
  contains !--------------------------------------------------------------------
 
  !> Calculate the stomach volume (ml) based on the fish body mass (g)
  !! using Salgado, A., Salgado, I., & Valdebenito, I. (2018). A direct and
  !! straightforward method for measurement real maximum fish stomach volume
  !! to improve aquaculture feeding research. Latin American Journal of
  !! Aquatic Research, 46(5), 880–889.
  !! https://doi.org/10.3856/vol46-issue5-fulltext-3
  elemental function stomach_vol( mass ) result (stomach)
    !> fish body mass
    real(srp), intent(in) :: mass
    real(srp) :: stomach
 
    stomach = 0.0007_srp * mass ** 1.829_srp
 
  end function stomach_vol
 
  !> Initialise an empty ::fish object to an empty zero state.
  elemental subroutine fish_init_new_zero(this, id_num)
    class(fish_skel), intent(inout) :: this
    !> Optional ID number of the fish, default is commondata::unknown
    integer, optional, intent(in) :: id_num
    !> Optional appetite factor level of the fish. Default is
    !! commondata::missing
 
    if (present(id_num)) then
      this%idnum = id_num
    else
      this%idnum = unknown
    end if
 
    this%history%food_mass_absorb_curr = missing
    this%history%food_mass_absorb_cum = missing
    this%history%energy_balance_curr = missing
    this%history%body_mass_current = missing
 
    this%history%total_n_food_ingested = 0
    this%history%total_mass_food_ingested = 0.0_srp
    this%history%total_mass_evacuated = 0.0_srp
 
    this%body_mass0 = global_body_mass
    call add_to_history(this%history%body_mass_current,this%body_mass0)
 
  end subroutine fish_init_new_zero
 
  !> SMR calculation from Evans 1992 data, *scalar* version
  pure function fish_smr_temp_bodymass_evan1990_scalar(this, temperature,     &
                                        is_per_hour, time_step) result(smr_out)
    class(fish_skel), intent(in) :: this
    !> Ambient temperature ºC
    real(srp), intent(in) :: temperature
    !> Optional flag if calculation is done per hour (default FALSE,
    !! meaning per one time step, second)
    logical, optional, intent(in) :: is_per_hour
    !! Time step, optional parameter.
    integer(LONG), optional, intent(in) :: time_step
    !> Standard metabolic rate (SMR) in mg O2 / kg per unit time (s or h)
    real(srp) :: smr_out
 
    logical :: is_hour_loc
    integer(LONG) :: step_loc
 
    if (present(is_per_hour)) then
      is_hour_loc = is_per_hour
    else
      is_hour_loc = .false.
    end if
 
    if (present(time_step)) then
      step_loc = time_step
    else
      step_loc = global_time_step
    end if
 
    smr_out=fish_smr_temp_bodymass_evan1990_value( this%mass(),               &
                                                   temperature,               &
                                                   is_hour_loc, step_loc )
 
  end function fish_smr_temp_bodymass_evan1990_scalar
 
  !> SMR calculation from Evans 1992 data, *non-object scalar* version
  pure function fish_smr_temp_bodymass_evan1990_value(body_mass, temperature, &
                                        is_per_hour, time_step) result(smr_out)
    !> Body mass of the fish, kg
    real(srp), intent(in) :: body_mass
    !> Ambient temperature ºC
    real(srp), intent(in) :: temperature
    !> Optional flag if calculation is done per hour (default FALSE, per
    !! second)
    logical, optional, intent(in)   :: is_per_hour
    !! Time step, optional parameter.
    integer(LONG), optional, intent(in) :: time_step
    !> Standard metabolic rate (SMR) in mg O2 / kg per unit time (s or h)
    real(srp) :: smr_out
 
    logical :: is_hour_loc
    integer(LONG) :: step_loc
 
    if (present(is_per_hour)) then
      is_hour_loc = is_per_hour
    else
      is_hour_loc = .false.
    end if
 
    if (present(time_step)) then
      step_loc = time_step
    else
      step_loc = global_time_step
    end if
 
    if ( is_hour_loc ) then
      !                  mg O2                       kg
      smr_out = oxygen_rate_std(temperature)*g2kg(body_mass)
    else
      !                  mg O2                       kg           per s
      smr_out = oxygen_rate_std(temperature)*g2kg(body_mass) / real(hour, srp)
    end if
 
    if (is_stress_enable)                                                     &
          smr_out = smr_out +                                                 &
                    smr_out * global_stress_cost_smr * stress_suppress(step_loc)
 
  end function fish_smr_temp_bodymass_evan1990_value
 
  !> Calculate the active metabolic rate of fish in units of oxygen
  !! consumption **per hour**.
  !! @f[ AMR = a M^{b} U^{c} @f]
  !!
  !! See Ohlberger, J., Staaks, G., Van Dijk, P. L. M., & Hölker, F. (2005).
  !! Modelling energetic costs of fish swimming. Journal of Experimental
  !! Zoology Part A: Comparative Experimental Biology, 303, 657–664.
  elemental function amr_swim_basic(fish_mass,swimming_speed,time_unit)             &
                                                            result(amr_h_out)
    !> Fish mass, g
    real(srp), intent(in) :: fish_mass
    !> Fish swimming speed cm/s
    real(srp), intent(in) :: swimming_speed
    !> Optional time unit, e.g. commondata::minute or commondata::hour
    !! The default value is s 1/60 of commondata:::minute
    integer, intent(in), optional :: time_unit
    !> Active metabolic rate, mg O2 per `time_unit`
    !! (equation default commondata::hour)
    real(srp) :: amr_h_out
 
    integer :: time_unit_loc
 
    ! default parameters of the AMR model
    ! TODO: Set new parameters
    real(srp), parameter :: a = 0.24_srp
    real(srp), parameter :: b = 0.80_srp
    real(srp), parameter :: c = 0.60_srp
 
    if (present(time_unit)) then
      time_unit_loc = time_unit
    else
      time_unit_loc = 1
    end if
 
    associate( m => fish_mass, u => swimming_speed, amr => amr_h_out )
 
      amr = a * m**b * u**c
 
    end associate
 
    ! We need to adjust the AMR according to the time_unit as the equation
    ! (with the empirical parameters) calculates per hour
    amr_h_out = time_unit_loc * amr_h_out/hour
 
  end function amr_swim_basic
 
  !> Calculate the active metabolic rate. This can be based on
  !! ::amr_swim_basic() or def_salmon::salmon_oxygen_consumption().
  elemental function fish_amr_locomotion_energy_cost(this,                    &
                    temperature, is_exclude_smr, is_per_hour) result (amr_out)
    class(fish), intent(in) :: this
    !> Ambient temperature, C
    real(srp), optional, intent(in) :: temperature
    !> Opional flag to force exclude SMR (if SMR is included into the equation
    !! making this function essentially SMR + AMR function).
    !! Note that SMR and AMR are added in the_fish::fish::energy_update().
    !! The default value is FALSE.
    logical, optional, intent(in) :: is_exclude_smr
    !> Optional flag if calculation is done per hour (default FALSE,
    !! meaning per one time step, second)
    logical, optional, intent(in) :: is_per_hour
    !> Active metabolic rate (AMR), mg O2 / s
    real(srp) :: amr_out
 
    logical :: is_exclude_smr_loc, is_per_hour_loc
    real(srp) :: current_activity, temp_loc
 
    if (present(is_exclude_smr)) then
      is_exclude_smr_loc = is_exclude_smr
    else
      is_exclude_smr_loc = .false.
    end if
 
    if (present(is_per_hour)) then
      is_per_hour_loc = is_per_hour
    else
       is_per_hour_loc = .false.
    end if
 
    if (present(temperature)) then
      temp_loc = temperature
    else
      temp_loc = global_temperature
    end if
 
    current_activity = this%activity_baseline() * this%activity_appetite()
 
    if (is_per_hour_loc) then
      amr_out = salmon_oxygen_consumption(                                    &
                              g2kg(this%mass()), temp_loc, current_activity,  &
                              time_unit=hour )
    else
      amr_out = salmon_oxygen_consumption(                                    &
                              g2kg(this%mass()), temp_loc, current_activity )
    end if
 
    if (is_exclude_smr_loc) amr_out = amr_out - this%smr(temp_loc)
 
  end function fish_amr_locomotion_energy_cost
 
  !> Calculate the whole body energy content.
  elemental function fish_whole_body_energy_trout(this) result (energy_out)
    class(fish_skel), intent(in) :: this
    !> Energy content of the body, kJ per whole fish
    real(srp) :: energy_out
    !                                    %%                       kJ/g
    real(srp), parameter :: protein_cont=0.15_srp, protein_energy=24.0_srp,   &
                            fat_cont    =0.23_srp, fat_energy    =35.0_srp
 
    real(srp) :: protein_fish, fat_fish ! kJ/fish
 
    !> @f$ E = M P_c P_e + M F_c F_e @f$
    protein_fish = this%mass() * protein_cont * protein_energy
    fat_fish     = this%mass() * fat_cont     * fat_energy
 
    energy_out = protein_fish + fat_fish
 
  end function fish_whole_body_energy_trout
 
  !> Calculate the body mass from whole body energy.
  !! @note Note that this procedure is the reverse of
  !!       the_fish::fish_whole_body_energy_trout().
  elemental function energy2mass( body_energy ) result (mass_get)
    !> Whole body energy of the fish, kJ
    real(srp), intent(in) :: body_energy
    !> Body mass, g
    real(srp) :: mass_get
    !                                    %%                       kJ/g
    real(srp), parameter :: protein_cont=0.15_srp, protein_energy=24.0_srp,   &
                            fat_cont    =0.23_srp, fat_energy    =35.0_srp
 
    !> @f$ M = \frac {E} { P_c P_e + F_c F_e } @f$
    mass_get = body_energy /                                                  &
                    ( protein_cont * protein_energy + fat_cont * fat_energy )
 
  end function energy2mass
 
  !> Initialise an empty food item as ingested by the fish. Cleanup and set
  !! zero state.
  !! @warning This procedure is by default used as a **destructor**, it sets
  !!          the food item  an undefined state (missing, unknown).
  elemental subroutine food_item_eaten_set_zero_init_state(this,              &
                                            mass, start_mass, delay)
    class(food_item_eaten), intent(inout) :: this
    !> Optional mass of the food item, default is set to
    !! commondata::global_food_item_mass
    real(srp), optional, intent(in) :: mass
 
    !> Optional starting reference mass of the food item.
    !! @note Not used here
    real(srp), optional, intent(in) :: start_mass
 
    !> Optional ingestion delay, default is set from the parameter
    !! commondata::global_ingestion_delay
    integer, optional, intent(in) :: delay
 
    if (present(mass)) then
      call this%init( mass )
    else
      call  this%init( missing )
    end if
 
    this%time_in_process = 0
 
    if (present(delay)) then
      this%ingestion_delay = delay
    else
      this%ingestion_delay = global_ingestion_delay_min * minute
    end if
 
  end subroutine food_item_eaten_set_zero_init_state
 
  !> Initialise an empty food item as ingested by the fish. Cleanup and set
  !! zero state.
  !! @warning This procedure is by default used as a **destructor**, it sets
  !!          the food item  an undefined state (missing, unknown).
  elemental subroutine food_item_eaten_mg_set_zero_init_state(this,           &
                                            mass, start_mass, delay)
    class(food_item_eaten_mg), intent(inout) :: this
    !> Optional mass of the food item, default is set to
    !! commondata::global_food_item_mass
    real(srp), optional, intent(in) :: mass
 
    !> Optional starting reference mass of the food item.
    real(srp), optional, intent(in) :: start_mass
 
    !> Optional ingestion delay, default is set from the parameter
    !! commondata::global_ingestion_delay
    integer, optional, intent(in) :: delay
 
    ! Local copies of optionals.
    real(srp) :: mass_loc
    integer :: delay_loc
 
    if (present(mass)) then
      mass_loc = mass
    else
      mass_loc = missing
    end if
 
    if (present(delay)) then
      delay_loc = delay
    else
      delay_loc = global_digestion_delay_min * minute
    end if
 
    ! initialise the stomach components
    call food_item_eaten_set_zero_init_state(this,mass_loc,missing,delay_loc)
    !call this%zero(mass_loc, MISSING, delay_loc) ! segmentation fault
 
    if (present(start_mass)) then
      this%mass_0 = start_mass
    else
      this%mass_0 = this%mass
    end if
 
    this%absorped = 0.0_srp
 
    this%time_in_mg = 0
 
    this%time_absorbed = 0
 
  end subroutine food_item_eaten_mg_set_zero_init_state
 
  !> Initialise a single food item as it is ingested by the fish.
  !! Note that the ingested food item is copied from the external food item
  !! environ::food_item that comes as the input argument.
  elemental subroutine food_item_eaten_eat_init_at_ingest(this, food_item_in, &
                                                                  ingest_delay)
    class(food_item_eaten), intent(inout) :: this
    !> The food item object (class) being ingested by the fish.
    class(food_item), intent(in) :: food_item_in
    !> Optional ingestion delay, if not provided the default value is
    !! commondata::global_ingestion_delay.
    !! @note Ingestion delay is assumed a property of the fish
    !!       gastrointestinal system rather than the food item.
    integer, optional, intent(in) :: ingest_delay
 
    !> - copy all data components from the input argument `food_item_in`
    !!   to `this`. Note that class assignment operator overloading is used
    !!   for `=` assignment with backend defined in
    !!   environ::food_item_copy_assign_to().
    !    @note The code is equivalent to ` this%mass = food_item_in%mass `
    !          see notes on ::copy().
    call this%copy( food_item_in )
 
    !> - The time counter of the food item within the gastrointestinal tract
    !!   of the fish is set to zero.
    !! .
    this%time_in_process = 0
 
    if (present(ingest_delay)) then
      this%ingestion_delay = ingest_delay
    else
      this%ingestion_delay = global_ingestion_delay_min * minute
    end if
 
  end subroutine food_item_eaten_eat_init_at_ingest
 
  !> This procedure updates the time the food item spent in the
  !! **stomach** of the fish.
  elemental subroutine food_item_eaten_update_age_in_stomach(this, t_increment)
    class(food_item_eaten), intent(inout) :: this
    !> Optional time step increment, default value is 1.
    integer, optional, intent(in) :: t_increment
 
    integer :: inc_loc
 
    if (present(t_increment)) then
      inc_loc = t_increment
    else
      inc_loc = 1
    end if
 
    this%time_in_process = this%time_in_process + inc_loc
 
  end subroutine food_item_eaten_update_age_in_stomach
 
  !> Check if the food item mass shrunk to near-zero and destruct it.
  elemental subroutine food_item_eaten_check_destruct( this,                  &
                                    is_destructed, tolerance, max_time )
    class(food_item_eaten), intent(inout) :: this
    !> Optional logical indicator that is set to TRUE if the food item was
    !! actually destructed.
    logical, optional, intent(out) :: is_destructed
    !> Optional numerical tolerance, the smallest value of the food item mass
    !! that counts as non-zero. Default is set by commondata::min_vol_toler.
    real(srp), optional, intent(in) :: tolerance
 
    !> Maximum time parameter, not used here, for compatibility with the
    !! respective midgut method
    integer, optional, intent(in) :: max_time
 
    ! Local copy of optional
    real(srp) :: toler_loc
 
    if (present(tolerance)) then
      toler_loc = tolerance
    else
      toler_loc = min_mass_toler
    end if
 
    if ( this%mass < toler_loc ) then
      call this%zero()
      if (present(is_destructed)) is_destructed = .true.
    else
      if (present(is_destructed)) is_destructed = .false.
    end if
 
  end subroutine food_item_eaten_check_destruct
 
  !> Check if the food item has been processed for more than
  !! commondata::global_maximum_duration_midgut s time in the midgut.
  elemental subroutine food_item_eaten_mg_check_destruct( this,               &
                                        is_destructed, tolerance, max_time )
    class(food_item_eaten_mg), intent(inout) :: this
    !> Optional logical indicator that is set to TRUE if the food item was
    !! actually destructed.
    logical, optional, intent(out) :: is_destructed
    !> Optional numerical tolerance parameter. Note that it is not used here
    !! but needed for compatibility with the food_item_eaten::zero() method
    !! implemented for the base type the_fish::food_item_eaten.
    real(srp), optional, intent(in) :: tolerance
 
    !> Maximum time processed in the midgut, by default all food items that
    !! have been in the process of absorption for more than this time will
    !! destructed from midgut and evacuated.
    integer, optional, intent(in) :: max_time
 
    integer :: max_time_loc
 
    if (present(max_time)) then
      max_time_loc = max_time
    else
      max_time_loc = global_maximum_duration_midgut_min * minute
    end if
 
    if ( this%is_evacuate(max_time_loc) ) then
      !> The food item is destructed with the food_item_eaten_mg::zero()
      !! function.
      call this%zero()
      if (present(is_destructed)) is_destructed = .true.
    else
      if (present(is_destructed)) is_destructed = .false.
    end if
 
  end subroutine food_item_eaten_mg_check_destruct
 
  !> Check if the food item in the midgut has spent for more than the maximum
  !! time in the midgut after full absorption and should be destructed and
  !! evacuated (TRUE).
  elemental function food_item_eaten_mg_exceed_max_is_to_destruct(this,       &
                                                  max_time) result(is_destruct)
    class(food_item_eaten_mg), intent(in) :: this
    !> Maximum time processed in the midgut, by default all food items that
    !! have been in the process of absorption for more than this time will
    !! destructed from midgut and evacuated.
    !! Default value is based on commondata::global_maximum_duration_midgut_min.
    integer, optional, intent(in) :: max_time
 
    logical :: is_destruct
 
    integer :: max_time_loc
 
    if (present(max_time)) then
      max_time_loc = max_time
    else
      max_time_loc = global_maximum_duration_midgut_min * minute
    end if
 
    if ( this%time_absorbed >= max_time_loc ) then
      is_destruct = .true.
    else
      is_destruct = .false.
    end if
 
  end function food_item_eaten_mg_exceed_max_is_to_destruct
 
  !> Check if the food item has been absorbed to the maximum possible mass
  elemental function food_item_eaten_is_absorbed_full(this, max_absorb)       &
                                                            result(is_absorbed)
    class(food_item_eaten_mg), intent(in) :: this
 
    real(srp), optional, intent(in) :: max_absorb
 
    logical :: is_absorbed
 
    real(srp) :: max_absorb_loc
 
    check_invalid: if (this%mass_0==missing.or.this%absorped==0.0_srp) then
      is_absorbed = .false.
      return
    end if check_invalid
 
    if (present(max_absorb)) then
      max_absorb_loc = max_absorb
    else
      max_absorb_loc = this%mass_0 * global_absorption_ratio
    end if
 
    if ( this%absorped <= max_absorb_loc - min_mass_toler ) then
      is_absorbed = .false.
    else
      is_absorbed = .true.
    end if
 
  end function food_item_eaten_is_absorbed_full
 
  !-----------------------------------------------------------------------------
  !> This function defines the mass of a single food item as it passes
  !! through the fish stomach.
  elemental function st_food_item_mass ( time_in_process, dry_mass,           &
                    ingestion_delay, water_uptake, rescale ) result (volume_out)
    !> The time since the food item has been ingested by the fish
    integer, intent(in) :: time_in_process
    !> Optional dry mass of the food item, if absent the default value
    !! commondata::global_food_item_mass is used.
    real(srp), optional, intent(in) :: dry_mass
    !> Ingestion delay: the time duration from ingestion of the food item
    !! until it starts transport into midgut. Water uptake occurs during
    !! the ingestion delay. Default is commondata::global_ingestion_delay.
    integer, optional, intent(in) :: ingestion_delay
    !> The relative mass of water added to the food item after ingestion
    !! during the ingestion delay. Default is commondata::global_water_uptake.
    real(srp), optional, intent(in) :: water_uptake
    !> Optional time rescaling parameter for the stomach transport array
    real(srp), optional, intent(in) :: rescale
    !> Returns the mass of the food item.
    real(srp) :: volume_out
 
    ! Local copies of optional arguments
    real(srp) :: dry_mass_loc, water_uptake_loc, rescale_loc
    integer :: ingestion_delay_loc
 
    ! maximum mass of the food item after water uptake
    real(srp) :: max_vol
 
    ! Process optional parameters, default values are set is absent.
    if (present(dry_mass)) then
      dry_mass_loc = dry_mass
    else
      dry_mass_loc = global_food_item_mass
    end if
 
    if (present(ingestion_delay)) then
      ingestion_delay_loc = ingestion_delay
    else
      ingestion_delay_loc = global_ingestion_delay_min * minute
    end if
 
    if (present(water_uptake)) then
      water_uptake_loc = water_uptake
    else
      water_uptake_loc = global_water_uptake
    end if
 
    if (present(rescale)) then
      rescale_loc = rescale
    else
      rescale_loc = 1.0_srp
    end if
 
    !> ### Implementation notes ###
    !> In stomach, the feed that is ingested is passed through two stages:
    !!  - ingestion delay (commondata::global_ingestion_delay) when water
    !!    uptake occurs and the mass of the food item **increases**:
    !!  - passage to midgut where the mass of the food item **decreases**.
    !!  .
    !!
    !> The food item is subjected to water uptake. The pattern of water
    !! uptake over time, the maximum the_fish::food_item_eaten::ingestion_delay
    !! seconds, is defined by a commondata::logistic() equation:
    !! @f[ v_t = c_{max} + \frac{c_{max}-c_{0}}{1+a \cdot e^{-r \cdot x}} @f]
    !! where @f$ c_{0} @f$ is the dry mass of the food item, as defined by
    !! the parameter commondata::global_food_item_mass, @f$ c_{max} @f$
    !! is the asymptotic maximum mass of the food item after water intake,
    !! @f$ a @f$ and @f$ r @f$ are logistic function parameters (see
    !! commondata::logistic()).
    !!
    !> - The asymptotic maximum mass of the food item after water uptake
    !!   @f$ c_{max} @f$ is defined by
    !!   @f[ c_{max} = c_{0} + c_{0} \cdot u @f]
    !!   where @f$ u @f$ is the water uptake parameter
    !!   commondata::global_water_uptake.
    !! .
    max_vol = dry_mass_loc + dry_mass_loc * water_uptake_loc
 
    !> The value of the function is calculated as follows.
    !> - If the time the food item spend in the stomach is smaller than the
    !!   ingestion delay (commondata::global_ingestion_delay), its mass
    !!   **increases** through water uptake up to the `max_vol` value
    !!   @f$ c_{max} @f$.
    if ( time_in_process <= 0 ) then
      volume_out = global_food_item_mass
    else if ( time_in_process < ingestion_delay_loc ) then
      associate( t => real(time_in_process, srp),                            &
                  cm => max_vol, c0 => dry_mass_loc,                          &
                  a => global_water_uptake_a,                                 &
                  r => global_water_uptake_r )
        volume_out = c0 + logistic( t, a, cm-c0, r )
      end associate
    !> - If the time the food item spent in stomach is longer than the
    !!   ingestion delay, its mass **decreases** (from the maximum mass
    !!   after water intake) as determined by the interpolation function that
    !!   is determined by the grid:
    !!     - abscissa: commondata::global_transport_pattern_t;
    !!     - ordinate: commondata::global_transport_pattern_r.
    !!     .
    !!   given the time after water intake completed (ingestion delay
    !!   commondata::global_ingestion_delay).
    !! .
    else
      volume_out =                                                            &
         max_vol *                                                            &
            cspline_scalar( real(global_transport_pattern_t, srp)*rescale_loc,&
                            global_transport_pattern_r,                       &
                            real(time_in_process - ingestion_delay_loc, srp) )
    end if
 
  end function st_food_item_mass
 
  !> Initialise an empty fish stomach object.
  !! @warning Note that the fish itself is not initialised here, it should be
  !!          done with the ::fish::new() method.
  elemental subroutine fish_stomach_init(this, mass)
    class(stomach), intent(inout) :: this
 
    real(srp), optional, intent(in) :: mass
 
    if(present(mass)) then
      this%mass_stomach = mass
    else
      this%mass_stomach = global_stomach_mass
    end if
 
    this%n_food_items_stomach = 0
 
    !> Note that at the init, all food items in the stomach are destructed
    !! with the ::food_item_eaten::zero() method.
    call this%food_items_stomach%zero()
 
  end subroutine fish_stomach_init
 
  !> Update one time step of the fish stomach life cycle
  impure subroutine fish_stomach_update_step(this)
    class(stomach), intent(inout) :: this
 
    ! Note: Fixed array sized to object is not allowed in elemental procedure:
    ! `logical, dimension( this%n_food_items_stomach ) :: is_destruct`
    logical, allocatable, dimension(:) :: is_destruct
 
    debug_check_max_stom: if (is_debug) then
      if (this%n_food_items_stomach >                                         &
                               max_food_items_index * debug_warn_level )  then
        write(*,*) "DEBUG_CHECK_MAX_STOM: N of food items in stomach ",       &
                   this%n_food_items_stomach," exceeds ",                     &
                   nint(max_food_items_index * debug_warn_level)
      end if
      if (this%n_food_items_stomach > max_food_items_index)  then
        write(*,*) "DEBUG_CHECK_MAX_STOM: N of food items in stomach ",       &
                   this%n_food_items_stomach," exceeds MAX_FOOD_ITEMS_INDEX=",&
                   max_food_items_index
        stop
      end if
    end if debug_check_max_stom
 
    allocate( is_destruct(this%n_food_items_stomach) )
 
    is_destruct = .false.
 
    !> ### Implementation notes ###
    associate( stom_foods                                                    &
                => this%food_items_stomach(1:this%n_food_items_stomach) )
 
      !> First, the time counter for each of the existing (non-zero)
      !! `n_food_items_stomach` food item (`time_in_process`) is updated, by
      !! default incremented by 1.
      call stom_foods%time_update()
 
      !> Second, the mass of each food item in stomach is recalculated
      !! according to its individual `time_in_process` counter.
      stom_foods%mass =  st_food_item_mass( stom_foods%time_in_process, rescale=stomach_adjust() )
 
      !> Check if the food items shrunk below the minimum size and should then
      !! disappear (be destructed). The `is_destruct` boolean indicators show
      !! how many food items should be destructed
      call stom_foods%check_disappear(is_destructed=is_destruct)
      this%n_food_items_stomach = this%n_food_items_stomach - count(is_destruct)
 
    end associate
 
    !> In the debug mode, also check and fix counter for negative values.
    if (is_debug) then
      if ( this%n_food_items_stomach < 0 ) this%n_food_items_stomach = 0
      if ( this%n_food_items_stomach >= max_food_items_index )                &
                            this%n_food_items_stomach = max_food_items_index
    end if
 
  end subroutine fish_stomach_update_step
 
  !> This procedure shows what happens when the fish eats a single food item.
  elemental subroutine fish_ingest_food_item(this, food_item_ingested)
    class(fish), intent(inout) :: this
    !> An environ::food_item class object that represents the food item
    !! being ingested.
    class(food_item), intent(in) :: food_item_ingested
 
    type(food_item_eaten) :: new_food_item_in
 
    !> - Individual totals are updated: N of food items ingested and the
    !!   total mass of food ingested.
    call add_to_history( this%history%total_n_food_ingested,                  &
                         last(this%history%total_n_food_ingested) + 1 )
 
    call add_to_history( this%history%total_mass_food_ingested,               &
                         last(this%history%total_mass_food_ingested) +        &
                                                    food_item_ingested%mass )
 
    !> - The newly ingested food item object is initialised, data structure
    !!   is copied from the input argument food_item object. Ingestion delay
    !!   is set default.
    call new_food_item_in%init_ingest( food_item_ingested )
 
    !> - The food item is then added to the stack array of the fish
    !!   stomach. It takes the first place shifting all other food items in
    !!   sequence.
    this%food_items_stomach(2:this%n_food_items_stomach+1) =                  &
                            this%food_items_stomach(1:this%n_food_items_stomach)
 
    this%food_items_stomach(1) = new_food_item_in
 
    !> - The food items counter is incremented.
    this%n_food_items_stomach = this%n_food_items_stomach + 1
 
    !> - Finally, whenever the food is eaten, initialise it in the midgut.
    call this%add_midgut( food_item_ingested )
 
  end subroutine fish_ingest_food_item
 
  !> Calculate the total mass of food in the fish stomach.
  elemental function fish_stomach_total_mass_food(this, tolerance)            &
                                                        result(total_mass)
    class(stomach), intent(in) :: this
    !> Optional numerical tolerance, the smallest value of the food item mass
    !! that counts as non-zero. Default is set by commondata::min_vol_toler.
    real(srp), optional, intent(in) :: tolerance
    !> Total mass of food in the stomach.
    real(srp) :: total_mass
 
    ! Local copy of optional
    real(srp) :: toler_loc
 
    if (present(tolerance)) then
      toler_loc = tolerance
    else
      toler_loc = min_mass_toler
    end if
 
    !> The total value is calculated as an array sum to the total number of
    !! food items being kept in the stomach (`this%n_food_items_stomach`).
    associate( foods_stomach                                                 &
                    => this%food_items_stomach(1:this%n_food_items_stomach) )
      total_mass = sum( foods_stomach%mass,                                  &
                            foods_stomach%mass > toler_loc .and.              &
                            foods_stomach%mass /= missing )
    end associate
 
  end function fish_stomach_total_mass_food
 
  !> Calculate the fish appetite factor based on the relative stomach fullness.
  !! Appetite factor is equivalent to the probability of consuming any food
  !! item.
  elemental function fish_appetite_factor_stomach(this, fcapacity)           &
                                                              result(appetite)
    class(stomach), intent(in) :: this
    real(srp), optional, intent(in) :: fcapacity
    real(srp) :: appetite
 
    real(srp) :: fcapacity_loc
 
    real(srp) :: mass_relative
 
    if (present(fcapacity)) then
      fcapacity_loc = fcapacity
    else
      fcapacity_loc = global_stomach_mass
    end if
 
    !> ### Implementation notes ###
    !! The appetite factor ranges between (0 and 1) and is based on the
    !! relative filling capacity. Here relative filling capacity @f$ c_r @f$
    !! is calculated as the ratio of the total mass of food in the stomach
    !! @f$ c_{max} @f$ to the maximum filling capacity @f$ c_t @f$
    !! @f[ c_r = \frac{c_t}{c_{max}} @f]
    mass_relative = this%st_food_mass()/fcapacity_loc
 
    !> Given the relative filling capacity @f$ c_r @f$, appetite factor
    !! that ranges from 0 to 1 is calculated as logistic function
    !! ::appetite_func().
    !!
    !! The @f$ a @f$ and @f$ r @f$ logistic parameters are defined by the
    !! commondata::global_appetite_logist_a and
    !! commondata::global_appetite_logist_r parameters.
    !!
    !! *gnuplot code:*
    !! @verbatim
    !!    set xrange [0:1]
    !!    plot 1 - 1/(1+500000*exp(-20*x))
    !! @endverbatim
    appetite = appetite_func( mass_relative, global_appetite_logist_a,        &
                                            global_appetite_logist_r )
 
  end function fish_appetite_factor_stomach
 
 
  !> Initialise an empty fish midgut object.
  !! @warning Note that the fish itself is not initialised here, it should be
  !!          done with the ::fish::new() method.
  elemental subroutine fish_midgut_init(this, mass)
    class(midgut), intent(inout) :: this
 
    real(srp), optional, intent(in) :: mass
 
    if(present(mass)) then
      this%mass_midgut = mass
    else
      this%mass_midgut = global_midgut_mass
    end if
 
    this%n_food_items_midgut = 0
 
    !> Note that at the init, all food items in the stomach are destructed
    !! with the ::food_item_eaten::zero() method.
    call this%food_items_midgut%zero()
 
  end subroutine fish_midgut_init
 
  !> This procedure processes the food item in the **midgut** of the fish.
  impure subroutine fish_midgut_update_from_stomach(this,dry_mass,t_increment)
    class(midgut), intent(inout) :: this
    !> Optional dry mass of the food item, if absent the default value
    !! commondata::global_food_item_mass is used.
    real(srp), optional, intent(in) :: dry_mass
    !> Optional time step increment, default value is 1.
    integer, optional, intent(in) :: t_increment
 
    ! Local copies of optional arguments
    real(srp) :: dry_mass_loc
    integer :: inc_loc
 
    integer :: i
 
    ! Absorption variables: absorption increment and the mass of the food
    ! prior to absorption
    real(srp) :: absorbed_i , mass_0
 
    ! Minimum mass of a food item at transferred to the midgut when the
    ! timer counting the duration of the time spent in the midgut starts.
    ! See commondata::min_vol_toler for details.
    real(srp), parameter :: threshold_min_vol = min_mass_toler
 
    ! An array of food item volumes that transits to midgut from stomach at
    ! the current time step, the delta
    real(srp), dimension(MAX_FOOD_ITEMS_INDEX) :: delta_vol
 
    ! Note: Fixed array sized to object is not allowed in elemental procedure:
    ! `logical, dimension( this%n_food_items_stomach ) :: is_destruct`
    logical, allocatable, dimension(:) :: is_destruct
 
    ! This debugging check is conducted to make sure the maximum number
    ! of food items in the gastrointestinal system is never exceeded.
    ! Whenever it exceeds commondata::debug_warn_level of the maximum, a
    ! warning message is issued.
    debug_check_max_midgut: if (is_debug) then
      if (this%n_food_items_midgut >                                          &
                               max_food_items_index * debug_warn_level )  then
        write(*,"(a,i4,a,i4)")                                                &
                   "DEBUG_CHECK_MAX_MIDGUT: N of food items in midgut ",      &
                   this%n_food_items_midgut," exceeds ",                      &
                   nint(max_food_items_index * debug_warn_level)
      end if
      if (this%n_food_items_midgut > max_food_items_index)  then
        write(*,"(a,i4,a,i4)")                                                &
                   "DEBUG_CHECK_MAX_MIDGUT: N of food items in midgut ",      &
                   this%n_food_items_midgut, " exceeds MAX_FOOD_ITEMS_INDEX=",&
                   max_food_items_index
        stop
      end if
    end if debug_check_max_midgut
 
    allocate( is_destruct(this%n_food_items_midgut) )
    is_destruct = .false.
    delta_vol = missing
 
    ! Process optional parameters, default values are set is absent.
    if (present(dry_mass)) then
      dry_mass_loc = dry_mass
    else
      dry_mass_loc = global_food_item_mass
    end if
 
    if (present(t_increment)) then
      inc_loc = t_increment
    else
      inc_loc = 1
    end if
 
    !> ### Implementation notes ###
    !> #### Transfer to midgut ####
    !> Each food item smoothly transfers from the stomach to the midgut, which
    !! is implemented as removing a small proportion if its mass with
    !! simultaneous addition of the same mass to the linked midgut data
    !! component.
    !!
    !! First, for each valid food item in the stomach, calculate the difference
    !! between mass of the food item in stomach at the previous time step
    !! *t-1* and the current step *t*: @f$ \Delta v_{i} @f$ (`delta_vol`).
    !! @f[ \Delta v_{i} = S_{i-1} - S_{i}  @f]
    !! where @f$ S_{i} @f$ is the mass of the *i*-th food item in the stomach.
    !! This difference is equal to the mass increment that transfers to
    !! the midgut.
    !!
    !! This increment is set to zero for all food items that stayed in the
    !! stomach for less than the `ingestion_delay` duration (this period is
    !! water uptake, no transfer, see stomach::st_step()).
    associate(stom_foods =>                                                  &
                         this%food_items_stomach(1:this%n_food_items_stomach))
 
      where ( stom_foods%time_in_process > stom_foods%ingestion_delay .and.   &
              stom_foods%time_in_process /= unknown  )
        delta_vol(1:this%n_food_items_stomach) =                              &
            st_food_item_mass(stom_foods%time_in_process-1, dry_mass_loc, rescale=stomach_adjust() ) -   &
            st_food_item_mass(stom_foods%time_in_process,   dry_mass_loc, rescale=stomach_adjust() )
      elsewhere
        delta_vol(1:this%n_food_items_stomach) = 0.0_srp
      end where
 
    end associate
 
    !> The increment @f$ \Delta v_{i} @f$ (`delta_vol`) is here also set
    !! to zero for all food items that have transferred fully from the
    !! stomach to midgut.
    delta_vol(this%n_food_items_stomach+1:this%n_food_items_midgut) = 0.0_srp
 
    !> Second, the mass of each food item already in the midgut is
    !! augmented by the respective increment @f$ \Delta v_{i} @f$.
    !! @f[ V_{i} + \Delta v_{i}  @f]
    !! The age of each food item (`time_in_process` data component) in the
    !! midgut is also incremented by one (i.e. discretization period, s).
    associate( midgut_foods =>                                               &
                         this%food_items_midgut(1:this%n_food_items_midgut) )
 
      midgut_foods%mass = midgut_foods%mass +                                 &
                                        delta_vol(1:this%n_food_items_midgut)
      midgut_foods%time_in_process = midgut_foods%time_in_process + inc_loc
 
      !> The timer for the duration spent in the midgut is also incremented.
      where ( midgut_foods%mass >= threshold_min_vol .and.                    &
                                                        delta_vol > 0.0_srp )
        midgut_foods%time_in_mg = midgut_foods%time_in_mg + 1
      end where
 
    end associate
 
    !> #### Absorption ####
    !> The **absorption** processis this: certain proportion of the mass of
    !! each food item in the midgut is subtracted following the equation:
    !! @f[ V_{i} - V_{i} \frac{r_{max} \sum V{i}}{K_{MM} + \sum V{i}} @f]
    !! where @f$ r_{max} @f$ and @f$ K_{MM} @f$ are the Michaelis-Menten
    !! equation parameters: the maximum rate and the K constant
    !! (see commondata::michaelis_menten() for more details). Here
    !! @f$ \sum V{i} @f$ is the total mass of food in the midgut.
    !!
    !! Only the food items that stay in the midgut for more than the
    !! duration defined by commondata::global_digestion_delay_min are subject
    !! to absorption. Also, there is a maximum limit on the absorption set
    !! by the commondata::global_absorption_ratio model parameter.  This limit
    !! is checked in the food_item_eaten_mg::is_absorbed_full() method.
    !!
    !! @note Note that each food item in the midgut is subjected to
    !!       Michaelis-Menten absorption depending on the **total** mass
    !!       of food in the midgut @f$ \sum V{i} @f$.
    do_absorp: do i=this%n_food_items_stomach+1, this%n_food_items_midgut
 
      associate( midgut_foods => this%food_items_midgut(i),                  &
                  r => this%mg_food_mass() / this%mass_midgut,                &
                  rmax =>  global_mid_gut_mm_r_max,                           &
                  km => global_mid_gut_mm_k_m )
 
        ! Time counter in midgut is incremented for all absorping food items
        midgut_foods%time_in_mg = midgut_foods%time_in_mg + 1
 
        debug_check_full_midgut_delta: if (is_debug) then
          if ( delta_vol(i) > 0.0_srp  ) then
            write(*,*) "DEBUG_CHECK_FULL_MIDGUT_DELTA: `delta_vol` nonzero [",&
                        i, "], ", delta_vol(i), ", certainly eror"
          end if
        end if debug_check_full_midgut_delta
 
        if ( (midgut_foods%time_in_mg > midgut_foods%ingestion_delay) .and.   &
             (.not. midgut_foods%is_absorbed_full()) ) then
 
          mass_0 = midgut_foods%mass
          absorbed_i = mass_0 * michaelis_menten( r, rmax, km )
 
          midgut_foods%mass = midgut_foods%mass - absorbed_i *                &
                                        midgut_temp_factor(global_temperature)
 
          if ( midgut_foods%mass < midgut_foods%mass_0 -                      &
                    midgut_foods%mass_0 * global_absorption_ratio  ) then
 
            midgut_foods%mass = midgut_foods%mass_0 -                         &
              midgut_foods%mass_0 * global_absorption_ratio
 
            midgut_foods%absorped = midgut_foods%mass_0*global_absorption_ratio
 
          else
 
            midgut_foods%absorped =                                           &
                        min( midgut_foods%mass_0 * global_absorption_ratio,   &
                                    midgut_foods%absorped +                   &
                                            ( mass_0 - midgut_foods%mass ) )
          end if
 
        end if
 
      end associate
 
    end do do_absorp
 
    debug_digest_checks: if (is_debug) then
      associate( midgut_foods =>                                             &
                      this%food_items_midgut(1:this%n_food_items_midgut) )
        if (any(midgut_foods%ingestion_delay/=global_digestion_delay_min      &
                                                            * minute)) then
          write(*,*) "DEBUG_DIGEST_CHECKS: midgut digestion delay mismatch: ",&
                     global_digestion_delay_min, midgut_foods%ingestion_delay
          stop 255
        end if
        if ( any(midgut_foods%absorped >                                     &
                       midgut_foods%mass_0 * global_absorption_ratio ) ) then
          write(*,*) "DEBUG_DIGEST_CHECKS: exceed absorption MAX ratio at ",  &
                    global_time_step , " step, count ",                       &
                    count(midgut_foods%absorped >                             &
                        midgut_foods%mass_0 * global_absorption_ratio ),      &
                    ", total Ns ", this%n_food_items_stomach, "/",            &
                    this%n_food_items_midgut, " stomach/midgut"
        end if
      end associate
    end if debug_digest_checks
 
    !> Finally, the food_item_eaten_mg::check_disappear() method is called for
    !! each food item in the midgut that has fully transferred from
    !! the stomach. This checks if the item has been processed in the midgut
    !! or more than commondata::global_maximum_duration_midgut_min. Those food
    !! items are deleted from the midgut and are excreted.
    associate( foods_midgut_not_stom =>                                      &
                  this%food_items_midgut(                                     &
                      this%n_food_items_stomach+1:this%n_food_items_midgut) )
 
      !> The time counter for the time since full absorption (`time_absorbed`
      !! data component) also updates (incremented by 1) for all food items
      !! that are fully absorbed.
      where ( foods_midgut_not_stom%is_absorbed_full() )
        foods_midgut_not_stom%time_absorbed =                                 &
                                    foods_midgut_not_stom%time_absorbed + 1
      end where
 
      debug_digest_disappear: if ( is_debug ) then
        if ( any( foods_midgut_not_stom%time_absorbed >                       &
                         global_maximum_duration_midgut_min * minute ) ) then
          write(*,*) "DEBUG_DIGEST_DISAPPEAR: must remove ",                  &
                      count( foods_midgut_not_stom%time_absorbed >            &
                                global_maximum_duration_midgut_min * minute ),&
                     " at time ", global_time_step
        end if
      end if debug_digest_disappear
 
      !> Update the individual totals for the mass of food evacuated.
      call add_to_history( this%history%total_mass_evacuated,                 &
                           last(this%history%total_mass_evacuated) +          &
                                                        this%evacuated_step() )
 
 
      !> Check if the food items must be evacuated from the midgut, disappear
      call foods_midgut_not_stom%check_disappear( is_destructed=is_destruct )
 
      !> The count of the food items in the midgut is updated after reduction
      !! at the previous step.
      this%n_food_items_midgut = this%n_food_items_midgut-count( is_destruct )
 
      debug_check_subzero_food: if (is_debug) then
        if (  any( foods_midgut_not_stom%mass < 0.00000001_srp .and.          &
                   foods_midgut_not_stom%mass /= missing ) ) then
           write(*,*) "DEBUG_CHECK_SUBZERO_FOOD: midgut food items are " //   &
                      "suspiciously small"
           write(*,*) "    ", foods_midgut_not_stom%mass
        end if
 
      end if debug_check_subzero_food
 
    end associate
 
    ! Calculate the cumulative history of absorption in midgut
    call this%history_update()
 
  end subroutine fish_midgut_update_from_stomach
 
  !> Calculate and update the current energy balance of the fish. Note that
  !! the energy balance is kept in the history array
  !! the_fish::fish_skel::history::energy_balance_curr.
  elemental subroutine fish_energy_balance_update (this, temperature, start_value)
    class(fish), intent(inout) :: this
    !> Optional current ambient temperature, if the parameter absent, use
    !! the global default commondata::global_temperature.
    real(srp), optional, intent(in) :: temperature
    !> Optional start value of the energy balance, if the parameter is absent,
    !! the nefault value oz zero is used
    real(srp), optional, intent(in) :: start_value
 
    integer :: hist_length, absorb_last_item
    real(srp) :: absorp_increment
 
    real(srp) :: temp_loc, start_loc
 
    real(srp), parameter :: balance_start = 0.0_srp ! starting energy balance
    real(srp) :: balance_current
 
    if (present(temperature)) then
      temp_loc = temperature
    else
      temp_loc = global_temperature
    end if
 
    if (present(start_value)) then
      start_loc = start_value
    else
      start_loc = balance_start
    end if
 
    balance_current = temp_loc
 
    hist_length = size(this%history%food_mass_absorb_cum) ! HISTORY_SIZE
 
    ! The history is initialized to commondata::missing values, so it is easy
    ! to determine what is the last valid item in the history stack array.
    absorb_last_item = this%absorp_last_i()
 
    !> ### Implementation notes
    !!
    !! @f[ E_{i+1} = E_{i} + F(\Delta a) - E_{SMR} - E_{AMR} , @f]
    !! where @f$ E_{i} @f$ is the energy balance of the fish at time @f$ i @f$
    !! (history array: the_fish::history_array_fish::energy_balance_curr),
    !! @f$ \Delta a @f$ is the energy that came from increment in the food
    !! absorption @f$ a @f$,  @f$ E_{SMR} @f$ is the energetic equivalent
    !! of the standard metabolic rate based on the_fish::fish_skel::smr(),
    !! and @f$ E_{AMR} @f$ is the energetic equivalent of the active metabolic
    !! rate bsed on the_fish::fish_skel::amr().
    if ( absorb_last_item == 1 ) then
      call add_to_history(this%history%energy_balance_curr, start_loc -       &
                                energy_o2(o2vol(mg2g( this%smr(temp_loc) ))) )
    else
      ! Cumulative absorption increment is obtained from history stack
      absorp_increment = this%history%food_mass_absorb_cum(hist_length) -     &
                              this%history%food_mass_absorb_cum(hist_length-1)
 
      balance_current = this%history%energy_balance_curr(hist_length)         &
                         + feed_energy(absorp_increment)                      &
                         - this%uptake_energy(temp_loc)
 
      call add_to_history( this%history%energy_balance_curr, balance_current )
    end if
 
  end subroutine fish_energy_balance_update
 
  !> Calculate urinal (UE) and branchial energy (ZE) as a proportion of
  !! the SMR
  elemental function fish_energy_branchial_step(this, temperature, is_oxygen) &
                                                                  result(ue_ze)
    class(fish), intent(in) :: this
    !> Optional current ambient temperature, if the parameter absent, use
    !! the global default commondata::global_temperature.
    real(srp), optional, intent(in) :: temperature
    !> Optional logical flag that defines that oxygen consumption (mass, g)
    !! rather than energy is output, default is FALSE (energy is calculated)
    !! @note Note that oxygen consumption in case of `is_oxygen=TRUE` is
    !!       output as mass in g, NOT volume.
    logical, optional, intent(in) :: is_oxygen
    real(srp) :: ue_ze
 
    real(srp) :: temp_loc
    logical :: set_is_oxygen
 
    if (present(temperature)) then
      temp_loc = temperature
    else
      temp_loc = global_temperature
    end if
 
    if (present(is_oxygen)) then
      set_is_oxygen = is_oxygen
    else
      set_is_oxygen = .false.
    end if
 
    if (set_is_oxygen) then
      ue_ze = mg2g(this%smr(temp_loc)) * global_ue_ze_factor
    else
      if (global_is_ue_ze_fixed_rate)  then
        ! Test fixed N mass based on Forsberg 1996: 36 mg TAN kg fish day
        ! Forsberg, O.I., 1997. The impact of varying feeding regimes on
        ! oxygen consumption and excretion of carbon dioxide and nitrogen in
        ! post-smolt Atlantic salmon Salmo salar L. Aquac Research 28, 29–41.
        ! https://doi.org/10.1046/j.1365-2109.1997.00826.x
        !ue_ze=
        !  branchial_nitrogen_to_energy(0.036_SRP*g2kg(this%mass())/(HOUR*24))
        ue_ze =                                                               &
          branchial_nitrogen_to_energy(                                       &
                             nitrogen_mass(                                   &
                                from_micro(global_ue_ze_ammonia_excretion) )  &
                                                            *this%mass()/hour )
      else
        ue_ze = energy_o2(o2vol(mg2g(this%smr(temp_loc)))) * global_ue_ze_factor
      end if
    end if
 
  end function fish_energy_branchial_step
 
  !> Calculate the mass of ammonia, mol to g
  elemental function nitrogen_mass(mol, is_ammonia) result(gram)
    real(srp), parameter :: ammonia_molar_mass=17.03052_srp   ! 17.03052 g/mol
    real(srp), parameter :: nitrogen_molar_mass=14.006747_srp ! 14.006747 g/mol
    !> Amount in mol
    real(srp), intent(in) :: mol
    !> Optional flag indicating that
    logical, intent(in), optional :: is_ammonia
    real(srp) :: gram
    logical :: use_ammonia
 
    if (present(is_ammonia)) then
      use_ammonia = is_ammonia
    else
      use_ammonia = .false.
    end if
 
    if (use_ammonia) then
      gram = mol * ammonia_molar_mass
    else
      gram = mol * nitrogen_molar_mass
    end if
 
  end function nitrogen_mass
 
  !> Calculate the branchal and urinary energy loss from branchial and
  !! urinary nitrogen loss.
  !! See Bureau, D.P., Kaushik, S.J., Cho, C.Y. (2002) Bioenergetics.
  !! In Hardy, R.W., Halver, J.E., eds. Fish Nutrition, 3 ed
  !! Academic Press, San Diego, CA, USA (pp 1–59):
  !!
  !! <<Because most nitrogen losses are as ammonia, and the difference in
  !! the amount of energy loss per gram of nitrogen between ammonia and
  !! urea is small, it has been proposed that the loss of1gofnitrogen by
  !! fish under normal conditions can be equated to an energy loss
  !! of 24.9 kJ.>> (p.21)
  !!
  !! see also
  !! Saravanan, S., Schrama, J. W., Figueiredo-Silva, A. C.,
  !! Kaushik, S. J., Verreth, J. A. J., & Geurden, I. (2012). Constraints on
  !! Energy Intake in Fish: The Link between Diet Composition, Energy
  !! Metabolism, and Energy Intake in Rainbow Trout. PLoS ONE, 7(4), e34743.
  !! https://doi.org/10.1371/journal.pone.0034743
  elemental function branchial_nitrogen_to_energy(excreted_nitrogen, is_kg)   &
                                                                    result(bue)
    !> Branchial and urinary nitrogen loss, g
    real(srp), intent(in) :: excreted_nitrogen
    !> Optional flag indicating that excreted nitrogen loss is given in kg
    logical, optional, intent(in) :: is_kg
    !> Energy loss kJ
    real(srp) :: bue
    logical :: is_kg_loc
    !> Energy equivalent kJ of 1 g excreted nitrogen NH3-N
    real(srp), parameter :: energ_g_nitro = 24.85_srp
 
    if (present(is_kg)) then
      is_kg_loc = is_kg
    else
      is_kg_loc = .false.
    end if
 
    bue = excreted_nitrogen * energ_g_nitro
 
    if (is_kg_loc) bue = bue / 1000.0_srp
 
  end function branchial_nitrogen_to_energy
 
  !> Calculate *total energy, kJ consumption* of the fish depending on the
  !! temperature E = SMR + AMR + (UE+ZE) + SDA
  !! Public interface: the_fish::uptake_energy().
  !!
  !! @note This function should be adapted together with the oxygen
  !!       consumption code ::fish_oxygen_consumption_step().
  elemental function fish_energy_consumption_step(this, temperature, sda_off) &
                                                                result(kj_out)
    class(fish), intent(in) :: this
    real(srp), optional, intent(in)  :: temperature
    !> Optional logical switch that sets SDA component OFF, by default FALSE,
    !! intended for use mainly for debugging and testing.
    logical, optional, intent(in) :: sda_off
    real(srp) :: kj_out
 
    real(srp) :: temp_loc
    logical :: sda_off_debug
 
    if (present(temperature)) then
      temp_loc = temperature
    else
      temp_loc = global_temperature
    end if
 
    if (present(sda_off)) then
      sda_off_debug = sda_off
    else
      sda_off_debug = .false.
    end if
 
    !! @f[ e = E_{SMR} + E_{AMR} + (E_{UE}+E_{ZE}) + E_{SDA} , @f]
    kj_out = energy_o2(o2vol(mg2g(this%smr(temp_loc)))) +             & ! SMR
             energy_o2(o2vol(mg2g(this%amr(temp_loc)))) +             & ! AMR
             this%ue_ze(temp_loc)                                       ! UE+ZE
 
    if (.NOT. sda_off_debug) then
      kj_out = kj_out +                                               &
             energy_o2(o2vol(mg2g(this%smr(temp_loc)))) *             & ! SDA
                                    this%sda_fact(global_sda_absorp_rate_max, &
                                                  global_sda_factor_max)
    end if
 
  end function fish_energy_consumption_step
 
  !> Calculate *total oxygen consumption* (l) of the fish depending on the
  !! temperature O2 = O2_SMR + O2_AMR + O2_{UE+ZE} O2_SDA
  !! @note This function should be adapted together with the energy
  !!       consumption code ::fish_energy_consumption_step().
  elemental function fish_oxygen_consumption_step(this, temperature,          &
                                            is_mass, sda_off) result (o2_out)
    class(fish), intent(in) :: this
    real(srp), optional, intent(in)  :: temperature
    !> Optional logical switch, if TRUE, determining that oxygen consumption is
    !! calculated as mass (g), the default is volume (l)
    logical, optional, intent(in) :: is_mass
    !> Optional logical switch that sets SDA component OFF, by default FALSE,
    !! intended for use mainly for debugging and testing.
    logical, optional, intent(in) :: sda_off
    real(srp) :: o2_out
 
    ! Is UE+ZE equivalent O2 included in consumption, default is NO
    ! The value TRUE only makes sense if Global_IS_UE_ZE_Fixed_Rate is FALSE
    logical, parameter :: is_include_ue_ze = .false.
 
    real(srp) :: temp_loc
    logical :: is_mass_g, sda_off_debug
 
    if (present(temperature)) then
      temp_loc = temperature
    else
      temp_loc = global_temperature
    end if
 
    if (present(sda_off)) then
      sda_off_debug = sda_off
    else
      sda_off_debug = .false.
    end if
 
    if (present(is_mass)) then
      is_mass_g = is_mass
    else
      is_mass_g = .false.
    end if
 
    if (is_mass_g) then
      !! @f[ e = E_{SMR} + E_{AMR} + E_{UE+ZE} + E_{SDA} , @f]
      o2_out =            mg2g(this%smr(temp_loc)) +                  & ! SMR
                          mg2g(this%amr(temp_loc))                      ! AMR
 
      if (is_include_ue_ze)                                                   &
            o2_out=o2_out+this%ue_ze(temp_loc, is_oxygen=.true.)        ! UE+ZE
      if (.NOT. sda_off_debug) then
        o2_out = o2_out + mg2g(this%smr(temp_loc)) *                  & ! SDA
                                    this%sda_fact(global_sda_absorp_rate_max, &
                                                  global_sda_factor_max)
      end if
    else
      !! @f[ e = E_{SMR} + E_{AMR} + E_{UE+ZE} + E_{SDA} , @f]
      o2_out =            o2vol(mg2g(this%smr(temp_loc))) +           & ! SMR
                          o2vol(mg2g(this%amr(temp_loc)))               ! AMR
      if (is_include_ue_ze)                                                   &
            o2_out=o2_out+o2vol(this%ue_ze(temp_loc, is_oxygen=.true.)) ! UE+ZE
      if (.NOT. sda_off_debug) then
        o2_out = o2_out + o2vol(mg2g(this%smr(temp_loc))) *           & ! SDA
                                    this%sda_fact(global_sda_absorp_rate_max, &
                                                  global_sda_factor_max)
      end if
    end if
 
  end function fish_oxygen_consumption_step
 
  !> Calculate the relative SDA (specific dynamic action) factor for SMR unit
  !! of energy uptake. SDA depends on the absorption rate up th the
  !! maxmum cap value max_sda that is achieved at the absorption rate `max_ap`
  !!
  !! The relative SDA factor is based on a linear relationship defined by
  !! the maximum absorption rate `max_ap` and max SDA factor `max_sda`.
  !!
  !! @verbatim
  !! max_sda +       * 2,0 :: SDA = 2.0 x SMR
  !!         |     * .                      y = k x + b
  !!         |   *   .                          ?     =0
  !!         | *     .
  !!       0 +-------+
  !!         0      max_ap = 0.0007
  !! @endverbatim
  elemental function fish_energy_sda_factor(this, max_ap, max_sda, is_limit)  &
                                                                result(ap_fact)
    use base_utils, only : solve_line
    class(fish), intent(in) :: this
    !> Maximum (absolute) absorption rate when the maximum value of SDA
    !! (`max_sda`) is achieved
    real(srp), intent(in)   :: max_ap
    !> Maximum value of SDA achieved at the maximum absorption rate (`max_ap`)
    real(srp), intent(in)   :: max_sda
    logical, optional, intent(in) :: is_limit
    real(srp) :: ap_fact
    real(srp), parameter :: max_ultimate = 24.0_srp  ! maximum limit for SDA
 
    real(srp) :: k, b
 
    call solve_line(k, b, 0.0_srp, 0.0_srp, max_ap, max_sda)
 
    ap_fact = k * this%absorption_rate() + b
 
    if (present(is_limit)) then
      if (is_limit)  ap_fact = within(ap_fact, 0.0_srp, max_sda)
    else
      ap_fact = within(ap_fact, 0.0_srp, max_ultimate)
    end if
 
  end function fish_energy_sda_factor
 
  !> Grow the fish body mass: increment mass based on energy.
  elemental subroutine fish_grow_increment_energy(this)
    class(fish), intent(inout) :: this
 
    real(srp) :: updated_mass, last_ebal, prev_ebal
 
    last_ebal = missing
    prev_ebal = missing
 
    ! Note that the first time step is covered by the_fish::fish:skel:new()
    ! method
    if (global_time_step > 1) then
 
      last_ebal = last(this%history%energy_balance_curr)
      prev_ebal = this%history%energy_balance_curr(                           &
                                    size(this%history%energy_balance_curr)-1 )
 
      !> ### Implementation details ###
      !! The current body mass is incremented as
      !! @f$ M_{i+1} = M_{i} + \Delta M  @f$ by the mass equivalent
      !! @f$ \Delta M  @f$ of the energy increment for the time step
      !! calculated using the_fish::energy2mass().
      updated_mass = last(this%history%body_mass_current) +                   &
                                    energy2mass( last_ebal - prev_ebal )
 
      call add_to_history( this%history%body_mass_current, updated_mass )
 
    end if
 
  end subroutine fish_grow_increment_energy
 
  !> Calculate the increment of the total mass of food evacuated from the
  !! midgut after full absorption and evacuation delay (see
  !! commondata:global_maximum_duration_midgut_min) that occurs at
  !! one time step.
  !! @warning Note that this method should be called before
  !!     the_fish::food_item_eaten_mg::check_disappear() so that all food
  !!     itemsthat should be evacuated have not yet been evacuated.
  elemental function fish_midgut_total_mass_evacuated_step(this,max_time)     &
                                                            result(volume_out)
    class(midgut), intent(in) :: this
    real(srp) :: volume_out
    !> Maximum time processed in the midgut, by default all food items that
    !! have been in the process of absorption for more than this time will
    !! destructed from midgut and evacuated.
    !! Default value is based on commondata::global_maximum_duration_midgut_min.
    integer, optional, intent(in) :: max_time
 
    integer ::  i
    integer :: max_time_loc
 
    if (present(max_time)) then
      max_time_loc = max_time
    else
      max_time_loc = global_maximum_duration_midgut_min * minute
    end if
 
    volume_out = 0.0_srp
 
    associate( foods_midgut_not_stom =>                                      &
                  this%food_items_midgut(                                     &
                      this%n_food_items_stomach+1:this%n_food_items_midgut) )
 
      do i = 1, size(foods_midgut_not_stom)
        if ( foods_midgut_not_stom(i)%is_evacuate(max_time_loc) )  then
          volume_out = volume_out + foods_midgut_not_stom(i)%mass_0 -         &
                                    foods_midgut_not_stom(i)%absorped
        end if
      end do
 
    end associate
 
  end function fish_midgut_total_mass_evacuated_step
 
  !> Calculate the number of the total mass of food evacuated from the
  !! midgut after full absorption and evacuation delay (see
  !! commondata:global_maximum_duration_midgut_min) that occurs at
  !! one time step.
  !! @warning Note that this method should be called before
  !!     the_fish::food_item_eaten_mg::check_disappear() so that all food
  !!     itemsthat should be evacuated have not yet been evacuated.
  elemental function fish_midgut_n_evacuated_step(this, max_time) result(n_out)
    class(midgut), intent(in) :: this
    integer :: n_out
    !> Maximum time processed in the midgut, by default all food items that
    !! have been in the process of absorption for more than this time will
    !! destructed from midgut and evacuated.
    !! Default value is based on commondata::global_maximum_duration_midgut_min.
    integer, optional, intent(in) :: max_time
 
    integer :: max_time_loc
 
    if (present(max_time)) then
      max_time_loc = max_time
    else
      max_time_loc = global_maximum_duration_midgut_min * minute
    end if
 
    associate( foods_midgut_not_stom =>                                      &
                  this%food_items_midgut(                                     &
                      this%n_food_items_stomach+1:this%n_food_items_midgut) )
      n_out = count( foods_midgut_not_stom%is_evacuate(max_time_loc) )
    end associate
 
  end function fish_midgut_n_evacuated_step
 
  !> Update the history stack arrays for the fish, ::midgut class.
  elemental subroutine fish_midgut_history_update(this)
    class(midgut), intent(inout) :: this
 
      integer :: absorb_last_item
      real(srp) :: absorp_prev, absorp_total, absorp_diff, absorp_cumulate
 
      ! The history is initialized to commondata::missing values, so it is easy
      ! to determine what is the last valid item in the history stack array.
      absorb_last_item = this%absorp_last_i()
 
      ! If the whole history stack is composed of commondata::missing, this
      ! means it is the first time step, then all historical values are zero.
      ! Alternatively, if all history stack items are non-commondata::missing,
      ! the last value is the latest history.
      if ( absorb_last_item == 0 ) then
        absorp_prev = 0.0_srp
        absorp_total = 0.0_srp
        absorp_cumulate = 0.0_srp
      else
        absorp_prev = this%history%food_mass_absorb_cum(absorb_last_item)
        absorp_total =                                                        &
            sum( this%food_items_midgut(1:this%n_food_items_midgut)%absorped, &
                 this%food_items_midgut(1:this%n_food_items_midgut)%absorped  &
                                                                  /= missing )
        ! Note: It is important that absorp_diff subtracts current
        !       instantaneous absorption value, not cumulative value
        !       from the history stack array. Otherwise cumulative
        !       absorption is wrongly calculated.
        absorp_diff = absorp_total -                                          &
                this%history%food_mass_absorb_curr( absorb_last_item )
        if ( absorp_diff > 0.0_srp ) then
          absorp_cumulate = absorp_prev + absorp_diff
        else
          absorp_cumulate = absorp_prev
        end if
      end if
 
      call add_to_history( this%history%food_mass_absorb_curr, absorp_total )
      call add_to_history( this%history%food_mass_absorb_cum, absorp_cumulate )
 
  end subroutine fish_midgut_history_update
 
  !> Determine the index of the latest absorption item in the history
  !! stack array. If the whole array is commondata::missing, this means the
  !! history stack is empty, then the function returns zero. Note that in
  !! most cases the function should return commondata::history_size as the
  !! history stack is filled with previous values.
  elemental function midgut_absorp_history_last_past_idx(this) result(inum)
    class(midgut), intent(in) :: this
    integer :: inum
    integer :: absorb_last_item
 
      ! The history is initialized to commondata::missing values, so it is easy
      ! to determine what is the last valid item in the history stack array.
      absorb_last_item = count( this%history%food_mass_absorb_cum == missing )
 
      inum = history_size - absorb_last_item
 
  end function midgut_absorp_history_last_past_idx
 
  !> Calculate the total mass of food in the fish midgut.
  elemental function fish_midgut_total_mass_food(this, tolerance)           &
                                                            result(total_mass)
    class(midgut), intent(in) :: this
    !> Optional numerical tolerance, the smallest value of the food item mass
    !! that counts as non-zero. Default is set by commondata::min_vol_toler.
    real(srp), optional, intent(in) :: tolerance
    !> Total mass of food in the midgut.
    real(srp) :: total_mass
 
    ! Local copy of optional
    real(srp) :: toler_loc
 
    if (present(tolerance)) then
      toler_loc = tolerance
    else
      toler_loc = min_mass_toler
    end if
 
    !> The total value is calculated as an array sum to the total number of
    !! food items being kept in the stomach (`this%n_food_items_stomach`).
    associate( foods_midgut                                                  &
                       => this%food_items_midgut(1:this%n_food_items_midgut)  )
      total_mass =                                                            &
        sum( foods_midgut%mass,                                              &
                  foods_midgut%mass > toler_loc .and.                         &
                  foods_midgut%mass /= missing )
    end associate
 
  end function fish_midgut_total_mass_food
 
  !> Calculate the fish appetite factor based on the relative midgut fullness.
  !! Appetite factor is equivalent to the probability of consuming any food
  !! item.
  !! @note Note that the parameters of the appetite function are identical for
  !!       both the **stomach** and **midgut** appetite factors:
  !!        - stomach::appetite_stomach
  !!        - midgut::
  !!        .
  elemental function fish_appetite_factor_midgut(this, fcapacity)           &
                                                              result(appetite)
    class(midgut), intent(in) :: this
    real(srp), optional, intent(in) :: fcapacity
    real(srp) :: appetite
 
    real(srp) :: fcapacity_loc
 
    real(srp) :: mass_relative, correction_fac
 
    if (present(fcapacity)) then
      fcapacity_loc = fcapacity
    else
      fcapacity_loc = this%mass_midgut
    end if
 
    !> ### Implementation notes ###
    !! The appetite factor ranges between (0 and 1) and is based on the
    !! relative filling capacity. Here relative filling capacity @f$ c_r @f$
    !! is calculated as the ratio of the total mass of food in the midgut
    !! @f$ c_{max} @f$ to the maximum filling capacity @f$ c_t @f$
    !! @f[ c_r = \frac{c_t}{c_{max}} @f]
    mass_relative = this%mg_food_mass()/fcapacity_loc
 
    !> Correction factor is introduced to take account of the fact that
    !! the capacity of the midgut is bigger than stomach but basically the
    !! same amount of food is coming from stomach (in case of one meal).
    correction_fac = this%mass_midgut / this%mass_stomach
 
    !> Given the relative filling capacity @f$ c_r @f$, appetite factor
    !! that ranges from 0 to 1 is calculated as logistic function
    !! ::appetite_func().
    !!
    !! The @f$ a @f$ and @f$ r @f$ logistic parameters are defined by the
    !! commondata::global_appetite_logist_a and
    !! commondata::global_appetite_logist_r parameters.
    !!
    !! *gnuplot code:*
    !! @verbatim
    !!    set xrange [0:1]
    !!    plot 1 - 1/(1+500000*exp(-20*x))
    !! @endverbatim
    appetite = appetite_func( mass_relative * correction_fac,                 &
                                            global_appetite_logist_a,         &
                                            global_appetite_logist_r )
 
  end function fish_appetite_factor_midgut
 
  !> Calculate instantaneous absorption rate  based on the current absorption
  !! data and previous absorption saved in Output Arrays
  elemental function fish_midgut_absorption_rate(this) result (out_rate)
    class(midgut), intent(in) :: this
    real(srp) :: out_rate
 
    if( global_time_step == 1 .or. this%absorp_last_i()<history_size/2 ) then
      out_rate = 0.0_srp
      return
    end if
 
    ! It is assumed that the history stack already includes the current
    ! absorption value as the last value, so previous is taken as n-1.
    ! History stack array for absorption is updated in ::mg_step().
    out_rate = this%history%food_mass_absorb_cum(  this%absorp_last_i() ) -   &
                  this%history%food_mass_absorb_cum(  this%absorp_last_i()-1  )
 
  end function fish_midgut_absorption_rate
 
  !> Appetite factor multiplier based on the total absorption rate.
  elemental function fish_midgut_absorption_factor(this) result (out_factor)
    class(midgut), intent(in) :: this
    real(srp) :: out_factor
 
    real(srp) :: absorp_rate
 
    real(srp), dimension(*), parameter :: grid_test_x = [0.0_srp, 0.00003_srp],&
                                          grid_test_y = [0.0_srp, 1.0_srp]
 
 
    absorp_rate = this%absorption_rate()
 
    if ( absorp_rate >= grid_test_x(1) .and.                                  &
         absorp_rate <= grid_test_x(size(grid_test_x)) ) then
      out_factor = cspline_scalar( grid_test_x, grid_test_y, absorp_rate )
    else if ( absorp_rate < grid_test_x(1) ) then
      out_factor = grid_test_y(1)
    else
      out_factor = grid_test_y(size(grid_test_y))
    end if
 
  end function fish_midgut_absorption_factor
 
  !> Shift food items in the midgut after a new food item is ingested and
  !! therefore food items are shifted in stomach adding the new item.
  elemental subroutine fish_midgut_shift_items_after_new_ingest( this,        &
                                                                 food_item_pass)
    class(midgut), intent(inout) :: this
    !> Optional food item that is ingested by the fish (environ::food_item) or
    !! any extension class. This food mass is used only to set the reference
    !! start mass of the food item.
    class(food_item), optional, intent(in) :: food_item_pass
 
    ! Local copy of the start mass of the food item
    real(srp) :: food_item_pass_start_mass_loc
 
    if (present(food_item_pass)) then
      select type ( food_item_pass )
        type is (food_item)
          food_item_pass_start_mass_loc = food_item_pass%mass
        type is (food_item_eaten)
          food_item_pass_start_mass_loc = food_item_pass%mass
        type is (food_item_eaten_mg)
          food_item_pass_start_mass_loc = food_item_pass%mass_0
        class default
          food_item_pass_start_mass_loc = food_item_pass%mass
      end select
    else
      food_item_pass_start_mass_loc = global_food_item_mass
    end if
 
    !> - The food item is then added to the stack array of the fish
    !!   stomach. It takes the first place shifting all other food items in
    !!   sequence.
    this%food_items_midgut(2:this%n_food_items_midgut+1) =                    &
                            this%food_items_midgut(1:this%n_food_items_midgut)
 
    !> - The first food item in the stack is then obtained from the input
    !!   parameter **food_item_pass**. But note that at this stage, its actual
    !!   environ::food_item::mass is set to zero, the starting mass
    !!   food_item_eaten::mass_0 is obtained from the input item's actual
    !!   mass. The food_item_eaten::time_in_process is set to zero. All
    !!   this is done via the food_item_eaten::zero() method.
    !! .
    call this%food_items_midgut(1)%zero(                                      &
                                    mass=0.0_srp,                           &
                                    start_mass=food_item_pass_start_mass_loc )
 
    this%n_food_items_midgut = this%n_food_items_midgut + 1
 
  end subroutine fish_midgut_shift_items_after_new_ingest
 
  !> Computational backend for the appetite function.
  !! Appetite is calculated as
  !! @f[ \alpha = 1 - \frac{1}{1+a \cdot e^{-r \cdot c_r}} @f]
  !! where @f$ c_r @f$ is the relative filling capacity (range 0 to 1) of
  !! the fish stomach, @f$ a @f$ and @f$ r @f$ are logistic function
  !! parameters.
  !!
  !! *gnuplot code:*
  !! @verbatim
  !!   set xrange [0:1]
  !!   plot 1-1/(1+500000*exp(-20*x))
  !! @endverbatim
  !!
  !! @note Note that m_ifcmd::save_appetite_function_data() saves the appetite
  !!       factor data for plotting, testing and analysis.
  elemental function appetite_func( mass_relative, A, R ) result(value_out)
    !> The mass of food in the stomach relative to the filling capacity.
    real(srp), intent(in) :: mass_relative
    !> Parameters of the logistic function @f$ a @f$.
    real(srp), intent(in) :: a
    !> Parameters of the logistic function @f$ r @f$.
    real(srp), intent(in) :: r
    !> Returns the appetite factor value that ranges from 0 to 1.
    real(srp) ::  value_out
 
    value_out =  1.0_srp - logistic( mass_relative, a, 1.0_srp, r )
 
  end function appetite_func
 
  !> Calculate the overall appetite of the fish, which defines the
  !! probability to ingest any food item.
  elemental function fish_appetite_probability_ingestion(this, time_step)     &
                                                              result(appetite)
    class(fish), intent(in) :: this
    !> Time step, optional parameter
    integer(LONG), optional, intent(in) :: time_step
    real(srp) :: appetite
    integer(LONG) :: step_loc
 
    !> Fish appetite at night is near-zero
    if ( .not. is_day() .and. global_appetite_fish_night >= 0.0_srp ) then
      appetite = global_appetite_fish_night
      return
    end if
 
    if(present(time_step)) then
      step_loc = time_step
    else
      step_loc = global_time_step
    end if
 
    !> The overall appetite is a function of individual appetite factors:
    !!  - appetite factor for stomach the_fish::stomach::appetite_stomach()
    !!  - appetite factor for midgut the_fish::midgut::appetite_midgut()
    !!  - appetite factor for energy balance the_fish::fish::appetite_energy()
    !!  .
    !!
    !!
    !! *Version 1* multiplication rule: @f$ A = A_{stomach} A_{midgut} @f$
    !! probability of a combination of events events.
    !appetite = this%appetite_stomach() * this%appetite_midgut()
    !----------------------
    !> *Version 2* maximum/competition rule:
    !! @f[
    !!  \left\{\begin{matrix} A = A_{stomach}, A_{stomach} < 0.2 \\
    !!      max( A_{stomach}, A_{midgut},A_{energy}  )
    !!  \end{matrix}\right.
    !! @f]
    if ( this%appetite_stomach() < global_appetite_stomach_threshold ) then
      appetite = this%appetite_stomach()
    else
      appetite = maxval( [ this%appetite_stomach(), this%appetite_midgut(),   &
                           this%appetite_energy() ] )
    end if
 
    ! Suppress appetite by stress (if enabled compile-time)
    if (is_stress_enable)                                                     &
                      appetite = appetite - appetite*stress_suppress(step_loc)
 
  end function fish_appetite_probability_ingestion
 
  !> Calculate suppression factor for appetite caused by stress intervention
  !! at time `time`: generic interface is ::stress_factor().
  !! This is the scalar version, vector -based version is
  !!::stress_suppress_factor_v().
  !!
  !! @verbatim
  !!   +**
  !! F |  ***
  !!   |    **
  !!   +-------*--
  !!           ^ time_step_incr
  !!          A-F
  !! @endverbatim
  pure function stress_suppress_factor_s(time_step_incr, is_force) result (fout)
    use base_utils, only : cspline_scalar
    !> Time step increment from the moment stress intervention has occurred
    integer, intent(in) :: time_step_incr
    !> optional flag to force calculation of the function even if stress
    !! is disabled
    logical, optional, intent(in) :: is_force
    real(srp) :: fout
 
    logical :: is_force_loc
 
    if (present(is_force)) then
      is_force_loc = is_force
    else
      is_force_loc = .false.
    end if
 
    if ( is_stress() .or. is_force_loc ) then
      fout = cspline_scalar( hour * global_stress_factor_hour,                &
                             global_stress_fact_suppress,                     &
                             real(time_step_incr, srp) )
    else
      fout = 0.0_srp
    end if
 
  end function stress_suppress_factor_s
 
  !> Vector-based version of the stress-based suppression function
  !! ::stress_factor(). There is a scalar version of the same function
  !! ::stress_suppress_factor_s.
  pure function stress_suppress_factor_v(time_step_incr, is_force) result (fout)
    use base_utils, only : cspline_vector
    !> Time step increment from the moment stress intervention has occurred
    integer, dimension(:), intent(in) :: time_step_incr
    !> optional flag to force calculation of the function even if stress
    !! is disabled
    logical, optional, intent(in) :: is_force
    real(srp), allocatable, dimension(:) :: fout
 
    logical :: is_force_loc
 
    if (present(is_force)) then
      is_force_loc = is_force
    else
      is_force_loc = .false.
    end if
 
    if (.NOT. allocated(fout)) allocate(fout(size(time_step_incr)))
 
    if ( is_stress() .or. is_force_loc ) then
      fout = cspline_vector( hour * global_stress_factor_hour,                &
                             global_stress_fact_suppress,                     &
                             real(time_step_incr, srp) )
    else
      fout = 0.0_srp
    end if
 
  end function stress_suppress_factor_v
 
  !> Calculate real stress suppression factor for a specific time step
  !! given the configured stress intervention pattern. This is a high-level
  !! interface to ::stress_factor().
  !! @verbatim
  !!             1**          2**
  !! F           |  ***       |  ***
  !!             |    **      |     **
  !!   --#-------+-------*----+-----#--*---
  !!     ^                          ^
  !!    F=0                   +<---># time_over
  !!                                ^
  !!                               A-F
  !! @endverbatim
  pure function stress_suppress(time_step) result (fout)
    !> Current time step
    integer(LONG), intent(in), optional :: time_step
    real(srp) :: fout
 
    integer :: i, last_stress, time_over
    integer(LONG) :: time_loc
 
    fout = 0.0_srp
 
    if (.not. is_stress() ) then
      return
    end if
 
    if (present(time_step)) then
        time_loc = time_step
    else
        time_loc = global_time_step
    end if
 
    last_stress = unknown
 
    ! get current stress intervention number
    do i = 1, size(global_stress_intervention_time)
      if ( time_loc >=  global_stress_intervention_time(i) ) then
        last_stress = i
      end if
    end do
 
    ! If the current time step is before the first stress intervention,
    ! zero suppression of appetite
    if (last_stress == unknown) then
      fout = 0.0_srp
      return
    end if
 
    ! time difference over the last stress intervention
    time_over = time_loc - global_stress_intervention_time(last_stress)
 
    if (time_over <= last(global_stress_factor_hour)*hour) then
      fout = stress_factor(time_over)
    else
      fout = 0.0_srp
    end if
 
    if (fout <= small ) fout = 0.0_srp
 
  end function stress_suppress
 
  !> Get the current body mass of the fish. Note that this is the shortcut to
  !! obtain the last value from the body mass history array
  !! the_fish::fish_skel::history::body_mass_current.
  elemental function fish_body_mass_current_get(this) result (current_mass_get)
    class(fish_skel), intent(in) :: this
    real(srp) :: current_mass_get
 
    current_mass_get = last( this%history%body_mass_current )
 
  end function fish_body_mass_current_get
 
  !> Calculate appetite factor that depends on the energy balance
  elemental function fish_appetite_factor_energy_balance(this)                &
                                                      result(appetite_factor)
    class(fish), intent(in) :: this
    real(srp) :: appetite_factor
 
    integer :: hist_len
    real(hrp) :: ebal_older, ebal_newer ! average past energy balance components
    real(hrp) :: ebal_diff
    integer, parameter :: block_size = 4
    real(hrp) :: delta_e, appetite_hg
 
    hist_len = size(this%history%energy_balance_curr) ! HISTORY_SIZE
 
    ! If simulation is just started and the history is empty, exit with zero
    if (global_time_step < block_size*2 ) then
      appetite_factor = 0.0_srp
      return
    end if
 
    associate( e_hist => this%history%energy_balance_curr )
 
      if (block_size < 4 .or. hist_len < block_size*2) then
        ebal_older = e_hist(hist_len-1)
        ebal_newer = e_hist(hist_len)
      else
        !array older ::
        ebal_older = average(e_hist(hist_len-2*block_size+1:hist_len-block_size), &
                                    undef_ret_null=.true.)
        !array newer ::
        ebal_newer = average(e_hist(hist_len-block_size+1:hist_len), &
                             undef_ret_null=.true.)
      end if
 
      ! Absolute difference in the past energy balance
      ebal_diff = ebal_newer - ebal_older
 
    end associate
 
    !> ### Implementation notes ###
    !!
    !! The energy component of the appetite is defined by the following
    !! formula:
    !! @f[ \alpha_E = \frac{1}{1 + e^{-r_E \cdot (- \Delta_E - b_E )}}  @f]
    !! @verbatim
    !!   gnuplot> set xrange [-1:1]
    !!   gnuplot> plot 1 / (1+2.718**(-40*(-x-0.2)) )
    !! @endverbatim
    associate( r_e => real(global_energy_appetite_rate,hrp),                 &
                b_e => real(global_energy_appetite_shift,hrp) )
 
      delta_e = ebal_diff/this%smr(global_temperature)
      ! Guard against extreme values of the basic difference:
      ! Note that the validity range values are chosen because the f(x)
      ! lie within [0,1] asympotically for all x<-1 and x>1, seegnuplot
      ! code
      if (delta_e < 10.0_hrp * tiny(1.0_srp)) delta_e = 0.0_hrp
      if (delta_e > 10.0_hrp) delta_e = 10.0_hrp
 
      appetite_hg = 1.0_hrp / (1.0_hrp + exp(-r_e * (-delta_e-b_e)))
      appetite_factor = real(appetite_hg, srp)
 
    end associate
 
  end function fish_appetite_factor_energy_balance
 
  !> This procedure determines if the fish decides to consume the food item
  !! and if yes, ingests it using the ::do_ingest() method.
  !! @important This version accepts **raw food item**.   Generic
  !!            interface linking both raw and pointer versions is ambiguous.
  !! @note This cannot be pure due to random operation
  impure subroutine fish_decide_eat_food_item(this, food_item_in, have_eaten )
    class(fish), intent(inout) :: this
    !> The food item that is available for decision and consumption.
    class(food_item), intent(in) :: food_item_in
    !> Optional indicator showing if the fish have actually eaten the food
    !! item (returns TRUE).
    logical, optional, intent(out) :: have_eaten
 
    if (present(have_eaten)) have_eaten = .false.
 
    if ( rand_uni() < this%appetite() ) then
      call this%do_ingest( food_item_in )
      if (present(have_eaten)) have_eaten = .true.
    end if
 
  end subroutine fish_decide_eat_food_item
 
  !> This procedure determines if the fish decides to consume the food item
  !! and if yes, ingests it using the ::do_ingest() method.
  !! @important This version accepts **pointer** to food item. Generic
  !!            interface linking both raw and pointer versions is ambiguous.
  !! @note This cannot be pure due to random operation
  impure subroutine fish_decide_eat_food_item_ptr(this,food_item_in,have_eaten)
    class(fish), intent(inout) :: this
    !> The food item that is available for decision and consumption.
    class(food_item), pointer, intent(in) :: food_item_in
    !> Optional indicator showing if the fish have actually eaten the food
    !! item (returns TRUE).
    logical, optional, intent(out) :: have_eaten
 
    if (present(have_eaten)) have_eaten = .false.
 
    if ( .NOT. associated(food_item_in) ) return
 
    if ( rand_uni() < this%appetite() ) then
      call this%do_ingest( food_item_in )
      if (present(have_eaten)) have_eaten = .true.
    end if
 
  end subroutine fish_decide_eat_food_item_ptr
 
  !> Calculate the baseline activity of the fish.
  !!
  !> Global baseline locomotor activity pattern  (swimming speed) for each
  !! time step consists of two components:
  !!  - baseline that differs at day ::global_baseline_activity_day and
  !!    night ::global_baseline_activity_night
  !!  - feeding activity increment defined by feeding activity factor
  !!    ::global_feeding_activity_factor
  !!  .
  !!
  !! This is illustrated by the following diagram:
  !!
  !! @verbatim
  !!                       +----+   feeding activity multiplier
  !!                       |    |
  !!                    +--+ - -+----------------+
  !!                    |                        |                     baseline
  !! -------------------+                        +-------------------> activity
  !! -night-------------+-day--------------------+-night------------->
  !! @endverbatim
  elemental function fish_activity_baseline_component(this, time_step)        &
                                                              result(act_out)
    class(fish), intent(in) :: this
    !> Optional time step of the model, default value is
    !! commondata::global_time_step
    integer(LONG), optional, intent(in) :: time_step
 
    real(srp) :: act_out
 
    integer(LONG) :: time_step_loc, i
 
    if (present(time_step)) then
      time_step_loc = time_step
    else
      time_step_loc = global_time_step
    end if
 
    if ( is_day(time_step_loc) ) then
      act_out = global_baseline_activity_day
    else
      act_out = global_baseline_activity_night
    end if
 
    if (is_stress_enable)                                                     &
            act_out = act_out -                                               &
                      act_out * global_stress_activity_decr *                 &
                                                stress_suppress(time_step_loc)
 
  end function fish_activity_baseline_component
 
  !> Calculate activity factor linked to appetite.
  !! Appetite factor has a linear relationship with the fish appetite
  !! !> @f$ l_{A} = k_{a} A + 1 @f$, here @f$ k_{a} + 1 @f$ is the maximum
  !! value of the appetite activity factor when appetite value is
  !! maximum @f$ A = 1 @f$. Note that the appetite activity factor works
  !! as a multiplier for the baseline actvity
  !! the_fish::fish::activity_baseline(), see
  !! ::fish_amr_locomotion_energy_cost().
  !!
  !! Gnuplot commands for plotting:
  !!
  !! @verbatim
  !!   set xrange [0:1]
  !!   set yrange [0:5]
  !!   plot 0.5*x + 1
  !! @endverbatim
  elemental function fish_activity_appetite_factor(this) result (act_out)
    class(fish), intent(in) :: this
    real(srp) :: act_out
 
    act_out = global_appetite_activity_factor * this%appetite() + 1.0_srp
 
  end function fish_activity_appetite_factor
 
  !> Calculate the adjustment factor for the transformation of the stomach
  !! transport pattern Time scale commondata::global_transport_pattern_t
  !! for a different temperature or fish size.
  elemental function stomach_adjust(target_emptying_time) result (adj_out)
    real(srp) :: adj_out
    !> Optional maximum stomach emptying time for a different target
    !! temperature, default value is calculated by ::stomach_emptying_time().
    integer, optional, intent(in) :: target_emptying_time
    real(srp), parameter :: min_diff = 0.00001
    integer :: i
 
    integer :: target_empty_loc
 
    if (present(target_emptying_time)) then
      target_empty_loc = target_emptying_time
    else
      target_empty_loc =                                                      &
                    stomach_emptying_time(global_body_mass, global_temperature)
    end if
 
    do i=1, global_transport_dim
      if ( global_transport_pattern_r(i) < min_diff ) then
        adj_out =                                                             &
          real(target_empty_loc,srp)/real(Global_Transport_Pattern_T(i),srp)
        exit
      end if
    end do
  end function stomach_adjust
 
  !> Calculate the stomach emptying time for a fish with specific body mass
  !! at a specific temperature. Calculation is based on spline interpolation
  !! of experimental data.
  elemental function stomach_emptying_time(fish_mass, temperature)            &
                                                        result(empty_time_raw)
    !> Optional fish mass, default value is commondata::global_body_mass
    real(srp), optional, intent(in) :: fish_mass
    !> Optional temperature, default value is commondata::global_temperature
    real(srp), optional, intent(in) :: temperature
    integer :: empty_time_raw
 
    integer :: temp_size, i
    real(srp) :: mass_loc, temperature_loc
 
    real(srp),                                                                &
      dimension(size(Global_Stomach_Emptying_Pattern%grid_temperature)) :: x_t
 
    if (present(fish_mass)) then
      mass_loc = fish_mass
    else
      mass_loc = global_body_mass
    end if
 
    if (present(temperature)) then
      temperature_loc = temperature
    else
      temperature_loc = global_temperature
    end if
 
    temp_size = size(global_stomach_emptying_pattern%grid_temperature)
 
    do i = 1, temp_size
      x_t(i) = cspline_scalar(                                               &
                global_stomach_emptying_pattern%grid_mass,                    &
                real(global_stomach_emptying_pattern%emptying_time(:,i),srp), &
                mass_loc )
    end do
 
    empty_time_raw = nint( cspline_scalar(                                    &
                            global_stomach_emptying_pattern%grid_temperature, &
                            x_t, temperature_loc  ) )
 
  end function stomach_emptying_time
 
  !> Object-oriented frontend for ::stomach_emptying_time(): Calculate the
  !! stomach emptying time for a fish with specific body mass at a specific
  !! temperature. Calculation is based on spline interpolation of experimental
  !! data.
  elemental function fish_stomach_emptying_time(this, temperature)            &
                                                        result(empty_time_raw)
    class(stomach), intent(in) :: this
    integer :: empty_time_raw
    !> Optional temperature, if absent use commondata::global_temperature
    !! parameter.
    real(srp), optional, intent(in) :: temperature
 
    real(srp) :: temp_loc
 
    if (present(temperature)) then
      temp_loc = temperature
    else
      temp_loc = global_temperature
    end if
 
    empty_time_raw = stomach_emptying_time(this%mass(), temperature)
 
  end function fish_stomach_emptying_time
 
  !> Get the stomach emptying time base matrix from a CSV input data file.
  !!
  !! File structure:
  !! @verbatim
  !!          |    50.00   100.00   300.00   500.00
  !! ---------+------------------------------------
  !!    5.00  |    24.00    35.00    75.00   100.00
  !!   10.00  |    15.00    20.00    50.00    75.00
  !!   15.00  |     9.00    15.00    40.00    50.00
  !!   20.00  |     5.00    10.00    30.00    35.00
  !!   22.00  |     4.00     8.00    29.00    33.00
  !! @endverbatim
  impure function get_stomach_emptying_matrix_csv(csv_file) result (pattern_out)
    use csv_io , only: csv_matrix_read
    !> Input file name
    character(len=*), intent(in) :: csv_file
    type(stomach_emptying_pattern) :: pattern_out
 
    real(srp),  allocatable, dimension(:,:) :: raw_data_csv
    logical :: is_success
 
    raw_data_csv = transpose(                                                 &
                  csv_matrix_read( csv_file, csv_file_status=is_success,      &
                                    missing_code=missing ) )
 
    ! First row is mass grid
    pattern_out%grid_mass = raw_data_csv(2:,1)
    ! First column is temperature grid
    pattern_out%grid_temperature = raw_data_csv(1,2:)
    ! Remaining data is stomach emptying time
    pattern_out%emptying_time = nint(raw_data_csv(2:,2:)) * hour
 
  end function get_stomach_emptying_matrix_csv
 
 
end module the_fish