#include "fabm_driver.h"

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                        !
!                                                                        !
!     ERGOM model CGT version ported to FABM by Hagen Radtke             !
!     (Leibniz Institute for Baltic Sea Research Warnemünde - IOW)._rk      !
!     Code automatically generated by Code Generation Tool Version <cgt_version> !
!     See www.ergom.net                                                  !
!                                                                        !
!                                                                        !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!-----------------------------------------------------------------------
!BOP
!
! !MODULE: fabm\_iow\_ergom --- Biogeochemical model ERGOM
!
! !INTERFACE:
   module fabm_iow_ergom
!
! !DESCRIPTION:
! Model: ERGOM IOW Baltic Version 1.0_rk \\
! Ecological ReGional Ocean Model \\
! Author(s): Dr. Thomas Neumann, Dr. Martin Schmidt, Dr. Hagen Radtke \\
! Contact: thomas.neumann@io-warnemuende.de \\
! Code generated: 2014-Jan-09 13:14
!
! !USES:
   use fabm_types

   implicit none

   private
!
! !PUBLIC MEMBER FUNCTIONS:
   public type_iow_ergom
!
! !REVISION HISTORY:
!  Author(s):
!
! !PUBLIC_DERIVED_TYPES:
   type,extends(type_base_model) :: type_iow_ergom
! Variable identifiers
      type (type_state_variable_id)        :: id_t_n2
      type (type_state_variable_id)        :: id_t_o2
      type (type_state_variable_id)        :: id_t_nh4
      type (type_state_variable_id)        :: id_t_no3
      type (type_state_variable_id)        :: id_t_po4
      type (type_state_variable_id)        :: id_t_spp
      type (type_state_variable_id)        :: id_t_zoo
      type (type_state_variable_id)        :: id_t_h2s
      type (type_state_variable_id)        :: id_t_sul
      type (type_state_variable_id)        :: id_t_ipw
      type (type_state_variable_id)        :: id_t_lpp
      type (type_state_variable_id)        :: id_t_cya
      type (type_state_variable_id)        :: id_t_det
      type (type_state_variable_id)        :: id_t_dic
      type (type_bottom_state_variable_id) :: id_t_sed
      type (type_bottom_state_variable_id) :: id_t_ips
      type (type_dependency_id)            :: id_par,id_temp,id_salt
      type (type_horizontal_dependency_id) :: id_I_0,id_wind,id_taub
      type (type_horizontal_dependency_id) :: id_depth, id_lat, id_lon
!      type (type_diagnostic_variable_id)   :: id_dPAR,id_GPP,id_NCP,id_PPR,id_NPR

! Model parameters
      real(rk) :: critical_stress
      real(rk) :: cya0
      real(rk) :: din_min_lpp
      real(rk) :: din_min_spp
      real(rk) :: dip_min_cya
      real(rk) :: epsilon
      real(rk) :: food_min_zoo
      real(rk) :: gamma0
      real(rk) :: gamma1
      real(rk) :: h2s_min_po4_liber
      real(rk) :: ips_threshold
      real(rk) :: ips_cl
      real(rk) :: k_h2s_no3
      real(rk) :: k_h2s_o2
      real(rk) :: k_sul_no3
      real(rk) :: k_sul_o2
      real(rk) :: light_opt_cya
      real(rk) :: light_opt_lpp
      real(rk) :: light_opt_spp
      real(rk) :: lpp0
      real(rk) :: no3_min_sed_denit
      real(rk) :: no3_min_det_denit
      real(rk) :: o2_min_det_resp
      real(rk) :: o2_min_nit
      real(rk) :: o2_min_po4_retent
      real(rk) :: o2_min_sed_resp
      real(rk) :: q10_det_rec
      real(rk) :: q10_h2s
      real(rk) :: q10_nit
      real(rk) :: q10_sed_rec
      real(rk) :: r_biores
      real(rk) :: r_cya_assim
      real(rk) :: r_det_rec
      real(rk) :: r_ips_burial
      real(rk) :: r_ips_ero
      real(rk) :: r_ips_liber
      real(rk) :: r_lpp_assim
      real(rk) :: r_nh4_nitrif
      real(rk) :: r_pp_mort
      real(rk) :: r_pp_resp
      real(rk) :: r_sed_ero
      real(rk) :: r_sed_rec
      real(rk) :: r_spp_assim
      real(rk) :: r_zoo_graz
      real(rk) :: r_zoo_mort
      real(rk) :: r_zoo_resp
      real(rk) :: rfr_c
      real(rk) :: rfr_h
      real(rk) :: rfr_o
      real(rk) :: rfr_p
      real(rk) :: sali_max_cya
      real(rk) :: sali_min_cya
      real(rk) :: sed_max
      real(rk) :: spp0
      real(rk) :: temp_min_cya
      real(rk) :: temp_min_spp
      real(rk) :: temp_opt_zoo
      real(rk) :: w_cya
      real(rk) :: w_det
      real(rk) :: w_det_sedi
      real(rk) :: w_ipw
      real(rk) :: w_ipw_sedi
      real(rk) :: w_lpp
      real(rk) :: w_n2_stf
      real(rk) :: w_o2_stf
      real(rk) :: zoo0
      real(rk) :: zoo_cl
      real(rk) :: stf_no3
      real(rk) :: stf_nh4
      real(rk) :: stf_po4

! model time step may is stored here preliminarily
      real(rk) :: cgt_timestep

      contains

!     Model procedures
      procedure :: initialize
      procedure :: do
      procedure :: do_bottom
      procedure :: get_light_extinction
      procedure :: do_surface

  end type
!EOP
!-----------------------------------------------------------------------

   CONTAINS

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the ergom model
!
! !INTERFACE:
pure function theta(x)
!
! !DESCRIPTION:
!   The heavyside step function
!
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   real(rk), intent(in) :: x
   real(rk)             :: theta
!
!EOP
!-----------------------------------------------------------------------
!BOC
     if (x .gt. 0.0_rk) then
       theta = 1.0_rk
     else
       theta = 0.0_rk
     end if
   END function theta
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the ergom model
!
! !INTERFACE:
subroutine initialize(self,configunit)
!
! !DESCRIPTION:
!   Here, the ergom namelist is read and the variables exported by the model are registered with FABM
!
! !USES
!   List any modules used by this routine.
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   class (type_iow_ergom), intent(inout), target :: self
   integer,                intent(in)            :: configunit
!
! !LOCAL VARIABLES:
   real(rk)           :: critical_stress=0.016_rk    ! critical shear stress for sediment erosion [N/m2]
   real(rk)           :: cya0           =4.5E-9_rk   ! seed concentration for diazotroph cyanobacteria [mol/kg]
   real(rk)           :: din_min_lpp    =1.125E-6_rk ! DIN half saturation constant for large-cell phytoplankton growth [mol/kg]
   real(rk)           :: din_min_spp    =4.5E-7_rk   ! DIN half saturation constant for small-cell phytoplankton growth [mol/kg]
   real(rk)           :: dip_min_cya    =1.125E-7_rk ! DIP half saturation constant for diazotroph cyanobacteria growth [mol/kg]
   real(rk)           :: epsilon        =4.5E-17_rk  ! no division by 0
   real(rk)           :: food_min_zoo   =4.108E-6_rk ! Ivlev phytoplankton concentration for zooplankton grazing [mol/kg]
   real(rk)           :: gamma0         =0.18_rk     ! light attenuation parameter (opacity of clear water) [1/m]
   real(rk)           :: gamma1         =30.0_rk     ! light attenuation parameter (opacity of POM containing chlorophyll) [m**2/mol]
   real(rk)           :: h2s_min_po4_liber=1.0E-6_rk   ! minimum h2s concentration for liberation of iron phosphate from the sediment [mol/kg]
   real(rk)           :: ips_threshold  =0.1_rk      ! threshold for increased PO4 burial [mol/m**2]
   real(rk)           :: ips_cl         =0.02025_rk  ! iron phosphate in sediment closure parameter [mol/m2]
   real(rk)           :: k_h2s_no3      =800000.0_rk ! reaction constant h2s oxidation with no3 [kg/mol/d]
   real(rk)           :: k_h2s_o2       =800000.0_rk ! reaction constant h2s oxidation with o2 [kg/mol/d]
   real(rk)           :: k_sul_no3      =20000.0_rk  ! reaction constant sulfur oxidation with no3 [kg/mol/d]
   real(rk)           :: k_sul_o2       =20000.0_rk  ! reaction constant sulfur oxidation with o2 [kg/mol/d]
   real(rk)           :: light_opt_cya  =50.0_rk     ! optimal light for diazotroph cyanobacteria growth [W/m²]
   real(rk)           :: light_opt_lpp  =35.0_rk     ! optimal light for large-cell phytoplankton growth [W/m²]
   real(rk)           :: light_opt_spp  =50.0_rk     ! optimal light for small-cell phytoplankton growth [W/m²]
   real(rk)           :: lpp0           =4.5E-9_rk   ! seed concentration for large-cell phytoplankton [mol/kg]
   real(rk)           :: no3_min_sed_denit=1.423E-7_rk ! nitrate half-saturation concentration for denitrification in the water column [mol/kg]
   real(rk)           :: no3_min_det_denit=1.0E-9_rk   ! minimum no3 concentration for recycling of detritus using nitrate (denitrification)
   real(rk)           :: o2_min_det_resp=1.0E-6_rk   ! oxygen half-saturation constant for detritus recycling [mol/kg]
   real(rk)           :: o2_min_nit     =3.75E-6_rk  ! oxygen half-saturation constant for nitrification [mol/kg]
   real(rk)           :: o2_min_po4_retent=3.75E-5_rk  ! oxygen half-saturation concentration for retension of phosphate during sediment denitrification [mol/kg]
   real(rk)           :: o2_min_sed_resp=6.4952E-5_rk ! oxygen half-saturation constant for recycling of sediment detritus using oxygen [mol/kg]
   real(rk)           :: q10_det_rec    =0.15_rk     ! temperature dependence of recycling according to van t Hoffs rule [1/K], equivalent to Q10=4.5_rk
   real(rk)           :: q10_h2s        =0.0693_rk   ! temperature dependence of oxidation of h2s and sul according to van t Hoffs rule [1/K], equivalent to Q10=2.0_rk
   real(rk)           :: q10_nit        =0.11_rk     ! temperature dependence of nitrification according to van t Hoffs rule [1/K], equivalent to Q10=3.0_rk
   real(rk)           :: q10_sed_rec    =0.175_rk    ! temperature dependence of detritus recycling in the sediment according to van t Hoffs rule [1/K], equivalent to Q10=5.8_rk
   real(rk)           :: r_biores       =0.03_rk     ! bio-resuspension rate [1/d]
   real(rk)           :: r_cya_assim    =0.75_rk     ! maximum rate for nutrient uptake of diazotroph cyanobacteria [1/d]
   real(rk)           :: r_det_rec      =0.003_rk    ! recycling rate (detritus to ammonium) at 0°C [1/d]
   real(rk)           :: r_ips_burial   =0.007_rk    ! final burial rate for PO4 [1/d]
   real(rk)           :: r_ips_ero      =6.0_rk      ! erosion rate for iron PO4 [1/d]
   real(rk)           :: r_ips_liber    =0.1_rk      ! PO4 liberation rate under anoxic conditions [1/d]
   real(rk)           :: r_lpp_assim    =1.3_rk      ! maximum rate for nutrient uptake of large-cell phytoplankton [1/d]
   real(rk)           :: r_nh4_nitrif   =0.1_rk      ! nitrification rate at 0°C [1/d]
   real(rk)           :: r_pp_mort      =0.02_rk     ! mortality rate of phytoplankton [1/d]
   real(rk)           :: r_pp_resp      =0.01_rk     ! respiration rate of phytoplankton to ammonium [1/d]
   real(rk)           :: r_sed_ero      =6.0_rk      ! maximum sediment detritus erosion rate [1/d]
   real(rk)           :: r_sed_rec      =0.002_rk    ! maximum recycling rate of sediment to ammonium [1/d]
   real(rk)           :: r_spp_assim    =0.4_rk      ! maximum rate for nutrient uptake of small-cell phytoplankton [1/d]
   real(rk)           :: r_zoo_graz     =0.5_rk      ! maximum zooplankton grazing rate [1/d]
   real(rk)           :: r_zoo_mort     =0.03_rk     ! mortality rate of zooplankton [1/d]
   real(rk)           :: r_zoo_resp     =0.01_rk     ! respiration rate of zooplankton [1/d]
   real(rk)           :: rfr_c          =6.625_rk    ! redfield ratio C/N
   real(rk)           :: rfr_h          =16.4375_rk  ! redfield ratio H/N
   real(rk)           :: rfr_o          =6.875_rk    ! redfield ratio O/N
   real(rk)           :: rfr_p          =0.0625_rk   ! redfield ratio P/N
   real(rk)           :: sali_max_cya   =10.0_rk     ! upper salinity limit - diazotroph cyanobacteria [psu]
   real(rk)           :: sali_min_cya   =1.0_rk      ! lower salinity limit - diazotroph cyanobacteria [psu]
   real(rk)           :: sed_max        =4.5_rk      ! maximum sediment detritus concentration that feels erosion [mol/m**2]
   real(rk)           :: spp0           =4.5E-9_rk   ! seed concentration for small-cell phytoplankton [mol/kg]
   real(rk)           :: temp_min_cya   =13.5_rk     ! lower temperature limit - diazotroph cyanobacteria [°C]
   real(rk)           :: temp_min_spp   =10.0_rk     ! lower temperature limit - small-cell phytoplankton [°C]
   real(rk)           :: temp_opt_zoo   =20.0_rk     ! optimal temperature for zooplankton grazing [°C]
   real(rk)           :: w_cya          =0.1_rk      ! vertical speed of diazotroph cyanobacteria [m/d]
   real(rk)           :: w_det          =-4.5_rk     ! vertical speed of detritus [m/d]
   real(rk)           :: w_det_sedi     =-2.25_rk    ! sedimentation speed (negative for downward) [m/d]
   real(rk)           :: w_ipw          =-1.0_rk     ! vertical speed of suspended iron PO4 [m/d]
   real(rk)           :: w_ipw_sedi     =-0.5_rk     ! sedimentation speed for iron PO4 [m/d]
   real(rk)           :: w_lpp          =-0.5_rk     ! vertical speed of large-cell phytoplankton [m/d]
   real(rk)           :: w_n2_stf       =5.0_rk      ! piston velocity for n2 surface flux [m/d]
   real(rk)           :: w_o2_stf       =5.0_rk      ! piston velocity for oxygen surface flux [m/d]
   real(rk)           :: zoo0           =4.5E-9_rk   ! seed concentration for zooplankton [mol/kg]
   real(rk)           :: zoo_cl         =9.0E-8_rk   ! zooplankton closure parameter [mol/kg]
   real(rk)           :: stf_no3        =0.0_rk      ! surface tracer flux for nitrate [mol/m2/d]
   real(rk)           :: stf_nh4        =0.0_rk      ! surface tracer flux for ammonium [mol/m2/d]
   real(rk)           :: stf_po4        =0.0_rk      ! surface tracer flux for phosphate [mol/m2/d]
   real(rk)           :: t_n2_initial ! dissolved molecular nitrogen
   real(rk)           :: t_o2_initial ! dissolved oxygen
   real(rk)           :: t_nh4_initial ! ammonium
   real(rk)           :: t_no3_initial ! nitrate
   real(rk)           :: t_po4_initial ! phosphate
   real(rk)           :: t_spp_initial ! small-cell phytoplankton
   real(rk)           :: t_zoo_initial ! zooplankton
   real(rk)           :: t_h2s_initial ! hydrogen sulfide
   real(rk)           :: t_sul_initial ! sulfur
   real(rk)           :: t_ipw_initial ! suspended iron phosphate
   real(rk)           :: t_lpp_initial ! large-cell phytoplankton
   real(rk)           :: t_cya_initial ! diazotroph cyanobacteria
   real(rk)           :: t_det_initial ! detritus
   real(rk)           :: t_sed_initial ! sediment detritus
   real(rk)           :: t_ips_initial ! iron phosphate in sediment
   real(rk)           :: cgt_timestep=1.0_rk  ! ecosystem model time step [d]
   real(rk), parameter:: secs_pr_day=86400.0_rk

   character(len=64)  :: dic_variable=''     ! Link variable to carbon  model

   namelist /iow_ergom/ &
                        critical_stress, &
                        cya0           , &
                        din_min_lpp    , &
                        din_min_spp    , &
                        dip_min_cya    , &
                        epsilon        , &
                        food_min_zoo   , &
                        gamma0         , &
                        gamma1         , &
                        h2s_min_po4_liber, &
                        ips_threshold  , &
                        ips_cl         , &
                        k_h2s_no3      , &
                        k_h2s_o2       , &
                        k_sul_no3      , &
                        k_sul_o2       , &
                        light_opt_cya  , &
                        light_opt_lpp  , &
                        light_opt_spp  , &
                        lpp0           , &
                        no3_min_sed_denit, &
                        no3_min_det_denit, &
                        o2_min_det_resp, &
                        o2_min_nit     , &
                        o2_min_po4_retent, &
                        o2_min_sed_resp, &
                        q10_det_rec    , &
                        q10_h2s        , &
                        q10_nit        , &
                        q10_sed_rec    , &
                        r_biores       , &
                        r_cya_assim    , &
                        r_det_rec      , &
                        r_ips_burial   , &
                        r_ips_ero      , &
                        r_ips_liber    , &
                        r_lpp_assim    , &
                        r_nh4_nitrif   , &
                        r_pp_mort      , &
                        r_pp_resp      , &
                        r_sed_ero      , &
                        r_sed_rec      , &
                        r_spp_assim    , &
                        r_zoo_graz     , &
                        r_zoo_mort     , &
                        r_zoo_resp     , &
                        rfr_c          , &
                        rfr_h          , &
                        rfr_o          , &
                        rfr_p          , &
                        sali_max_cya   , &
                        sali_min_cya   , &
                        sed_max        , &
                        spp0           , &
                        temp_min_cya   , &
                        temp_min_spp   , &
                        temp_opt_zoo   , &
                        w_cya          , &
                        w_det          , &
                        w_det_sedi     , &
                        w_ipw          , &
                        w_ipw_sedi     , &
                        w_lpp          , &
                        w_n2_stf       , &
                        w_o2_stf       , &
                        zoo0           , &
                        zoo_cl         , &
                        stf_no3        , &
                        stf_nh4        , &
                        stf_po4        , &
                        t_n2_initial, &
                        t_o2_initial, &
                        t_nh4_initial, &
                        t_no3_initial, &
                        t_po4_initial, &
                        t_spp_initial, &
                        t_zoo_initial, &
                        t_h2s_initial, &
                        t_sul_initial, &
                        t_ipw_initial, &
                        t_lpp_initial, &
                        t_cya_initial, &
                        t_det_initial, &
                        t_sed_initial, &
                        t_ips_initial, &
                        cgt_timestep,  &
                        dic_variable
!EOP
!-----------------------------------------------------------------------
!BOC
   ! Read the namelist
   if (configunit>0) read(configunit,nml=iow_ergom,err=99,end=100)

!  Store parameter values in our own derived type
   self%critical_stress   = critical_stress
   self%cya0              = cya0
   self%din_min_lpp       = din_min_lpp
   self%din_min_spp       = din_min_spp
   self%dip_min_cya       = dip_min_cya
   self%epsilon           = epsilon
   self%food_min_zoo      = food_min_zoo
   self%gamma0            = gamma0
   self%gamma1            = gamma1
   self%h2s_min_po4_liber   = h2s_min_po4_liber
   self%ips_threshold     = ips_threshold
   self%ips_cl            = ips_cl
   self%k_h2s_no3         = k_h2s_no3
   self%k_h2s_o2          = k_h2s_o2
   self%k_sul_no3         = k_sul_no3
   self%k_sul_o2          = k_sul_o2
   self%light_opt_cya     = light_opt_cya
   self%light_opt_lpp     = light_opt_lpp
   self%light_opt_spp     = light_opt_spp
   self%lpp0              = lpp0
   self%no3_min_sed_denit   = no3_min_sed_denit
   self%no3_min_det_denit   = no3_min_det_denit
   self%o2_min_det_resp   = o2_min_det_resp
   self%o2_min_nit        = o2_min_nit
   self%o2_min_po4_retent   = o2_min_po4_retent
   self%o2_min_sed_resp   = o2_min_sed_resp
   self%q10_det_rec       = q10_det_rec
   self%q10_h2s           = q10_h2s
   self%q10_nit           = q10_nit
   self%q10_sed_rec       = q10_sed_rec
   self%r_biores          = r_biores
   self%r_cya_assim       = r_cya_assim
   self%r_det_rec         = r_det_rec
   self%r_ips_burial      = r_ips_burial
   self%r_ips_ero         = r_ips_ero
   self%r_ips_liber       = r_ips_liber
   self%r_lpp_assim       = r_lpp_assim
   self%r_nh4_nitrif      = r_nh4_nitrif
   self%r_pp_mort         = r_pp_mort
   self%r_pp_resp         = r_pp_resp
   self%r_sed_ero         = r_sed_ero
   self%r_sed_rec         = r_sed_rec
   self%r_spp_assim       = r_spp_assim
   self%r_zoo_graz        = r_zoo_graz
   self%r_zoo_mort        = r_zoo_mort
   self%r_zoo_resp        = r_zoo_resp
   self%rfr_c             = rfr_c
   self%rfr_h             = rfr_h
   self%rfr_o             = rfr_o
   self%rfr_p             = rfr_p
   self%sali_max_cya      = sali_max_cya
   self%sali_min_cya      = sali_min_cya
   self%sed_max           = sed_max
   self%spp0              = spp0
   self%temp_min_cya      = temp_min_cya
   self%temp_min_spp      = temp_min_spp
   self%temp_opt_zoo      = temp_opt_zoo
   self%w_cya             = w_cya
   self%w_det             = w_det
   self%w_det_sedi        = w_det_sedi
   self%w_ipw             = w_ipw
   self%w_ipw_sedi        = w_ipw_sedi
   self%w_lpp             = w_lpp
   self%w_n2_stf          = w_n2_stf
   self%w_o2_stf          = w_o2_stf
   self%zoo0              = zoo0
   self%zoo_cl            = zoo_cl
   self%stf_no3           = stf_no3
   self%stf_nh4           = stf_nh4
   self%stf_po4           = stf_po4

   call self%get_parameter(self%cgt_timestep,'cgt_timestep','d','timestep',default=cgt_timestep)


!  Register state variables
   call self%register_state_variable(self%id_t_n2           ,'t_n2','mol/kg', &
         'dissolved molecular nitrogen', t_n2_initial &
         , minimum=_ZERO_ &
         )
   call self%register_state_variable(self%id_t_o2           ,'t_o2','mol/kg', &
         'dissolved oxygen', t_o2_initial &
         , minimum=_ZERO_ &
         , no_river_dilution=.true. &
         )
   call self%register_state_variable(self%id_t_nh4          ,'t_nh4','mol/kg', &
         'ammonium', t_nh4_initial &
         , minimum=_ZERO_ &
         , no_river_dilution=.true. &
         )
   call self%register_state_variable(self%id_t_no3          ,'t_no3','mol/kg', &
         'nitrate', t_no3_initial &
         , minimum=_ZERO_ &
         , no_river_dilution=.true. &
         )
   call self%register_state_variable(self%id_t_po4          ,'t_po4','mol/kg', &
         'phosphate', t_po4_initial &
         , minimum=_ZERO_ &
         , no_river_dilution=.true. &
         )
   call self%register_state_variable(self%id_t_spp          ,'t_spp','mol/kg', &
         'small-cell phytoplankton', t_spp_initial &
         , minimum=_ZERO_ &
         )
   call self%register_state_variable(self%id_t_zoo          ,'t_zoo','mol/kg', &
         'zooplankton', t_zoo_initial &
         , minimum=_ZERO_ &
         )
   call self%register_state_variable(self%id_t_h2s          ,'t_h2s','mol/kg', &
         'hydrogen sulfide', t_h2s_initial &
         , minimum=_ZERO_ &
         )
   call self%register_state_variable(self%id_t_sul          ,'t_sul','mol/kg', &
         'sulfur', t_sul_initial &
         , minimum=_ZERO_ &
         )
   call self%register_state_variable(self%id_t_ipw          ,'t_ipw','mol/kg', &
         'suspended iron phosphate', t_ipw_initial &
         , minimum=_ZERO_ &
         , vertical_movement=w_ipw/secs_pr_day &
         )
   call self%register_state_variable(self%id_t_lpp          ,'t_lpp','mol/kg', &
         'large-cell phytoplankton', t_lpp_initial &
         , minimum=_ZERO_ &
         , vertical_movement=w_lpp/secs_pr_day &
         )
   call self%register_state_variable(self%id_t_cya          ,'t_cya','mol/kg', &
         'diazotroph cyanobacteria', t_cya_initial &
         , minimum=_ZERO_ &
         , vertical_movement=w_cya/secs_pr_day &
         )
   call self%register_state_variable(self%id_t_det          ,'t_det','mol/kg', &
         'detritus', t_det_initial &
         , minimum=_ZERO_ &
         , vertical_movement=w_det/secs_pr_day &
         , no_river_dilution=.true. &
         )
   call self%register_state_variable(self%id_t_sed          ,'t_sed','mol/m**2','sediment detritus', &
         0.0_rk &
         , minimum=_ZERO_ &
         )
   call self%register_state_variable(self%id_t_ips          ,'t_ips','mol/m**2','iron phosphate in sediment', &
         0.0_rk &
         , minimum=_ZERO_ &
         )

   ! Register the contribution of all state variables to total nitrogen
!   call self%add_to_aggregate_variable(standard_variables%total_nitrogen,self%id_t_spp)
!   call self%add_to_aggregate_variable(standard_variables%total_nitrogen,self%id_t_lpp)
!   call self%add_to_aggregate_variable(standard_variables%total_nitrogen,self%id_t_zoo)
!   call self%add_to_aggregate_variable(standard_variables%total_nitrogen,self%id_t_det)
!   call self%add_to_aggregate_variable(standard_variables%total_nitrogen,self%id_t_cya)

! Register diagnostic variables
!   call self%register_diagnostic_variable(self%id_dPAR,'PAR','W/m**2','photosynthetically active radiation',   &
!         time_treatment=time_treatment_averaged)
!   call self%register_diagnostic_variable(self%id_GPP,'GPP','mmol/m**3','gross primary production',            &
!         time_treatment=time_treatment_step_integrated)
!   call self%register_diagnostic_variable(self%id_NCP,'NCP','mmol/m**3','net community production',            &
!         time_treatment=time_treatment_step_integrated)
!   call self%register_diagnostic_variable(self%id_PPR,'PPR','mmol/m**3/d','gross primary production rate',     &
!         time_treatment=time_treatment_averaged)
!   call self%register_diagnostic_variable(self%id_NPR,'NPR','mmol/m**3/d','net community production rate',     &
!         time_treatment=time_treatment_averaged)


   ! Register link to external DIC pool, if DIC variable name is provided in namelist.
   call self%register_state_dependency(self%id_t_dic,'dic','mmol/m**3','total dissolved inorganic carbon',required=.false.)
   if (dic_variable/='') call self%request_coupling(self%id_t_dic,dic_variable)

! Register environmental dependencies
   call self%register_dependency(self%id_par,standard_variables%downwelling_photosynthetic_radiative_flux)
   call self%register_dependency(self%id_temp,standard_variables%temperature)
   call self%register_dependency(self%id_salt,standard_variables%practical_salinity)
   call self%register_dependency(self%id_I_0,standard_variables%surface_downwelling_photosynthetic_radiative_flux)
   call self%register_dependency(self%id_wind,standard_variables%wind_speed)
   call self%register_dependency(self%id_depth,standard_variables%bottom_depth )
   call self%register_dependency(self%id_taub,standard_variables%bottom_stress)
   call self%register_dependency(self%id_lat,standard_variables%latitude)
   call self%register_dependency(self%id_lon,standard_variables%longitude)

   return

 99 call self%fatal_error('iow_ergom_init','Error reading namelist iow_ergom.')

100 call self%fatal_error('iow_ergom_init','Namelist iow_ergom was not found.')

   END subroutine initialize
!EOC


!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Right hand sides of ergom model
!
! !INTERFACE:
   subroutine do(self,_ARGUMENTS_DO_)
!!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   class(type_iow_ergom), INTENT(IN) :: self
  _DECLARE_ARGUMENTS_DO_
!
!
! !LOCAL VARIABLES:
    real(rk)        :: critical_stress
    real(rk)        :: cya0
    real(rk)        :: din_min_lpp
    real(rk)        :: din_min_spp
    real(rk)        :: dip_min_cya
    real(rk)        :: epsilon
    real(rk)        :: food_min_zoo
    real(rk)        :: gamma0
    real(rk)        :: gamma1
    real(rk)        :: h2s_min_po4_liber
    real(rk)        :: ips_threshold
    real(rk)        :: ips_cl
    real(rk)        :: k_h2s_no3
    real(rk)        :: k_h2s_o2
    real(rk)        :: k_sul_no3
    real(rk)        :: k_sul_o2
    real(rk)        :: light_opt_cya
    real(rk)        :: light_opt_lpp
    real(rk)        :: light_opt_spp
    real(rk)        :: lpp0
    real(rk)        :: no3_min_sed_denit
    real(rk)        :: no3_min_det_denit
    real(rk)        :: o2_min_det_resp
    real(rk)        :: o2_min_nit
    real(rk)        :: o2_min_po4_retent
    real(rk)        :: o2_min_sed_resp
    real(rk)        :: q10_det_rec
    real(rk)        :: q10_h2s
    real(rk)        :: q10_nit
    real(rk)        :: q10_sed_rec
    real(rk)        :: r_biores
    real(rk)        :: r_cya_assim
    real(rk)        :: r_det_rec
    real(rk)        :: r_ips_burial
    real(rk)        :: r_ips_ero
    real(rk)        :: r_ips_liber
    real(rk)        :: r_lpp_assim
    real(rk)        :: r_nh4_nitrif
    real(rk)        :: r_pp_mort
    real(rk)        :: r_pp_resp
    real(rk)        :: r_sed_ero
    real(rk)        :: r_sed_rec
    real(rk)        :: r_spp_assim
    real(rk)        :: r_zoo_graz
    real(rk)        :: r_zoo_mort
    real(rk)        :: r_zoo_resp
    real(rk)        :: rfr_c
    real(rk)        :: rfr_h
    real(rk)        :: rfr_o
    real(rk)        :: rfr_p
    real(rk)        :: sali_max_cya
    real(rk)        :: sali_min_cya
    real(rk)        :: sed_max
    real(rk)        :: spp0
    real(rk)        :: temp_min_cya
    real(rk)        :: temp_min_spp
    real(rk)        :: temp_opt_zoo
    real(rk)        :: w_cya
    real(rk)        :: w_det
    real(rk)        :: w_det_sedi
    real(rk)        :: w_ipw
    real(rk)        :: w_ipw_sedi
    real(rk)        :: w_lpp
    real(rk)        :: w_n2_stf
    real(rk)        :: w_o2_stf
    real(rk)        :: zoo0
    real(rk)        :: zoo_cl
    real(rk)        :: stf_no3
    real(rk)        :: stf_nh4
    real(rk)        :: stf_po4
    real(rk)        :: t_n2
    real(rk)        :: t_o2
    real(rk)        :: t_nh4
    real(rk)        :: t_no3
    real(rk)        :: t_po4
    real(rk)        :: t_spp
    real(rk)        :: t_zoo
    real(rk)        :: t_h2s
    real(rk)        :: t_sul
    real(rk)        :: t_ipw
    real(rk)        :: t_lpp
    real(rk)        :: t_cya
    real(rk)        :: t_det
    real(rk)        :: temp_k
    real(rk)        :: o2_sat
    real(rk)        :: n2_sat
    real(rk)        :: solubility_n2
    real(rk)        :: temp_sq
    real(rk)        :: zoo_eff
    real(rk)        :: din
    real(rk)        :: din_sq
    real(rk)        :: po4_sq
    real(rk)        :: pp
    real(rk)        :: lpp_plus_lpp0
    real(rk)        :: spp_plus_spp0
    real(rk)        :: cya_plus_cya0
    real(rk)        :: food_zoo
    real(rk)        :: lim_light_lpp
    real(rk)        :: lim_light_spp
    real(rk)        :: lim_light_cya
    real(rk)        :: lr_assim_lpp
    real(rk)        :: lr_assim_spp
    real(rk)        :: lr_assim_cya
    real(rk)        :: lr_graz_zoo
    real(rk)        :: frac_po4retent
    real(rk)        :: lim_t_n2_7
    real(rk)        :: lim_t_o2_0
    real(rk)        :: lim_t_o2_2
    real(rk)        :: lim_t_o2_4
    real(rk)        :: lim_t_o2_6
    real(rk)        :: lim_t_nh4_10
    real(rk)        :: lim_t_no3_1
    real(rk)        :: lim_t_no3_3
    real(rk)        :: lim_t_no3_8
    real(rk)        :: lim_t_po4_9
    real(rk)        :: lim_t_spp_12
    real(rk)        :: lim_t_zoo_14
    real(rk)        :: lim_t_h2s_5
    real(rk)        :: lim_t_h2s_18
    real(rk)        :: lim_t_sul_19
    real(rk)        :: lim_t_ipw_20
    real(rk)        :: lim_t_lpp_11
    real(rk)        :: lim_t_cya_13
    real(rk)        :: lim_t_det_15
    real(rk)        :: p_no3_assim_lpp
    real(rk)        :: p_nh4_assim_lpp
    real(rk)        :: p_no3_assim_spp
    real(rk)        :: p_nh4_assim_spp
    real(rk)        :: p_n2_assim_cya
    real(rk)        :: p_lpp_graz_zoo
    real(rk)        :: p_spp_graz_zoo
    real(rk)        :: p_cya_graz_zoo
    real(rk)        :: p_lpp_resp_nh4
    real(rk)        :: p_spp_resp_nh4
    real(rk)        :: p_cya_resp_nh4
    real(rk)        :: p_zoo_resp_nh4
    real(rk)        :: p_lpp_mort_det
    real(rk)        :: p_spp_mort_det
    real(rk)        :: p_cya_mort_det
    real(rk)        :: p_zoo_mort_det
    real(rk)        :: p_nh4_nit_no3
    real(rk)        :: p_det_resp_nh4
    real(rk)        :: p_det_denit_nh4
    real(rk)        :: p_det_sulf_nh4
    real(rk)        :: p_h2s_oxo2_sul
    real(rk)        :: p_h2s_oxno3_sul
    real(rk)        :: p_sul_oxo2_so4
    real(rk)        :: p_sul_oxno3_so4
    real(rk)        :: cgt_temp
    real(rk)        :: cgt_sali
    real(rk)        :: cgt_light
    real(rk)        :: cgt_density
    real(rk)        :: cgt_longitude
    real(rk)        :: cgt_latitude
    real(rk)        :: cgt_timestep
    real(rk)        :: secs_pr_day=86400.0_rk
    real(rk)        :: temp1, temp2, temp3, temp4, temp5, &
                       temp6, temp7, temp8, temp9

!EOP
!-----------------------------------------------------------------------
!BOC
!   Store values of constants in local variables
    critical_stress = self%critical_stress
    cya0            = self%cya0
    din_min_lpp     = self%din_min_lpp
    din_min_spp     = self%din_min_spp
    dip_min_cya     = self%dip_min_cya
    epsilon         = self%epsilon
    food_min_zoo    = self%food_min_zoo
    gamma0          = self%gamma0
    gamma1          = self%gamma1
    h2s_min_po4_liber = self%h2s_min_po4_liber
    ips_threshold   = self%ips_threshold
    ips_cl          = self%ips_cl
    k_h2s_no3       = self%k_h2s_no3
    k_h2s_o2        = self%k_h2s_o2
    k_sul_no3       = self%k_sul_no3
    k_sul_o2        = self%k_sul_o2
    light_opt_cya   = self%light_opt_cya
    light_opt_lpp   = self%light_opt_lpp
    light_opt_spp   = self%light_opt_spp
    lpp0            = self%lpp0
    no3_min_sed_denit = self%no3_min_sed_denit
    no3_min_det_denit = self%no3_min_det_denit
    o2_min_det_resp = self%o2_min_det_resp
    o2_min_nit      = self%o2_min_nit
    o2_min_po4_retent = self%o2_min_po4_retent
    o2_min_sed_resp = self%o2_min_sed_resp
    q10_det_rec     = self%q10_det_rec
    q10_h2s         = self%q10_h2s
    q10_nit         = self%q10_nit
    q10_sed_rec     = self%q10_sed_rec
    r_biores        = self%r_biores
    r_cya_assim     = self%r_cya_assim
    r_det_rec       = self%r_det_rec
    r_ips_burial    = self%r_ips_burial
    r_ips_ero       = self%r_ips_ero
    r_ips_liber     = self%r_ips_liber
    r_lpp_assim     = self%r_lpp_assim
    r_nh4_nitrif    = self%r_nh4_nitrif
    r_pp_mort       = self%r_pp_mort
    r_pp_resp       = self%r_pp_resp
    r_sed_ero       = self%r_sed_ero
    r_sed_rec       = self%r_sed_rec
    r_spp_assim     = self%r_spp_assim
    r_zoo_graz      = self%r_zoo_graz
    r_zoo_mort      = self%r_zoo_mort
    r_zoo_resp      = self%r_zoo_resp
    rfr_c           = self%rfr_c
    rfr_h           = self%rfr_h
    rfr_o           = self%rfr_o
    rfr_p           = self%rfr_p
    sali_max_cya    = self%sali_max_cya
    sali_min_cya    = self%sali_min_cya
    sed_max         = self%sed_max
    spp0            = self%spp0
    temp_min_cya    = self%temp_min_cya
    temp_min_spp    = self%temp_min_spp
    temp_opt_zoo    = self%temp_opt_zoo
    w_cya           = self%w_cya
    w_det           = self%w_det
    w_det_sedi      = self%w_det_sedi
    w_ipw           = self%w_ipw
    w_ipw_sedi      = self%w_ipw_sedi
    w_lpp           = self%w_lpp
    w_n2_stf        = self%w_n2_stf
    w_o2_stf        = self%w_o2_stf
    zoo0            = self%zoo0
    zoo_cl          = self%zoo_cl
    stf_no3         = self%stf_no3
    stf_nh4         = self%stf_nh4
    stf_po4         = self%stf_po4

!   Enter spatial_loops (if any)
    _LOOP_BEGIN_

!   Retrieve current (local) state variable values
    _GET_(self%id_t_n2           ,t_n2           ) !dissolved molecular nitrogen
    _GET_(self%id_t_o2           ,t_o2           ) !dissolved oxygen
    _GET_(self%id_t_nh4          ,t_nh4          ) !ammonium
    _GET_(self%id_t_no3          ,t_no3          ) !nitrate
    _GET_(self%id_t_po4          ,t_po4          ) !phosphate
    _GET_(self%id_t_spp          ,t_spp          ) !small-cell phytoplankton
    _GET_(self%id_t_zoo          ,t_zoo          ) !zooplankton
    _GET_(self%id_t_h2s          ,t_h2s          ) !hydrogen sulfide
    _GET_(self%id_t_sul          ,t_sul          ) !sulfur
    _GET_(self%id_t_ipw          ,t_ipw          ) !suspended iron phosphate
    _GET_(self%id_t_lpp          ,t_lpp          ) !large-cell phytoplankton
    _GET_(self%id_t_cya          ,t_cya          ) !diazotroph cyanobacteria
    _GET_(self%id_t_det          ,t_det          ) !detritus

!   Retrieve current environmental conditions
    _GET_HORIZONTAL_(self%id_lon,cgt_longitude)! Longitude [degE]
    _GET_HORIZONTAL_(self%id_lat,cgt_latitude) ! Latitude  [degN]
    _GET_   (self%id_par,cgt_light)    ! local photosynthetically active radiation
    _GET_   (self%id_temp,cgt_temp)    ! local water temperature
    _GET_   (self%id_salt,cgt_sali)    ! local water salinity
    cgt_density = 1035.0               ! Boussinesq density
    cgt_timestep  = self%cgt_timestep  ! ecosystem model time step [d]

!   Calculation of auxiliary variables
    temp_k          = cgt_temp + 273.15_rk
    o2_sat          = (10.18e0_rk+((5.306e-3_rk-4.8725e-5_rk*cgt_temp)*cgt_temp-0.2785e0_rk)*cgt_temp+cgt_sali*((2.2258e-3_rk+(4.39e-7_rk*cgt_temp-4.645e-5_rk)*cgt_temp)*cgt_temp-6.33e-2_rk))*44.66e0_rk*1e-6_rk
    temp1 = log((298.15_rk-cgt_temp)/(273.15_rk+cgt_temp))
    temp2 = temp1*temp1
    temp3 = temp2*temp1
    n2_sat          = 1e-6_rk*exp(6.42931_rk + 2.92704_rk*temp1 + 4.32531_rk*temp2 + 4.69149_rk*temp3 + cgt_sali*(0.0_rk -7.44129e-3_rk - 8.02566e-3_rk*temp1 - 1.46775e-2_rk*temp2))
    solubility_n2   = n2_sat * 1.263925e-5_rk
    temp_sq         = max(0.0_rk,cgt_temp)*max(0.0_rk,cgt_temp)
    zoo_eff         = t_zoo*t_zoo/zoo_cl
    din             = t_no3+t_nh4
    din_sq          = din*din
    po4_sq          = t_po4*t_po4
    pp              = t_lpp+t_spp+t_cya
    lpp_plus_lpp0   = t_lpp+lpp0
    spp_plus_spp0   = t_spp+spp0
    cya_plus_cya0   = t_cya+cya0
    food_zoo        = t_lpp+t_spp+0.5_rk*t_cya
    lim_light_lpp   = cgt_light/max(cgt_light/2.0_rk,light_opt_lpp)*exp(1-cgt_light/max(cgt_light/2.0_rk,light_opt_lpp))
    temp1 = max(cgt_light/2.0_rk,light_opt_spp)
    lim_light_spp   = cgt_light/temp1*exp(1-cgt_light/temp1)
    lim_light_cya   = cgt_light/max(cgt_light/2.0_rk,light_opt_lpp)*exp(1-cgt_light/max(cgt_light/2.0_rk,light_opt_lpp))
    lr_assim_lpp    = r_lpp_assim*theta(t_o2-2*t_h2s)*min(din_sq/(din_sq+din_min_lpp*din_min_lpp),min(po4_sq/(po4_sq+din_min_lpp*din_min_lpp*rfr_p*rfr_p),lim_light_lpp))
    lr_assim_spp    = r_spp_assim*theta(t_o2-2*t_h2s)*min(din_sq/(din_sq+din_min_spp*din_min_spp),min(po4_sq/(po4_sq+din_min_spp*din_min_spp*rfr_p*rfr_p),lim_light_spp))*(1+temp_sq/(temp_sq+temp_min_spp*temp_min_spp))
    lr_assim_cya    = r_cya_assim*theta(t_o2-2*t_h2s)*min(po4_sq/(po4_sq+dip_min_cya*dip_min_cya),lim_light_spp)*(1/(1+exp(temp_min_cya-cgt_temp)))*(1/(1+exp(cgt_sali-sali_max_cya)))*(1/(1+exp(sali_min_cya-cgt_sali)))
    lr_graz_zoo     = r_zoo_graz*(1-exp(-food_zoo*food_zoo/(food_min_zoo*food_min_zoo)))*theta(t_o2-2*t_h2s)*(1.0_rk+temp_sq/(temp_opt_zoo*temp_opt_zoo)*exp(2.0_rk-cgt_temp*2.0_rk/temp_opt_zoo))
    frac_po4retent  = 0.15_rk + 0.3_rk*theta(cgt_latitude-60.75_rk)

!   Calculation of process limitation factors
             lim_t_n2_7           = theta(t_n2-0.0_rk)
             lim_t_o2_0           = 1.0_rk-exp(-t_o2/o2_min_det_resp)
             lim_t_o2_2           = theta(t_o2-0.0_rk)
             lim_t_o2_4           = t_o2*t_o2/(t_o2*t_o2+o2_min_po4_retent*o2_min_po4_retent)
             lim_t_o2_6           = t_o2*t_o2/(t_o2*t_o2+o2_min_sed_resp*o2_min_sed_resp)
             lim_t_nh4_10         = theta(t_nh4-0.0_rk)
             lim_t_no3_1          = 1.0_rk-exp(-t_no3/no3_min_det_denit)
             lim_t_no3_3          = t_no3*t_no3/(t_no3*t_no3+no3_min_sed_denit*no3_min_sed_denit)
             lim_t_no3_8          = theta(t_no3-0.0_rk)
             lim_t_po4_9          = theta(t_po4-0.0_rk)
             lim_t_spp_12         = theta(t_spp-0.0_rk)
             lim_t_zoo_14         = theta(t_zoo-0.0_rk)
             lim_t_h2s_5          = theta(t_h2s-h2s_min_po4_liber)
             lim_t_h2s_18         = theta(t_h2s-0.0_rk)
             lim_t_sul_19         = theta(t_sul-0.0_rk)
             lim_t_ipw_20         = theta(t_ipw-0.0_rk)
             lim_t_lpp_11         = theta(t_lpp-0.0_rk)
             lim_t_cya_13         = theta(t_cya-0.0_rk)
             lim_t_det_15         = theta(t_det-0.0_rk)

!   Calculation of process rates
    p_no3_assim_lpp = (lpp_plus_lpp0*lr_assim_lpp*t_no3/(din+epsilon))*lim_t_no3_8*lim_t_po4_9
    p_nh4_assim_lpp = (lpp_plus_lpp0*lr_assim_lpp*t_nh4/(din+epsilon))*lim_t_po4_9*lim_t_nh4_10
    p_no3_assim_spp = (spp_plus_spp0*lr_assim_spp*t_no3/(din+epsilon))*lim_t_no3_8*lim_t_po4_9
    p_nh4_assim_spp = (spp_plus_spp0*lr_assim_spp*t_nh4/(din+epsilon))*lim_t_nh4_10*lim_t_po4_9
    p_n2_assim_cya  = (cya_plus_cya0*lr_assim_cya)*lim_t_n2_7*lim_t_po4_9
    p_lpp_graz_zoo  = ((t_zoo+zoo0)*lr_graz_zoo*t_lpp/max(food_zoo,epsilon))*lim_t_lpp_11
    p_spp_graz_zoo  = ((t_zoo+zoo0)*lr_graz_zoo*t_spp/max(food_zoo,epsilon))*lim_t_spp_12
    p_cya_graz_zoo  = ((t_zoo+zoo0)*lr_graz_zoo*(0.5_rk*t_cya)/max(food_zoo,epsilon))*lim_t_cya_13
    p_lpp_resp_nh4  = (t_lpp*r_pp_resp)*lim_t_lpp_11*lim_t_o2_2
    p_spp_resp_nh4  = (t_spp*r_pp_resp)*lim_t_o2_2*lim_t_spp_12
    p_cya_resp_nh4  = (t_cya*r_pp_resp)*lim_t_o2_2*lim_t_cya_13
    p_zoo_resp_nh4  = (zoo_eff*r_zoo_resp)*lim_t_zoo_14*lim_t_o2_2
    p_lpp_mort_det  = (t_lpp*r_pp_mort*(1+9*theta(5.0e-6_rk-t_o2)))*lim_t_lpp_11
    p_spp_mort_det  = (t_spp*r_pp_mort*(1+9*theta(5.0e-6_rk-t_o2)))*lim_t_spp_12
    p_cya_mort_det  = (t_cya*r_pp_mort*(1+9*theta(5.0e-6_rk-t_o2)))*lim_t_cya_13
    p_zoo_mort_det  = (zoo_eff*r_zoo_mort*(1+9*theta(5.0e-6_rk-t_o2)))*lim_t_zoo_14
    p_nh4_nit_no3   = (t_nh4*r_nh4_nitrif*exp(q10_nit*cgt_temp))*lim_t_o2_2*lim_t_nh4_10
    p_det_resp_nh4  = (t_det*r_det_rec*exp(q10_det_rec*cgt_temp))*lim_t_o2_0*lim_t_det_15
    p_det_denit_nh4 = (t_det*r_det_rec*exp(q10_det_rec*cgt_temp))*(1.0_rk-lim_t_o2_0)*lim_t_no3_1*lim_t_det_15
    p_det_sulf_nh4  = (t_det*r_det_rec*exp(q10_det_rec*cgt_temp))*(1.0_rk-lim_t_o2_0)*(1.0_rk-lim_t_no3_1)*lim_t_det_15
    p_h2s_oxo2_sul  = (t_h2s*t_o2*k_h2s_o2*exp(q10_h2s*cgt_temp))*lim_t_o2_2*lim_t_h2s_18
    p_h2s_oxno3_sul = (t_h2s*t_no3*k_h2s_no3*exp(q10_h2s*cgt_temp))*lim_t_h2s_18*lim_t_no3_8
    p_sul_oxo2_so4  = (t_sul*t_o2*k_sul_o2*exp(q10_h2s*cgt_temp))*lim_t_sul_19*lim_t_o2_2
    p_sul_oxno3_so4 = (t_sul*t_no3*k_sul_no3*exp(q10_h2s*cgt_temp))*lim_t_sul_19*lim_t_no3_8

!   Calculation of late auxiliary variables

!  Applying time tendencies for tracers
  _SET_ODE_(self%id_t_n2           ,0.0_rk       + (p_det_denit_nh4)*(2.65_rk)    /secs_pr_day       + (p_h2s_oxno3_sul)*(0.2_rk)     /secs_pr_day       + (p_sul_oxno3_so4)*(0.6_rk)     /secs_pr_day       - (p_n2_assim_cya)*(0.5_rk)      /secs_pr_day   )
  _SET_ODE_(self%id_t_o2           ,0.0_rk       + (p_no3_assim_lpp)*(8.625_rk)   /secs_pr_day       + (p_nh4_assim_lpp)*(6.625_rk)   /secs_pr_day       + (p_no3_assim_spp)*(8.625_rk)   /secs_pr_day       + (p_nh4_assim_spp)*(6.625_rk)   /secs_pr_day       + (p_n2_assim_cya)*(7.375_rk)    /secs_pr_day       - (p_lpp_resp_nh4)*(6.625_rk)    /secs_pr_day       - (p_spp_resp_nh4)*(6.625_rk)    /secs_pr_day       - (p_cya_resp_nh4)*(6.625_rk)    /secs_pr_day       - (p_zoo_resp_nh4)*(6.625_rk)    /secs_pr_day       - (p_nh4_nit_no3)*(2)         /secs_pr_day       - (p_det_resp_nh4)*(6.625_rk)    /secs_pr_day       - (p_h2s_oxo2_sul)*(0.5_rk)      /secs_pr_day       - (p_sul_oxo2_so4)*(1.5_rk)      /secs_pr_day   )
  _SET_ODE_(self%id_t_nh4          ,0.0_rk       + p_lpp_resp_nh4              /secs_pr_day       + p_spp_resp_nh4              /secs_pr_day       + p_cya_resp_nh4              /secs_pr_day       + p_zoo_resp_nh4              /secs_pr_day       + p_det_resp_nh4              /secs_pr_day       + p_det_denit_nh4             /secs_pr_day       + p_det_sulf_nh4              /secs_pr_day       - p_nh4_assim_lpp             /secs_pr_day       - p_nh4_assim_spp             /secs_pr_day       - p_nh4_nit_no3               /secs_pr_day   )
  _SET_ODE_(self%id_t_no3          ,0.0_rk       + p_nh4_nit_no3               /secs_pr_day       - p_no3_assim_lpp             /secs_pr_day       - p_no3_assim_spp             /secs_pr_day       - (p_det_denit_nh4)*(5.3_rk)     /secs_pr_day       - (p_h2s_oxno3_sul)*(0.4_rk)     /secs_pr_day       - (p_sul_oxno3_so4)*(1.2_rk)     /secs_pr_day   )
  _SET_ODE_(self%id_t_po4          ,0.0_rk       + (p_lpp_resp_nh4)*(rfr_p)    /secs_pr_day       + (p_spp_resp_nh4)*(rfr_p)    /secs_pr_day       + (p_cya_resp_nh4)*(rfr_p)    /secs_pr_day       + (p_zoo_resp_nh4)*(rfr_p)    /secs_pr_day       + (p_det_resp_nh4)*(rfr_p)    /secs_pr_day       + (p_det_denit_nh4)*(rfr_p)   /secs_pr_day       + (p_det_sulf_nh4)*(rfr_p)    /secs_pr_day       - (p_no3_assim_lpp)*(rfr_p)   /secs_pr_day       - (p_nh4_assim_lpp)*(rfr_p)   /secs_pr_day       - (p_no3_assim_spp)*(rfr_p)   /secs_pr_day       - (p_nh4_assim_spp)*(rfr_p)   /secs_pr_day       - (p_n2_assim_cya)*(rfr_p)    /secs_pr_day   )
  _SET_ODE_(self%id_t_spp          ,0.0_rk       + p_no3_assim_spp             /secs_pr_day       + p_nh4_assim_spp             /secs_pr_day       - p_spp_graz_zoo              /secs_pr_day       - p_spp_resp_nh4              /secs_pr_day       - p_spp_mort_det              /secs_pr_day   )
  _SET_ODE_(self%id_t_zoo          ,0.0_rk       + p_lpp_graz_zoo              /secs_pr_day       + p_spp_graz_zoo              /secs_pr_day       + p_cya_graz_zoo              /secs_pr_day       - p_zoo_resp_nh4              /secs_pr_day       - p_zoo_mort_det              /secs_pr_day   )
  _SET_ODE_(self%id_t_h2s          ,0.0_rk       + (p_det_sulf_nh4)*(3.3125_rk)   /secs_pr_day       - p_h2s_oxo2_sul              /secs_pr_day       - p_h2s_oxno3_sul             /secs_pr_day   )
  _SET_ODE_(self%id_t_sul          ,0.0_rk       + p_h2s_oxo2_sul              /secs_pr_day       + p_h2s_oxno3_sul             /secs_pr_day       - p_sul_oxo2_so4              /secs_pr_day       - p_sul_oxno3_so4             /secs_pr_day   )
  _SET_ODE_(self%id_t_ipw          ,0.0_rk   )
  _SET_ODE_(self%id_t_lpp          ,0.0_rk       + p_no3_assim_lpp             /secs_pr_day       + p_nh4_assim_lpp             /secs_pr_day       - p_lpp_graz_zoo              /secs_pr_day       - p_lpp_resp_nh4              /secs_pr_day       - p_lpp_mort_det              /secs_pr_day   )
  _SET_ODE_(self%id_t_cya          ,0.0_rk       + p_n2_assim_cya              /secs_pr_day       - p_cya_graz_zoo              /secs_pr_day       - p_cya_resp_nh4              /secs_pr_day       - p_cya_mort_det              /secs_pr_day   )
  _SET_ODE_(self%id_t_det          ,0.0_rk       + p_lpp_mort_det              /secs_pr_day       + p_spp_mort_det              /secs_pr_day       + p_cya_mort_det              /secs_pr_day       + p_zoo_mort_det              /secs_pr_day       - p_det_resp_nh4              /secs_pr_day       - p_det_denit_nh4             /secs_pr_day       - p_det_sulf_nh4              /secs_pr_day   )

    if (_AVAILABLE_(self%id_t_dic)) _SET_ODE_(self%id_t_dic, rfr_c*(r_pp_mort * food_zoo + r_zoo_resp * zoo_eff + r_sed_rec * t_det - r_det_rec * spp_plus_spp0 - r_lpp_assim * lpp_plus_lpp0 - r_cya_assim * cya_plus_cya0)/secs_pr_day)

! Export diagnostic variables
!   _SET_DIAGNOSTIC_(self%id_dPAR,par)
!   _SET_DIAGNOSTIC_(self%id_GPP ,r1/(p1+self%p10)+r2/(p2+self%p20)+r3/(p2+self%p30))
!   _SET_DIAGNOSTIC_(self%id_NCP ,r1/(p1+self%p10)+r2/(p2+self%p20)+r3/(p2+self%p30) - self%lpa*(p1+p2+p3))
!   _SET_DIAGNOSTIC_(self%id_PPR ,(r1/(p1+self%p10)+r2/(p2+self%p20)+r3/(p2+self%p30))*secs_pr_day)
!   _SET_DIAGNOSTIC_(self%id_NPR ,(r1/(p1+self%p10)+r2/(p2+self%p20)+r3/(p2+self%p30) - self%lpa*(p1+p2+p3))*secs_pr_day)
!   Leave spatial loops (if any)
   _LOOP_END_

   END subroutine do
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Get the light extinction coefficient due to biogeochemical
! variables
!
! !DESCRIPTION:

! !INTERFACE:
   subroutine get_light_extinction(self,_ARGUMENTS_GET_EXTINCTION_)
!
! !INPUT PARAMETERS:
   class (type_iow_ergom), intent(in) :: self
   _DECLARE_ARGUMENTS_GET_EXTINCTION_
!
! !REVISION HISTORY:
!  Original author(s): Jorn Bruggeman
!
! !LOCAL VARIABLES:
   real(rk)                     :: gamma1
   real(rk)                     :: t_spp
   real(rk)                     :: t_lpp
   real(rk)                     :: t_cya
   real(rk)                     :: t_det

!EOP
!-----------------------------------------------------------------------
!BOC
   ! Store values of constants in local variables
   gamma1          = self%gamma1

   ! Enter spatial loops (if any)
   _LOOP_BEGIN_

   ! Retrieve current (local) state variable values.

   _GET_(self%id_t_spp          ,t_spp          ) ! small-cell phytoplankton
   _GET_(self%id_t_lpp          ,t_lpp          ) ! large-cell phytoplankton
   _GET_(self%id_t_cya          ,t_cya          ) ! diazotroph cyanobacteria
   _GET_(self%id_t_det          ,t_det          ) ! detritus

   ! Self-shading with explicit contribution from background phytoplankton concentration.
   _SET_EXTINCTION_(0.0_rk                     + t_spp           * gamma1                     + t_lpp           * gamma1                     + t_cya           * gamma1                     + t_det           * gamma1                    )

   ! Leave spatial loops (if any)
   _LOOP_END_

   end subroutine get_light_extinction
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE:
!
! !INTERFACE:
   subroutine do_bottom(self,_ARGUMENTS_DO_BOTTOM_)
!
! !DESCRIPTION:
! Calculating the benthic fluxes
!
   implicit none

! !INPUT PARAMETERS:
   class (type_iow_ergom), intent(in) :: self
   _DECLARE_ARGUMENTS_DO_BOTTOM_
!
! !REVISION HISTORY:
!  Original author(s):
!
! !LOCAL VARIABLES:
   real(rk)                   :: critical_stress
   real(rk)                   :: cya0
   real(rk)                   :: din_min_lpp
   real(rk)                   :: din_min_spp
   real(rk)                   :: dip_min_cya
   real(rk)                   :: epsilon
   real(rk)                   :: food_min_zoo
   real(rk)                   :: gamma0
   real(rk)                   :: gamma1
   real(rk)                   :: h2s_min_po4_liber
   real(rk)                   :: ips_threshold
   real(rk)                   :: ips_cl
   real(rk)                   :: k_h2s_no3
   real(rk)                   :: k_h2s_o2
   real(rk)                   :: k_sul_no3
   real(rk)                   :: k_sul_o2
   real(rk)                   :: light_opt_cya
   real(rk)                   :: light_opt_lpp
   real(rk)                   :: light_opt_spp
   real(rk)                   :: lpp0
   real(rk)                   :: no3_min_sed_denit
   real(rk)                   :: no3_min_det_denit
   real(rk)                   :: o2_min_det_resp
   real(rk)                   :: o2_min_nit
   real(rk)                   :: o2_min_po4_retent
   real(rk)                   :: o2_min_sed_resp
   real(rk)                   :: q10_det_rec
   real(rk)                   :: q10_h2s
   real(rk)                   :: q10_nit
   real(rk)                   :: q10_sed_rec
   real(rk)                   :: r_biores
   real(rk)                   :: r_cya_assim
   real(rk)                   :: r_det_rec
   real(rk)                   :: r_ips_burial
   real(rk)                   :: r_ips_ero
   real(rk)                   :: r_ips_liber
   real(rk)                   :: r_lpp_assim
   real(rk)                   :: r_nh4_nitrif
   real(rk)                   :: r_pp_mort
   real(rk)                   :: r_pp_resp
   real(rk)                   :: r_sed_ero
   real(rk)                   :: r_sed_rec
   real(rk)                   :: r_spp_assim
   real(rk)                   :: r_zoo_graz
   real(rk)                   :: r_zoo_mort
   real(rk)                   :: r_zoo_resp
   real(rk)                   :: rfr_c
   real(rk)                   :: rfr_h
   real(rk)                   :: rfr_o
   real(rk)                   :: rfr_p
   real(rk)                   :: sali_max_cya
   real(rk)                   :: sali_min_cya
   real(rk)                   :: sed_max
   real(rk)                   :: spp0
   real(rk)                   :: temp_min_cya
   real(rk)                   :: temp_min_spp
   real(rk)                   :: temp_opt_zoo
   real(rk)                   :: w_cya
   real(rk)                   :: w_det
   real(rk)                   :: w_det_sedi
   real(rk)                   :: w_ipw
   real(rk)                   :: w_ipw_sedi
   real(rk)                   :: w_lpp
   real(rk)                   :: w_n2_stf
   real(rk)                   :: w_o2_stf
   real(rk)                   :: zoo0
   real(rk)                   :: zoo_cl
   real(rk)                   :: stf_no3
   real(rk)                   :: stf_nh4
   real(rk)                   :: stf_po4
   real(rk)                   :: t_n2
   real(rk)                   :: t_o2
   real(rk)                   :: t_nh4
   real(rk)                   :: t_no3
   real(rk)                   :: t_po4
   real(rk)                   :: t_spp
   real(rk)                   :: t_zoo
   real(rk)                   :: t_h2s
   real(rk)                   :: t_sul
   real(rk)                   :: t_ipw
   real(rk)                   :: t_lpp
   real(rk)                   :: t_cya
   real(rk)                   :: t_det
   real(rk)                   :: t_sed
   real(rk)                   :: t_ips
   real(rk)                   :: temp_k
   real(rk)                   :: o2_sat
   real(rk)                   :: n2_sat
   real(rk)                   :: solubility_n2
   real(rk)                   :: temp_sq
   real(rk)                   :: zoo_eff
   real(rk)                   :: din
   real(rk)                   :: din_sq
   real(rk)                   :: po4_sq
   real(rk)                   :: pp
   real(rk)                   :: lpp_plus_lpp0
   real(rk)                   :: spp_plus_spp0
   real(rk)                   :: cya_plus_cya0
   real(rk)                   :: food_zoo
   real(rk)                   :: lim_light_lpp
   real(rk)                   :: lim_light_spp
   real(rk)                   :: lim_light_cya
   real(rk)                   :: lr_assim_lpp
   real(rk)                   :: lr_assim_spp
   real(rk)                   :: lr_assim_cya
   real(rk)                   :: lr_graz_zoo
   real(rk)                   :: frac_po4retent
   real(rk)                   :: frac_denit_sed
   real(rk)                   :: sed_active
   real(rk)                   :: lr_sed_rec
   real(rk)                   :: ips_eff
   real(rk)                   :: erosion_is_active
    real(rk)        :: lim_t_n2_7
    real(rk)        :: lim_t_o2_0
    real(rk)        :: lim_t_o2_2
    real(rk)        :: lim_t_o2_4
    real(rk)        :: lim_t_o2_6
    real(rk)        :: lim_t_nh4_10
    real(rk)        :: lim_t_no3_1
    real(rk)        :: lim_t_no3_3
    real(rk)        :: lim_t_no3_8
    real(rk)        :: lim_t_po4_9
    real(rk)        :: lim_t_spp_12
    real(rk)        :: lim_t_zoo_14
    real(rk)        :: lim_t_h2s_5
    real(rk)        :: lim_t_h2s_18
    real(rk)        :: lim_t_sul_19
    real(rk)        :: lim_t_ipw_20
    real(rk)        :: lim_t_lpp_11
    real(rk)        :: lim_t_cya_13
    real(rk)        :: lim_t_det_15
    real(rk)        :: lim_t_sed_16
    real(rk)        :: lim_t_ips_17
   real(rk)                   :: p_sed_resp_nh4
   real(rk)                   :: p_nh4_nitdenit_n2
   real(rk)                   :: p_sed_denit_nh4
   real(rk)                   :: p_sed_sulf_nh4
   real(rk)                   :: p_po4_retent_ips
   real(rk)                   :: p_ips_liber_po4
   real(rk)                   :: p_det_sedi_sed
   real(rk)                   :: p_ipw_sedi_ips
   real(rk)                   :: p_sed_ero_det
   real(rk)                   :: p_ips_ero_ipw
   real(rk)                   :: p_sed_biores_det
   real(rk)                   :: p_ips_biores_ipw
   real(rk)                   :: p_sed_burial
   real(rk)                   :: p_ips_burial
   real(rk)                   :: cgt_temp
   real(rk)                   :: cgt_sali
   real(rk)                   :: cgt_light
   real(rk)                   :: cgt_density
   real(rk)                   :: cgt_longitude
   real(rk)                   :: cgt_latitude
   real(rk)                   :: cgt_current_wave_stress
   real(rk)                   :: cgt_bottomdepth
   real(rk)                   :: cgt_timestep
   real(rk)                   :: temp1, temp2, temp3, temp4, temp5, &
                                 temp6, temp7, temp8, temp9
   real(rk), parameter:: secs_pr_day=86400.0_rk


!EOP
!-----------------------------------------------------------------------
!BOC

   ! Store values of constants in local variables
      critical_stress = self%critical_stress
      cya0            = self%cya0
      din_min_lpp     = self%din_min_lpp
      din_min_spp     = self%din_min_spp
      dip_min_cya     = self%dip_min_cya
      epsilon         = self%epsilon
      food_min_zoo    = self%food_min_zoo
      gamma0          = self%gamma0
      gamma1          = self%gamma1
      h2s_min_po4_liber = self%h2s_min_po4_liber
      ips_threshold   = self%ips_threshold
      ips_cl          = self%ips_cl
      k_h2s_no3       = self%k_h2s_no3
      k_h2s_o2        = self%k_h2s_o2
      k_sul_no3       = self%k_sul_no3
      k_sul_o2        = self%k_sul_o2
      light_opt_cya   = self%light_opt_cya
      light_opt_lpp   = self%light_opt_lpp
      light_opt_spp   = self%light_opt_spp
      lpp0            = self%lpp0
      no3_min_sed_denit = self%no3_min_sed_denit
      no3_min_det_denit = self%no3_min_det_denit
      o2_min_det_resp = self%o2_min_det_resp
      o2_min_nit      = self%o2_min_nit
      o2_min_po4_retent = self%o2_min_po4_retent
      o2_min_sed_resp = self%o2_min_sed_resp
      q10_det_rec     = self%q10_det_rec
      q10_h2s         = self%q10_h2s
      q10_nit         = self%q10_nit
      q10_sed_rec     = self%q10_sed_rec
      r_biores        = self%r_biores
      r_cya_assim     = self%r_cya_assim
      r_det_rec       = self%r_det_rec
      r_ips_burial    = self%r_ips_burial
      r_ips_ero       = self%r_ips_ero
      r_ips_liber     = self%r_ips_liber
      r_lpp_assim     = self%r_lpp_assim
      r_nh4_nitrif    = self%r_nh4_nitrif
      r_pp_mort       = self%r_pp_mort
      r_pp_resp       = self%r_pp_resp
      r_sed_ero       = self%r_sed_ero
      r_sed_rec       = self%r_sed_rec
      r_spp_assim     = self%r_spp_assim
      r_zoo_graz      = self%r_zoo_graz
      r_zoo_mort      = self%r_zoo_mort
      r_zoo_resp      = self%r_zoo_resp
      rfr_c           = self%rfr_c
      rfr_h           = self%rfr_h
      rfr_o           = self%rfr_o
      rfr_p           = self%rfr_p
      sali_max_cya    = self%sali_max_cya
      sali_min_cya    = self%sali_min_cya
      sed_max         = self%sed_max
      spp0            = self%spp0
      temp_min_cya    = self%temp_min_cya
      temp_min_spp    = self%temp_min_spp
      temp_opt_zoo    = self%temp_opt_zoo
      w_cya           = self%w_cya
      w_det           = self%w_det
      w_det_sedi      = self%w_det_sedi
      w_ipw           = self%w_ipw
      w_ipw_sedi      = self%w_ipw_sedi
      w_lpp           = self%w_lpp
      w_n2_stf        = self%w_n2_stf
      w_o2_stf        = self%w_o2_stf
      zoo0            = self%zoo0
      zoo_cl          = self%zoo_cl
      stf_no3         = self%stf_no3
      stf_nh4         = self%stf_nh4
      stf_po4         = self%stf_po4

   ! Enter spatial loops over the horizontal domain (if any)._rk
   _HORIZONTAL_LOOP_BEGIN_

   ! Retrieve abiotic parameters
      _GET_HORIZONTAL_(self%id_taub,cgt_current_wave_stress)
      _GET_HORIZONTAL_(self%id_depth,cgt_bottomdepth) !local water depth
      _GET_HORIZONTAL_(self%id_lon,cgt_longitude)     ! Longitude [degE]
      _GET_HORIZONTAL_(self%id_lat,cgt_latitude)      ! Latitude  [degN]
      _GET_   (self%id_par,cgt_light)    ! local photosynthetically active radiation
      _GET_   (self%id_temp,cgt_temp)    ! local water temperature
      _GET_   (self%id_salt,cgt_sali)    ! local water salinity
      cgt_density   = 1035.0_rk          ! Boussinesq density
      cgt_timestep  = self%cgt_timestep  ! ecosystem model time step [d]

   ! Retrieve current (local) state variable values.
      _GET_(self%id_t_n2           ,t_n2           )
      _GET_(self%id_t_o2           ,t_o2           )
      _GET_(self%id_t_nh4          ,t_nh4          )
      _GET_(self%id_t_no3          ,t_no3          )
      _GET_(self%id_t_po4          ,t_po4          )
      _GET_(self%id_t_spp          ,t_spp          )
      _GET_(self%id_t_zoo          ,t_zoo          )
      _GET_(self%id_t_h2s          ,t_h2s          )
      _GET_(self%id_t_sul          ,t_sul          )
      _GET_(self%id_t_ipw          ,t_ipw          )
      _GET_(self%id_t_lpp          ,t_lpp          )
      _GET_(self%id_t_cya          ,t_cya          )
      _GET_(self%id_t_det          ,t_det          )
      _GET_HORIZONTAL_(self%id_t_sed          ,t_sed          )
      _GET_HORIZONTAL_(self%id_t_ips          ,t_ips          )

   ! Calculate values of auxiliary variables
       temp_k          = cgt_temp + 273.15_rk
       o2_sat          = (10.18e0_rk+((5.306e-3_rk-4.8725e-5_rk*cgt_temp)*cgt_temp-0.2785e0_rk)*cgt_temp+cgt_sali*((2.2258e-3_rk+(4.39e-7_rk*cgt_temp-4.645e-5_rk)*cgt_temp)*cgt_temp-6.33e-2_rk))*44.66e0_rk*1e-6_rk
       temp1 = log((298.15_rk-cgt_temp)/(273.15_rk+cgt_temp))
       temp2 = temp1*temp1
       temp3 = temp2*temp1
       n2_sat          = 1e-6_rk*exp(6.42931_rk + 2.92704_rk*temp1 + 4.32531_rk*temp2 + 4.69149_rk*temp3 + cgt_sali*(0.0_rk -7.44129e-3_rk - 8.02566e-3_rk*temp1 - 1.46775e-2_rk*temp2))
       solubility_n2   = n2_sat * 1.263925e-5_rk
       temp_sq         = max(0.0_rk,cgt_temp)*max(0.0_rk,cgt_temp)
       zoo_eff         = t_zoo*t_zoo/zoo_cl
       din             = t_no3+t_nh4
       din_sq          = din*din
       po4_sq          = t_po4*t_po4
       pp              = t_lpp+t_spp+t_cya
       lpp_plus_lpp0   = t_lpp+lpp0
       spp_plus_spp0   = t_spp+spp0
       cya_plus_cya0   = t_cya+cya0
       food_zoo        = t_lpp+t_spp+0.5_rk*t_cya
       lim_light_lpp   = cgt_light/max(cgt_light/2.0_rk,light_opt_lpp)*exp(1-cgt_light/max(cgt_light/2.0_rk,light_opt_lpp))
       temp1 = max(cgt_light/2.0_rk,light_opt_spp)
       lim_light_spp   = cgt_light/temp1*exp(1-cgt_light/temp1)
       lim_light_cya   = cgt_light/max(cgt_light/2.0_rk,light_opt_lpp)*exp(1-cgt_light/max(cgt_light/2.0_rk,light_opt_lpp))
       lr_assim_lpp    = r_lpp_assim*theta(t_o2-2*t_h2s)*min(din_sq/(din_sq+din_min_lpp*din_min_lpp),min(po4_sq/(po4_sq+din_min_lpp*din_min_lpp*rfr_p*rfr_p),lim_light_lpp))
       lr_assim_spp    = r_spp_assim*theta(t_o2-2*t_h2s)*min(din_sq/(din_sq+din_min_spp*din_min_spp),min(po4_sq/(po4_sq+din_min_spp*din_min_spp*rfr_p*rfr_p),lim_light_spp))*(1+temp_sq/(temp_sq+temp_min_spp*temp_min_spp))
       lr_assim_cya    = r_cya_assim*theta(t_o2-2*t_h2s)*min(po4_sq/(po4_sq+dip_min_cya*dip_min_cya),lim_light_spp)*(1/(1+exp(temp_min_cya-cgt_temp)))*(1/(1+exp(cgt_sali-sali_max_cya)))*(1/(1+exp(sali_min_cya-cgt_sali)))
       lr_graz_zoo     = r_zoo_graz*(1-exp(-food_zoo*food_zoo/(food_min_zoo*food_min_zoo)))*theta(t_o2-2*t_h2s)*(1.0_rk+temp_sq/(temp_opt_zoo*temp_opt_zoo)*exp(2.0_rk-cgt_temp*2.0_rk/temp_opt_zoo))
       frac_po4retent  = 0.15_rk + 0.3_rk*theta(cgt_latitude-60.75_rk)
       frac_denit_sed  = (1-theta(cgt_bottomdepth-6.8_rk))*0.85_rk + theta(cgt_bottomdepth-6.8_rk)*(1-theta(cgt_bottomdepth-9.8_rk))*(0.85_rk-(0.1_rk/3.0_rk)*(cgt_bottomdepth-6.8_rk)) + theta(cgt_bottomdepth-9.8_rk)*(1-theta(cgt_bottomdepth-92.7_rk))*0.75_rk + theta(cgt_bottomdepth-92.7_rk)*(1-theta(cgt_bottomdepth-167.5_rk))*(0.75_rk-(0.15_rk/75.0_rk)*(cgt_bottomdepth-92.7_rk)) + theta(cgt_bottomdepth-167.5_rk)*0.6_rk
       sed_active      = max(0.0_rk,min(t_sed,sed_max/4.5_rk))
       lr_sed_rec      = r_sed_rec*exp(q10_sed_rec*cgt_temp)*(1.0_rk-0.9_rk*theta(2*t_h2s-t_o2))
       ips_eff         = max(t_ips,t_ips+t_ips*(t_ips-ips_threshold)/ips_cl)
       erosion_is_active = theta(cgt_current_wave_stress - critical_stress)

   ! Calculate process limitation functions
       lim_t_n2_7           = theta(t_n2-0.0_rk)
       lim_t_o2_0           = 1.0_rk-exp(-t_o2/o2_min_det_resp)
       lim_t_o2_2           = theta(t_o2-0.0_rk)
       lim_t_o2_4           = t_o2*t_o2/(t_o2*t_o2+o2_min_po4_retent*o2_min_po4_retent)
       lim_t_o2_6           = t_o2*t_o2/(t_o2*t_o2+o2_min_sed_resp*o2_min_sed_resp)
       lim_t_nh4_10         = theta(t_nh4-0.0_rk)
       lim_t_no3_1          = 1.0_rk-exp(-t_no3/no3_min_det_denit)
       lim_t_no3_3          = t_no3*t_no3/(t_no3*t_no3+no3_min_sed_denit*no3_min_sed_denit)
       lim_t_no3_8          = theta(t_no3-0.0_rk)
       lim_t_po4_9          = theta(t_po4-0.0_rk)
       lim_t_spp_12         = theta(t_spp-0.0_rk)
       lim_t_zoo_14         = theta(t_zoo-0.0_rk)
       lim_t_h2s_5          = theta(t_h2s-h2s_min_po4_liber)
       lim_t_h2s_18         = theta(t_h2s-0.0_rk)
       lim_t_sul_19         = theta(t_sul-0.0_rk)
       lim_t_ipw_20         = theta(t_ipw-0.0_rk)
       lim_t_lpp_11         = theta(t_lpp-0.0_rk)
       lim_t_cya_13         = theta(t_cya-0.0_rk)
       lim_t_det_15         = theta(t_det-0.0_rk)
       lim_t_sed_16         = theta(t_sed-0.0_rk)
       lim_t_ips_17         = theta(t_ips-0.0_rk)

   ! Calculate process rates
       p_sed_resp_nh4  = (lr_sed_rec*sed_active)*lim_t_sed_16*lim_t_o2_2
       p_nh4_nitdenit_n2 = (frac_denit_sed*p_sed_resp_nh4*theta(t_o2-5.0e-6_rk))*lim_t_o2_2*lim_t_nh4_10
       p_sed_denit_nh4 = (lr_sed_rec*sed_active)*(1.0_rk-lim_t_o2_2)*lim_t_no3_3*lim_t_sed_16
       p_sed_sulf_nh4  = (lr_sed_rec*sed_active)*(1.0_rk-lim_t_o2_2)*(1.0_rk-lim_t_no3_3)*lim_t_sed_16
       p_po4_retent_ips = (p_sed_resp_nh4*frac_po4retent)*lim_t_o2_4*lim_t_po4_9
       p_ips_liber_po4 = (t_ips*r_ips_liber)*lim_t_h2s_5*lim_t_ips_17
       p_det_sedi_sed  = ((1.0_rk-erosion_is_active)*(0.0_rk-w_det_sedi)*t_det*cgt_density)*lim_t_det_15
       p_ipw_sedi_ips  = ((1.0_rk-erosion_is_active)*(0.0_rk-w_ipw_sedi)*t_ipw*cgt_density)*lim_t_ipw_20
       p_sed_ero_det   = (erosion_is_active*r_sed_ero*sed_active)*lim_t_sed_16
       p_ips_ero_ipw   = (erosion_is_active*r_ips_ero*t_ips)*lim_t_ips_17
       p_sed_biores_det = (r_biores*sed_active)*lim_t_o2_6*lim_t_sed_16
       p_ips_biores_ipw = (r_biores*t_ips)*lim_t_o2_6*lim_t_ips_17
       p_sed_burial    = ((t_sed-sed_active)/cgt_timestep)*lim_t_sed_16
       p_ips_burial    = (ips_eff*r_ips_burial*sed_active/sed_max)*lim_t_ips_17

   ! Apply time tendencies to tracers
      _SET_ODE_BEN_(self%id_t_sed          ,0.0_rk                     + p_det_sedi_sed              /secs_pr_day                     - p_sed_resp_nh4              /secs_pr_day                     - p_sed_denit_nh4             /secs_pr_day                     - p_sed_sulf_nh4              /secs_pr_day                     - p_sed_ero_det               /secs_pr_day                     - p_sed_biores_det            /secs_pr_day                     - p_sed_burial                /secs_pr_day                    )
      _SET_ODE_BEN_(self%id_t_ips          ,0.0_rk                     + (p_po4_retent_ips)*(rfr_p)  /secs_pr_day                     + p_ipw_sedi_ips              /secs_pr_day                     - p_ips_liber_po4             /secs_pr_day                     - p_ips_ero_ipw               /secs_pr_day                     - p_ips_biores_ipw            /secs_pr_day                     - p_ips_burial                /secs_pr_day                    )
      _SET_BOTTOM_EXCHANGE_(self%id_t_n2           ,0.0_rk                     + (p_nh4_nitdenit_n2)*(0.5_rk)   /(1000.0_rk*cgt_density*secs_pr_day)                     + (p_sed_denit_nh4)*(2.65_rk)    /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_BOTTOM_EXCHANGE_(self%id_t_o2           ,0.0_rk                     - (p_sed_resp_nh4)*(6.625_rk)    /(1000.0_rk*cgt_density*secs_pr_day)                     - (p_nh4_nitdenit_n2)*(0.75_rk)  /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_BOTTOM_EXCHANGE_(self%id_t_nh4          ,0.0_rk                     + p_sed_resp_nh4              /(1000.0_rk*cgt_density*secs_pr_day)                     + p_sed_denit_nh4             /(1000.0_rk*cgt_density*secs_pr_day)                     + p_sed_sulf_nh4              /(1000.0_rk*cgt_density*secs_pr_day)                     - p_nh4_nitdenit_n2           /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_BOTTOM_EXCHANGE_(self%id_t_no3          ,0.0_rk                     - (p_sed_denit_nh4)*(5.3_rk)     /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_BOTTOM_EXCHANGE_(self%id_t_po4          ,0.0_rk                     + (p_sed_resp_nh4)*(rfr_p)    /(1000.0_rk*cgt_density*secs_pr_day)                     + (p_sed_denit_nh4)*(rfr_p)   /(1000.0_rk*cgt_density*secs_pr_day)                     + (p_sed_sulf_nh4)*(rfr_p)    /(1000.0_rk*cgt_density*secs_pr_day)                     + p_ips_liber_po4             /(1000.0_rk*cgt_density*secs_pr_day)                     - (p_po4_retent_ips)*(rfr_p)  /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_BOTTOM_EXCHANGE_(self%id_t_h2s          ,0.0_rk                     + (p_sed_sulf_nh4)*(3.3125_rk)   /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_BOTTOM_EXCHANGE_(self%id_t_ipw          ,0.0_rk                     + p_ips_ero_ipw               /(1000.0_rk*cgt_density*secs_pr_day)                     + p_ips_biores_ipw            /(1000.0_rk*cgt_density*secs_pr_day)                     - p_ipw_sedi_ips              /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_BOTTOM_EXCHANGE_(self%id_t_det          ,0.0_rk                     + p_sed_ero_det               /(1000.0_rk*cgt_density*secs_pr_day)                     + p_sed_biores_det            /(1000.0_rk*cgt_density*secs_pr_day)                     - p_det_sedi_sed              /(1000.0_rk*cgt_density*secs_pr_day)                            )

    if (_AVAILABLE_(self%id_t_dic)) _SET_BOTTOM_EXCHANGE_(self%id_t_dic, rfr_c*lr_sed_rec/(1000.0_rk*cgt_density*secs_pr_day))

   ! Leave spatial loops over the horizontal domain (if any)._rk
   _HORIZONTAL_LOOP_END_

   end subroutine do_bottom
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Surface fluxes for the ergom model
!
! !INTERFACE:

   subroutine do_surface(self,_ARGUMENTS_DO_SURFACE_)
   class (type_iow_ergom),intent(in) :: self
   _DECLARE_ARGUMENTS_DO_SURFACE_
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
! !LOCAL VARIABLES:
   real(rk)                   :: critical_stress
   real(rk)                   :: cya0
   real(rk)                   :: din_min_lpp
   real(rk)                   :: din_min_spp
   real(rk)                   :: dip_min_cya
   real(rk)                   :: epsilon
   real(rk)                   :: food_min_zoo
   real(rk)                   :: gamma0
   real(rk)                   :: gamma1
   real(rk)                   :: h2s_min_po4_liber
   real(rk)                   :: ips_threshold
   real(rk)                   :: ips_cl
   real(rk)                   :: k_h2s_no3
   real(rk)                   :: k_h2s_o2
   real(rk)                   :: k_sul_no3
   real(rk)                   :: k_sul_o2
   real(rk)                   :: light_opt_cya
   real(rk)                   :: light_opt_lpp
   real(rk)                   :: light_opt_spp
   real(rk)                   :: lpp0
   real(rk)                   :: no3_min_sed_denit
   real(rk)                   :: no3_min_det_denit
   real(rk)                   :: o2_min_det_resp
   real(rk)                   :: o2_min_nit
   real(rk)                   :: o2_min_po4_retent
   real(rk)                   :: o2_min_sed_resp
   real(rk)                   :: q10_det_rec
   real(rk)                   :: q10_h2s
   real(rk)                   :: q10_nit
   real(rk)                   :: q10_sed_rec
   real(rk)                   :: r_biores
   real(rk)                   :: r_cya_assim
   real(rk)                   :: r_det_rec
   real(rk)                   :: r_ips_burial
   real(rk)                   :: r_ips_ero
   real(rk)                   :: r_ips_liber
   real(rk)                   :: r_lpp_assim
   real(rk)                   :: r_nh4_nitrif
   real(rk)                   :: r_pp_mort
   real(rk)                   :: r_pp_resp
   real(rk)                   :: r_sed_ero
   real(rk)                   :: r_sed_rec
   real(rk)                   :: r_spp_assim
   real(rk)                   :: r_zoo_graz
   real(rk)                   :: r_zoo_mort
   real(rk)                   :: r_zoo_resp
   real(rk)                   :: rfr_c
   real(rk)                   :: rfr_h
   real(rk)                   :: rfr_o
   real(rk)                   :: rfr_p
   real(rk)                   :: sali_max_cya
   real(rk)                   :: sali_min_cya
   real(rk)                   :: sed_max
   real(rk)                   :: spp0
   real(rk)                   :: temp_min_cya
   real(rk)                   :: temp_min_spp
   real(rk)                   :: temp_opt_zoo
   real(rk)                   :: w_cya
   real(rk)                   :: w_det
   real(rk)                   :: w_det_sedi
   real(rk)                   :: w_ipw
   real(rk)                   :: w_ipw_sedi
   real(rk)                   :: w_lpp
   real(rk)                   :: w_n2_stf
   real(rk)                   :: w_o2_stf
   real(rk)                   :: zoo0
   real(rk)                   :: zoo_cl
   real(rk)                   :: stf_no3
   real(rk)                   :: stf_nh4
   real(rk)                   :: stf_po4
   real(rk)                   :: t_n2
   real(rk)                   :: t_o2
   real(rk)                   :: t_nh4
   real(rk)                   :: t_no3
   real(rk)                   :: t_po4
   real(rk)                   :: t_spp
   real(rk)                   :: t_zoo
   real(rk)                   :: t_h2s
   real(rk)                   :: t_sul
   real(rk)                   :: t_ipw
   real(rk)                   :: t_lpp
   real(rk)                   :: t_cya
   real(rk)                   :: t_det
   real(rk)                   :: temp_k
   real(rk)                   :: o2_sat
   real(rk)                   :: n2_sat
   real(rk)                   :: solubility_n2
   real(rk)                   :: temp_sq
   real(rk)                   :: zoo_eff
   real(rk)                   :: din
   real(rk)                   :: din_sq
   real(rk)                   :: po4_sq
   real(rk)                   :: pp
   real(rk)                   :: lpp_plus_lpp0
   real(rk)                   :: spp_plus_spp0
   real(rk)                   :: cya_plus_cya0
   real(rk)                   :: food_zoo
   real(rk)                   :: lim_light_lpp
   real(rk)                   :: lim_light_spp
   real(rk)                   :: lim_light_cya
   real(rk)                   :: lr_assim_lpp
   real(rk)                   :: lr_assim_spp
   real(rk)                   :: lr_assim_cya
   real(rk)                   :: lr_graz_zoo
   real(rk)                   :: frac_po4retent
   real(rk)                   :: solubility_o2
   real(rk)                   :: schmidtnumber_o2
   real(rk)                   :: schmidtnumber_n2
   real(rk)                   :: lim_t_n2_7
   real(rk)                   :: lim_t_o2_0
   real(rk)                   :: lim_t_o2_2
   real(rk)                   :: lim_t_o2_4
   real(rk)                   :: lim_t_o2_6
   real(rk)                   :: lim_t_nh4_10
   real(rk)                   :: lim_t_no3_1
   real(rk)                   :: lim_t_no3_3
   real(rk)                   :: lim_t_no3_8
   real(rk)                   :: lim_t_po4_9
   real(rk)                   :: lim_t_spp_12
   real(rk)                   :: lim_t_zoo_14
   real(rk)                   :: lim_t_h2s_5
   real(rk)                   :: lim_t_h2s_18
   real(rk)                   :: lim_t_sul_19
   real(rk)                   :: lim_t_ipw_20
   real(rk)                   :: lim_t_lpp_11
   real(rk)                   :: lim_t_cya_13
   real(rk)                   :: lim_t_det_15
   real(rk)                   :: p_n2_stf_down
   real(rk)                   :: p_n2_stf_up
   real(rk)                   :: p_o2_stf_down
   real(rk)                   :: p_o2_stf_up
   real(rk)                   :: p_no3_stf_down
   real(rk)                   :: p_nh4_stf_down
   real(rk)                   :: p_po4_stf_down
   real(rk)                   :: cgt_temp
   real(rk)                   :: cgt_sali
   real(rk)                   :: cgt_light
   real(rk)                   :: cgt_density
   real(rk)                   :: cgt_longitude
   real(rk)                   :: cgt_latitude
   real(rk)                   :: cgt_current_wave_stress
   real(rk)                   :: cgt_timestep
   real(rk)                   :: temp1, temp2, temp3, temp4, temp5, &
                                 temp6, temp7, temp8, temp9
   real(rk), parameter:: secs_pr_day=86400.0_rk

!EOP
!-----------------------------------------------------------------------
!BOC
      ! Store values of constants in local variables
      critical_stress = self%critical_stress
      cya0            = self%cya0
      din_min_lpp     = self%din_min_lpp
      din_min_spp     = self%din_min_spp
      dip_min_cya     = self%dip_min_cya
      epsilon         = self%epsilon
      food_min_zoo    = self%food_min_zoo
      gamma0          = self%gamma0
      gamma1          = self%gamma1
      h2s_min_po4_liber = self%h2s_min_po4_liber
      ips_threshold   = self%ips_threshold
      ips_cl          = self%ips_cl
      k_h2s_no3       = self%k_h2s_no3
      k_h2s_o2        = self%k_h2s_o2
      k_sul_no3       = self%k_sul_no3
      k_sul_o2        = self%k_sul_o2
      light_opt_cya   = self%light_opt_cya
      light_opt_lpp   = self%light_opt_lpp
      light_opt_spp   = self%light_opt_spp
      lpp0            = self%lpp0
      no3_min_sed_denit = self%no3_min_sed_denit
      no3_min_det_denit = self%no3_min_det_denit
      o2_min_det_resp = self%o2_min_det_resp
      o2_min_nit      = self%o2_min_nit
      o2_min_po4_retent = self%o2_min_po4_retent
      o2_min_sed_resp = self%o2_min_sed_resp
      q10_det_rec     = self%q10_det_rec
      q10_h2s         = self%q10_h2s
      q10_nit         = self%q10_nit
      q10_sed_rec     = self%q10_sed_rec
      r_biores        = self%r_biores
      r_cya_assim     = self%r_cya_assim
      r_det_rec       = self%r_det_rec
      r_ips_burial    = self%r_ips_burial
      r_ips_ero       = self%r_ips_ero
      r_ips_liber     = self%r_ips_liber
      r_lpp_assim     = self%r_lpp_assim
      r_nh4_nitrif    = self%r_nh4_nitrif
      r_pp_mort       = self%r_pp_mort
      r_pp_resp       = self%r_pp_resp
      r_sed_ero       = self%r_sed_ero
      r_sed_rec       = self%r_sed_rec
      r_spp_assim     = self%r_spp_assim
      r_zoo_graz      = self%r_zoo_graz
      r_zoo_mort      = self%r_zoo_mort
      r_zoo_resp      = self%r_zoo_resp
      rfr_c           = self%rfr_c
      rfr_h           = self%rfr_h
      rfr_o           = self%rfr_o
      rfr_p           = self%rfr_p
      sali_max_cya    = self%sali_max_cya
      sali_min_cya    = self%sali_min_cya
      sed_max         = self%sed_max
      spp0            = self%spp0
      temp_min_cya    = self%temp_min_cya
      temp_min_spp    = self%temp_min_spp
      temp_opt_zoo    = self%temp_opt_zoo
      w_cya           = self%w_cya
      w_det           = self%w_det
      w_det_sedi      = self%w_det_sedi
      w_ipw           = self%w_ipw
      w_ipw_sedi      = self%w_ipw_sedi
      w_lpp           = self%w_lpp
      w_n2_stf        = self%w_n2_stf
      w_o2_stf        = self%w_o2_stf
      zoo0            = self%zoo0
      zoo_cl          = self%zoo_cl
      stf_no3         = self%stf_no3
      stf_nh4         = self%stf_nh4
      stf_po4         = self%stf_po4

   ! Enter spatial loops over the horizontal domain (if any)._rk
   _HORIZONTAL_LOOP_BEGIN_

   ! Retrieve abiotic parameters
      _GET_HORIZONTAL_(self%id_taub,cgt_current_wave_stress)
      _GET_   (self%id_par,cgt_light) ! local photosynthetically active radiation
      _GET_   (self%id_temp,cgt_temp) !local water temperature
      _GET_   (self%id_salt,cgt_sali) !local water temperature
      cgt_density = 1035.0_rk
      cgt_longitude = 0.0_rk             ! Longitude [degE]
      cgt_latitude  = 0.0_rk             ! Latitude  [degN]
      cgt_timestep  = self%cgt_timestep ! ecosystem model time step [d]


   ! Retrieve current (local) state variable values.
      _GET_(self%id_t_n2           ,t_n2           )
      _GET_(self%id_t_o2           ,t_o2           )

   ! Calculate values of auxiliary variables
       temp_k          = cgt_temp + 273.15_rk
       o2_sat          = (10.18e0_rk+((5.306e-3_rk-4.8725e-5_rk*cgt_temp)*cgt_temp-0.2785e0_rk)*cgt_temp+cgt_sali*((2.2258e-3_rk+(4.39e-7_rk*cgt_temp-4.645e-5_rk)*cgt_temp)*cgt_temp-6.33e-2_rk))*44.66e0_rk*1e-6_rk
       temp1 = log((298.15_rk-cgt_temp)/(273.15_rk+cgt_temp))
       temp2 = temp1*temp1
       temp3 = temp2*temp1
       n2_sat          = 1e-6_rk*exp(6.42931_rk + 2.92704_rk*temp1 + 4.32531_rk*temp2 + 4.69149_rk*temp3 + cgt_sali*(0.0_rk -7.44129e-3_rk - 8.02566e-3_rk*temp1 - 1.46775e-2_rk*temp2))

   ! Calculate process limitation functions
       lim_t_n2_7           = theta(t_n2-0.0_rk)
       lim_t_o2_2           = theta(t_o2-0.0_rk)

   ! Calculate process rates
       p_n2_stf_down   = w_n2_stf*(n2_sat-t_n2)*theta(n2_sat-t_n2)*cgt_density
       p_n2_stf_up     = (w_n2_stf*(t_n2-n2_sat)*theta(t_n2-n2_sat)*cgt_density)*lim_t_n2_7
       p_o2_stf_down   = w_o2_stf*(o2_sat-t_o2)*theta(o2_sat-t_o2)*cgt_density
       p_o2_stf_up     = (w_o2_stf*(t_o2-o2_sat)*theta(t_o2-o2_sat)*cgt_density)*lim_t_o2_2
       p_no3_stf_down  = stf_no3
       p_nh4_stf_down  = stf_nh4
       p_po4_stf_down  = stf_po4

   ! Apply time tendencies to tracers
      _SET_SURFACE_EXCHANGE_(self%id_t_n2           ,0.0_rk                     + p_n2_stf_down               /(1000.0_rk*cgt_density*secs_pr_day)                     - p_n2_stf_up                 /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_SURFACE_EXCHANGE_(self%id_t_o2           ,0.0_rk                     + p_o2_stf_down               /(1000.0_rk*cgt_density*secs_pr_day)                     - p_o2_stf_up                 /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_SURFACE_EXCHANGE_(self%id_t_nh4          ,0.0_rk                     + p_nh4_stf_down              /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_SURFACE_EXCHANGE_(self%id_t_no3          ,0.0_rk                     + p_no3_stf_down              /(1000.0_rk*cgt_density*secs_pr_day)                            )
      _SET_SURFACE_EXCHANGE_(self%id_t_po4          ,0.0_rk                     + p_po4_stf_down              /(1000.0_rk*cgt_density*secs_pr_day)                            )

   ! Leave spatial loops over the horizontal domain (if any)._rk
   _HORIZONTAL_LOOP_END_

   end subroutine do_surface
!EOC


!-----------------------------------------------------------------------

  END MODULE fabm_iow_ergom

!-----------------------------------------------------------------------
! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
!-----------------------------------------------------------------------
