! This model code was generated by the Code Generation Tool CGT
! www.ergom.net
!
! Model Version: <version>, code generation date: <now> 
!----------------------------------------------------------------
! <CONTACT EMAIL="thomas.neumann@io-warnemuende.de"> T. Neumann
! </CONTACT>
! <CONTACT EMAIL="hagen.radtke@io-warnemuende.de"> H. Radtke
! </CONTACT>
! <CONTACT EMAIL="<contact>"> <author> 
! </CONTACT>
!
! <REVIEWER EMAIL="none yet">
! </REVIEWER>
!<OVERVIEW>
! Ecological Regional Ocean Model - ERGOM
!</OVERVIEW>
!
!<DESCRIPTION>
!       ERGOM is developed for marine ecosystems in which bentho-pelagic
!       fluxes and changing redox conditions are important.It describes
!       the nitrogen and phosphorus cycle and partly the sulfor cycle.
!       Included is a carbon module and the oxygen dynamics.
!</DESCRIPTION>
!
! <INFO>
! http://www.ergom.net
! <REFERENCE>
! T. Neumann, 2000: Towards a 3D-ecosystem model of the Baltic Sea
! Journal of Marine Systems, 25, 405-419
! </REFERENCE>
! <REFERENCE>
! T. Neumann, W. Fennel, Ch. Kremp 2002: Experimental Simulations with
! an Ecosystem Model of the Baltic Sea: A Nutrient Load Reduction Experiment,
! Global Biogeochemical Cycles 16, No 3, 7-1 – 7-19
! </REFERENCE>
! <REFERENCE>
! W. Fennel, Th. Neumann, 2004: Introduction to the Modelling of Marine Ecosystems,
! Elsevier Oceanography Series 72, Elsevier
! </REFERENCE>
! <REFERENCE>
! Radtke, H., T. Neumann, M. Voss, and W. Fennel (2012), Modeling pathways of
! riverine nitrogen and phosphorus in the Baltic Sea,
! J. Geophys. Res., 117, C09024, doi:10.1029/2012JC008119
! </REFERENCE>
! <REFERENCE>
! I. Kuznetsov, T. Neumann, Simulation of carbon dynamics in the Baltic Sea with a 3D model,
! Journal of Marine Systems, Available online 22 October 2012,
! ISSN 0924-7963, doi: 10.1016/j.jmarsys.2012.10.011.
! </REFERENCE>
! <DEVELOPER_NOTES>
! Ergom will be continuously developed.
! </DEVELOPER_NOTES>
! </INFO>
!
!----------------------------------------------------------------

#define salt4yellowSubs

module generic_ERGOM

  use coupler_types_mod,   only: coupler_2d_bc_type
  use field_manager_mod,   only: fm_string_len
  use fms_mod,             only: write_version_number, open_namelist_file, check_nml_error, close_file  
  use mpp_mod,             only: mpp_error, NOTE, WARNING, FATAL, stdout, stdlog
  use mpp_mod,             only: mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_ROUTINE
  use time_manager_mod,    only: time_type, get_date, days_in_year, time_type_to_real
  use fm_util_mod,         only: fm_util_start_namelist, fm_util_end_namelist
  use diag_manager_mod,    only: register_diag_field, send_data

  use g_tracer_utils, only : g_tracer_type,g_tracer_start_param_list,g_tracer_end_param_list
  use g_tracer_utils, only : g_tracer_add,g_tracer_add_param, g_tracer_set_files
  use g_tracer_utils, only : g_tracer_set_values,g_tracer_get_pointer,g_tracer_get_common
  use g_tracer_utils, only : g_tracer_coupler_set,g_tracer_coupler_get
  use g_tracer_utils, only : g_tracer_send_diag, g_tracer_get_values
  use g_tracer_utils, only : g_diag_type, g_diag_field_add

  use data_override_mod, only: data_override

  implicit none ; private

  character(len=fm_string_len), parameter :: mod_name       = 'generic_ERGOM'
  character(len=fm_string_len), parameter :: package_name   = 'generic_ergom'

  character(len=128) :: version=&
         '$Id: generic_ERGOM.F90,v <version> <now> $'
  character (len=128) :: tagname = &
         '$Name: IOW $'


  public do_generic_ERGOM
  public generic_ERGOM_register
  public generic_ERGOM_init
  public generic_ERGOM_update_from_coupler
  public generic_ERGOM_update_from_source
  public generic_ERGOM_update_from_bottom
  public generic_ERGOM_set_boundary_values
  public generic_ERGOM_register_diag
  public generic_ERGOM_end

  !
  !This type contains all the parameters and arrays used in this module.
  !
  !Note that there is no programatic reason for treating
  !the following as a type. These are the parameters used only in this module.
  !It suffices for varables to be a declared at the top of the module.
  !nnz: Find out about the timing overhead for using type%x rather than x

  !An auxiliary type for storing varible names
  type, public :: vardesc
     character(len=fm_string_len) :: name     ! variable name in a NetCDF file.
     character(len=fm_string_len) :: longname ! long name of that variable.
     character(len=1)  :: hor_grid            ! hor. grid:  u, v, h, q, or 1.
     character(len=1)  :: z_grid              ! vert. grid:  L, i, or 1.
     character(len=1)  :: t_grid              ! time description: s, a, m, or 1.
     character(len=fm_string_len) :: units    ! dimensions of the variable.
     character(len=1)  :: mem_size            ! size in memory: d or f.
  end type vardesc

  type generic_ERGOM_type

<tracers solubility/=0>
     real, ALLOCATABLE, dimension(:,:) :: &
          <gasName>_<trimName>_csurf               ! surface concentration of <gasName> [mol/m3]
</tracers>
<tracers solubility/=0; childOf=none>
     real, ALLOCATABLE, dimension(:,:) :: &
          <gasName>_alpha, &                       ! surface solubility of <gasName>    [mol/m3/Pa]
          <gasName>_sc_no                          ! Schmidt number of <gasName>        [1]
</tracers>


<processes vertLoc=SED>
     real, ALLOCATABLE, dimension(:,:) :: &
          <name> , &                      ! <description>
          saved_rate_<name>               ! for Runge-Kutta method
</processes>
<processes vertLoc=SUR>
     real, ALLOCATABLE, dimension(:,:) :: &
          <name> , &                      ! <description>
          saved_rate_<name>               ! for Runge-Kutta method
</processes>
<auxiliaries vertLoc=SED>
     real, ALLOCATABLE, dimension(:,:) :: &
          <name>                          ! <description>
</auxiliaries>
<auxiliaries vertLoc=SUR>
     real, ALLOCATABLE, dimension(:,:) :: &
          <name>                          ! <description>
</auxiliaries>
<auxiliaries isZIntegral=1>
     real, ALLOCATABLE, dimension(:,:) :: &
          zintegralarray_<name>           ! place to accumulate <description> over k layers
</auxiliaries>
<tracers vertLoc=FIS>
     real, ALLOCATABLE, dimension(:,:) :: &
          cumulated_change_of_<name>      ! place to accumulate changes in <description> over k layers
</tracers>

     real, ALLOCATABLE, dimension(:,:,:) :: &
<constants dependsOn=xyz>
          <name> , &                      ! <description>
</constants>
<constants dependsOn=xyzt>
          <name> , &                      ! <description>
</constants>
<auxiliaries vertLoc=WAT>
          <name> , &                      ! <description>
</auxiliaries>
<processes vertLoc=WAT>
          <name> , &                      ! <description>
          saved_rate_<name> , &           ! for Runge-Kutta method
</processes>
<celements isAging/=0>
          old_<aged>  , &  ! temporary storage of age concentration of <description>
</celements>
          temp3darray , &                 ! to temporarily store sediment tracers to extract fluffy layer
          vertspeedarray , &                 ! to temporarily store sediment tracers to extract fluffy layer
          irr_inst , &                    ! instantanuous radiation [W/m2]
          irr_zw , &                      ! instantanuous radiation at bottom of cell [W/m2]
          bio_opacity                     ! opacity of biological tracers [1/m]
<tracers vertLoc=WAT>
     real, dimension(:,:,:,:), pointer :: &
          <name>          ! <description>
     real, dimension(:,:), allocatable :: &
          btf_<name>      ! bottom flux of tracer <name> [mol/m2/s] if implicit_bottomflux=.true.
</tracers>
!! TODO: workaround; these tracers are not automatically declared
<celements>
     ! ersterer wird oben als Pointer allokiert
     real, dimension(:,:,:), allocatable   :: &
          sed_<total>       ! total tagged element concentration in the sediment, porewater
</celements>

<tracers vertLoc/=WAT>
     real, dimension(:,:), allocatable :: &
          <name>          ! fluffy layer concentration [mol/m2] of <description>
</tracers>
<tracers vertLoc=WAT; vertSpeed/=0>
     real, dimension(:,:,:), pointer :: &
          vmove_<name> , & ! vertical movement [m/s]
          vdiff_<name>     ! vertical diffusion [m2/s]
</tracers>
<tracers vertLoc=WAT>
     real, dimension(:,:,:), allocatable :: &
          tracer_array_intermediate_<name>, & ! for Runge-Kutta time step
          patankar_modification_<name>        ! for Runge-Kutta Patankar
</tracers>
<tracers vertLoc/=WAT>
     real, dimension(:,:), allocatable ::   &
          tracer_array_intermediate_<name>, & ! for Runge-Kutta time step
          patankar_modification_<name>        ! for Runge-Kutta Patankar
</tracers>

     real    :: gamma_1                  ! constant for opacity contribution of chlorophyll [kg/µg/m]
     real    :: rho_0                    ! water density for sediment module [kg/m3]
     integer :: &
<processes>
          id_<name> , &                  ! diag id for <description>
</processes>
<constants dependsOn=xyz>
          id_<name> , &                  ! diag id for <description>
</constants>
<constants dependsOn=xyzt>
          id_<name> , &                  ! diag id for <description>
</constants>
<auxiliaries>
          id_<name> , &                  ! diag id for <description>
</auxiliaries>
<tracers vertLoc=WAT; atmosDep/=0>
          id_dep_wet_<name>     = -1,  & ! <description> atmospheric deposition
</tracers>
<tracers vertLoc=WAT; riverDep/=0>
          id_runoff_flux_<name> = -1,  & ! <description> runoff flux to the ocean
          id_diffusive_flux_<name> = -1,  & ! <description> diffusive flux from land to the ocean
</tracers>
          id_irr_inst   = -1             ! instantaneous light [W/m2]
     ! SEDIMENT MODEL FIELDS FOLLOW
<tracers vertLoc=WAT; isInPorewater=1>
     real, dimension(:,:,:), allocatable :: &
          sed_<name>                          , & ! porewater concentration of <description>
          tracer_array_intermediate_sed_<name>, & ! for Runge-Kutta time step
          patankar_modification_sed_<name>        ! for Runge-Kutta Patankar
</tracers>
<tracers vertLoc=SED>
     real, dimension(:,:,:), allocatable :: &
          sed_<name>                          , & ! sediment concentration of <description>
          tracer_array_intermediate_sed_<name>, & ! for Runge-Kutta time step
          patankar_modification_sed_<name>        ! for Runge-Kutta Patankar
</tracers>
<auxiliaries>
     real, ALLOCATABLE, dimension(:,:,:) :: &
          sed_<name>            ! array to be allocated if sedimentary values are written to output
</auxiliaries>
<processes vertLoc=WAT; isInPorewater=1>
     real, ALLOCATABLE, dimension(:,:,:) :: &
          sed_<name> , &        ! array to be allocated if sedimentary values are written to output
          saved_rate_sed_<name> ! for Runge-Kutta method
</processes>
<tracers vertLoc=WAT; isInPorewater=1>
     real, ALLOCATABLE, dimension(:,:,:) :: &
          moldiff_<name>        ! molecular diffusivity [m2/s]
</tracers>
<tracers vertLoc=SED>
     real, ALLOCATABLE, dimension(:,:,:) :: &
          moldiff_<name>        ! molecular diffusivity [m2/s]
</tracers>
<processes vertLoc=SED; isInPorewater=1>
     real, ALLOCATABLE, dimension(:,:,:) :: &
          sed_<name> , &        ! array to be allocated if sedimentary values are written to output
          saved_rate_sed_<name> ! for Runge-Kutta method
</processes>
<auxiliaries>
     integer :: id_sed_<name>   ! diag id for sediment values of <description>
</auxiliaries>
<processes vertLoc=WAT; isInPorewater=1>
     integer :: id_sed_<name>   ! diag id for sediment values of <description>
</processes>
<processes vertLoc=SED; isInPorewater=1>
     integer :: id_sed_<name>   ! diag id for sediment values of <description>
</processes>

     real, dimension(:,:,:), allocatable :: & ! sediment model fields follow
          sed_cellheights           , &  ! cell heights in the sediment [m]
          sed_diffusivity_porewater , &  ! bioturbation-generated diffusivity for pore water species [m2/s]
          sed_diffusivity_solids    , &  ! bioturbation-generated diffusivity for solid species [m2/s]
          sed_2d_params             , &  ! sed_inert_deposition[m/s], temperature[degC], salinity[g/kg], bottomdepth[m], fluffy_layer_thickness[m], diffusive layer thickness[m]
          sed_inert_ratio                ! volume fraction of bioinert material in the sediments [1]
     ! END OF SEDIMENT MODEL FIELDS
     character(len=fm_string_len) :: ice_restart_file
     character(len=fm_string_len) :: ocean_restart_file, IC_file
  end type generic_ERGOM_type

  type tracer_2d
    character(len=fm_string_len) ::  name               = 'none'  ! name of the 2d tracer
    character(len=fm_string_len) ::  longname           = 'none'  ! longname of the 2d tracer
    character(len=fm_string_len) ::  units              = 'none'  ! units of the 2d tracer
    character(len=fm_string_len) ::  name_of_3d_tracer  = 'none'  ! name of the 3d tracer
    integer                      ::  layer_in_3d_tracer = 0       ! z-level inside the 3d tracer
    integer                      ::  diag_field         = -1      ! handler returned from register_diag_field
    real, dimension(:,:), pointer::  p_field                      ! pointer to the 2d field stored
    logical                      ::  field_assigned     = .false. ! whether a 2d field has been assigned
  end type tracer_2d

  real, parameter :: missing_value1=-1.0e+20

<constants dependsOn=none>
  real :: <name> = <value> ! <description>
</constants>

  type(generic_ERGOM_type), save   :: ergom
  type(tracer_2d),            ALLOCATABLE, dimension(:), save :: tracers_2d

  !The following logical for using this module is overwritten
  ! by generic_tracer_nml namelist
  logical, save :: do_generic_ERGOM = .false.

  real :: gamma_1 = 0.0150943    ! constant for opacity contribution of chlorophyll [kg/µg/m]
                                 !
                                 ! This is used to estimate a chlorophyll concentration
                                 ! from the bio_opacity of all tracers, assuming all bio_opacity is
                                 ! caused by chlorophyll.
                                 ! This chl concentration is used to provide it for
                                 ! ocean_shortwave_gfdl.F90 which does not want opacity but
                                 ! chl concentration from the biomodel and calculates its own
                                 ! opacity from it

  integer :: NUM_SEDIMENT_LAYERS = 1   ! Default setting is 1 which means fluff layer only.
                                       ! Setting it to a larger value implies the use of a sediment model.
  integer :: NUM_VMOVE_STEPS       = 1
  integer :: RUNGE_KUTTA_DEPTH     = 1
  integer :: SPLITTING_FACTOR_ERGOM= 1
  integer :: n
  integer :: vlev_sed = 2                ! number of 2d tracers that may be stored in one diagnostic 3d tracer

  logical :: genus_style_par = .false.       ! calculate photosynthetically active radiation as in GENUS project
  logical :: implicit_bottomfluxes  = .true. ! tracer concentration changes due to water-sediment interaction are
                                             !   applied outside this module
  logical :: implicit_surfacefluxes = .true. ! tracer concentration changes due to gas exchange in the surface are
                                             !   applied outside this module
  logical :: implicit_movement = .true.      ! vertical migration, sinking and diffusion is applied outside this
                                             !   module

  ! identification numbers for mpp clocks
  integer :: id_source, id_init, id_alloc, id_preloop, id_mainloop, id_vertmig, id_output
  integer :: id_sedsource, id_sedalloc, id_sedpreloop, id_sedmainloop, id_sedvertmig, id_sedoutput

  namelist /ergom_nml/  &
! constants for biomodel
<constants dependsOn=none>
   <name>       , & ! <description>
</constants>
   NUM_SEDIMENT_LAYERS , & ! Default value is 1 which means one fluffy layer only.
                           ! Setting it to a value larger than 1 means a sediment model will be used.
   NUM_VMOVE_STEPS     , &
   RUNGE_KUTTA_DEPTH   , &
   SPLITTING_FACTOR_ERGOM , & ! how many ocean model time steps are one ERGOM time step
   gamma_1                , & ! constant for opacity contribution of chlorophyll [kg/µg/m]
   genus_style_par        , & ! calculate light intensity as it was done in the GENUS model
   implicit_bottomfluxes  , & ! tracer concentration changes due to water-sediment interaction are
                              !   applied outside this module
   implicit_surfacefluxes , & ! tracer concentration changes due to gas exchange in the surface are
                              !   applied outside this module
   implicit_movement      , & ! vertical migration, sinking and diffusion is applied outside this
                              !   module by an implicit method
   vlev_sed

contains

  subroutine generic_ERGOM_register(tracer_list)
    type(g_tracer_type), pointer :: tracer_list

    character(len=fm_string_len), parameter :: sub_name = 'generic_ERGOM_register'
    character(len=fm_string_len) :: errorstring
    integer :: ioun, io_status, ierr, i, j
    logical :: found

    ! provide for namelist over-ride of defaults
    ioun = open_namelist_file()
    read  (ioun, ergom_nml,iostat=io_status)
    write (stdout(),'(/)')
    write (stdout(), ergom_nml)
    write (stdlog(), ergom_nml)
    ierr = check_nml_error(io_status,'ergom_nml')
    call close_file (ioun)

    !Specify all prognostic and diagnostic tracers of this modules.
    call user_add_tracers(tracer_list)

  end subroutine generic_ERGOM_register

  ! <SUBROUTINE NAME="generic_ERGOM_init">
  !  <OVERVIEW>
  !   Initialize the generic ERGOM module
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   This subroutine:
  !       Adds all the CFC Tracers to the list of generic Tracers passed to it
  !       via utility subroutine g_tracer_add().
  !       Adds all the parameters used by this module via utility subroutine g_tracer_add_param().
  !       Allocates all work arrays used in the module.
  !  </DESCRIPTION>
  !  <TEMPLATE>
  !   call generic_ERGOM_init(tracer_list)
  !  </TEMPLATE>
  !  <IN NAME="tracer_list" TYPE="type(g_tracer_type), pointer">
  !   Pointer to the head of generic tracer list.
  !  </IN>
  ! </SUBROUTINE>

  subroutine generic_ERGOM_init(tracer_list)
    type(g_tracer_type), pointer :: tracer_list

    character(len=fm_string_len), parameter :: sub_name = 'generic_ERGOM_init'
    integer :: ioun, io_status, ierr

    id_init     = mpp_clock_id('(ERGOM init) '            ,grain=CLOCK_ROUTINE)
    id_alloc    = mpp_clock_id('(ERGOM allocate) '        ,grain=CLOCK_ROUTINE)
    id_source   = mpp_clock_id('(ERGOM source) '          ,grain=CLOCK_ROUTINE)
    id_preloop  = mpp_clock_id('(ERGOM source pre-loop) ' ,grain=CLOCK_ROUTINE)
    id_mainloop = mpp_clock_id('(ERGOM source main loop) ',grain=CLOCK_ROUTINE)
    id_vertmig  = mpp_clock_id('(ERGOM source vmove) '    ,grain=CLOCK_ROUTINE)
    id_output   = mpp_clock_id('(ERGOM source output) '   ,grain=CLOCK_ROUTINE)

    call mpp_clock_begin(id_init)
    call write_version_number( version, tagname )

    ! provide for namelist over-ride of defaults
    ioun = open_namelist_file()
    read  (ioun, ergom_nml,iostat=io_status)
    write (stdout(),'(/)')
    write (stdout(), ergom_nml)
    write (stdlog(), ergom_nml)
    ierr = check_nml_error(io_status,'ergom_nml')
    call close_file (ioun)

    if (NUM_SEDIMENT_LAYERS .gt. 1) then
       id_sedalloc    = mpp_clock_id('(ERGOM SED allocate) '        ,grain=CLOCK_ROUTINE)
       id_sedsource   = mpp_clock_id('(ERGOM SED source) '          ,grain=CLOCK_ROUTINE)
       id_sedpreloop  = mpp_clock_id('(ERGOM SED source pre-loop) ' ,grain=CLOCK_ROUTINE)
       id_sedmainloop = mpp_clock_id('(ERGOM SED source main loop) ',grain=CLOCK_ROUTINE)
       id_sedvertmig  = mpp_clock_id('(ERGOM SED source vmove) '    ,grain=CLOCK_ROUTINE)
       id_sedoutput   = mpp_clock_id('(ERGOM SED source output) '   ,grain=CLOCK_ROUTINE)
    endif

    !Specify and initialize all parameters used by this package
    call user_add_params

    !Allocate and initiate all the private work arrays used by this module.
    call user_allocate_arrays

    ! now print a summary of all parameters
    write (stdout(),'(/)')
    write (stdout(),*) 'Summary of the ERGOM model setup'
    write (stdout(),*) 'CONSTANTS:'
<constants dependsOn=none>
    write (stdout(),'((a), e13.6)')'    <name>		        : ', <name>
    write (stdout(),'(a)')         '        = <description>'
</constants>
<constants dependsOn=xyz>
    write (stdout(),'((a))')'    <name>		        : depends on x,y,z'
    write (stdout(),'(a)')         '        = <description>'
</constants>
<constants dependsOn=xyzt>
    write (stdout(),'((a))')'    <name>		        : depends on x,y,z,t'
    write (stdout(),'(a)')         '        = <description>'
</constants>
    write (stdout(),'(/)')
    write (stdout(),*) 'TRACERS (3d, no vertical movement):'
<tracers vertLoc=WAT; vertSpeed=0>
    write (stdout(),'(a)')         '    <name>'
    write (stdout(),'((a), e13.6)')'        opacity        [m2/mol]: ',<opacity>
</tracers>
    write (stdout(),'(/)')
    write (stdout(),*) 'TRACERS (3d, vertical movement):'
<tracers vertLoc=WAT; vertSpeed/=0>
    write (stdout(),'(a)')         '    <name>  (<description>)'
    write (stdout(),'((a), e13.6)')'        vertical speed [m /day]: '//'<vertSpeed>=',<vertSpeedValue>
    write (stdout(),'((a), e13.6)')'        opacity        [m2/mol]: ',<opacity>
</tracers>
    write (stdout(),'(/)')
    write (stdout(),*) 'TRACERS (2d, at the surface):'
<tracers vertLoc=SUR>
    write (stdout(),'(a)')         '    <name>  (<description>)'
    write (stdout(),'((a), e13.6)')'        opacity        [m2/mol]: ',<opacity>
</tracers>
    write (stdout(),'(/)')
    write (stdout(),*) 'TRACERS (2d, at the bottom):'
<tracers vertLoc=SED>
    write (stdout(),'(a)')         '    <name>  (<description>)'
</tracers>
    write (stdout(),'(/)')
write (stdout(),*) 'TRACERS (pseudo-3d (fish)):'
<tracers vertLoc=FIS>
    write (stdout(),'(a)')         '    <name>  (<description>)'
</tracers>
    write (stdout(),'(/)')
    write (stdout(),*) 'PROCESSES in the surface layer:'
<processes vertLoc=SUR>
    write (stdout(),'(a)')         '    <name>  (<description>)'
</processes>
    write (stdout(),'(/)')
    write (stdout(),*) 'PROCESSES in the water column:'
<processes vertLoc=WAT>
    write (stdout(),'(a)')         '    <name>  (<description>)'
</processes>
    write (stdout(),'(/)')
    write (stdout(),*) 'PROCESSES in the bottom layer:'
<processes vertLoc=SED>
    write (stdout(),'(a)')         '    <name>  (<description>)'
</processes>
    write (stdout(),'(/)')
    if (implicit_surfacefluxes) then
       write (stdout(),'(a)') 'Using implicit surface fluxes for gas exchange'
       write (stdout(),'(a)') '   TRACERS with solubility and Schmidt number:'
<tracers vertLoc=WAT; solubility/=0>
       write (stdout(),'(a)')    '       <name>'
</tracers>
    else
       write (stdout(),'(a)') 'Using explicit surface fluxes for gas exchange'
       write (stdout(),'(a)') '   (see processes in the surface layer)'
    endif
    if (implicit_movement) then
       write (stdout(),'(a)') 'Using implicit scheme for vertical migration and additional diffusion'
    else
       write (stdout(),'(a)') 'Using explicit scheme for vertical migration and additional diffusion'
       write (stdout(),'(a,I2)') '   Number of vertical movement steps per bio-timestep: ', NUM_VMOVE_STEPS
    endif
    if (implicit_bottomfluxes) then
       write (stdout(),'(a)') 'Using implicit bottom fluxes'
       write (stdout(),'(a)') '   (Update of bottom cell concentration is done via providing bottom fluxes) '
    else
       write (stdout(),'(a)') 'Using explicit bottom fluxes'
       write (stdout(),'(a)') '   (Update of bottom cell concentration is done in this module) '
    endif
    if (genus_style_par) then
       write (stdout(),'(a)') 'Using GENUS model style vertical light extinction'
    else
       write (stdout(),'(a)') 'Using classical ERGOM style vertical light extinction'
    endif
    write (stdout(),'(a,I2)') 'Depth of Runge-Kutta scheme (1=Euler-Forward) : ', RUNGE_KUTTA_DEPTH
    write (stdout(),'(a,I2)') 'split factor (ERGOM timestep/ocean timestep)  :' , SPLITTING_FACTOR_ERGOM
    write (stdout(),'(a,I2)') 'Number of vertical sediment layers            : ', NUM_SEDIMENT_LAYERS
    if (NUM_SEDIMENT_LAYERS .gt. 1) then
       write (stdout(),'(a)') '  => vertically resolved sediment model is used '
       write (stdout(),'(a)') '     TRACERS in pore water: '
<tracers vertLoc=WAT; isInPorewater=1>
    write (stdout(),'(a)')    '       <name>'
</tracers>
    else
       write (stdout(),'(a)') '  => no vertically resolved sediment model is used '
    endif

    call mpp_clock_end(id_init)

  end subroutine generic_ERGOM_init

  subroutine user_allocate_arrays
    integer :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau, n, i
    character(len=fm_string_len) :: mystring

    call mpp_clock_begin(id_alloc)

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau)
    !Allocate all the private arrays.

    if (implicit_surfacefluxes) then
<tracers solubility/=0>
       allocate(ergom%<gasName>_<trimName>_csurf(isd:ied,jsd:jed)); ergom%<gasName>_<trimName>_csurf=0.0
</tracers>
<tracers solubility/=0; childOf=none>
       allocate(ergom%<gasName>_alpha(isd:ied,jsd:jed)); ergom%<gasName>_alpha=0.0
       allocate(ergom%<gasName>_sc_no(isd:ied,jsd:jed)); ergom%<gasName>_sc_no=0.0
</tracers>
    endif

<tracers vertLoc/=WAT>
    allocate(ergom%<name>(isd:ied,jsd:jed)); ergom%<name>=0.0
</tracers>
<tracers vertLoc=FIS>
    call user_2d_tracer_assign_array('<trimName>', ergom%<name>)
</tracers>
    if (NUM_SEDIMENT_LAYERS .eq. 1) then
<tracers vertLoc=SED>
       call user_2d_tracer_assign_array('<trimName>', ergom%<name>)
</tracers>
    else
       allocate(ergom%temp3darray(isd:ied,jsd:jed,nk)); ergom%temp3darray=0.0
       allocate(ergom%vertspeedarray(isd:ied,jsd:jed,nk)); ergom%vertspeedarray=0.0
<tracers vertLoc=WAT; isInPorewater=1>
       allocate(ergom%sed_<name>(isd:ied,jsd:jed,nk)); ergom%sed_<name>=0.0
</tracers>
<tracers vertLoc=WAT; isInPorewater=1>
       allocate(ergom%moldiff_<name>(isd:ied,jsd:jed,nk)); ergom%moldiff_<name>=0.0
</tracers>
<tracers vertLoc=SED>
       allocate(ergom%moldiff_<name>(isd:ied,jsd:jed,nk)); ergom%moldiff_<name>=0.0
</tracers>
<tracers vertLoc=SED>
       allocate(ergom%sed_<name>(isd:ied,jsd:jed,nk)); ergom%sed_<name>=0.0
</tracers>
!! TODO: workaround; these tracers are not automatically allocated
<celements>
       allocate(ergom%sed_<total>(isd:ied,jsd:jed,nk)); ergom%sed_<total>=0.0
</celements>
       allocate(ergom%sed_cellheights          (isd:ied,jsd:jed,nk)); ergom%sed_cellheights=0.0
       allocate(ergom%sed_diffusivity_porewater(isd:ied,jsd:jed,nk)); ergom%sed_diffusivity_porewater=0.0
       allocate(ergom%sed_diffusivity_solids   (isd:ied,jsd:jed,nk)); ergom%sed_diffusivity_solids=0.0
       allocate(ergom%sed_2d_params            (isd:ied,jsd:jed,nk)); ergom%sed_2d_params=0.0
       allocate(ergom%sed_inert_ratio          (isd:ied,jsd:jed,nk)); ergom%sed_inert_ratio=0.0
    endif

    if (implicit_bottomfluxes) then
<tracers vertLoc=WAT>
      allocate(ergom%btf_<name>(isd:ied,jsd:jed)); ergom%btf_<name>=0.0
</tracers>
    endif

<auxiliaries isZIntegral=1>
    allocate(ergom%zintegralarray_<name>(isd:ied,jsd:jed)); ergom%zintegralarray_<name>=0.0
</auxiliaries>
<tracers vertLoc=FIS>
    allocate(ergom%cumulated_change_of_<name>(isd:ied,jsd:jed)); ergom%cumulated_change_of_<name>=0.0
</tracers>
    allocate(ergom%irr_inst(isd:ied,jsd:jed,nk)); ergom%irr_inst=0.0
    allocate(ergom%irr_zw(isd:ied,jsd:jed,nk)); ergom%irr_zw=0.0
    allocate(ergom%bio_opacity(isd:ied,jsd:jed,nk)); ergom%bio_opacity=0.0
<constants dependsOn=xyz>
    allocate(ergom%<name>(isd:ied,jsd:jed,nk)); ergom%<name>=0.0
</constants>
<constants dependsOn=xyzt>
    allocate(ergom%<name>(isd:ied,jsd:jed,nk)); ergom%<name>=0.0
</constants>
<auxiliaries vertLoc=WAT; isUsedElsewhere=1>
    allocate(ergom%<name>(isd:ied,jsd:jed,nk)); ergom%<name>=0.0
</auxiliaries>
<celements isAging/=0>
    allocate(ergom%old_<aged>(isd:ied,jsd:jed,nk)); ergom%old_<aged>=0.0
</celements>

    !allocations only for Runge-Kutta method
    if (RUNGE_KUTTA_DEPTH .eq. 2) then
<tracers vertLoc=WAT>
       allocate(ergom%tracer_array_intermediate_<name>(isd:ied,jsd:jed,nk)) ; ergom%tracer_array_intermediate_<name>=0.0
       allocate(ergom%patankar_modification_<name>(isd:ied,jsd:jed,nk))     ; ergom%patankar_modification_<name>    =0.0
</tracers>
<tracers vertLoc/=WAT>
       allocate(ergom%tracer_array_intermediate_<name>(isd:ied,jsd:jed))    ; ergom%tracer_array_intermediate_<name>=0.0
       allocate(ergom%patankar_modification_<name>(isd:ied,jsd:jed))        ; ergom%patankar_modification_<name>    =0.0
</tracers>
<processes vertLoc=WAT>
       allocate(ergom%saved_rate_<name>(isd:ied,jsd:jed,nk)) ; ergom%saved_rate_<name>=0.0
</processes>
<processes vertLoc/=WAT>
       allocate(ergom%saved_rate_<name>(isd:ied,jsd:jed)) ;    ergom%saved_rate_<name>=0.0
</processes>
       if (NUM_SEDIMENT_LAYERS .gt. 1) then
<tracers vertLoc=WAT; isInPorewater=1>
          allocate(ergom%tracer_array_intermediate_sed_<name>(isd:ied,jsd:jed,nk)) ; ergom%tracer_array_intermediate_sed_<name>=0.0
          allocate(ergom%patankar_modification_sed_<name>(isd:ied,jsd:jed,nk))     ; ergom%patankar_modification_sed_<name>    =0.0
</tracers>
<tracers vertLoc=SED>
          allocate(ergom%tracer_array_intermediate_sed_<name>(isd:ied,jsd:jed,nk)) ; ergom%tracer_array_intermediate_sed_<name>=0.0
          allocate(ergom%patankar_modification_sed_<name>(isd:ied,jsd:jed,nk))     ; ergom%patankar_modification_sed_<name>    =0.0
</tracers>
<processes vertLoc=WAT; isInPorewater=1>
          allocate(ergom%saved_rate_sed_<name>(isd:ied,jsd:jed,nk)) ; ergom%saved_rate_sed_<name>=0.0
</processes>
<processes vertLoc=SED; isInPorewater=1>
          allocate(ergom%saved_rate_sed_<name>(isd:ied,jsd:jed,nk)) ; ergom%saved_rate_sed_<name>=0.0
</processes>
       endif
    endif

    call mpp_clock_end(id_alloc)

  end subroutine user_allocate_arrays

  subroutine generic_ERGOM_register_diag(diag_list)
    type(g_diag_type), pointer :: diag_list
    type(vardesc)  :: vardesc_temp
    integer        :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau, axes(3)
    type(time_type):: init_time

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,axes=axes,init_time=init_time)

<processes vertLoc=WAT>
    vardesc_temp = vardesc("<name>","<description>",'h','1','s','mol kg-1 d-1','f')
    ergom%id_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:3),&
         init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
    if (ergom%id_<name> .gt. 0) then
      allocate(ergom%<name>(isd:ied,jsd:jed,nk))
      ergom%<name> = 0.0
    endif
</processes>

<processes vertLoc/=WAT>
    vardesc_temp = vardesc("<name>","<description>",'h','1','s','mol m-2 d-1','f')
    ergom%id_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:2),&
         init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
    if (ergom%id_<name> .gt. 0) then
      allocate(ergom%<name>(isd:ied,jsd:jed))
      ergom%<name> = 0.0
    endif
</processes>

    if (NUM_SEDIMENT_LAYERS .gt. 1) then
<processes vertLoc=WAT; isInPorewater=1>
       vardesc_temp = vardesc("sed_<name>","<description>",'h','1','s','mol kg-1 d-1','f')
       ergom%id_sed_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:3),&
            init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
       if (ergom%id_sed_<name> .gt. 0) then
         allocate(ergom%sed_<name>(isd:ied,jsd:jed,nk))
         ergom%<name> = 0.0
       endif
</processes>

<processes vertLoc=SED; isInPorewater=1>
       vardesc_temp = vardesc("sed_<name>","<description>",'h','1','s','mol m-2 d-1','f')
       ergom%id_sed_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:3),&
            init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
       if (ergom%id_<name> .gt. 0) then
         allocate(ergom%sed_<name>(isd:ied,jsd:jed,nk))
         ergom%<name> = 0.0
       endif
</processes>
    endif

<auxiliaries vertLoc=WAT; isUsedElsewhere=0>
    vardesc_temp = vardesc("<name>","<description>",'h','1','s','1','f')
    ergom%id_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:3),&
         init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
    if (ergom%id_<name> .gt. 0) then
      allocate(ergom%<name>(isd:ied,jsd:jed,nk))
      ergom%<name> = 0.0
    endif
</auxiliaries>

<auxiliaries vertLoc/=WAT>
    vardesc_temp = vardesc("<name>","<description>",'h','1','s','1','f')
    ergom%id_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:2),&
         init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
    if (ergom%id_<name> .gt. 0) then
      allocate(ergom%<name>(isd:ied,jsd:jed))
      ergom%<name> = 0.0
    endif
</auxiliaries>

    if (NUM_SEDIMENT_LAYERS .gt. 1) then
<auxiliaries vertLoc=WAT; isUsedElsewhere=0>
       vardesc_temp = vardesc("sed_<name>","<description>",'h','1','s','1','f')
       ergom%id_sed_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:3),&
            init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
       if (ergom%id_sed_<name> .gt. 0) then
         allocate(ergom%sed_<name>(isd:ied,jsd:jed,nk))
         ergom%sed_<name> = 0.0
       endif
</auxiliaries>

<auxiliaries vertLoc/=WAT>
       vardesc_temp = vardesc("sed_<name>","<description>",'h','1','s','1','f')
       ergom%id_sed_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:3),&
            init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
       if (ergom%id_sed_<name> .gt. 0) then
         allocate(ergom%sed_<name>(isd:ied,jsd:jed,nk))
         ergom%sed_<name> = 0.0
       endif
</auxiliaries>
    endif

<tracers vertLoc=WAT; atmosDep/=0>
    vardesc_temp = vardesc("dep_wet_<name>","Atmospheric deposition of <description> to the ocean",'h','1','s','mol m-2 s-1','f')
    ergom%id_dep_wet_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:2),&
         init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)

</tracers>
<tracers vertLoc=WAT; riverDep/=0>
    vardesc_temp = vardesc("runoff_flux_<name>","<description> runoff flux to the ocean",'h','1','s','mol m-2 s-1','f')
    ergom%id_runoff_flux_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:2),&
         init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
    vardesc_temp = vardesc("diffusive_flux_<name>","<description> diffusive flux from land to ocean",'h','1','s','mol m-2 s-1','f')
    ergom%id_diffusive_flux_<name> = register_diag_field(package_name, vardesc_temp%name, axes(1:2),&
         init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)

</tracers>

    vardesc_temp = vardesc("irr_inst","Instantaneous Light",'h','L','s','W m-2','f')
    ergom%id_irr_inst = register_diag_field(package_name, vardesc_temp%name, axes(1:3),&
         init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)

    call user_register_2d_tracers

  end subroutine generic_ERGOM_register_diag

  !
  !   This is an internal sub, not a public interface.
  !   Add all the parameters to be used in this module.
  !
  subroutine user_add_params

  !Specify all parameters used in this modules.
  !=============================================================================

  !Add the known experimental parameters used for calculations
  !in this module.
  !All the g_tracer_add_param calls must happen between
  !g_tracer_start_param_list and g_tracer_end_param_list  calls.
  !This implementation enables runtime overwrite via field_table.

    integer                      :: n
    character(len=fm_string_len) :: mystring
    call g_tracer_start_param_list(package_name)

    call g_tracer_add_param('gamma_1'  , ergom%gamma_1       , gamma_1   )
    call g_tracer_add_param('RHO_0'    , ergom%rho_0         , 1035.0   )

    call g_tracer_end_param_list(package_name)

  end subroutine user_add_params

  !
  !   This is an internal sub, not a public interface.
  !   Add all the tracers to be used in this module.
  !
  subroutine user_add_tracers(tracer_list)
    type(g_tracer_type), pointer :: tracer_list

    character(len=fm_string_len), parameter :: sub_name = 'user_add_tracers'
    integer                                 :: i, n
    character(len=fm_string_len)            :: mystring

    call g_tracer_start_param_list(package_name)
    call g_tracer_add_param('ice_restart_file'   , ergom%ice_restart_file   , 'ice_ergom.res.nc')
    call g_tracer_add_param('ocean_restart_file' , ergom%ocean_restart_file , 'ocean_ergom.res.nc' )
    call g_tracer_add_param('IC_file'            , ergom%IC_file       , '')
    call g_tracer_end_param_list(package_name)

    ! Set Restart files
    call g_tracer_set_files(ice_restart_file=ergom%ice_restart_file, ocean_restart_file=ergom%ocean_restart_file )

    !=====================================================
    !Specify all prognostic tracers of this modules.
    !=====================================================
    !User adds one call for each prognostic tracer below!
    !User should specify if fluxes must be extracted from boundary
    !by passing one or more of the following methods as .true.
    !and provide the corresponding parameters array
    !methods: flux_gas,flux_runoff,flux_wetdep,flux_drydep
    !
    !prog_tracers: nitrogen
    !diag_tracers: none
    !

    call user_init_2d_tracer_list

<backwardTracers vertLoc=WAT; vertSpeed/=0>
    call g_tracer_add(tracer_list,package_name,&
       name       = '<name>',          &
       longname   = '<description>',   &
  <if useInitValue/=0>
       const_init_value = <initValue>, &
  </if>
       flux_runoff= (<riverDep> > 0),  &
  <if solubility=0>
         flux_gas   = .false.,         &
  </if>
  <if solubility/=0>
         flux_gas   = implicit_surfacefluxes,        &
         flux_gas_name  = '<trimname>_flux',          &
         flux_gas_type  = 'air_sea_gas_flux_generic',&
         flux_gas_param = (/ 6.67e-07, 1.0 /),       &
         flux_gas_restart_file  = 'ocean_ergom_airsea_flux.res.nc', &
         flux_gas_molwt = <molarMass>, &
  </if>
       flux_wetdep= (<atmosDep> > 0),  &
       flux_drydep= (<riverDep> > 0),  &
       flux_param = (/ 1.0 /),         &
       units      = 'mol/kg',          &
       prog       = .true. ,           &
       flux_bottom= implicit_bottomfluxes , &
       move_vertical = .true. ,        &
       diff_vertical = .true.)
</tracers>
<backwardTracers vertLoc=WAT; vertSpeed=0>
  <if solubility/=0>
    ! a comment to "flux_gas_param":
    ! The first term of flux_gas_param [s/m] describes the relation between (10m wind speed)**2 and gas transfer velocity:
    ! flux       = (csat-c) * sqrt(660/schmidt_number) * param1 * u10**2
    ! [mol/m2/s] = [mol/m3] * [1]                      * [s/m]  * [m2/s2]
    ! param1 = 8.61e-7 according to eqn. 3 in: Wannikhof (1992): Relationship between Wind Speed and Gas Exchange over the Ocean, JGR, 1997,
    ! Wannikhof (1992): Relationship between Wind Speed and Gas Exchange over the Ocean, JGR, Vol. 97, page 7376
    ! more recent Wanninkhof estimetes propose piston velocity by 0.24*U10*U10 
  </if>
    call g_tracer_add(tracer_list,package_name,&
         name       = '<name>',    &
         longname   = '<description>',      &
  <if useInitValue/=0>
    <if initValue=0.0>
         const_init_value = 1.0e-20,     &
    </if>
    <if initValue/=0.0>
         const_init_value = <initValue>, &
    </if>
  </if>
         units      = 'mol/kg',          &
         prog       = .true.,            &
  <if solubility=0>
         flux_gas   = .false.,           &
  </if>
  <if solubility/=0>
         flux_gas   = implicit_surfacefluxes,        &
         flux_gas_name  = '<trimname>_flux',          &
         flux_gas_type  = 'air_sea_gas_flux_generic',&
         flux_gas_param = (/ 6.67e-07, 1.0 /),       &
         flux_gas_restart_file  = 'ocean_ergom_airsea_flux.res.nc', &
         flux_gas_molwt = <molarMass>, &
  </if>
         flux_bottom= implicit_bottomfluxes,           &
         flux_runoff= (<riverDep> > 0),  &
         flux_param = (/ 1.0 /),         &
         flux_wetdep= (<atmosDep> > 0),  &
         flux_drydep= (<riverDep> > 0),  &
         btm_reservoir=.false. )
</backwardTracers>

<tracers vertLoc=FIS>
    call user_add_2d_tracer(tracer_list,    &
                            name       = '<trimName>',  &
                            longname   = '<description>', &
                            units      = 'mol/m^2')
</tracers>
    if (NUM_SEDIMENT_LAYERS .eq. 1) then
       ! only one sediment layer exists (=fluffy layer), so try to store several sediment tracers in one
       ! 3-d diagnostic tracer
<tracers vertLoc=SED>
       call user_add_2d_tracer(tracer_list,    &
                               name       = '<trimName>',  &
                               longname   = '<description>', &
                               units      = 'mol/m^2')
</tracers>
    else
       ! several sediment layers exist (1=fluffy layer, 2..NUM_SEDIMENT_LAYERS=inside the sediment),
       ! so every sediment tracer gets its own 3-d diagnostic tracer
<tracers vertLoc=SED>
       call g_tracer_add(tracer_list,package_name,&
                         name       = '<name>', &
                         longname   = '<description>', &
                         units      = 'mol m-2',     &
                         prog       = .false.                )
</tracers>
<tracers vertLoc=WAT; isInPorewater=1>
       call g_tracer_add(tracer_list,package_name,&
                         name       = 'sed_<name>', &
                         longname   = 'pore water concentration of <description>', &
                         units      = 'mol kg-1',     &
                         prog       = .false.                )
</tracers>
       ! In this case, we need additional information about the sediment
       call g_tracer_add(tracer_list,package_name,&
                         name       = 'sed_cellheights', &
                         longname   = 'cell heights in the sediment', &
                         units      = 'm',     &
                         prog       = .false.                )
       call g_tracer_add(tracer_list,package_name,&
                         name       = 'sed_diffusivity_porewater', &
                         longname   = 'bioturbation-generated diffusivity for pore water species', &
                         units      = 'm2/s',     &
                         prog       = .false.                )
       call g_tracer_add(tracer_list,package_name,&
                         name       = 'sed_diffusivity_solids', &
                         longname   = 'bioturbation-generated diffusivity for solid species', &
                         units      = 'm2/s',     &
                         prog       = .false.                )
       call g_tracer_add(tracer_list,package_name,&
                         name       = 'sed_2d_params', &
                         longname   = 'sed_inert_deposition, temperature, salinity, bottomdepth, fluffy_layer_thickness, diffusive_layer_thickness', &
                         units      = 'm/s',     &
                         prog       = .false.            )
       call g_tracer_add(tracer_list,package_name,&
                         name       = 'sed_inert_ratio', &
                         longname   = 'volume fraction of bioinert material in the sediments', &
                         units      = '1.0',     &
                         prog       = .false.            )
    endif

<constants dependsOn=xyz>
   call g_tracer_add(tracer_list,package_name,&
         name       = '<name>', &
         longname   = '<description>', &
         units      = '1',     &
         prog       = .false.,       &
         const_init_value = 0.0          )

</constants>
<constants dependsOn=xyzt>
   call g_tracer_add(tracer_list,package_name,&
         name       = '<name>', &
         longname   = '<description>', &
         units      = '1',     &
         prog       = .false.,       &
         const_init_value = 0.0          )

</constants>

<auxiliaries vertLoc=WAT; isUsedElsewhere=1>
   call g_tracer_add(tracer_list,package_name,&
         name       = '<name>', &
         longname   = '<description>', &
         units      = '1',     &
         prog       = .false.,       &
         const_init_value = 0.0          )

</auxiliaries>

   call g_tracer_add(tracer_list,package_name,&
         name       = 'chl',         &
         longname   = 'Chlorophyll', &
         units      = 'ug kg-1',     &
         prog       = .false.,       &
         const_init_value = 0.08           )

    call g_tracer_add(tracer_list,package_name,&
         name       = 'opacity_bio', &
         longname   = 'Opacity of water ingredients', &
         units      = 'm-1',     &
         prog       = .false.,       &
         const_init_value = 1.0e-20           )

    call g_tracer_add(tracer_list,package_name,&
         name       = 'opacity_water', &
         longname   = 'Opacity of pure water', &
         units      = 'm-1',     &
         prog       = .false.,       &
         const_init_value = gamma0           )

  end subroutine user_add_tracers

  !
  !   This is an internal sub, not a public interface.
  !   initialize the list of 2d tracers with a length of zero
  !
  subroutine user_init_2d_tracer_list
    character(len=fm_string_len), parameter :: sub_name = 'user_init_2d_tracer_list'
    allocate(tracers_2d(0))
  end subroutine user_init_2d_tracer_list
  !
  !   This is an internal sub, not a public interface.
  !   Insert a 2d tracer into the 2d tracer list
  !   2d tracers are stored  together in diagnostic 3d tracers.
  !   Those are registered for output and in the restart list.
  !
  subroutine user_add_2d_tracer(tracer_list,name,longname,units)
    type(g_tracer_type), pointer :: tracer_list
    character(len=*), intent(in) :: name, longname, units

    character(len=fm_string_len), parameter :: sub_name = 'user_add_2d_tracer'
    integer                      :: m, n, dummy
    character(len=fm_string_len) :: temp_string
    type(tracer_2d), ALLOCATABLE, dimension(:) :: temp_tracers_2d

    n  = size(tracers_2d)
    allocate(temp_tracers_2d(n))
    do m=1,n
      temp_tracers_2d(m)=tracers_2d(m)
    end do
    deallocate(tracers_2d)
    allocate(tracers_2d(n+1))
    do m=1,n
      tracers_2d(m)=temp_tracers_2d(m)
    end do
    deallocate(temp_tracers_2d)

    tracers_2d(n+1)%name     = trim(adjustl(name))
    tracers_2d(n+1)%longname = trim(adjustl(longname))
    tracers_2d(n+1)%units    = trim(adjustl(units))
    tracers_2d(n+1)%layer_in_3d_tracer = mod(n,vlev_sed)+1
    tracers_2d(n+1)%field_assigned     = .false.  ! whether a 2d field has been assigned

    m = n/vlev_sed+1
    write( temp_string, '(i4)' )  m; temp_string = adjustl(temp_string)
    tracers_2d(n+1)%name_of_3d_tracer = 'tracer_2d_'//trim(temp_string)
    if (mod(n,vlev_sed) .eq. 0) then
      call g_tracer_add(tracer_list,package_name,                              &
         name       = 'tracer_2d_'//trim(temp_string),                         &
         longname   = 'Tracer '//trim(temp_string)//' containing 2d variables',&
         units      = 'none',                                                  &
         prog       = .false.)
    endif
  end subroutine user_add_2d_tracer

  !   This is an internal sub, not a public interface.
  !   Register the 2d tracer for output via a 2d output field
  subroutine user_register_2d_tracers
    real,parameter :: missing_value1=-1.0e+20
    type(vardesc)  :: vardesc_temp
    integer        :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau, axes(3)
    type(time_type):: init_time
    integer        :: m,n

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,axes=axes,init_time=init_time)

    n=size(tracers_2d)
    do m=1,n
      vardesc_temp = vardesc(                                                                   &
      tracers_2d(m)%name, tracers_2d(m)%longname,'h','L','s',tracers_2d(m)%units,'f')
      tracers_2d(m)%diag_field = register_diag_field(package_name, vardesc_temp%name, axes(1:2),&
	      init_time, vardesc_temp%longname,vardesc_temp%units, missing_value = missing_value1)
    end do
  end subroutine user_register_2d_tracers

  !   This is an internal sub, not a public interface.
  !   Let the p_field pointer of the 2d tracer point to the data array (for read and write)
  subroutine user_2d_tracer_assign_array(name,array)
    character(len=*), intent(in)           :: name
    real,dimension(:,:),target             :: array
    integer                                :: n
    logical                                :: found_tracer
    character(len=fm_string_len)           :: errorstring

    found_tracer=.false.
    do n=1,size(tracers_2d)
      if (trim(adjustl(name)) .eq. trim(adjustl(tracers_2d(n)%name))) then
        tracers_2d(n)%p_field => array
        tracers_2d(n)%field_assigned = .true.
        found_tracer=.true.
      end if
    end do
    if (.not. found_tracer) then
       write(errorstring, '(a)') &
       'array assigned to tracer '//trim(name)// &
       ', but that tracer was not added by user_add_2d_tracer'
       call  mpp_error(FATAL, errorstring)
    end if
  end subroutine user_2d_tracer_assign_array

  !   This is an internal sub, not a public interface.
  !   Load all data stored in the (3d tracers containing 2d tracer data) into their temporary 2d arrays
   subroutine user_get_2d_tracer_values(tracer_list,isd,ied,jsd,jed,nk)
    type(g_tracer_type),    pointer      :: tracer_list
    integer,                  intent(in) :: isd,ied,jsd,jed,nk
    real,dimension(:,:,:),allocatable    :: a3d
    integer                              :: n
    character(len=fm_string_len)         :: loaded_3d_tracer
    allocate(a3d(isd:ied,jsd:jed,1:nk))
    loaded_3d_tracer=''
    do n=1,size(tracers_2d)
      if (tracers_2d(n)%field_assigned) then
        if (trim(adjustl(tracers_2d(n)%name_of_3d_tracer)) .ne. trim(adjustl(loaded_3d_tracer))) then
          call g_tracer_get_values(     &
	  tracer_list,tracers_2d(n)%name_of_3d_tracer,'field',a3d,isd,jsd,ntau=1,positive=.true.)
          loaded_3d_tracer=tracers_2d(n)%name_of_3d_tracer
        end if
        tracers_2d(n)%p_field = a3d(:,:,tracers_2d(n)%layer_in_3d_tracer)
      end if
    end do
    deallocate(a3d)
  end subroutine user_get_2d_tracer_values

  !   This is an internal sub, not a public interface.
  !   Save all 2d tracer data from their temporary 2d arrays into the 3d tracers containing 2d tracer data
  subroutine user_set_2d_tracer_values(tracer_list ,model_time)
    type(g_tracer_type),      pointer         :: tracer_list
    type(time_type),          intent(in)      :: model_time
    real, dimension(:,:,:),   pointer         :: grid_tmask

    integer                                   :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau
    real,dimension(:,:,:),allocatable         :: a3d
    integer                                   :: n, layer
    character(len=fm_string_len)              :: loaded_3d_tracer
    logical                                   :: used

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,grid_tmask=grid_tmask)

    allocate(a3d(isd:ied,jsd:jed,1:nk))
    loaded_3d_tracer=''
    do n=1,size(tracers_2d)
      if (tracers_2d(n)%field_assigned) then
        if (trim(adjustl(tracers_2d(n)%name_of_3d_tracer)) .ne. trim(adjustl(loaded_3d_tracer))) then
          ! a new 3d tracer starts -> first save the temporarily stored data into the old 3d tracer
          if (loaded_3d_tracer .ne. '') then
            call g_tracer_set_values(tracer_list,loaded_3d_tracer,'field',a3d,isd,jsd,ntau=1)
          end if
          ! then, nullify the temporary tracer
          a3d=0.0
          loaded_3d_tracer=tracers_2d(n)%name_of_3d_tracer
        end if
        ! save the data to the temporary tracer
        layer = tracers_2d(n)%layer_in_3d_tracer
        a3d(:,:,layer)  = tracers_2d(n)%p_field
        ! save the values to the diagnostic field
        used = send_data(tracers_2d(n)%diag_field, a3d(:,:,layer),      &
                  model_time, rmask = grid_tmask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
      end if
    end do
    !now, save the last tracer
    if (loaded_3d_tracer .ne. '') then
      call g_tracer_set_values(tracer_list,loaded_3d_tracer,'field',a3d,isd,jsd,ntau=1)
    end if
    deallocate(a3d)
  end subroutine user_set_2d_tracer_values

  ! <SUBROUTINE NAME="generic_ERGOM_update_from_coupler">
  !  <OVERVIEW>
  !   Modify the values obtained from the coupler if necessary.
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Currently an empty stub for CFCs.
  !   Some tracer fields need to be modified after values are obtained from the coupler.
  !   This subroutine is the place for specific tracer manipulations.
  !  </DESCRIPTION>
  !  <TEMPLATE>
  !   call generic_ERGOM_update_from_coupler(tracer_list)
  !  </TEMPLATE>
  !  <IN NAME="tracer_list" TYPE="type(g_tracer_type), pointer">
  !   Pointer to the head of generic tracer list.
  !  </IN>
  ! </SUBROUTINE>
  subroutine generic_ERGOM_update_from_coupler(tracer_list)
    type(g_tracer_type), pointer :: tracer_list
    character(len=fm_string_len), parameter :: sub_name = 'generic_ERGOM_update_from_coupler'
    !
    ! Nothing to be done for ERGOM
    !
    return
  end subroutine generic_ERGOM_update_from_coupler

  ! <SUBROUTINE NAME="generic_ERGOM_update_from_bottom">
  !  <OVERVIEW>
  !   Set values of bottom fluxes and reservoirs
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Some tracers have bottom fluxes and reservoirs.
  !   This subroutine is the place for specific tracer manipulations.
  !  </DESCRIPTION>
  !  <TEMPLATE>
  !   call generic_ERGOM_update_from_bottom(tracer_list,dt, tau, model_time)
  !  </TEMPLATE>
  !  <IN NAME="tracer_list" TYPE="type(g_tracer_type), pointer">
  !   Pointer to the head of generic tracer list.
  !  </IN>
  !  <IN NAME="dt" TYPE="real">
  !   Time step increment
  !  </IN>
  !  <IN NAME="tau" TYPE="integer">
  !   Time step index to be used for %field
  !  </IN>
  ! </SUBROUTINE>
  subroutine generic_ERGOM_update_from_bottom(tracer_list, dt, tau, model_time)
    type(g_tracer_type), pointer :: tracer_list
    real,               intent(in) :: dt
    integer,            intent(in) :: tau
    type(time_type),    intent(in) :: model_time
    integer :: isc,iec, jsc,jec,isd,ied,jsd,jed,nk,ntau
    logical :: used
    real, dimension(:,:,:),pointer :: grid_tmask
    real, dimension(:,:,:,:),ALLOCATABLE :: temp_field
    !
    ! Nothing to be done for ERGOM
    !
   end subroutine generic_ERGOM_update_from_bottom

   real function theta(x)
      real, intent(in) :: x
      if (x > 0) then
         theta=1.0
      else
         theta=0.0
      endif
      return
   end function theta

   real function power(x,y)
      real, intent(in) :: x, y
      power = x**y
      return
   end function power


  ! <SUBROUTINE NAME="cgt_bio_timestep">
  !  <OVERVIEW>
  !   Update tracer concentration fields due to the source/sink contributions.
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Ecosystem model code automatically generated by Code Generation Tool - see www.ergom.net
  !  </DESCRIPTION>
  ! </SUBROUTINE>
  subroutine cgt_bio_timestep(tracer_list,Temp,Salt, rho_dzt,dzt,xt,yt,hblt_depth,&
       ilb,jlb,tau,dt,grid_dat,model_time,nbands,max_wavelength_band, sw_pen_band,opacity_band, &
       current_wave_stress,use_saved_rates,intermediate,diff_cbt)

    type(g_tracer_type),            pointer    :: tracer_list
    real, dimension(ilb:,jlb:,:),   intent(in) :: Temp,Salt,rho_dzt,dzt
    real, dimension(ilb:,jlb:),     intent(in) :: xt,yt
    real, dimension(ilb:,jlb:),     intent(in) :: hblt_depth
    integer,                        intent(in) :: ilb,jlb,tau
    real,                           intent(in) :: dt
    real, dimension(ilb:,jlb:),     intent(in) :: grid_dat
    type(time_type),                intent(in) :: model_time

    integer,                        intent(in) :: nbands
    real, dimension(:),             intent(in) :: max_wavelength_band
    real, dimension(:,ilb:,jlb:),   intent(in) :: sw_pen_band
    real, dimension(:,ilb:,jlb:,:), intent(in) :: opacity_band
    real, dimension(ilb:,jlb:),     intent(in) :: current_wave_stress
    logical,                        intent(in) :: use_saved_rates
    logical,                        intent(in) :: intermediate
    real, dimension(ilb:,jlb:,:,:),optional,intent(in) :: diff_cbt

    ! abiotic parameters for biological processes
    real :: cgt_temp                ! potential temperature     [Celsius]
    real :: cgt_sali                ! salinity                  [g/kg]
    real :: cgt_light               ! downward light flux (PAR) [W/m2]
    real :: cgt_cellheight          ! cell height               [m]
    real :: cgt_density             ! density                   [kg/m3]
    real :: cgt_bottomdepth         ! bottom depth              [m]
    real :: cgt_timestep            ! timestep in days          [days]
    real :: cgt_current_wave_stress ! combined bottom stress of waves and current [N/m2]
    real :: cgt_longitude           ! longitude of grid cell    [deg]
    real :: cgt_latitude            ! latitude of grid cell     [deg]
    real :: cgt_diffusivity         ! diffusivity at cell bottom for T [m**2/s]
    integer :: cgt_iteration        ! current number of iteration in the iterative loop
    integer :: cgt_year, cgt_dayofyear, cgt_mymonth, cgt_myday, cgt_myhour, cgt_myminute, cgt_mysecond ! some date/time variables
    real :: cgt_hour                ! fractional hour (0..23.9999)
    real :: cgt_in_sediment         ! whether we are below the sediment surface

    character(len=fm_string_len), parameter :: sub_name = 'cgt_bio_timestep'
    integer :: isc,iec, jsc,jec,isd,ied,jsd,jed,nk,ntau, i, j, k , kblt, n, m, nband
    real    :: light_itemp, light_ntemp
    real, dimension(:,:,:) ,pointer :: grid_tmask
    integer, dimension(:,:),pointer :: mask_coast,grid_kmt
    logical                         :: used

    integer                         :: number_of_loop

    real                            :: temp1
    real                            :: temp2
    real                            :: temp3
    real                            :: temp4
    real                            :: temp5
    real                            :: temp6
    real                            :: temp7
    real                            :: temp8
    real                            :: temp9
    real,  allocatable, dimension(:,:) :: cgt_bottomdepth_temp
    real,  allocatable, dimension(:,:) :: light1, light2
<tracers>
    real                            :: <name>   ! <description>
  <limitations>
    real                            :: <name>
  </limitations>
</tracers>

<constants dependsOn=xyz>
    real                            :: <name>   ! <description>
</constants>
<constants dependsOn=xyzt>
    real                            :: <name>   ! <description>
</constants>

<auxiliaries>
    real                            :: <name>   ! <description>
</auxiliaries>

<processes>
    real                            :: <name>   ! <description>
</processes>

<tracers vertLoc=WAT>
    real                            :: above_<name>   ! value one grid cell above
</tracers>
<tracers vertLoc=FIS>
    real                            :: cumulated_change_of_<name>
</tracers>
<processes isStiff/=0>
    real                            :: after_patankar_<name>
</processes>

<tracers vertLoc=WAT; riverDep/=0>
    real, dimension(:,:)   ,pointer :: runoff_flux_<name>
    real, dimension(:,:)   ,pointer :: diffusive_flux_<name>
</tracers>
<tracers vertLoc=WAT; atmosDep/=0>
    real, dimension(:,:)   ,pointer :: wet_<name>
</tracers>

    ! the following variables are used for the positive-definite scheme
    real :: timestep_fraction           ! ratio between allowed timestep and attempted timestep
    real :: timestep_fraction_new       ! ratio between timestep allowed by each tracer and attempted timestep
    real :: fraction_of_total_timestep  ! ratio between remaining timestep and original timestep
<tracers>
    real :: change_of_<name>            ! possible change during remaining timestep
    real :: change_btf_<name>           ! possible bottom tracer flux during remaining timestep
</tracers>
<processes>
    real :: total_rate_<name>           ! total process rate accumulated during all sub-timesteps
</processes>
    integer :: which_tracer_exhausted   ! stores the number of the exhausted tracer

    call mpp_clock_begin(id_source)
    call mpp_clock_begin(id_preloop)

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,&
         grid_tmask=grid_tmask,grid_mask_coast=mask_coast,grid_kmt=grid_kmt)

    !Get all 2d tracer values
    call user_get_2d_tracer_values(tracer_list,isd,ied,jsd,jed,nk)
    if (NUM_SEDIMENT_LAYERS .gt. 1) then
       !Get values for fluffy layer
<tracers vertLoc=SED>
       call g_tracer_get_values(tracer_list,'<name>'  ,'field',ergom%temp3darray,isd,jsd,ntau=1)
       ergom%<name> = ergom%temp3darray(:,:,1)
</tracers>
    endif

    !Get the pointers to the prognostic tracers after applying the physics.
<tracers vertLoc=WAT>
    call g_tracer_get_pointer(tracer_list,'<name>'  ,'field',ergom%<name>    )
</tracers>
<constants dependsOn=xyz>
    call g_tracer_get_values(tracer_list,'<name>'  ,'field',ergom%<name>,isd,jsd,ntau=1)
</constants>
<constants dependsOn=xyzt>
    call g_tracer_get_values(tracer_list,'<name>'  ,'field',ergom%<name>,isd,jsd,ntau=1)
    call data_override('OCN', '<trimName>', ergom%<trimName>(isc:iec,jsc:jec,1:nk), model_time)
</constants>
<auxiliaries vertLoc=WAT; isUsedElsewhere=1>
    call g_tracer_get_values(tracer_list,'<name>'  ,'field',ergom%<name>,isd,jsd,ntau=1)
</auxiliaries>

    ! get the pointers to vertical velocity and diffusivity fields
<tracers vertLoc=WAT; vertSpeed/=0>
    call g_tracer_get_pointer(tracer_list,'<name>','vmove',ergom%vmove_<name>)
    call g_tracer_get_pointer(tracer_list,'<name>','vdiff',ergom%vdiff_<name>)
</tracers>

! light calculation
! note: opacity_band already includes light path extension due to sun angle, if respective namelist option is set
!       e.g. ocean_shortwave_jerlov.F90
    allocate(light1(isc:iec,jsc:jec))
    allocate(light2(isc:iec,jsc:jec))

    nband = 1
    if (.NOT. genus_style_par) then
!       k = 1
!       do j = jsc, jec ; do i = isc, iec  !{
!         ergom%irr_inst(i,j,k) = 0.5 * sw_pen_band(nband,i,j) * exp(-opacity_band(nband,i,j,k)* dzt(i,j,k)*0.5)
!       enddo; enddo  !} i,j
!       do k = 2, nk ; do j = jsc, jec ; do i = isc, iec  !{
!         ergom%irr_inst(i,j,k) = ergom%irr_inst(i,j,k-1) * exp(-opacity_band(nband,i,j,k)* dzt(i,j,k))
!       enddo; enddo ; enddo  !} i,j,k
      k=1  
!      estimate light at the middle of the tracer layer, this might not be exactly where the tracer point is
      do j = jsc, jec ; do i = isc, iec
         light1(i,j) = 0.5*sw_pen_band(nband,i,j)
      enddo; enddo
      do k=1, nk
         do j = jsc, jec ; do i = isc, iec
            light2(i,j) = light1(i,j)*exp(-opacity_band(nband,i,j,k)* dzt(i,j,k))
            ergom%irr_inst(i,j,k) = light1(i,j) * exp(-opacity_band(nband,i,j,k)* 0.5*dzt(i,j,k))
            light1(i,j) = light2(i,j)
         enddo; enddo
      enddo
    else
       k = 1
       do j = jsc, jec ; do i = isc, iec  !{
         light_itemp = opacity_band(nband,i,j,k)* dzt(i,j,k)
         light_ntemp = exp(-light_itemp)
         ergom%irr_zw(i,j,k)   = 0.5 * sw_pen_band(nband,i,j) * light_ntemp
         ergom%irr_inst(i,j,k) = 0.5 * sw_pen_band(nband,i,j) * ( 1. -  light_ntemp ) / light_itemp
       enddo; enddo  !} i,j
       do k = 2, nk ; do j = jsc, jec ; do i = isc, iec  !{
         light_itemp = opacity_band(nband,i,j,k)* dzt(i,j,k) + 1.0e-30  ! some schemes give zero opa. in the nk-layer
         light_ntemp = exp(-light_itemp)
         ergom%irr_zw(i,j,k)   = ergom%irr_zw(i,j,k-1) * light_ntemp
         ergom%irr_inst(i,j,k) = ergom%irr_zw(i,j,k-1) * ( 1.-  light_ntemp ) / light_itemp
       enddo; enddo ; enddo  !} i,j,k
    endif
    deallocate(light1)
    deallocate(light2)

    !------------------------------------
    ! STEP 1: calculate total element concentrations
    !------------------------------------

    !------------------------------------
    ! STEP 1.1: calculate total element concentrations in water column
    !------------------------------------
    do k = 1, nk
       do j = jsc, jec
          do i = isc, iec
          <celements>
             ergom%<total>(i,j,k,tau) = &
             <containingTracers vertLoc=WAT>
                max(0.0,ergom%<ct>(i,j,k,tau))*<ctAmount> + &
             </containingTracers>
                0.0
          </celements>
          enddo
       enddo
    enddo

    !------------------------------------
    ! STEP 1.2: calculate total element concentrations in the fluff layer
    !------------------------------------
    do j = jsc, jec
       do i = isc, iec
       <celements>
          ergom%<totalBottom>(i,j) = &
          <containingTracers vertLoc=SED>
             ergom%<ct>(i,j)*<ctAmount> + &
          </containingTracers>
             0.0
       </celements>
       enddo
    enddo

    !------------------------------------
    ! STEP 2: calculate space-independent forcing variables
    !------------------------------------

    call get_date(model_time,cgt_year,cgt_mymonth,cgt_myday,cgt_myhour,cgt_myminute,cgt_mysecond)
    cgt_hour = cgt_myhour + cgt_myminute/60 + cgt_mysecond/3600
    cgt_myhour=days_in_year(model_time)
    if (cgt_myhour .eq. 366) then
       if (cgt_mymonth .eq.  1) cgt_dayofyear = cgt_myday
       if (cgt_mymonth .eq.  2) cgt_dayofyear = cgt_myday +  31
       if (cgt_mymonth .eq.  3) cgt_dayofyear = cgt_myday +  60
       if (cgt_mymonth .eq.  4) cgt_dayofyear = cgt_myday +  91
       if (cgt_mymonth .eq.  5) cgt_dayofyear = cgt_myday + 121
       if (cgt_mymonth .eq.  6) cgt_dayofyear = cgt_myday + 152
       if (cgt_mymonth .eq.  7) cgt_dayofyear = cgt_myday + 182
       if (cgt_mymonth .eq.  8) cgt_dayofyear = cgt_myday + 213
       if (cgt_mymonth .eq.  9) cgt_dayofyear = cgt_myday + 244
       if (cgt_mymonth .eq. 10) cgt_dayofyear = cgt_myday + 274
       if (cgt_mymonth .eq. 11) cgt_dayofyear = cgt_myday + 305
       if (cgt_mymonth .eq. 12) cgt_dayofyear = cgt_myday + 335
    else
       if (cgt_mymonth .eq.  1) cgt_dayofyear = cgt_myday
       if (cgt_mymonth .eq.  2) cgt_dayofyear = cgt_myday +  31
       if (cgt_mymonth .eq.  3) cgt_dayofyear = cgt_myday +  59
       if (cgt_mymonth .eq.  4) cgt_dayofyear = cgt_myday +  90
       if (cgt_mymonth .eq.  5) cgt_dayofyear = cgt_myday + 120
       if (cgt_mymonth .eq.  6) cgt_dayofyear = cgt_myday + 151
       if (cgt_mymonth .eq.  7) cgt_dayofyear = cgt_myday + 181
       if (cgt_mymonth .eq.  8) cgt_dayofyear = cgt_myday + 212
       if (cgt_mymonth .eq.  9) cgt_dayofyear = cgt_myday + 243
       if (cgt_mymonth .eq. 10) cgt_dayofyear = cgt_myday + 273
       if (cgt_mymonth .eq. 11) cgt_dayofyear = cgt_myday + 304
       if (cgt_mymonth .eq. 12) cgt_dayofyear = cgt_myday + 334
    endif

    !------------------------------------
    ! STEP 3: initialize isZIntegral=1 auxiliaries with zero
    !------------------------------------
    <auxiliaries isZIntegral=1; calcAfterProcesses=0>
       ergom%zintegralarray_<name> = 0.0
    </auxiliaries>

    !------------------------------------
    ! STEP 4: Pre-loop where the isZIntegral=1 auxiliaries are calculated
    !------------------------------------
    allocate(cgt_bottomdepth_temp(isc:iec,jsc:jec))
    cgt_bottomdepth_temp = 0.0
    do k = 1, nk
       do j = jsc, jec
          do i = isc, iec
             if (k .gt. grid_kmt(i,j)) cycle
             cgt_bottomdepth_temp(i,j) = cgt_bottomdepth_temp(i,j) + dzt(i,j,k)
             cgt_bottomdepth = cgt_bottomdepth_temp(i,j)
             !------------------------------------
             ! STEP 4.1: prepare abiotic parameters
             !------------------------------------
             cgt_temp       = temp(i,j,k)               ! potential temperature     [Celsius]
             cgt_sali       = salt(i,j,k)               ! salinity                  [g/kg]
             cgt_diffusivity= diff_cbt(i,j,k,1)         ! diffusivity               [m**2/s]
             cgt_light      = ergom%irr_inst(i,j,k)     ! light intensity           [W/m2]
             cgt_cellheight = dzt(i,j,k)                ! cell height               [m]
             cgt_density    = rho_dzt(i,j,k)/dzt(i,j,k) ! density [kg/m3]
             cgt_timestep   = dt/(24*3600)              ! timestep in days
             cgt_longitude  = xt(i,j)                   ! longitude [deg]
             cgt_latitude   = yt(i,j)                   ! latitude [deg]
             cgt_in_sediment = 0.0                      ! we are above sediment surface
             if (k == grid_kmt(i,j)) then
                cgt_current_wave_stress=current_wave_stress(i,j)
             endif

             !------------------------------------
             ! STEP 4.2: load tracer values
             !------------------------------------
<tracers vertLoc=WAT; calcBeforeZIntegral=1>
             <name> = ergom%<name>(i,j,k,tau) ! <description>
             if (k .eq. 1) then
                above_<name> = ergom%<name>(i,j,k,tau)
             else
                above_<name> = ergom%<name>(i,j,k-1,tau)
             endif
</tracers>
<tracers vertLoc=FIS; calcBeforeZIntegral=1>
             <name> = ergom%<name>(i,j) ! <description>
</tracers>
<constants dependsOn=xyz>
             <name> = ergom%<name>(i,j,k) ! <description>
</constants>
<constants dependsOn=xyzt>
             <name> = ergom%<name>(i,j,k) ! <description>
</constants>

<auxiliaries vertLoc=WAT; isUsedElsewhere=1; calcBeforeZIntegral=1>
             <name> = ergom%<name>(i,j,k) ! <description>
</auxiliaries>

<tracers vertLoc=WAT; isPositive=1; calcBeforeZIntegral=1>
             <name>       = max(<name>,0.0)
             above_<name> = max(above_<name>,0.0)
</tracers>
<tracers vertLoc=FIS; isPositive=1; calcBeforeZIntegral=1>
             <name>       = max(<name>,0.0)
</tracers>

             if (k == grid_kmt(i,j)) then
<tracers vertLoc=SED>
                <name> = ergom%<name>(i,j) ! <description>
</tracers>

<tracers vertLoc=SED; isPositive=1>
                <name> = max(<name>,0.0)
</tracers>
             endif

             !------------------------------------
             ! STEP 4.3: calculate auxiliaries
             !------------------------------------

<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=1; calcBeforeZIntegral=1>
             ! <description> :
             <name> = (above_<formula>-<formula>)/cgt_cellheight

</auxiliaries>

             ! initialize auxiliaries for iterative loop
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=0; calcBeforeZIntegral=1; iterations/=0>
             <name> = <iterInit>
</auxiliaries>

             ! iterative loop follows
             do cgt_iteration=1,<maxIterations>
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=0; calcBeforeZIntegral=1; iterations/=0>
                if (cgt_iteration .le. <iterations>) then
                   ! <description> :
                   temp1 = <temp1>
                   temp2 = <temp2>
                   temp3 = <temp3>
                   temp4 = <temp4>
                   temp5 = <temp5>
                   temp6 = <temp6>
                   temp7 = <temp7>
                   temp8 = <temp8>
                   temp9 = <temp9>
                   <name> = <formula>
                endif
</auxiliaries>
             enddo

             ! normal auxiliaries (not iterative)
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=0; calcBeforeZIntegral=1; iterations=0>
             ! <description> :
             temp1 = <temp1>
             temp2 = <temp2>
             temp3 = <temp3>
             temp4 = <temp4>
             temp5 = <temp5>
             temp6 = <temp6>
             temp7 = <temp7>
             temp8 = <temp8>
             temp9 = <temp9>
             <name> = <formula>

</auxiliaries>

             !------------------------------------
             ! STEP 4.4: add contribution of the current layer k to the zIntegral
             !------------------------------------
<auxiliaries isZIntegral=1; calcAfterProcesses=0>
             ergom%zintegralarray_<name>(i,j) = ergom%zintegralarray_<name>(i,j) + (<formula>)*cgt_cellheight*cgt_density
</auxiliaries>

             !------------------------------------
             ! STEP 4.5: pre-calculate auxiliaries with vertLoc=SED or vertLoc=SUR
             !------------------------------------

             if (k == grid_kmt(i,j)) then
                ! initialize auxiliaries for iterative loop
<auxiliaries vertLoc=SED; calcAfterProcesses=0; isZIntegral=0; calcBeforeZIntegral=1; iterations/=0>
                <name> = <iterInit>
</auxiliaries>

                ! iterative loop follows
                do cgt_iteration=1,<maxIterations>
<auxiliaries vertLoc=SED; calcAfterProcesses=0; isZIntegral=0; calcBeforeZIntegral=1; iterations/=0>
                   if (cgt_iteration .le. <iterations>) then
                      ! <description> :
                      temp1 = <temp1>
                      temp2 = <temp2>
                      temp3 = <temp3>
                      temp4 = <temp4>
                      temp5 = <temp5>
                      temp6 = <temp6>
                      temp7 = <temp7>
                      temp8 = <temp8>
                      temp9 = <temp9>
                      <name> = <formula>
                   endif
</auxiliaries>
                enddo

                ! normal auxiliaries (not iterative)
<auxiliaries vertLoc=SED; calcAfterProcesses=0; isZIntegral=0; calcBeforeZIntegral=1; iterations=0>
                ! <description> :
                temp1 = <temp1>
                temp2 = <temp2>
                temp3 = <temp3>
                temp4 = <temp4>
                temp5 = <temp5>
                temp6 = <temp6>
                temp7 = <temp7>
                temp8 = <temp8>
                temp9 = <temp9>
                <name> = <formula>

</auxiliaries>
             endif

             if (k == 1) then
                ! initialize auxiliaries for iterative loop
<auxiliaries vertLoc=SUR; calcAfterProcesses=0; calcBeforeZIntegral=1; iterations/=0>
                <name> = <iterInit>
</auxiliaries>

                ! iterative loop follows
                do cgt_iteration=1,<maxIterations>
<auxiliaries vertLoc=SUR; calcAfterProcesses=0; calcBeforeZIntegral=1; iterations/=0>
                   if (cgt_iteration .le. <iterations>) then
                      ! <description> :
                      temp1 = <temp1>
                      temp2 = <temp2>
                      temp3 = <temp3>
                      temp4 = <temp4>
                      temp5 = <temp5>
                      temp6 = <temp6>
                      temp7 = <temp7>
                      temp8 = <temp8>
                      temp9 = <temp9>
                      <name> = <formula>
                   endif
</auxiliaries>
                enddo

                ! normal auxiliaries (not iterative)
<auxiliaries vertLoc=SUR; calcAfterProcesses=0; calcBeforeZIntegral=1; iterations=0>
                ! <description> :
                temp1 = <temp1>
                temp2 = <temp2>
                temp3 = <temp3>
                temp4 = <temp4>
                temp5 = <temp5>
                temp6 = <temp6>
                temp7 = <temp7>
                temp8 = <temp8>
                temp9 = <temp9>
                <name> = <formula>

</auxiliaries>
             endif

          enddo
       enddo
    enddo

    call mpp_clock_end(id_preloop)
    call mpp_clock_begin(id_mainloop)

    !------------------------------------
    ! STEP 5: initializations before the main loop
    !------------------------------------

    !------------------------------------
    ! STEP 5.1: initialize isZIntegral=1 auxiliaries which are calculated after the process rates with zero
    !------------------------------------

<auxiliaries isZIntegral=1; calcAfterProcesses=1>
       ergom%zintegralarray_<name> = 0.0
</auxiliaries>

    !------------------------------------
    ! STEP 5.2: initialize the arrays which cumulate the change of vertLoc=FIS tracers with zero
    !------------------------------------
<tracers vertLoc=FIS>
    ergom%cumulated_change_of_<name> = 0.0
</tracers>

    !------------------------------------
    ! STEP 6: Main loop
    !------------------------------------

    cgt_bottomdepth_temp = 0.0
    do k = 1, nk
       do j = jsc, jec
          do i = isc, iec
             if (k .gt. grid_kmt(i,j)) cycle
             cgt_bottomdepth_temp(i,j) = cgt_bottomdepth_temp(i,j) + dzt(i,j,k)
             cgt_bottomdepth = cgt_bottomdepth_temp(i,j)
             !------------------------------------
             ! STEP 6.1: prepare abiotic parameters
             !------------------------------------
             cgt_temp       = temp(i,j,k)               ! potential temperature     [Celsius]
             cgt_sali       = salt(i,j,k)               ! salinity                  [g/kg]
             cgt_diffusivity= diff_cbt(i,j,k,1)         ! diffusivity               [m**2/s]
             cgt_light      = ergom%irr_inst(i,j,k)     ! light intensity           [W/m2]
             cgt_cellheight = dzt(i,j,k)                ! cell height               [m]
             cgt_density    = rho_dzt(i,j,k)/dzt(i,j,k) ! density [kg/m3]
             cgt_timestep   = dt/(24*3600)              ! timestep in days
             cgt_longitude  = xt(i,j)                   ! longitude [deg]
             cgt_latitude   = yt(i,j)                   ! latitude [deg]
             cgt_in_sediment = 0.0                      ! we are above sediment surface
             if (k == grid_kmt(i,j)) then
                cgt_current_wave_stress=current_wave_stress(i,j)
             endif

             !------------------------------------
             ! STEP 6.2: load tracer values
             !------------------------------------
<tracers vertLoc=WAT>
             <name> = ergom%<name>(i,j,k,tau) ! <description>
             if (k .eq. 1) then
                above_<name> = ergom%<name>(i,j,k,tau)
             else
                above_<name> = ergom%<name>(i,j,k-1,tau)
             endif
</tracers>
<tracers vertLoc=FIS>
             <name> = ergom%<name>(i,j) ! <description>
</tracers>
<constants dependsOn=xyz>
             <name> = ergom%<name>(i,j,k) ! <description>
</constants>
<constants dependsOn=xyzt>
             <name> = ergom%<name>(i,j,k) ! <description>
</constants>

<auxiliaries vertLoc=WAT; isUsedElsewhere=1>
             <name> = ergom%<name>(i,j,k) ! <description>
</auxiliaries>

<tracers vertLoc=WAT; isPositive=1>
             <name>       = max(<name>,0.0)
             above_<name> = max(above_<name>,0.0)
</tracers>
<tracers vertLoc=FIS; isPositive=1>
             <name>       = max(<name>,0.0)
</tracers>

             if (k == grid_kmt(i,j)) then
<tracers vertLoc=SED>
                <name> = ergom%<name>(i,j) ! <description>
</tracers>

<tracers vertLoc=SED; isPositive=1>
                <name> = max(<name>,0.0)
</tracers>
             endif

             !------------------------------------
             ! STEP 6.3: opacity calculation
             !------------------------------------
             ! some tracers contribute to the opacity
             ! opacity in m**2/mol, concentration in mol/kg -> result in m**2/kg
             ergom%bio_opacity(i,j,k) = (0.0 &
             <tracers vertLoc=WAT; opacity/=0>
                + <opacity> * <name> &
             </tracers>
                )

             ! now, convert it from m**2/kg to 1/m
# if defined salt4yellowSubs
             ! in case salinity is used to parametrize yellow substances
             !  note: cgt_sali is reduced since MOM overestimates salinity
             ergom%bio_opacity(i,j,k) = ergom%bio_opacity(i,j,k) * ergom%rho_0 &
             + 0.2787 * max(0.3,cgt_sali-1.0)**(-0.627)
# else
             ergom%bio_opacity(i,j,k) = ergom%bio_opacity(i,j,k) * ergom%rho_0
# endif

             !------------------------------------
             ! STEP 6.4: calculate auxiliaries
             !------------------------------------
             !------------------------------------
             ! STEP 6.4.1: get isZIntegral values from pre-loop
             !------------------------------------
<auxiliaries isZIntegral=1; calcAfterProcesses=0>
             ! <description> :
             <name> = ergom%zintegralarray_<name>(i,j)
</auxiliaries>

             !------------------------------------
             ! STEP 6.4.2: calculate isZGradient auxiliaries
             !------------------------------------
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=1>
             ! <description> :
             <name> = (above_<formula>-<formula>)/cgt_cellheight

</auxiliaries>

             !------------------------------------
             ! STEP 6.4.3: vertLoc=WAT auxiliaries calclated iteratively
             !------------------------------------

             ! initialize auxiliaries for iterative loop
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=0; iterations/=0>
             <name> = <iterInit>
</auxiliaries>

             ! iterative loop follows
             do cgt_iteration=1,<maxIterations>
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=0; iterations/=0>
                if (cgt_iteration .le. <iterations>) then
                   ! <description> :
                   temp1 = <temp1>
                   temp2 = <temp2>
                   temp3 = <temp3>
                   temp4 = <temp4>
                   temp5 = <temp5>
                   temp6 = <temp6>
                   temp7 = <temp7>
                   temp8 = <temp8>
                   temp9 = <temp9>
                   <name> = <formula>
                endif
</auxiliaries>
             enddo

             !------------------------------------
             ! STEP 6.4.4: vertLoc=WAT auxiliaries calclated normally (not iteratively)
             !------------------------------------
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=0; iterations=0>
             ! <description> :
             temp1 = <temp1>
             temp2 = <temp2>
             temp3 = <temp3>
             temp4 = <temp4>
             temp5 = <temp5>
             temp6 = <temp6>
             temp7 = <temp7>
             temp8 = <temp8>
             temp9 = <temp9>
             <name> = <formula>

</auxiliaries>

             !------------------------------------
             ! STEP 6.4.5: vertLoc=SED auxiliaries calclated iteratively
             !------------------------------------
             if (k == grid_kmt(i,j)) then
                ! initialize auxiliaries for iterative loop
<auxiliaries vertLoc=SED; calcAfterProcesses=0; isZIntegral=0; iterations/=0>
                <name> = <iterInit>
</auxiliaries>

                ! iterative loop follows
                do cgt_iteration=1,<maxIterations>
<auxiliaries vertLoc=SED; calcAfterProcesses=0; isZIntegral=0; iterations/=0>
                   if (cgt_iteration .l. <iterations>) then
                      ! <description> :
                      temp1 = <temp1>
                      temp2 = <temp2>
                      temp3 = <temp3>
                      temp4 = <temp4>
                      temp5 = <temp5>
                      temp6 = <temp6>
                      temp7 = <temp7>
                      temp8 = <temp8>
                      temp9 = <temp9>
                      <name> = <formula>
                   endif
</auxiliaries>
                enddo

             !------------------------------------
             ! STEP 6.4.6: vertLoc=SED auxiliaries calclated normally (not iteratively)
             !------------------------------------

<auxiliaries vertLoc=SED; calcAfterProcesses=0; isZIntegral=0; iterations=0>
                ! <description> :
                temp1 = <temp1>
                temp2 = <temp2>
                temp3 = <temp3>
                temp4 = <temp4>
                temp5 = <temp5>
                temp6 = <temp6>
                temp7 = <temp7>
                temp8 = <temp8>
                temp9 = <temp9>
                <name> = <formula>

</auxiliaries>
             endif

             !------------------------------------
             ! STEP 6.4.7: vertLoc=SUR auxiliaries calclated iteratively
             !------------------------------------

             if (k == 1) then
                ! initialize auxiliaries for iterative loop
<auxiliaries vertLoc=SUR; calcAfterProcesses=0; iterations/=0>
                <name> = <iterInit>
</auxiliaries>

                ! iterative loop follows
                do cgt_iteration=1,<maxIterations>
<auxiliaries vertLoc=SUR; calcAfterProcesses=0; iterations/=0>
                   if (cgt_iteration .le. <iterations>) then
                      ! <description> :
                      temp1 = <temp1>
                      temp2 = <temp2>
                      temp3 = <temp3>
                      temp4 = <temp4>
                      temp5 = <temp5>
                      temp6 = <temp6>
                      temp7 = <temp7>
                      temp8 = <temp8>
                      temp9 = <temp9>
                      <name> = <formula>
                   endif
</auxiliaries>
                enddo

             !------------------------------------
             ! STEP 6.4.8: vertLoc=SUR auxiliaries calclated normally (not iteratively)
             !------------------------------------
<auxiliaries vertLoc=SUR; calcAfterProcesses=0; iterations=0>
                ! <description> :
                temp1 = <temp1>
                temp2 = <temp2>
                temp3 = <temp3>
                temp4 = <temp4>
                temp5 = <temp5>
                temp6 = <temp6>
                temp7 = <temp7>
                temp8 = <temp8>
                temp9 = <temp9>
                <name> = <formula>

</auxiliaries>
             endif


             !------------------------------------
             ! STEP 6.5: output of auxiliaries
             !------------------------------------
             if (.not. intermediate) then

<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isUsedElsewhere=0>
                if (ergom%id_<name> .gt. 0) then
                  ergom%<name>(i,j,k) = <name>
                endif
</auxiliaries>
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isUsedElsewhere=1>
                ergom%<name>(i,j,k) = <name>
</auxiliaries>

                if (k == grid_kmt(i,j)) then
<auxiliaries vertLoc=SED; calcAfterProcesses=0>
                   if (ergom%id_<name> .gt. 0) then
                      ergom%<name>(i,j) = <name>
                   endif
</auxiliaries>
                endif
                if (k == 1) then
<auxiliaries vertLoc=SUR; calcAfterProcesses=0>
                   if (ergom%id_<name> .gt. 0) then
                      ergom%<name>(i,j) = <name>
                   endif
</auxiliaries>
                endif
             endif

             !------------------------------------
             ! STEP 6.6: calculate process limitations
             !------------------------------------
<tracers vertLoc=WAT>
  <limitations>
             <name> = <formula>
  </limitations>
</tracers>

             if (k == grid_kmt(i,j)) then
<tracers vertLoc=SED>
  <limitations>
                <name> = <formula>
  </limitations>
</tracers>
             endif

             if (k == 1) then
<tracers vertLoc=SUR>
  <limitations>
                <name> = <formula>
  </limitations>
</tracers>
             endif

             !------------------------------------
             ! STEP 6.7
             !
             !-- POSITIVE-DEFINITE SCHEME --------
             !-- means the following steps will be repeated as often as nessecary
             !------------------------------------
             number_of_loop = 1
             fraction_of_total_timestep = 1.0   ! how much of the original timestep is remaining
<processes>
             total_rate_<name>          = 0.0
</processes>
             if (implicit_bottomfluxes) then
<tracers vertLoc=WAT>
                ergom%btf_<name>(i,j) = 0.0
</tracers>
             endif

             do while (cgt_timestep .gt. 0.0)

                !------------------------------------
                ! STEP 6.7.1: calculate process rates
                !------------------------------------
<processes vertLoc=WAT>
                ! <description> :
                <name> = <turnover>
                <name> = max(<name>,0.0)

</processes>

                if (k == grid_kmt(i,j)) then
<processes vertLoc=SED>
                   ! <description> :
                   <name> = <turnover>
                   <name> = max(<name>,0.0)

</processes>
                endif

                if (k == 1) then
<processes vertLoc=SUR>
                   ! <description> :
                   <name> = <turnover>
                   <name> = max(<name>,0.0)

</processes>
                endif

                !------------------------------------
                ! STEP 6.7.2: save the process rates if intermediate time step (Runge-Kutta)
                !------------------------------------
                if ((number_of_loop == 1) .and. intermediate) then
<processes vertLoc=WAT>
                   ergom%saved_rate_<name>(i,j,k)=ergom%saved_rate_<name>(i,j,k)+<name>
</processes>
                   if (k == grid_kmt(i,j)) then
<processes vertLoc=SED>
                      ergom%saved_rate_<name>(i,j)=ergom%saved_rate_<name>(i,j)+<name>
</processes>
                   endif
                   if (k == 1) then
<processes vertLoc=SUR>
                      ergom%saved_rate_<name>(i,j)=ergom%saved_rate_<name>(i,j)+<name>
</processes>
                   endif
                endif

                !------------------------------------
                ! STEP 6.7.3: use saved rates if they apply (Runge-Kutta)
                !------------------------------------
                if (use_saved_rates) then
<processes vertLoc=WAT>
                   if ((<name> > 0.0) .and. (ergom%saved_rate_<name>(i,j,k) >= 0.0)) then
                      <name> = ergom%saved_rate_<name>(i,j,k)
                   endif
</processes>
<processes vertLoc=WAT; isStiff/=0>
                   if ((<name> > 0.0) .and. (ergom%saved_rate_<name>(i,j,k) >= 0.0)) then
                      <name> = <name> * patankar_modification_<stiffTracer>(i,j,k)
                   endif
</processes>
                   if (k == grid_kmt(i,j)) then
<processes vertLoc=SED>
                      if ((<name> > 0.0) .and. (ergom%saved_rate_<name>(i,j) >= 0.0)) then
                         <name> = ergom%saved_rate_<name>(i,j)
                      endif
</processes>
<processes vertLoc=SED; isStiff/=0>
                      if ((<name> > 0.0) .and. (ergom%saved_rate_<name>(i,j) >= 0.0)) then
                         <name> = <name> * patankar_modification_<stiffTracer>(i,j,k)
                      endif
</processes>
                   endif
                   if (k == 1) then
<processes vertLoc=SUR>
                      if ((<name> > 0.0) .and. (ergom%saved_rate_<name>(i,j) >= 0.0)) then
                         <name> = ergom%saved_rate_<name>(i,j)
                      endif
</processes>
<processes vertLoc=SUR; isStiff/=0>
                      if ((<name> > 0.0) .and. (ergom%saved_rate_<name>(i,j) >= 0.0)) then
                         <name> = <name> * patankar_modification_<stiffTracer>(i,j,k)
                      endif
</processes>
                   endif
                endif

                !------------------------------------
                ! STEP 6.7.4: apply Patankar limitations
                !------------------------------------

<processes vertLoc=WAT; isStiff/=0>
                ! <description> :
                after_patankar_<name> = <name> * <stiffFactor>

</processes>
                if (k == grid_kmt(i,j)) then
<processes vertLoc=SED; isStiff/=0>
                   ! <description> :
                   after_patankar_<name> = <name> * <stiffFactor>

</processes>
                endif
                if (k == 1) then
<processes vertLoc=SUR; isStiff/=0>
                   ! <description> :
                   after_patankar_<name> = <name> * <stiffFactor>

</processes>
                endif

<processes vertLoc=WAT; isStiff/=0>
                ! <description> :
                <name> = after_patankar_<name>

</processes>
                if (k == grid_kmt(i,j)) then
<processes vertLoc=SED; isStiff/=0>
                   ! <description> :
                   <name> = after_patankar_<name>

</processes>
                endif
                if (k == 1) then
<processes vertLoc=SUR; isStiff/=0>
                   ! <description> :
                   <name> = after_patankar_<name>

</processes>
                endif

                !------------------------------------
                ! STEP 6.7.5: calculate possible euler-forward change (in a full timestep)
                !------------------------------------
<tracers>
                change_of_<name> = 0.0
</tracers>
                if (implicit_bottomfluxes) then
<tracers vertLoc=WAT>
                   change_btf_<name> = 0.0
</tracers>
                endif


<tracers vertLoc=WAT; hasTimeTendenciesVertLoc=WAT>

                change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                <timeTendencies vertLoc=WAT>
                   <timeTendency> & ! <description>
                </timeTendencies>
                )
</tracers>
<tracers vertLoc=FIS; hasTimeTendenciesVertLoc=WAT>

                change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                <timeTendencies vertLoc=WAT>
                   <timeTendency> & ! <description>
                </timeTendencies>
                )
</tracers>

                if (k == 1) then
<tracers vertLoc=WAT; hasTimeTendenciesVertLoc=SUR>

                   change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                   <timeTendencies vertLoc=SUR>
                      <timeTendency> & ! <description>
                   </timeTendencies>
                   )
</tracers>
               endif

                if (k == grid_kmt(i,j)) then
                   if (implicit_bottomfluxes) then
<tracers vertLoc=WAT; hasTimeTendenciesVertLoc=SED>

                      change_btf_<name> = change_btf_<name> + cgt_timestep*(0.0 &
                      <timeTendencies vertLoc=SED>
                         <timeTendency2d> & ! <description>
                      </timeTendencies>
                      )
</tracers>
                   else
<tracers vertLoc=WAT; hasTimeTendenciesVertLoc=SED>

                      change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                      <timeTendencies vertLoc=SED>
                         <timeTendency> & ! <description>
                      </timeTendencies>
                      )
</tracers>
                   endif
<tracers vertLoc=SED; hasTimeTendencies>

                   change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                   <timeTendencies>
                      <timeTendency> & ! <description>
                   </timeTendencies>
                   )
</tracers>
<tracers vertLoc=FIS; hasTimeTendenciesVertLoc=SED>

                   change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                   <timeTendencies vertLoc=SED>
                      <timeTendency> & ! <description>
                   </timeTendencies>
                   )
</tracers>

                endif

                !------------------------------------
                ! STEP 6.7.6: calculate maximum fraction of the timestep before some tracer gets exhausted
                !------------------------------------

                timestep_fraction = 1.0
                which_tracer_exhausted = -1

                ! find the tracer which is exhausted after the shortest period of time

                ! in the water column
<tracers vertLoc=WAT; isPositive=1>

                ! check if tracer <name> was exhausted from the beginning and is still consumed
                if ((ergom%<name>(i,j,k,tau) .le. 0.0) .and. (change_of_<name>  .lt. 0.0)) then
                   timestep_fraction = 0.0
                   which_tracer_exhausted = <numTracer>
                endif
                ! check if tracer <name> was present, but got exhausted
                if ((ergom%<name>(i,j,k,tau) .gt. 0.0) .and. (ergom%<name>(i,j,k,tau) + change_of_<name> .lt. 0.0)) then
                   timestep_fraction_new = ergom%<name>(i,j,k,tau) / (0.0 - change_of_<name>)
                   if (timestep_fraction_new .le. timestep_fraction) then
                      which_tracer_exhausted = <numTracer>
                      timestep_fraction = timestep_fraction_new
                   endif
                endif
</tracers>

                ! in the bottom layer
                if (k == grid_kmt(i,j)) then
<tracers vertLoc=SED; isPositive=1>

                   ! check if tracer <name> was exhausted from the beginning and is still consumed
                   if ((ergom%<name>(i,j) .le. 0.0) .and. (change_of_<name> .lt. 0.0)) then
                      timestep_fraction = 0.0
                      which_tracer_exhausted = <numTracer>
                   endif
                   ! check if tracer <name> was present, but got exhausted
                   if ((ergom%<name>(i,j) .gt. 0.0) .and. (ergom%<name>(i,j) + change_of_<name> .lt. 0.0)) then
                      timestep_fraction_new = ergom%<name>(i,j) / (0.0 - change_of_<name>)
                      if (timestep_fraction_new .le. timestep_fraction) then
                         which_tracer_exhausted = <numTracer>
                         timestep_fraction = timestep_fraction_new
                      endif
                   endif
</tracers>
                endif

                !------------------------------------
                ! STEP 6.7.7: update the limitations: processes limited by this tracer become zero in the future
                !------------------------------------

<tracers isPositive=1; vertLoc/=FIS>
                if (<numTracer> .eq. which_tracer_exhausted) then
                  <limitations>
                   <name> = 0.0
                  </limitations>
                endif
</tracers>

                !------------------------------------
                ! STEP 6.7.8: apply a Euler-forward timestep with the fraction of the time
                !------------------------------------
                ! in the water column
<tracers vertLoc=WAT>

                ! tracer <name> (<description>):
                ergom%<name>(i,j,k,tau) = ergom%<name>(i,j,k,tau) + change_of_<name> * timestep_fraction
</tracers>
                if (implicit_bottomfluxes) then
<tracers vertLoc=WAT>

                ! tracer <name> (<description>):
                ergom%btf_<name>(i,j) = ergom%btf_<name>(i,j) + change_btf_<name> * timestep_fraction
</tracers>
                endif
<tracers vertLoc=FIS>

                ! tracer <name> (<description>)
                ergom%cumulated_change_of_<name>(i,j) = ergom%cumulated_change_of_<name>(i,j) + change_of_<name> * timestep_fraction
</tracers>


                ! in the bottom layer
                if (k == grid_kmt(i,j)) then
<tracers vertLoc=SED>

                   ! tracer <name> (<description>)
                   ergom%<name>(i,j) = ergom%<name>(i,j) + change_of_<name> * timestep_fraction
</tracers>
                endif

                !------------------------------------
                ! STEP 6.7.9: save process rates to obtain total rate over full timestep in the end
                !------------------------------------
                if (.not. intermediate) then
<processes vertLoc=WAT>
                   total_rate_<name> = total_rate_<name> + <name> * timestep_fraction * fraction_of_total_timestep
</processes>
                   if (k == grid_kmt(i,j)) then
<processes vertLoc=SED>
                      total_rate_<name> = total_rate_<name> + <name> * timestep_fraction * fraction_of_total_timestep
</processes>
                   endif
                   if (k == 1) then
<processes vertLoc=SUR>
                      total_rate_<name> = total_rate_<name> + <name> * timestep_fraction * fraction_of_total_timestep
</processes>
                   endif
                endif

                !------------------------------------
                ! STEP 6.7.10: set timestep to remaining timestep only
                !------------------------------------
                cgt_timestep = cgt_timestep * (1.0 - timestep_fraction)                         ! remaining timestep
                fraction_of_total_timestep = fraction_of_total_timestep * (1.0 - timestep_fraction) ! how much of the original timestep is remaining

                if (number_of_loop .gt. 100) then
                   print*,"LOOP: ",number_of_loop
                   print*,"i=",i,", j=",j,", k=",k
                 <tracers isPositive=1>
                   if (<numTracer> .eq. which_tracer_exhausted) then
                     print*,"limiting tracer: <name>"
                     print*,"  <name> = ",<name>
                     print*,"  change_of_<name> = ",change_of_<name>
                   endif
                 </tracers>
                   print*,"tracers:"
                 <tracers>
                   print*,"  <name> = ",<name>
                 </tracers>
                   print*,"changes:"
                 <tracers>
                   print*,"  change_of_<name> = ",change_of_<name>
                 </tracers>
                   print*,"limitations:"
                 <tracers>
                  <limitations>
                   print*,"  <name> = ",<name>
                  </limitations>
                 </tracers>
                   print*,"cgt_timestep: ",cgt_timestep
                   print*,"timestep_fraction: ",timestep_fraction
                   print*,"  "
                   if (number_of_loop .gt. 200) then
                      stop
                   endif
                endif
                number_of_loop=number_of_loop+1

             enddo
             !------------------------------------
             !-- END OF POSITIVE-DEFINITE SCHEME -
             !------------------------------------

             !------------------------------------
             ! STEP 6.8: cumulate change of vertLoc=FIS tracers over the k levels
             !------------------------------------
             if (k == grid_kmt(i,j)) then
<tracers vertLoc=FIS>
                ! apply change of <description>:
                ergom%<name>(i,j) = ergom%<name>(i,j) + cumulated_change_of_<name>
</tracers>
             endif

             if (.not. intermediate) then

                !------------------------------------
                ! STEP 6.9: output of new tracer concentrations
                !------------------------------------

                ! not required in MOM5

                !------------------------------------
                ! STEP 6.10: get total process rates over full time step
                !------------------------------------

<processes vertLoc=WAT>
                <name> = total_rate_<name>
</processes>
                if (k == grid_kmt(i,j)) then
<processes vertLoc=SED>
                   <name> = total_rate_<name>
</processes>
                endif
                if (k == 1) then
<processes vertLoc=SUR>
                   <name> = total_rate_<name>
</processes>
                endif

                !------------------------------------
                ! STEP 6.11: output of process rates
                !------------------------------------
<processes vertLoc=WAT>
                if (ergom%id_<name> .gt. 0) then
                  ergom%<name>(i,j,k) = total_rate_<name>
                endif
</processes>
                if (k == grid_kmt(i,j)) then
<processes vertLoc=SED>
                   if (ergom%id_<name> .gt. 0) then
                      ergom%<name>(i,j) = total_rate_<name>
                   endif
</processes>
                endif
                if (k == 1) then
<processes vertLoc=SUR>
                   if (ergom%id_<name> .gt. 0) then
                      ergom%<name>(i,j) = total_rate_<name>
                   endif
</processes>
                endif

                !------------------------------------
                ! STEP 6.12: calculate "late" auxiliaries
                !------------------------------------
                !------------------------------------
                ! STEP 6.12.1: calculate all but the isZIntegral auxiliaries
                !------------------------------------
<auxiliaries vertLoc=WAT; calcAfterProcesses=1; isZGradient=0>
                ! <description> :
                temp1 = <temp1>
                temp2 = <temp2>
                temp3 = <temp3>
                temp4 = <temp4>
                temp5 = <temp5>
                temp6 = <temp6>
                temp7 = <temp7>
                temp8 = <temp8>
                temp9 = <temp9>
                <name> = <formula>

</auxiliaries>

                if (k == grid_kmt(i,j)) then
<auxiliaries vertLoc=SED; calcAfterProcesses=1; isZIntegral=0>
                   ! <description> :
                   temp1 = <temp1>
                   temp2 = <temp2>
                   temp3 = <temp3>
                   temp4 = <temp4>
                   temp5 = <temp5>
                   temp6 = <temp6>
                   temp7 = <temp7>
                   temp8 = <temp8>
                   temp9 = <temp9>
                   <name> = <formula>

</auxiliaries>
                endif

                if (k == 1) then
<auxiliaries vertLoc=SUR; calcAfterProcesses=1>
                   ! <description> :
                   temp1 = <temp1>
                   temp2 = <temp2>
                   temp3 = <temp3>
                   temp4 = <temp4>
                   temp5 = <temp5>
                   temp6 = <temp6>
                   temp7 = <temp7>
                   temp8 = <temp8>
                   temp9 = <temp9>
                   <name> = <formula>

</auxiliaries>
                endif

                !------------------------------------
                ! STEP 6.12.2: cumulate the isZIntegral auxiliaries over the k levels
                !------------------------------------
<auxiliaries vertLoc=SED; calcAfterProcesses=1; isZIntegral=1>
                ergom%zintegralarray_<name>(i,j) = ergom%zintegralarray_<name>(i,j) + (<formula>)*cgt_cellheight*cgt_density
</auxiliaries>

                !------------------------------------
                ! STEP 6.13: output of "late" auxiliaries
                !------------------------------------
<auxiliaries vertLoc=WAT; calcAfterProcesses=1; Isusedelsewhere=0>
                if (ergom%id_<name> .gt. 0) then
                  ergom%<name>(i,j,k) = <name>
                endif
</auxiliaries>
<auxiliaries vertLoc=WAT; calcAfterProcesses=1; isUsedElsewhere=1>
                ergom%<name>(i,j,k) = <name>
</auxiliaries>
                if (k == grid_kmt(i,j)) then
<auxiliaries vertLoc=SED; calcAfterProcesses=1; isZIntegral=0>
                   if (ergom%id_<name> .gt. 0) then
                      ergom%<name>(i,j) = <name>
                   endif
</auxiliaries>
<auxiliaries vertLoc=SED; calcAfterProcesses=1; isZIntegral=1>
                   if (ergom%id_<name> .gt. 0) then
                      ergom%<name>(i,j) = ergom%zintegralarray_<name>(i,j)
                   endif
</auxiliaries>
                endif
                if (k == 1) then
<auxiliaries vertLoc=SUR; calcAfterProcesses=1>
                   if (ergom%id_<name> .gt. 0) then
                      ergom%<name>(i,j) = <name>
                   endif
</auxiliaries>
                endif

                !---------------------------------------
                ! STEP 6.14: passing vertical velocity and diffusivity to the coupler
                !---------------------------------------

                <tracers vertLoc=WAT; vertSpeed/=0>
                  ergom%vmove_<name>(i,j,k) = (<vertSpeed>)/(24*3600)
                  ergom%vdiff_<name>(i,j,k) = <vertDiff>
                </tracers>

                !---------------------------------------
                ! STEP 6.15: output of parameters for surface flux of gasses
                !---------------------------------------

                ! Solubility has to be converted from mol/kg/Pa to mol/m3/Pa
                ! Surface concentration has to be converted from mol/kg to mol/m3

                if ((k == 1) .and. (implicit_surfacefluxes)) then
                  <tracers solubility/=0; childOf=none>
                   ergom%<gasName>_<trimName>_csurf(i,j) = <gasName> * ergom%Rho_0  ! surface concentration of dissolved <gasName> [mol/m3]
                  </tracers>
                  <tracers solubility/=0; childOf/=none>
                   ergom%<gasName>_<trimName>_csurf(i,j) = <gasName> * min(1.0,<trimName>/max(<childOf>,0.0000000001)) * ergom%Rho_0  ! surface concentration of dissolved <gasName> [mol/m3]
                  </tracers>
                  <tracers solubility/=0; childOf=none>
                   ergom%<gasName>_alpha(i,j) = <solubility> * ergom%Rho_0  ! solubility of <gasName> [mol/m3/Pa]
                   ergom%<gasName>_sc_no(i,j) = <schmidtNumber>  ! Schmidt number of <gasName> [1]
                  </tracers>
                endif

             endif

          enddo
       enddo
    enddo
    deallocate(cgt_bottomdepth_temp)

    !-----------------------------------
    ! biological timestep has ended
    !-----------------------------------

    call mpp_clock_end(id_mainloop)
    if (.not. intermediate) then

       call mpp_clock_begin(id_vertmig)

       !-----------------------------------
       ! STEP 7: send surface concentration, solubility and Schmidt number to the coupler for gas exchange
       !-----------------------------------

       if (implicit_surfacefluxes) then
          <tracers solubility/=0>
          call g_tracer_set_values(tracer_list,'<name>','alpha',ergom%<gasName>_alpha ,isd,jsd)
          call g_tracer_set_values(tracer_list,'<name>','csurf',ergom%<gasName>_<trimName>_csurf ,isd,jsd)
          call g_tracer_set_values(tracer_list,'<name>','sc_no',ergom%<gasName>_sc_no ,isd,jsd)
          </tracers>
       endif

       !-----------------------------------
       ! STEP 8: vertical movement of tracers
       !-----------------------------------
       !-----------------------------------
       ! STEP 8.1: calculate total element concentrations in the water column
       !-----------------------------------
       do k = 1, nk
          do j = jsc, jec
             do i = isc, iec
                if (k .gt. grid_kmt(i,j)) cycle
             <celements>
                ergom%<total>(i,j,k,tau) = &
                <containingTracers vertLoc=WAT>
                   max(0.0,ergom%<ct>(i,j,k,tau))*<ctAmount> + &
                </containingTracers>
                   0.0
             </celements>
             enddo
          enddo
       enddo

       if (.not. implicit_movement) then

          !-----------------------------------
          ! STEP 8.2: do several vertical movement steps
          !-----------------------------------
          do m = 1, NUM_VMOVE_STEPS
             !-----------------------------------
             ! STEP 8.2.1: store age concentrations in an array which is not modified by movement
             !-----------------------------------
             <celements isAging/=0>
               ergom%old_<aged> = ergom%<aged>(:,:,:,tau)
             </celements>
             !-----------------------------------
             ! STEP 8.2.2: move the age concentrations of the marked elements
             !-----------------------------------
             <tracers childOf/=none; vertLoc=WAT; vertSpeed/=0; hasCeAged>
                call generic_ERGOM_vmove(ergom%vmove_<name>, &
                                         ergom%<ceAgedName>(:,:,:,tau), ergom%old_<ceAgedName>(:,:,:,tau), &
                                         ergom%<name>(:,:,:,tau)*<ceAmount>, ergom%<ceTotalName>(:,:,:,tau), &
                                         dzt, dt/NUM_VMOVE_STEPS, isd, ied, jsd, jed)

             </tracers>
             !-----------------------------------
             ! STEP 8.2.3: move the tracers (including tagged tracers) themselves
             !-----------------------------------
             <tracers vertLoc=WAT; vertSpeed/=0>
                call generic_ERGOM_vmove(ergom%vmove_<name>, &
                                         ergom%<name>(:,:,:,tau), ergom%<name>(:,:,:,tau), &
                                         ergom%<name>(:,:,:,tau), ergom%<name>(:,:,:,tau), &
                                         dzt, dt/NUM_VMOVE_STEPS, isd, ied, jsd, jed)
             </tracers>

             !-----------------------------------
             ! STEP 8.2.4: calculate the new total marked element concentrations
             !-----------------------------------
             do k = 1, nk
                do j = jsc, jec
                   do i = isc, iec
                      if (k .gt. grid_kmt(i,j)) cycle
                   <celements>
                      ergom%<total>(i,j,k,tau) = &
                      <containingTracers vertLoc=WAT>
                         max(0.0,ergom%<ct>(i,j,k,tau))*<ctAmount> + &
                      </containingTracers>
                         0.0
                   </celements>
                   enddo
                enddo
             enddo
          enddo

          !-----------------------------------
          ! STEP 9: vertical diffusion of tracers
          !-----------------------------------
          !-----------------------------------
          ! STEP 9.1: calculate total element concentrations in the water column
          !-----------------------------------
          ! not necessary, as it was done in the movement steps before
          !-----------------------------------
          ! STEP 9.2: do several vertical diffusion steps
          !-----------------------------------
          do m = 1, NUM_VMOVE_STEPS
             !-----------------------------------
             ! STEP 9.2.1: store age concentrations in an array which is not modified by diffusion
             !-----------------------------------
             <celements isAging/=0>
               ergom%old_<aged> = ergom%<aged>(:,:,:,tau)
             </celements>
             !-----------------------------------
             ! STEP 9.2.2: diffuse the age concentrations of the marked elements
             !-----------------------------------
             <tracers childOf/=none; vertLoc=WAT; vertSpeed/=0; hasCeAged>
                call generic_ERGOM_vdiff(ergom%vdiff_<name>, &
                                         ergom%<ceAgedName>(:,:,:,tau), ergom%old_<ceAgedName>(:,:,:,tau),&
                                         ergom%<name>(:,:,:,tau)*<ceAmount>, ergom%<ceTotalName>(:,:,:,tau), &
                                         dzt, dt/NUM_VMOVE_STEPS, isd, ied, jsd, jed)

             </tracers>

             !-----------------------------------
             ! STEP 9.2.3: move the tracers (including tagged tracers) themselves
             !-----------------------------------

             <tracers vertLoc=WAT; vertSpeed/=0>
                call generic_ERGOM_vdiff(ergom%vdiff_<name>, &
                                         ergom%<name>(:,:,:,tau), ergom%<name>(:,:,:,tau), &
                                         ergom%<name>(:,:,:,tau), ergom%<name>(:,:,:,tau), &
                                         dzt, dt/NUM_VMOVE_STEPS, isd, ied, jsd, jed)
             </tracers>

             !-----------------------------------
             ! STEP 9.2.4: calculate the new total marked element concentrations
             !-----------------------------------
             do k = 1, nk
                do j = jsc, jec
                   do i = isc, iec
                      if (k .gt. grid_kmt(i,j)) cycle
                   <celements>
                      ergom%<total>(i,j,k,tau) = &
                      <containingTracers vertLoc=WAT>
                         max(0.0,ergom%<ct>(i,j,k,tau))*<ctAmount> + &
                      </containingTracers>
                         0.0
                   </celements>
                   enddo
                enddo
             enddo
          enddo

          ! switch off movement now to avoid that implicit routine does additional movement/diffusion
       <tracers vertLoc=WAT; vertSpeed/=0>
          ergom%vmove_<name>(i,j,k) = 0.0
          ergom%vdiff_<name>(i,j,k) = 0.0
       </tracers>
       endif

       !-----------------------------------
       ! STEP 10: calculate the new total marked element concentrations at the bottom
       !-----------------------------------

       do j = jsc, jec
          do i = isc, iec
             if (1 .gt. grid_kmt(i,j)) cycle
          <celements>
             ergom%<totalBottom>(i,j) = &
             <containingTracers vertLoc=SED>
                max(0.0,ergom%<ct>(i,j))*<ctAmount> + &
             </containingTracers>
                0.0
          </celements>
          enddo
       enddo

       call mpp_clock_end(id_vertmig)
       call mpp_clock_begin(id_output)

       !-----------------------------------
       ! STEP 11: save the new values to their designated places
       !-----------------------------------

       call user_set_2d_tracer_values(tracer_list, model_time)  ! save all 2d tracers
       if (implicit_bottomfluxes) then
<tracers vertLoc=WAT>
          call g_tracer_set_values(tracer_list,'<name>'  ,'btf',0.0-ergom%btf_<name>/(24.0*3600.0),isd,jsd)
          ! bottom flux is directed into the sediment and in [mol/m2/s]
</tracers>
       endif
       if (NUM_SEDIMENT_LAYERS .gt. 1) then
          ! store some values for use in the sediment model
          ! store fluffy layer concentration
<tracers vertLoc=SED>
          call g_tracer_get_values(tracer_list,'<name>'  ,'field',ergom%temp3darray,isd,jsd,ntau=1)
          ergom%temp3darray(:,:,1) = ergom%<name>
          call g_tracer_set_values(tracer_list,'<name>'  ,'field',ergom%temp3darray,isd,jsd,ntau=1)
</tracers>
          ! store bottom cell concentrations of pore water species
<tracers vertLoc=WAT; isInPorewater=1>
          call g_tracer_get_values(tracer_list,'sed_<name>'  ,'field',ergom%temp3darray,isd,jsd,ntau=1)
          do j = jsc, jec
             do i = isc, iec
                if (1 .gt. grid_kmt(i,j)) cycle
                ergom%temp3darray(i,j,1) = ergom%<name>(i,j,max(1,grid_kmt(i,j)),tau)
             enddo
          enddo
          call g_tracer_set_values(tracer_list,'sed_<name>'  ,'field',ergom%temp3darray,isd,jsd,ntau=1)
</tracers>
          ! store height of bottom cell
          call g_tracer_get_values(tracer_list,'sed_cellheights'  ,'field',ergom%temp3darray,isd,jsd,ntau=1)
          do j = jsc, jec
             do i = isc, iec
                if (1 .gt. grid_kmt(i,j)) cycle
                ergom%temp3darray(i,j,1) = dzt(i,j,max(1,grid_kmt(i,j)))
             enddo
          enddo
          call g_tracer_set_values(tracer_list,'sed_cellheights'  ,'field',ergom%temp3darray,isd,jsd,ntau=1)
          ! store temperature and salinity
          call g_tracer_get_values(tracer_list,'sed_2d_params'  ,'field',ergom%temp3darray,isd,jsd,ntau=1)
          do j = jsc, jec
             do i = isc, iec
                if (1 .gt. grid_kmt(i,j)) cycle
                ergom%temp3darray(i,j,2) = temp(i,j,max(1,grid_kmt(i,j)))
                ergom%temp3darray(i,j,3) = salt(i,j,max(1,grid_kmt(i,j)))
             enddo
          enddo
          call g_tracer_set_values(tracer_list,'sed_2d_params'  ,'field',ergom%temp3darray,isd,jsd,ntau=1)
       endif
       call g_tracer_set_values(tracer_list,'chl',     'field',ergom%bio_opacity/gamma_1, isd,jsd,ntau=1)
       call g_tracer_set_values(tracer_list,'opacity_bio','field',ergom%bio_opacity, isd,jsd,ntau=1)
<constants dependsOn=xyz>
       call g_tracer_set_values(tracer_list,'<name>','field',ergom%<name>, isd,jsd,ntau=1)
</constants>
<constants dependsOn=xyzt>
       call g_tracer_set_values(tracer_list,'<name>','field',ergom%<name>, isd,jsd,ntau=1)
</constants>
<auxiliaries vertLoc=WAT; isUsedElsewhere=1>
       call g_tracer_set_values(tracer_list,'<name>','field',ergom%<name>, isd,jsd,ntau=1)
</auxiliaries>

!   Diagnostics

<tracers vertLoc=WAT; riverDep/=0>
       call g_tracer_get_pointer(tracer_list,'<name>','runoff_tracer_flux',runoff_flux_<name>)
       call g_tracer_get_pointer(tracer_list,'<name>','drydep',diffusive_flux_<name>)
</tracers>
<tracers vertLoc=WAT; atmosDep/=0>
       call g_tracer_get_pointer(tracer_list,'<name>','wetdep',wet_<name>)
</tracers>

<tracers vertLoc=WAT; riverDep/=0>
       if (ergom%id_runoff_flux_<name> .gt. 0)      &
          used = send_data(ergom%id_runoff_flux_<name>, runoff_flux_<name>,     &
          model_time, rmask = grid_tmask(:,:,1), &
          is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
       if (ergom%id_runoff_flux_<name> .gt. 0)      &
          used = send_data(ergom%id_diffusive_flux_<name>, diffusive_flux_<name>,     &
          model_time, rmask = grid_tmask(:,:,1), &
          is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
</tracers>

<tracers vertLoc=WAT; atmosDep/=0>
       if (ergom%id_dep_wet_<name> .gt. 0)     &
          used = send_data(ergom%id_dep_wet_<name>, wet_<name>,                 &
          model_time, rmask = grid_tmask(:,:,1), &
          is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
</tracers>

<auxiliaries vertLoc=WAT; isUsedElsewhere=0>
       if (ergom%id_<name> .gt. 0) &
          used = send_data(ergom%id_<name>, ergom%<name>,                 &
          model_time, rmask = grid_tmask(:,:,:), &
          is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
</auxiliaries>

<auxiliaries vertLoc/=WAT>
       if (ergom%id_<name> .gt. 0) &
          used = send_data(ergom%id_<name>, ergom%<name>,                 &
          model_time, rmask = grid_tmask(:,:,1), &
          is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
</auxiliaries>


<processes vertLoc=WAT>
       if (ergom%id_<name> .gt. 0) &
          used = send_data(ergom%id_<name>, ergom%<name>,                 &
          model_time, rmask = grid_tmask(:,:,:), &
          is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
</processes>

<processes vertLoc/=WAT>
       if (ergom%id_<name> .gt. 0) &
          used = send_data(ergom%id_<name>, ergom%<name>,                 &
          model_time, rmask = grid_tmask(:,:,1), &
          is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
</processes>

       if (ergom%id_irr_inst .gt. 0)           &
            used = send_data(ergom%id_irr_inst,     ergom%irr_inst,             &
            model_time, rmask = grid_tmask(:,:,:),&
            is_in=isc, js_in=jsc, ks_in=1,ie_in=iec, je_in=jec, ke_in=nk)

       call mpp_clock_end(id_output)

    endif

    call mpp_clock_end(id_source)

    return
  end subroutine cgt_bio_timestep

! <SUBROUTINE NAME="cgt_bio_timestep_sed">
  !  <OVERVIEW>
  !   Update tracer concentration fields in the sediment due to the source/sink contributions.
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Ecosystem model code automatically generated by Code Generation Tool - see www.ergom.net
  !  </DESCRIPTION>
  ! </SUBROUTINE>
  subroutine cgt_bio_timestep_sed(tracer_list,Temp,Salt, rho_dzt,dzt,xt,yt,hblt_depth,&
       ilb,jlb,tau,dt,grid_dat,model_time,nbands,max_wavelength_band, sw_pen_band,opacity_band, &
       current_wave_stress,use_saved_rates,intermediate,diff_cbt)

    type(g_tracer_type),            pointer    :: tracer_list
    real, dimension(ilb:,jlb:,:),   intent(in) :: Temp,Salt,rho_dzt,dzt
    real, dimension(ilb:,jlb:),     intent(in) :: xt,yt
    real, dimension(ilb:,jlb:),     intent(in) :: hblt_depth
    integer,                        intent(in) :: ilb,jlb,tau
    real,                           intent(in) :: dt
    real, dimension(ilb:,jlb:),     intent(in) :: grid_dat
    type(time_type),                intent(in) :: model_time

    integer,                        intent(in) :: nbands
    real, dimension(:),             intent(in) :: max_wavelength_band
    real, dimension(:,ilb:,jlb:),   intent(in) :: sw_pen_band
    real, dimension(:,ilb:,jlb:,:), intent(in) :: opacity_band
    real, dimension(ilb:,jlb:),     intent(in) :: current_wave_stress
    logical,                        intent(in) :: use_saved_rates
    logical,                        intent(in) :: intermediate
    real, dimension(ilb:,jlb:,:,:),optional,intent(in) :: diff_cbt

    ! abiotic parameters for biological processes
    real :: cgt_temp                ! potential temperature     [Celsius]
    real :: cgt_sali                ! salinity                  [g/kg]
    real :: cgt_light               ! downward light flux (PAR) [W/m2]
    real :: cgt_cellheight          ! cell height               [m]
    real :: cgt_density             ! density                   [kg/m3]
    real :: cgt_bottomdepth         ! bottom depth              [m]
    real :: cgt_timestep            ! timestep in days          [days]
    real :: cgt_current_wave_stress ! combined bottom stress of waves and current [N/m2]
    real :: cgt_longitude           ! longitude of grid cell    [deg]
    real :: cgt_latitude            ! latitude of grid cell     [deg]
    real :: cgt_diffusivity         ! diffusivity at cell bottom [m*2/s]
    integer :: cgt_iteration        ! current number of iteration in the iterative loop
    integer :: cgt_year, cgt_dayofyear, cgt_mymonth, cgt_myday, cgt_myhour, cgt_myminute, cgt_mysecond ! some date/time variables
    real :: cgt_hour                ! fractional hour (0..23.9999)
    real :: cgt_in_sediment         ! whether we are below the sediment surface
    real :: sed_porosity, sed_tortuosity ! porosity and tortuosity of the sediment

    character(len=fm_string_len), parameter :: sub_name = 'cgt_bio_timestep_sed'
    integer :: isc,iec, jsc,jec,isd,ied,jsd,jed,nk,ntau, i, j, k, kk , kblt, n, m, nband
    real    :: light_itemp, light_ntemp
    real, dimension(:,:,:) ,pointer :: grid_tmask
    integer, dimension(:,:),pointer :: mask_coast,grid_kmt
    logical                         :: used

    integer                         :: number_of_loop

    real                            :: temp1
    real                            :: temp2
    real                            :: temp3
    real                            :: temp4
    real                            :: temp5
    real                            :: temp6
    real                            :: temp7
    real                            :: temp8
    real                            :: temp9

<tracers>
    real                            :: <name>   ! <description>
  <limitations>
    real                            :: <name>
  </limitations>
</tracers>

<constants dependsOn=xyz>
    real                            :: <name>   ! <description>
</constants>
<constants dependsOn=xyzt>
    real                            :: <name>   ! <description>
</constants>

<auxiliaries>
    real                            :: <name>   ! <description>
</auxiliaries>

<processes>
    real                            :: <name>   ! <description>
</processes>

<processes isStiff/=0>
    real                            :: after_patankar_<name>
</processes>

    ! the following variables are used for the positive-definite scheme
    real :: timestep_fraction           ! ratio between allowed timestep and attempted timestep
    real :: timestep_fraction_new       ! ratio between timestep allowed by each tracer and attempted timestep
    real :: fraction_of_total_timestep  ! ratio between remaining timestep and original timestep
<tracers>
    real :: change_of_<name>            ! possible change during remaining timestep
</tracers>
<processes>
    real :: total_rate_<name>           ! total process rate accumulated during all sub-timesteps
</processes>
    integer :: which_tracer_exhausted   ! stores the number of the exhausted tracer

    call mpp_clock_begin(id_sedsource)
    call mpp_clock_begin(id_sedpreloop)

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,&
         grid_tmask=grid_tmask,grid_mask_coast=mask_coast,grid_kmt=grid_kmt)

    !Get the pointers to the prognostic tracers after applying the physics.
<tracers vertLoc=WAT; isInPorewater=1>
    call g_tracer_get_pointer(tracer_list,'<name>'  ,'field',ergom%<name>)
    call g_tracer_get_values(tracer_list,'sed_<name>'  ,'field',ergom%sed_<name>,isd,jsd,ntau=1)
    do j = jsc, jec
       do i = isc, iec
          if (1 .gt. grid_kmt(i,j)) cycle
          ergom%sed_<name>(i,j,1) = ergom%<name>(i,j,grid_kmt(i,j),tau)
       enddo
    enddo
</tracers>
<tracers vertLoc=SED>
    call g_tracer_get_values(tracer_list,'<name>'  ,'field',ergom%sed_<name>,isd,jsd,ntau=1)
</tracers>

    !Get the values of the diagnostic tracers where important properties for the sediment are stored.
    call g_tracer_get_values(tracer_list,'sed_cellheights'            ,'field',ergom%sed_cellheights          ,isd,jsd,ntau=1 ) ! cell heights in the sediment [m]
    do j = jsc, jec
       do i = isc, iec
          if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
          ergom%sed_cellheights(i,j,1) = dzt(i,j,grid_kmt(i,j))
          ergom%sed_2d_params(i,j,2) = temp(i,j,grid_kmt(i,j))
          ergom%sed_2d_params(i,j,3) = salt(i,j,grid_kmt(i,j))
          ergom%sed_2d_params(i,j,4) = 0.0
          do k=1,grid_kmt(i,j)
             ergom%sed_2d_params(i,j,4) = ergom%sed_2d_params(i,j,4) + dzt(i,j,k)
          enddo
       enddo
    enddo
    call g_tracer_get_values(tracer_list,'sed_diffusivity_porewater'  ,'field',ergom%sed_diffusivity_porewater,isd,jsd,ntau=1 ) ! bioturbation-generated diffusivity for pore water species [m2/s]
    call g_tracer_get_values(tracer_list,'sed_diffusivity_solids'     ,'field',ergom%sed_diffusivity_solids   ,isd,jsd,ntau=1 ) ! bioturbation-generated diffusivity for solid species [m2/s]
    call g_tracer_get_values(tracer_list,'sed_2d_params'              ,'field',ergom%sed_2d_params            ,isd,jsd,ntau=1 ) ! sed_inert_deposition[m/s], temperature[degC], salinity[g/kg], bottomdepth[m]
    call g_tracer_get_values(tracer_list,'sed_inert_ratio'            ,'field',ergom%sed_inert_ratio          ,isd,jsd,ntau=1 ) ! volume fraction of bioinert material in the sediments [1]
    !Get the values of the spatially varying constants and elsewhere-used auxiliaries
<constants dependsOn=xyz>
    call g_tracer_get_values(tracer_list,'<name>'  ,'field',ergom%<name>,isd,jsd,ntau=1)
</constants>
<constants dependsOn=xyzt>
    call g_tracer_get_values(tracer_list,'<name>'  ,'field',ergom%<name>,isd,jsd,ntau=1)
    call data_override('OCN', '<trimName>', ergom%<trimName>(isc:iec,jsc:jec,1:nk), model_time)
</constants>
<auxiliaries vertLoc=WAT; isUsedElsewhere=1>
    call g_tracer_get_values(tracer_list,'<name>'  ,'field',ergom%<name>,isd,jsd,ntau=1)
</auxiliaries>


    !------------------------------------
    ! STEP 1: calculate total element concentrations
    !------------------------------------
    !------------------------------------
    ! STEP 1.1: calculate total element concentrations in pore water
    !------------------------------------
    do k = 2,NUM_SEDIMENT_LAYERS
       do j = jsc, jec
          do i = isc, iec
             if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
          <celements>
             ergom%sed_<total>(i,j,k) = &
             <containingTracers vertLoc=WAT; isInPorewater=1>
                max(0.0,ergom%sed_<ct>(i,j,k))*<ctAmount> + &
             </containingTracers>
                0.0
          </celements>
          enddo
       enddo
    enddo

    !------------------------------------
    ! STEP 1.2: calculate total element concentrations in solid phase
    !------------------------------------
    do k = 2,NUM_SEDIMENT_LAYERS
       do j = jsc, jec
          do i = isc, iec
             if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
             ! TODO: Modification by Daniel Neumann; because we are in the sediment model section
             !  from:
             !    ergom%<ct>(i,j)*<ctAmount> + &
             !  to:
             !    ergom%sed_<ct>(i,j,k)*<ctAmount> + &
          <celements>
             ergom%sed_<totalBottom>(i,j,k) = &
             <containingTracers vertLoc=SED>
                ergom%sed_<ct>(i,j,k)*<ctAmount> + &
             </containingTracers>
                0.0
          </celements>
          enddo
       enddo
    enddo

    !------------------------------------
    ! STEP 2: calculate space-independent forcing variables
    !------------------------------------

    call get_date(model_time,cgt_year,cgt_mymonth,cgt_myday,cgt_myhour,cgt_myminute,cgt_mysecond)
    cgt_hour = cgt_myhour + cgt_myminute/60 + cgt_mysecond/3600
    cgt_myhour=days_in_year(model_time)
    if (cgt_myhour .eq. 366) then
       if (cgt_mymonth .eq.  1) cgt_dayofyear = cgt_myday
       if (cgt_mymonth .eq.  2) cgt_dayofyear = cgt_myday +  31
       if (cgt_mymonth .eq.  3) cgt_dayofyear = cgt_myday +  60
       if (cgt_mymonth .eq.  4) cgt_dayofyear = cgt_myday +  91
       if (cgt_mymonth .eq.  5) cgt_dayofyear = cgt_myday + 121
       if (cgt_mymonth .eq.  6) cgt_dayofyear = cgt_myday + 152
       if (cgt_mymonth .eq.  7) cgt_dayofyear = cgt_myday + 182
       if (cgt_mymonth .eq.  8) cgt_dayofyear = cgt_myday + 213
       if (cgt_mymonth .eq.  9) cgt_dayofyear = cgt_myday + 244
       if (cgt_mymonth .eq. 10) cgt_dayofyear = cgt_myday + 274
       if (cgt_mymonth .eq. 11) cgt_dayofyear = cgt_myday + 305
       if (cgt_mymonth .eq. 12) cgt_dayofyear = cgt_myday + 335
    else
       if (cgt_mymonth .eq.  1) cgt_dayofyear = cgt_myday
       if (cgt_mymonth .eq.  2) cgt_dayofyear = cgt_myday +  31
       if (cgt_mymonth .eq.  3) cgt_dayofyear = cgt_myday +  59
       if (cgt_mymonth .eq.  4) cgt_dayofyear = cgt_myday +  90
       if (cgt_mymonth .eq.  5) cgt_dayofyear = cgt_myday + 120
       if (cgt_mymonth .eq.  6) cgt_dayofyear = cgt_myday + 151
       if (cgt_mymonth .eq.  7) cgt_dayofyear = cgt_myday + 181
       if (cgt_mymonth .eq.  8) cgt_dayofyear = cgt_myday + 212
       if (cgt_mymonth .eq.  9) cgt_dayofyear = cgt_myday + 243
       if (cgt_mymonth .eq. 10) cgt_dayofyear = cgt_myday + 273
       if (cgt_mymonth .eq. 11) cgt_dayofyear = cgt_myday + 304
       if (cgt_mymonth .eq. 12) cgt_dayofyear = cgt_myday + 334
    endif

    !------------------------------------
    ! STEP 3: initialize isZIntegral=1 auxiliaries with zero
    !------------------------------------
    ! not in sediment

    !------------------------------------
    ! STEP 4: Pre-loop where the isZIntegral=1 auxiliaries are calculated
    !------------------------------------
    ! not in sediment

    call mpp_clock_end(id_sedpreloop)
    call mpp_clock_begin(id_sedmainloop)

    !------------------------------------
    ! STEP 5: initializations before the main loop
    !------------------------------------

    !------------------------------------
    ! STEP 5.1: initialize isZIntegral=1 auxiliaries which are calculated after the process rates with zero
    !------------------------------------
    ! not in sediment

    !------------------------------------
    ! STEP 5.2: initialize the arrays which cumulate the change of vertLoc=FIS tracers with zero
    !------------------------------------
    ! not in sediment

    !------------------------------------
    ! STEP 6: Main loop
    !------------------------------------

    do k = 2, NUM_SEDIMENT_LAYERS
       do j = jsc, jec
          do i = isc, iec
             if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
             !------------------------------------
             ! STEP 6.1: prepare abiotic parameters
             !------------------------------------
             cgt_bottomdepth= ergom%sed_2d_params(i,j,4)
             do kk=2,k
                cgt_bottomdepth = cgt_bottomdepth + ergom%sed_cellheights(i,j,kk)
             enddo
             cgt_temp       = ergom%sed_2d_params(i,j,2)                                         ! potential temperature     [Celsius]
             cgt_sali       = ergom%sed_2d_params(i,j,3)                                         ! salinity                  [g/kg]
             cgt_diffusivity= 0.0
             cgt_light      = 0.0                                                                ! light intensity           [W/m2]
             cgt_cellheight = ergom%sed_cellheights(i,j,k)*(1.0-ergom%sed_inert_ratio(i,j,k))    ! cell height               [m]
             cgt_density    = 1035.0                                                             ! density [kg/m3]
             cgt_timestep   = dt/(24*3600)                                                       ! timestep in days
             cgt_longitude  = xt(i,j)                                                            ! longitude [deg]
             cgt_latitude   = yt(i,j)                                                            ! latitude [deg]
             cgt_in_sediment = 1.0                                                               ! we are below sediment surface
             cgt_current_wave_stress=current_wave_stress(i,j)

             !------------------------------------
             ! STEP 6.2: load tracer values
             !------------------------------------
<constants dependsOn=xyz>
             <name> = ergom%<name>(i,j,grid_kmt(i,j)) ! <description>
</constants>
<constants dependsOn=xyzt>
             <name> = ergom%<name>(i,j,grid_kmt(i,j)) ! <description>
</constants>
<auxiliaries vertLoc=WAT; isUsedElsewhere=1>
             <name> = ergom%<name>(i,j,grid_kmt(i,j)) ! <description>
</auxiliaries>
<tracers vertLoc=WAT; isInPorewater=1>
             <name> = ergom%sed_<name>(i,j,k) ! <description>
</tracers>
<tracers vertLoc=WAT; isInPorewater=0>
             <name> = 0.0
</tracers>
<tracers vertLoc=WAT; isPositive=1; isInPorewater=1>
             <name>       = max(<name>,0.0)
</tracers>

<tracers vertLoc=SED>
                <name> = ergom%sed_<name>(i,j,k) ! <description>
</tracers>
<tracers vertLoc=SED; isPositive=1>
                <name> = max(<name>,0.0)
</tracers>

             !------------------------------------
             ! STEP 6.3: opacity calculation
             !------------------------------------
             ! not in sediment

             !------------------------------------
             ! STEP 6.4: calculate auxiliaries
             !------------------------------------
             !------------------------------------
             ! STEP 6.4.1: get isZIntegral values from pre-loop
             !------------------------------------
             ! not in sediment
<auxiliaries isZIntegral=1; calcAfterProcesses=0>
             ! <description> :
             <name> = 0.0
</auxiliaries>

             !------------------------------------
             ! STEP 6.4.2: calculate isZGradient auxiliaries
             !------------------------------------
             ! not in sediment
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=1>
             ! <description> :
             <name> = 0.0

</auxiliaries>

             !------------------------------------
             ! STEP 6.4.3: vertLoc=WAT auxiliaries calclated iteratively
             !------------------------------------

             ! initialize auxiliaries for iterative loop
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=0; iterations/=0>
             <name> = <iterInit>
</auxiliaries>

             ! iterative loop follows
             do cgt_iteration=1,<maxIterations>
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=0; iterations/=0>
                if (cgt_iteration .le. <iterations>) then
                   ! <description> :
                   temp1 = <temp1>
                   temp2 = <temp2>
                   temp3 = <temp3>
                   temp4 = <temp4>
                   temp5 = <temp5>
                   temp6 = <temp6>
                   temp7 = <temp7>
                   temp8 = <temp8>
                   temp9 = <temp9>
                   <name> = <formula>
                endif
</auxiliaries>
             enddo

             !------------------------------------
             ! STEP 6.4.4: vertLoc=WAT auxiliaries calclated normally (not iteratively)
             !------------------------------------
<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isZGradient=0; iterations=0>
             ! <description> :
             temp1 = <temp1>
             temp2 = <temp2>
             temp3 = <temp3>
             temp4 = <temp4>
             temp5 = <temp5>
             temp6 = <temp6>
             temp7 = <temp7>
             temp8 = <temp8>
             temp9 = <temp9>
             <name> = <formula>

</auxiliaries>

             !------------------------------------
             ! STEP 6.4.5: vertLoc=SED auxiliaries calclated iteratively
             !------------------------------------
             ! initialize auxiliaries for iterative loop
<auxiliaries vertLoc=SED; calcAfterProcesses=0; isZIntegral=0; iterations/=0>
             <name> = <iterInit>
</auxiliaries>

             ! iterative loop follows
             do cgt_iteration=1,<maxIterations>
<auxiliaries vertLoc=SED; calcAfterProcesses=0; isZIntegral=0; iterations/=0>
                if (cgt_iteration .l. <iterations>) then
                   ! <description> :
                   temp1 = <temp1>
                   temp2 = <temp2>
                   temp3 = <temp3>
                   temp4 = <temp4>
                   temp5 = <temp5>
                   temp6 = <temp6>
                   temp7 = <temp7>
                   temp8 = <temp8>
                   temp9 = <temp9>
                   <name> = <formula>
                endif
</auxiliaries>
             enddo

             !------------------------------------
             ! STEP 6.4.6: vertLoc=SED auxiliaries calclated normally (not iteratively)
             !------------------------------------

<auxiliaries vertLoc=SED; calcAfterProcesses=0; isZIntegral=0; iterations=0>
             ! <description> :
             temp1 = <temp1>
             temp2 = <temp2>
             temp3 = <temp3>
             temp4 = <temp4>
             temp5 = <temp5>
             temp6 = <temp6>
             temp7 = <temp7>
             temp8 = <temp8>
             temp9 = <temp9>
             <name> = <formula>

</auxiliaries>

             !------------------------------------
             ! STEP 6.4.7: vertLoc=SUR auxiliaries calclated iteratively
             !------------------------------------
             ! not in sediment

             !------------------------------------
             ! STEP 6.4.8: vertLoc=SUR auxiliaries calclated normally (not iteratively)
             !------------------------------------
             ! not in sediment
<auxiliaries vertLoc=SUR; calcAfterProcesses=0>
                <name> = 0.0
</auxiliaries>

             !------------------------------------
             ! STEP 6.5: output of auxiliaries
             !------------------------------------
             if (.not. intermediate) then

<auxiliaries vertLoc=WAT; calcAfterProcesses=0; isUsedElsewhere=0>
                if (ergom%id_sed_<name> .gt. 0) then
                  ergom%sed_<name>(i,j,k) = <name>
                endif
</auxiliaries>

<auxiliaries vertLoc=SED; calcAfterProcesses=0>
                if (ergom%id_sed_<name> .gt. 0) then
                   ergom%sed_<name>(i,j,k) = <name>
                endif
</auxiliaries>
             endif

             !------------------------------------
             ! STEP 6.6: calculate process limitations
             !------------------------------------
<tracers vertLoc=WAT>
  <limitations>
             <name> = <formula>
  </limitations>
</tracers>

<tracers vertLoc=SED>
  <limitations>
             <name> = <formula>
  </limitations>
</tracers>

             !------------------------------------
             ! STEP 6.7
             !
             !-- POSITIVE-DEFINITE SCHEME --------
             !-- means the following steps will be repeated as often as nessecary
             !------------------------------------
             number_of_loop = 1
             fraction_of_total_timestep = 1.0   ! how much of the original timestep is remaining
<processes>
             total_rate_<name>          = 0.0
</processes>

             do while (cgt_timestep .gt. 0.0)

                !------------------------------------
                ! STEP 6.7.1: calculate process rates
                !------------------------------------
<processes vertLoc=WAT>
                ! <description> :
                <name> = <turnover>
                <name> = max(<name>,0.0)

</processes>

<processes vertLoc=SED>
                ! <description> :
                <name> = <turnover>
                <name> = max(<name>,0.0)

</processes>

                !------------------------------------
                ! STEP 6.7.2: save the process rates if intermediate time step (Runge-Kutta)
                !------------------------------------
                if ((number_of_loop == 1) .and. intermediate) then
<processes vertLoc=WAT; isInPorewater=1>
                   ergom%saved_rate_sed_<name>(i,j,k)=ergom%saved_rate_sed_<name>(i,j,k)+<name>
</processes>
<processes vertLoc=SED; isInPorewater=1>
                   ergom%saved_rate_sed_<name>(i,j,k)=ergom%saved_rate_sed_<name>(i,j,k)+<name>
</processes>
                endif

                !------------------------------------
                ! STEP 6.7.3: use saved rates if they apply (Runge-Kutta)
                !------------------------------------
                if (use_saved_rates) then
<processes vertLoc=WAT; isInPorewater=1>
                   if ((<name> > 0.0) .and. (ergom%saved_rate_sed_<name>(i,j,k) >= 0.0)) then
                      <name> = ergom%saved_rate_sed_<name>(i,j,k)
                   endif
</processes>
<processes vertLoc=WAT; isStiff/=0; isInPorewater=1>
                   if ((<name> > 0.0) .and. (ergom%saved_rate_sed_<name>(i,j,k) >= 0.0)) then
                      <name> = <name> * patankar_modification_<stiffTracer>(i,j,k)
                   endif
</processes>
<processes vertLoc=SED; isInPorewater=1>
                   if ((<name> > 0.0) .and. (ergom%saved_rate_sed_<name>(i,j,k) >= 0.0)) then
                      <name> = ergom%saved_rate_sed_<name>(i,j,k)
                   endif
</processes>
<processes vertLoc=SED; isStiff/=0; isInPorewater=1>
                   if ((<name> > 0.0) .and. (ergom%saved_rate_sed_<name>(i,j,k) >= 0.0)) then
                      <name> = <name> * patankar_modification_<stiffTracer>(i,j,k)
                   endif
</processes>
                endif

                !------------------------------------
                ! STEP 6.7.4: apply Patankar limitations
                !------------------------------------

<processes vertLoc=WAT; isStiff/=0; isInPorewater=1>
                ! <description> :
                after_patankar_<name> = <name> * <stiffFactor>

</processes>
<processes vertLoc=SED; isStiff/=0; isInPorewater=1>
                ! <description> :
                after_patankar_<name> = <name> * <stiffFactor>

</processes>

<processes vertLoc=WAT; isStiff/=0; isInPorewater=1>
                ! <description> :
                <name> = after_patankar_<name>

</processes>
<processes vertLoc=SED; isStiff/=0; isInPorewater=1>
                ! <description> :
                <name> = after_patankar_<name>

</processes>

                !------------------------------------
                ! STEP 6.7.5: calculate possible euler-forward change (in a full timestep)
                !------------------------------------
<tracers>
                change_of_<name> = 0.0
</tracers>

<tracers vertLoc=WAT; hasTimeTendenciesVertLoc=WAT>

                change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                <timeTendencies vertLoc=WAT>
                   <timeTendency> & ! <description>
                </timeTendencies>
                )
</tracers>
<tracers vertLoc=FIS; hasTimeTendenciesVertLoc=WAT>

                change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                <timeTendencies vertLoc=WAT>
                   <timeTendency> & ! <description>
                </timeTendencies>
                )
</tracers>
<tracers vertLoc=WAT; hasTimeTendenciesVertLoc=SED>

                change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                <timeTendencies vertLoc=SED>
                   <timeTendency> & ! <description>
                </timeTendencies>
                )
</tracers>

<tracers vertLoc=SED; hasTimeTendencies>

                change_of_<name> = change_of_<name> + cgt_timestep*(0.0 &
                <timeTendencies>
                   <timeTendency> & ! <description>
                </timeTendencies>
                )
</tracers>

                !------------------------------------
                ! STEP 6.7.6: calculate maximum fraction of the timestep before some tracer gets exhausted
                !------------------------------------

                timestep_fraction = 1.0
                which_tracer_exhausted = -1

                ! find the tracer which is exhausted after the shortest period of time

                ! in the water column
<tracers vertLoc=WAT; isInPorewater=1; isPositive=1>

                ! check if tracer <name> was exhausted from the beginning and is still consumed
                if ((ergom%sed_<name>(i,j,k) .le. 0.0) .and. (change_of_<name>  .lt. 0.0)) then
                   timestep_fraction = 0.0
                   which_tracer_exhausted = <numTracer>
                endif
                ! check if tracer <name> was present, but got exhausted
                if ((ergom%sed_<name>(i,j,k) .gt. 0.0) .and. (ergom%sed_<name>(i,j,k) + change_of_<name> .lt. 0.0)) then
                   timestep_fraction_new = ergom%sed_<name>(i,j,k) / (0.0 - change_of_<name>)
                   if (timestep_fraction_new .le. timestep_fraction) then
                      which_tracer_exhausted = <numTracer>
                      timestep_fraction = timestep_fraction_new
                   endif
                endif
</tracers>

                ! in the bottom layer
<tracers vertLoc=SED; isPositive=1>

                ! check if tracer <name> was exhausted from the beginning and is still consumed
                if ((ergom%sed_<name>(i,j,k) .le. 0.0) .and. (change_of_<name> .lt. 0.0)) then
                   timestep_fraction = 0.0
                   which_tracer_exhausted = <numTracer>
                endif
                ! check if tracer <name> was present, but got exhausted
                if ((ergom%sed_<name>(i,j,k) .gt. 0.0) .and. (ergom%sed_<name>(i,j,k) + change_of_<name> .lt. 0.0)) then
                   timestep_fraction_new = ergom%sed_<name>(i,j,k) / (0.0 - change_of_<name>)
                   if (timestep_fraction_new .le. timestep_fraction) then
                      which_tracer_exhausted = <numTracer>
                      timestep_fraction = timestep_fraction_new
                   endif
                endif
</tracers>

                !------------------------------------
                ! STEP 6.7.7: update the limitations: processes limited by this tracer become zero in the future
                !------------------------------------

<tracers isPositive=1; vertLoc/=FIS>
                if (<numTracer> .eq. which_tracer_exhausted) then
                  <limitations>
                   <name> = 0.0
                  </limitations>
                endif
</tracers>

                !------------------------------------
                ! STEP 6.7.8: apply a Euler-forward timestep with the fraction of the time
                !------------------------------------
                ! in the water column
<tracers vertLoc=WAT; isInPorewater=1>

                ! tracer <name> (<description>):
                ergom%sed_<name>(i,j,k) = ergom%sed_<name>(i,j,k) + change_of_<name> * timestep_fraction
</tracers>
<tracers vertLoc=SED>

                ! tracer <name> (<description>)
                ergom%sed_<name>(i,j,k) = ergom%sed_<name>(i,j,k) + change_of_<name> * timestep_fraction
</tracers>

                !------------------------------------
                ! STEP 6.7.9: save process rates to obtain total rate over full timestep in the end
                !------------------------------------
                if (.not. intermediate) then
<processes vertLoc=WAT; isInPorewater=1>
                   total_rate_<name> = total_rate_<name> + <name> * timestep_fraction * fraction_of_total_timestep
</processes>
<processes vertLoc=SED>
                   total_rate_<name> = total_rate_<name> + <name> * timestep_fraction * fraction_of_total_timestep
</processes>
                endif

                !------------------------------------
                ! STEP 6.7.10: set timestep to remaining timestep only
                !------------------------------------
                cgt_timestep = cgt_timestep * (1.0 - timestep_fraction)                         ! remaining timestep
                fraction_of_total_timestep = fraction_of_total_timestep * (1.0 - timestep_fraction) ! how much of the original timestep is remaining

                if (number_of_loop .gt. 100) then
                   print*,"LOOP: ",number_of_loop
                   print*,"i=",i,", j=",j,", k=",k
                 <tracers isPositive=1>
                   if (<numTracer> .eq. which_tracer_exhausted) then
                     print*,"limiting tracer: <name>"
                     print*,"  <name> = ",<name>
                     print*,"  change_of_<name> = ",change_of_<name>
                   endif
                 </tracers>
                   print*,"tracers:"
                 <tracers>
                   print*,"  <name> = ",<name>
                 </tracers>
                   print*,"changes:"
                 <tracers>
                   print*,"  change_of_<name> = ",change_of_<name>
                 </tracers>
                   print*,"limitations:"
                 <tracers>
                  <limitations>
                   print*,"  <name> = ",<name>
                  </limitations>
                 </tracers>
                   print*,"cgt_timestep: ",cgt_timestep
                   print*,"timestep_fraction: ",timestep_fraction
                   print*,"  "
                   if (number_of_loop .gt. 200) then
                      stop
                   endif
                endif
                number_of_loop=number_of_loop+1

             enddo
             !------------------------------------
             !-- END OF POSITIVE-DEFINITE SCHEME -
             !------------------------------------

             !------------------------------------
             ! STEP 6.8: cumulate change of vertLoc=FIS tracers over the k levels
             !------------------------------------
             ! not in sediment

             if (.not. intermediate) then

                !------------------------------------
                ! STEP 6.9: output of new tracer concentrations
                !------------------------------------

                ! not required in MOM5

                !------------------------------------
                ! STEP 6.10: get total process rates over full time step
                !------------------------------------

<processes vertLoc=WAT; isInPorewater=1>
                <name> = total_rate_<name>
</processes>
<processes vertLoc=SED; isInPorewater=1>
                <name> = total_rate_<name>
</processes>

                !------------------------------------
                ! STEP 6.11: output of process rates
                !------------------------------------
<processes vertLoc=WAT; isInPorewater=1>
                if (ergom%id_Sed_<Name> .gt. 0) then
                  ergom%sed_<name>(i,j,k) = total_rate_<name>
                endif
</processes>
<processes vertLoc=SED; isInPorewater=1>
                if (ergom%id_sed_<name> .gt. 0) then
                   ergom%sed_<name>(i,j,k) = total_rate_<name>
                endif
</processes>

                !------------------------------------
                ! STEP 6.12: calculate "late" auxiliaries
                !------------------------------------
                !------------------------------------
                ! STEP 6.12.1: calculate all but the isZIntegral auxiliaries
                !------------------------------------
<auxiliaries vertLoc=WAT; calcAfterProcesses=1; isZGradient=0>
                ! <description> :
                temp1 = <temp1>
                temp2 = <temp2>
                temp3 = <temp3>
                temp4 = <temp4>
                temp5 = <temp5>
                temp6 = <temp6>
                temp7 = <temp7>
                temp8 = <temp8>
                temp9 = <temp9>
                <name> = <formula>

</auxiliaries>

<auxiliaries vertLoc=SED; calcAfterProcesses=1; isZIntegral=0>
                   ! <description> :
                   temp1 = <temp1>
                   temp2 = <temp2>
                   temp3 = <temp3>
                   temp4 = <temp4>
                   temp5 = <temp5>
                   temp6 = <temp6>
                   temp7 = <temp7>
                   temp8 = <temp8>
                   temp9 = <temp9>
                   <name> = <formula>

</auxiliaries>

                !------------------------------------
                ! STEP 6.12.2: cumulate the isZIntegral auxiliaries over the k levels
                !------------------------------------
                ! not in sediment

                !------------------------------------
                ! STEP 6.13: output of "late" auxiliaries
                !------------------------------------
<auxiliaries vertLoc=WAT; calcAfterProcesses=1>
                if (ergom%id_sed_<name> .gt. 0) then
                  ergom%sed_<name>(i,j,k) = <name>
                endif
</auxiliaries>
<auxiliaries vertLoc=SED; calcAfterProcesses=1; isZIntegral=0>
                if (ergom%id_sed_<name> .gt. 0) then
                   ergom%sed_<name>(i,j,k) = <name>
                endif
</auxiliaries>
<auxiliaries vertLoc=SED; calcAfterProcesses=1; isZIntegral=1>
                if (ergom%id_sed_<name> .gt. 0) then
                   ergom%sed_<name>(i,j,k) = 0.0
                endif
</auxiliaries>

                !---------------------------------------
                ! STEP 6.14: passing vertical velocity and diffusivity to the coupler
                !---------------------------------------
                ! in sediment, we pass molecular diffusivity
<tracers vertLoc=WAT; molDiff/=0; isInPorewater=1>
                ergom%moldiff_<name>(i,j,k)=<molDiff> ! in m2/s
</tracers>

                !---------------------------------------
                ! STEP 6.15: output of parameters for surface flux of gasses
                !---------------------------------------
                ! not in sediment

             endif

          enddo
       enddo
    enddo

    !-----------------------------------
    ! biological timestep has ended
    !-----------------------------------

    call mpp_clock_end(id_sedmainloop)
    if (.not. intermediate) then

       call mpp_clock_begin(id_sedvertmig)

       !-----------------------------------
       ! STEP 7: send surface concentration, solubility and Schmidt number to the coupler for gas exchange
       !-----------------------------------

       ! in the sediment, "movement" is downward advection of vertLoc=SED tracers due to sediment accumulation

       !-----------------------------------
       ! STEP 8: vertical movement of tracers
       !-----------------------------------
       !-----------------------------------
       ! STEP 8.1: calculate total element concentrations in the pore water
       !-----------------------------------
       ! in sediment, we need to calculate this for the vertLoc=SED tracers (solid phase) instead, as they are moved
       do k = 2,NUM_SEDIMENT_LAYERS
          do j = jsc, jec
             do i = isc, iec
                if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
             ! TODO: Modification by Daniel Neumann; because we are in the sediment model section
             !  from:
             !    ergom%<ct>(i,j)*<ctAmount> + &
             !  to:
             !    ergom%sed_<ct>(i,j,k)*<ctAmount> + &
             <celements>
                ergom%sed_<totalBottom>(i,j,k) = &
                <containingTracers vertLoc=SED>
                   ergom%sed_<ct>(i,j,k)*<ctAmount> + &
                </containingTracers>
                   0.0
             </celements>
             enddo
          enddo
       enddo

       do k=2,NUM_SEDIMENT_LAYERS
          do j=jsc,jec
             do i=isc,iec
                ! first, calculate porosity
                sed_porosity = 1.0-ergom%sed_inert_ratio(i,j,k)
                ! second, save cell heights in new array
                ergom%temp3darray(i,j,k) = ergom%sed_cellheights(i,j,k)*max(ergom%sed_inert_ratio(i,j,k),epsilon)
                if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
                ! third, convert all tracer concentrations from mol/m2 to mol/kg in solid phase
<tracers vertLoc=SED>
                ergom%sed_<name>(i,j,k) = ergom%sed_<name>(i,j,k)/(ergom%sed_cellheights(i,j,k)*1035.0)/max(ergom%sed_inert_ratio(i,j,k),epsilon)
</tracers>
                ! fourth, store velocity in m/s
                ergom%vertspeedarray(i,j,k) = 0.0-ergom%sed_2d_params(i,j,1)
             enddo
          enddo
       enddo
       do j=jsc,jec
          do i=isc,iec
             if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
             ergom%vertspeedarray(i,j,1) = 0.0-ergom%sed_2d_params(i,j,1)
             ergom%temp3darray(i,j,1) = ergom%sed_2d_params(i,j,5)
<tracers vertLoc=SED>
             ergom%sed_<name>(i,j,1) = ergom%sed_<name>(i,j,1)/(ergom%sed_2d_params(i,j,5)*1035.0)
</tracers>
          enddo
       enddo

       !-----------------------------------
       ! STEP 8.2: do several vertical movement steps
       !-----------------------------------
       do m = 1, NUM_VMOVE_STEPS
          !-----------------------------------
          ! STEP 8.2.1: store age concentrations in an array which is not modified by movement
          !-----------------------------------
          <celements isAging/=0>
            ergom%old_<aged> = ergom%sed_<agedBottom>(:,:,:)
          </celements>
          !-----------------------------------
          ! STEP 8.2.2: move the age concentrations of the marked elements
          !-----------------------------------
          <tracers childOf/=none; vertLoc=SED; hasCeAged>
             call generic_ERGOM_sedmove(ergom%vertspeedarray, &
                                      ergom%sed_<ceAgedName>(:,:,:), ergom%old_sed_<ceAgedName>(:,:,:), &
                                      ergom%sed_<name>(:,:,:)*<ceAmount>, ergom%sed_<ceTotalName>(:,:,:), &
                                      ergom%temp3darray, ergom%temp3darray(:,:,1), dt/NUM_VMOVE_STEPS, &
                                      isd, ied, jsd, jed, NUM_SEDIMENT_LAYERS)
          </tracers>

          !-----------------------------------
          ! STEP 8.2.3: move the tracers (including tagged tracers) themselves
          !-----------------------------------
          <tracers vertLoc=SED>
             call generic_ERGOM_sedmove(ergom%vertspeedarray, &
                                      ergom%sed_<name>(:,:,:), ergom%sed_<name>(:,:,:), &
                                      ergom%sed_<name>(:,:,:), ergom%sed_<name>(:,:,:), &
                                      ergom%temp3darray, ergom%temp3darray(:,:,1), dt/NUM_VMOVE_STEPS, &
                                      isd, ied, jsd, jed, NUM_SEDIMENT_LAYERS)
           </tracers>

          !-----------------------------------
          ! STEP 8.2.4: calculate the new total marked element concentrations
          !-----------------------------------
          do k = 1, NUM_SEDIMENT_LAYERS
             do j = jsc, jec
                do i = isc, iec
                   if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
                ! TODO: Modification by Daniel Neumann; because we are in the sediment model section
                !  from:
                !    ergom%<ct>(i,j)*<ctAmount> + &
                !  to:
                !    ergom%sed_<ct>(i,j,k)*<ctAmount> + &
                <celements>
                   ergom%sed_<totalBottom>(i,j,k) = &
                   <containingTracers vertLoc=SED>
                      ergom%sed_<ct>(i,j,k)*<ctAmount> + &
                   </containingTracers>
                      0.0
                </celements>
                enddo
             enddo
          enddo
       enddo
       ! backtransformation to mol/m2
       do k=2,NUM_SEDIMENT_LAYERS
          do j=jsc,jec
             do i=isc,iec
                if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
<tracers vertLoc=SED>
                ergom%sed_<name>(i,j,k) = ergom%sed_<name>(i,j,k)*(ergom%sed_cellheights(i,j,k)*1035.0)*max(ergom%sed_inert_ratio(i,j,k),epsilon)
</tracers>
             enddo
          enddo
       enddo
       do j=jsc,jec
          do i=isc,iec
             if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
<tracers vertLoc=SED>
             ergom%sed_<name>(i,j,1) = ergom%sed_<name>(i,j,1)*(ergom%sed_2d_params(i,j,5)*1035.0)
</tracers>
          enddo
       enddo

       !-----------------------------------
       ! STEP 9: vertical diffusion of tracers
       !-----------------------------------
       ! first, calculate total diffusivity (molecular + bioturbation)
       do k=2,NUM_SEDIMENT_LAYERS
          do j=jsc,jec
             do i=isc,iec
                ! first, calculate porosity
                sed_porosity = 1.0-ergom%sed_inert_ratio(i,j,k)
                ! second, multiply cell heights by porosity
                ergom%temp3darray(i,j,k) = ergom%sed_cellheights(i,j,k)*max(sed_porosity,epsilon)
                if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
                ! third, calculate tortuosity from porosity
                sed_tortuosity = sqrt(1.0-2.02*log(sed_porosity))
                ! fourth, calculate effective molecular diffusivity
<tracers vertLoc=WAT; isInPorewater=1; molDiff/=0>
                ergom%moldiff_<name>(i,j,k) = (ergom%moldiff_<name>(i,j,k)/(sed_tortuosity*sed_tortuosity)+ ergom%sed_diffusivity_porewater(i,j,k))*power(sed_porosity,2.0)
</tracers>
<tracers vertLoc=WAT; isInPorewater=1; molDiff=0>
                ergom%moldiff_<name>(i,j,k) = ergom%sed_diffusivity_porewater(i,j,k)*power(sed_porosity,2.0)
</tracers>
             enddo
          enddo
       enddo
       do j=jsc,jec
          do i=isc,iec
             if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
             ergom%temp3darray(i,j,1) = dzt(i,j,grid_kmt(i,j))
          enddo
       enddo
       !-----------------------------------
       ! STEP 9.1: calculate total element concentrations in the water column
       !-----------------------------------
       do k = 2,NUM_SEDIMENT_LAYERS
          do j = jsc, jec
             do i = isc, iec
                if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
                ! TODO: Modification by Daniel Neumann; because we are in the sediment model section
                !  from:
                !    ergom%<ct>(i,j,k,tau)*<ctAmount> + &
                !  to:
                !    ergom%sed_<ct>(i,j,k)*<ctAmount> + &
             <celements>
                ergom%sed_<total>(i,j,k) = &
                <containingTracers vertLoc=WAT; isInPorewater=1>
                   ergom%sed_<ct>(i,j,k)*<ctAmount> + &
                </containingTracers>
                   0.0
             </celements>
             enddo
          enddo
       enddo
       !-----------------------------------
       ! STEP 9.2: do several vertical diffusion steps
       !-----------------------------------
       do m = 1, NUM_VMOVE_STEPS
          !-----------------------------------
          ! STEP 9.2.1: store age concentrations in an array which is not modified by diffusion
          !-----------------------------------
          <celements isAging/=0>
            ergom%old_<aged>(:,:,:,1) = ergom%sed_<aged>(:,:,:)
          </celements>
          !-----------------------------------
          ! STEP 9.2.2: diffuse the age concentrations of the marked elements
          !-----------------------------------
          <tracers childOf/=none; vertLoc=WAT; isInPorewater=1; hasCeAged>
             call generic_ERGOM_seddiff(ergom%moldiff_<name>, &
                                      ergom%sed_<ceAgedName>(:,:,:), ergom%old_<ceAgedName>(:,:,:,1),&
                                      ergom%sed_<name>(:,:,:)*<ceAmount>, ergom%sed_<ceTotalName>(:,:,:), &
                                      ergom%temp3darray(:,:,:), ergom%sed_2d_params(:,:,6), dt/NUM_VMOVE_STEPS, &
                                      isd, ied, jsd, jed, NUM_SEDIMENT_LAYERS)
           </tracers>

          !-----------------------------------
          ! STEP 9.2.3: move the tracers (including tagged tracers) themselves
          !-----------------------------------

          <tracers vertLoc=WAT; isInPorewater=1>
             call generic_ERGOM_seddiff(ergom%moldiff_<name>, &
                                      ergom%sed_<name>(:,:,:), ergom%sed_<name>(:,:,:), &
                                      ergom%sed_<name>(:,:,:), ergom%sed_<name>(:,:,:), &
                                      ergom%temp3darray(:,:,:), ergom%sed_2d_params(:,:,6), dt/NUM_VMOVE_STEPS, &
                                      isd, ied, jsd, jed, NUM_SEDIMENT_LAYERS)
          </tracers>

          !-----------------------------------
          ! STEP 9.2.4: calculate the new total marked element concentrations
          !-----------------------------------
          do k = 2,NUM_SEDIMENT_LAYERS
             do j = jsc, jec
                do i = isc, iec
                   if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
                <celements>
                   ergom%sed_<total>(i,j,k) = &
                   <containingTracers vertLoc=WAT; isInPorewater=1>
                      max(0.0,ergom%sed_<ct>(i,j,k))*<ctAmount> + &
                   </containingTracers>
                      0.0
                </celements>
                enddo
             enddo
          enddo
       enddo

       !-----------------------------------
       ! STEP 9a: vertical diffusion of solid phase tracers
       !-----------------------------------
       ! first, calculate total diffusivity (molecular + bioturbation)
       do k=2,NUM_SEDIMENT_LAYERS
          do j=jsc,jec
             do i=isc,iec
                ! first, calculate porosity
                sed_porosity = 1.0-ergom%sed_inert_ratio(i,j,k)
                ! second, save cell heights in new array
                ergom%temp3darray(i,j,k) = ergom%sed_cellheights(i,j,k)*max(ergom%sed_inert_ratio(i,j,k),epsilon)
                if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
                ! third, convert all tracer concentrations from mol/m2 to mol/kg in solid phase
<tracers vertLoc=SED>
                ergom%sed_<name>(i,j,k) = ergom%sed_<name>(i,j,k)/(ergom%sed_cellheights(i,j,k)*1035.0)/max(ergom%sed_inert_ratio(i,j,k),epsilon)
</tracers>
                ! fourth, store diffusivity
                if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
<tracers vertLoc=SED>
                ergom%moldiff_<name>(i,j,k) = ergom%sed_diffusivity_solids(i,j,k)*power(ergom%sed_inert_ratio(i,j,k),2.0)
</tracers>
             enddo
          enddo
       enddo
       do j=jsc,jec
          do i=isc,iec
             if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
             ergom%temp3darray(i,j,1) = ergom%sed_2d_params(i,j,5) ! fluffy layer thickness
<tracers vertLoc=SED>
             ergom%sed_<name>(i,j,1) = ergom%sed_<name>(i,j,1)/(ergom%sed_2d_params(i,j,5)*1035.0)
</tracers>
          enddo
       enddo

       !-----------------------------------
       ! STEP 9a.1: calculate total element concentrations in the solid phase
       !-----------------------------------
       ! not necessary, as it was done in the movement steps before
       !-----------------------------------
       ! STEP 9.2: do several vertical diffusion steps
       !-----------------------------------
       do m = 1, NUM_VMOVE_STEPS
          !-----------------------------------
          ! STEP 9.2.1: store age concentrations in an array which is not modified by diffusion
          !-----------------------------------
          <celements isAging/=0>
            ergom%old_<aged>(:,:,:,1) = ergom%sed_<aged>(:,:,:)
          </celements>
          !-----------------------------------
          ! STEP 9.2.2: diffuse the age concentrations of the marked elements
          !-----------------------------------
          <tracers childOf/=none; vertLoc=SED; hasCeAged>
             call generic_ERGOM_seddiff(ergom%moldiff_<name>, &
                                      ergom%sed_<ceAgedName>(:,:,:), ergom%old_<ceAgedName>(:,:,:,1),&
                                      ergom%sed_<name>(:,:,:)*<ceAmount>, ergom%sed_<ceTotalName>(:,:,:), &
                                      ergom%temp3darray(:,:,:), ergom%temp3darray(:,:,1)*0.5, dt/NUM_VMOVE_STEPS, &
                                      isd, ied, jsd, jed, NUM_SEDIMENT_LAYERS)
           </tracers>

          !-----------------------------------
          ! STEP 9.2.3: move the tracers (including tagged tracers) themselves
          !-----------------------------------

          <tracers vertLoc=SED>
             call generic_ERGOM_seddiff(ergom%moldiff_<name>, &
                                      ergom%sed_<name>(:,:,:), ergom%sed_<name>(:,:,:), &
                                      ergom%sed_<name>(:,:,:), ergom%sed_<name>(:,:,:), &
                                      ergom%temp3darray(:,:,:), ergom%temp3darray(:,:,1)*0.5, dt/NUM_VMOVE_STEPS, &
                                      isd, ied, jsd, jed, NUM_SEDIMENT_LAYERS)
          </tracers>

          !-----------------------------------
          ! STEP 9.2.4: calculate the new total marked element concentrations
          !-----------------------------------
          do k = 2,NUM_SEDIMENT_LAYERS
             do j = jsc, jec
                do i = isc, iec
                   if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
                <celements>
                   ergom%sed_<total>(i,j,k) = &
                   <containingTracers vertLoc=WAT; isInPorewater=1>
                      max(0.0,ergom%sed_<ct>(i,j,k))*<ctAmount> + &
                   </containingTracers>
                      0.0
                </celements>
                enddo
             enddo
          enddo
       enddo

       ! finally, sediment concentrations have to be converted backwards from mol/kg to mol/m2
       do k=2,NUM_SEDIMENT_LAYERS
          do j=jsc,jec
             do i=isc,iec
                if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
<tracers vertLoc=SED>
                ergom%sed_<name>(i,j,k) = ergom%sed_<name>(i,j,k)*(ergom%sed_cellheights(i,j,k)*1035.0)*max(ergom%sed_inert_ratio(i,j,k),epsilon)
</tracers>
             enddo
          enddo
       enddo
       do j=jsc,jec
          do i=isc,iec
             if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
<tracers vertLoc=SED>
             ergom%sed_<name>(i,j,1) = ergom%sed_<name>(i,j,1)*(ergom%sed_2d_params(i,j,5)*1035.0)
</tracers>
          enddo
       enddo


       !-----------------------------------
       ! STEP 10: calculate the new total marked element concentrations at the bottom
       !-----------------------------------
        do j = jsc, jec
           do i = isc, iec
           <celements>
              ergom%<totalBottom>(i,j) = &
              <containingTracers vertLoc=SED>
                 max(0.0,ergom%<ct>(i,j))*<ctAmount> + &
              </containingTracers>
                 0.0
           </celements>
           enddo
        enddo

       call mpp_clock_end(id_sedvertmig)
       call mpp_clock_begin(id_sedoutput)
!
       !-----------------------------------
       ! STEP 11: save the new values to their designated places
       !-----------------------------------
<tracers vertLoc=WAT; isInPorewater=1>
       call g_tracer_set_values(tracer_list,'sed_<name>'  ,'field',ergom%sed_<name>,isd,jsd,ntau=1)
</tracers>
<tracers vertLoc=SED>
       call g_tracer_set_values(tracer_list,'<name>'  ,'field',ergom%sed_<name>,isd,jsd,ntau=1)
</tracers>

       ! store bottom water values for use in water column model
<tracers vertLoc=WAT; isInPorewater=1>
       call g_tracer_get_pointer(tracer_list,'<name>'  ,'field',ergom%<name>)
       do j = jsc, jec
          do i = isc, iec
             if ((1 .gt. grid_kmt(i,j)) .or. (0.5e-9 .gt. ergom%sed_cellheights(i,j,2))) cycle
             ergom%<name>(i,j,grid_kmt(i,j),tau) = ergom%sed_<name>(i,j,1)
          enddo
       enddo
</tracers>

!   Diagnostics

<auxiliaries>
       if (ergom%id_sed_<name> .gt. 0) &
          used = send_data(ergom%id_sed_<name>, ergom%sed_<name>,                 &
          model_time, &
          is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
</auxiliaries>

<processes vertLoc=WAT; isInPorewater=1>
       if (ergom%id_sed_<name> .gt. 0) &
          used = send_data(ergom%id_sed_<name>, ergom%sed_<name>,                 &
          model_time, &
          is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
</processes>

<processes vertLoc=SED; isInPorewater=1>
       if (ergom%id_sed_<name> .gt. 0) &
          used = send_data(ergom%id_sed_<name>, ergom%sed_<name>,                 &
          model_time, &
          is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
</processes>

       call mpp_clock_end(id_sedoutput)

    endif

    call mpp_clock_end(id_sedsource)

    return
  end subroutine cgt_bio_timestep_sed

  ! <SUBROUTINE NAME="generic_ERGOM_update_from_source">
  !  <OVERVIEW>
  !   Update tracer concentration fields due to the source/sink contributions.
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Ecosystem model code automatically generated by Code Generation Tool - see www.ergom.net
  !  </DESCRIPTION>
  ! </SUBROUTINE>
  subroutine generic_ERGOM_update_from_source(tracer_list,Temp,Salt, rho_dzt,dzt,xt,yt,hblt_depth,&
       ilb,jlb,tau,dt,grid_dat,model_time,nbands,max_wavelength_band, sw_pen_band,opacity_band, &
       current_wave_stress,diff_cbt)

    type(g_tracer_type),            pointer    :: tracer_list
    real, dimension(ilb:,jlb:,:),   intent(in) :: Temp,Salt,rho_dzt,dzt
    real, dimension(ilb:,jlb:),     intent(in) :: xt,yt
    real, dimension(ilb:,jlb:),     intent(in) :: hblt_depth
    integer,                        intent(in) :: ilb,jlb,tau
    real,                           intent(in) :: dt
    real, dimension(ilb:,jlb:),     intent(in) :: grid_dat
    type(time_type),                intent(in) :: model_time

    integer,                        intent(in) :: nbands
    real, dimension(:),             intent(in) :: max_wavelength_band
    real, dimension(:,ilb:,jlb:),   intent(in) :: sw_pen_band
    real, dimension(:,ilb:,jlb:,:), intent(in) :: opacity_band
    real, dimension(ilb:,jlb:),     intent(in) :: current_wave_stress
    real, dimension(ilb:,jlb:,:,:),optional,intent(in) :: diff_cbt

    integer :: isc,iec, jsc,jec,isd,ied,jsd,jed,nk,ntau, i, j, k , kblt, n, m, nband
    real, dimension(:,:,:) ,pointer :: grid_tmask
    integer, dimension(:,:),pointer :: mask_coast,grid_kmt

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,&
                             grid_tmask=grid_tmask,grid_mask_coast=mask_coast,grid_kmt=grid_kmt)

    if (MOD(INT(time_type_to_real(model_time)/dt+0.5),splitting_factor_ergom) .eq. 0) then
       if (RUNGE_KUTTA_DEPTH .eq. 1) then ! positive Euler forward or positive Euler forward Patankar

          call cgt_bio_timestep(tracer_list,Temp,Salt, rho_dzt,dzt,xt,yt,hblt_depth,&
             ilb,jlb,tau,dt*splitting_factor_ergom,grid_dat,model_time,nbands,max_wavelength_band, sw_pen_band,opacity_band, &
             current_wave_stress,.false.,.false.,diff_cbt)
          if (NUM_SEDIMENT_LAYERS .gt. 1) then
             call cgt_bio_timestep_sed(tracer_list,Temp,Salt, rho_dzt,dzt,xt,yt,hblt_depth,&
                ilb,jlb,tau,dt*splitting_factor_ergom,grid_dat,model_time,nbands,max_wavelength_band, sw_pen_band,opacity_band, &
                current_wave_stress,.false.,.false.)
          endif

       else if (RUNGE_KUTTA_DEPTH .eq. 2) then ! positive Runge-Kutta or positive Runge-Kutta Patankar
          !Step 1: save initial tracer values
          call user_get_2d_tracer_values(tracer_list,isd,ied,jsd,jed,nk)
          !Get the pointers to the prognostic tracers after applying the physics.
<tracers vertLoc=WAT>
          call g_tracer_get_pointer(tracer_list,'<name>'  ,'field',ergom%<name>    )
</tracers>
<tracers vertLoc=WAT>
          ergom%tracer_array_intermediate_<name>=ergom%<name>(:,:,:,tau)
</tracers>
<tracers vertLoc/=WAT>
          ergom%tracer_array_intermediate_<name>=ergom%<name>(:,:)
</tracers>

          !Step 2: initialize saved process rates
<processes>
          ergom%saved_rate_<name>=0.0;
</processes>

          !Step 3: Do intermediate timestep
          call cgt_bio_timestep(tracer_list,Temp,Salt, rho_dzt,dzt,xt,yt,hblt_depth,&
             ilb,jlb,tau,dt*splitting_factor_ergom,grid_dat,model_time,nbands,max_wavelength_band, sw_pen_band,opacity_band, &
             current_wave_stress,.false.,.true.)

          !Step 4: Calculate the patankar modification factors
          do k=1,nk
             do j=jsc,jec
                do i=isc,iec
<tracers vertLoc=WAT>
                   ergom%patankar_modification_<name>(i,j,k) = ergom%tracer_array_intermediate_<name>(i,j,k)/max(ergom%<name>(i,j,k,tau),1.0e-30)
</tracers>
                enddo
             enddo
          enddo
          do j=jsc,jec
             do i=isc,iec
<tracers vertLoc/=WAT>
                ergom%patankar_modification_<name>(i,j) = ergom%tracer_array_intermediate_<name>(i,j)/max(ergom%<name>(i,j),1.0e-30)
</tracers>
             enddo
          enddo

          !Step 5: Do another intermediate timestep to get the new process rates
          call cgt_bio_timestep(tracer_list,Temp,Salt, rho_dzt,dzt,xt,yt,hblt_depth,&
             ilb,jlb,tau,dt*splitting_factor_ergom,grid_dat,model_time,nbands,max_wavelength_band, sw_pen_band,opacity_band, &
             current_wave_stress,.false.,.true.)

          !Step 6: calculate the average of the process rates
<processes>
          ergom%saved_rate_<name>=ergom%saved_rate_<name>/2.0
</processes>

          !Step 7: Reload original tracer values
<tracers vertLoc=WAT>
          ergom%<name>(:,:,:,tau)=ergom%tracer_array_intermediate_<name>
</tracers>
<tracers vertLoc/=WAT>
          ergom%<name>(:,:)=ergom%tracer_array_intermediate_<name>
</tracers>
          call user_set_2d_tracer_values(tracer_list, model_time)  ! save all 2d tracers

          !Step 8: Do the actual time step
          call cgt_bio_timestep(tracer_list,Temp,Salt, rho_dzt,dzt,xt,yt,hblt_depth,&
             ilb,jlb,tau,dt*splitting_factor_ergom,grid_dat,model_time,nbands,max_wavelength_band, sw_pen_band,opacity_band, &
             current_wave_stress,.true.,.false.)

       endif
    endif

    return
  end subroutine generic_ERGOM_update_from_source


  ! <SUBROUTINE NAME="generic_ERGOM_set_boundary_values">
  !  <OVERVIEW>
  !   Calculate and set coupler values at the surface / bottom
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Subroutine is empty as these values are calculated in the main routine
  !  </DESCRIPTION>
  !  <TEMPLATE>
  !   call generic_ERGOM_set_boundary_values(tracer_list,SST,SSS,rho,ilb,jlb,tau)
  !  </TEMPLATE>
  !  <IN NAME="tracer_list" TYPE="type(g_tracer_type), pointer">
  !   Pointer to the head of generic tracer list.
  !  </IN>
  !  <IN NAME="ilb,jlb" TYPE="integer">
  !   Lower bounds of x and y extents of input arrays on data domain
  !  </IN>
  !  <IN NAME="SST" TYPE="real, dimension(ilb:,jlb:)">
  !   Sea Surface Temperature
  !  </IN>
  !  <IN NAME="SSS" TYPE="real, dimension(ilb:,jlb:)">
  !   Sea Surface Salinity
  !  </IN>
  !  <IN NAME="rho" TYPE="real, dimension(ilb:,jlb:,:,:)">
  !   Ocean density
  !  </IN>
  !  <IN NAME="tau" TYPE="integer">
  !   Time step index of %field
  !  </IN>
  ! </SUBROUTINE>

  !User must provide the calculations for these boundary values.
  subroutine generic_ERGOM_set_boundary_values(tracer_list,SST,SSS,rho,ilb,jlb,taum1)
    type(g_tracer_type),          pointer    :: tracer_list
    real, dimension(ilb:,jlb:),     intent(in) :: SST, SSS
    real, dimension(ilb:,jlb:,:,:), intent(in) :: rho
    integer,                        intent(in) :: ilb,jlb,taum1
  end subroutine generic_ERGOM_set_boundary_values

  ! <SUBROUTINE NAME="generic_ERGOM_vmove">
  !  <OVERVIEW>
  !   Performs vertical movement (up or down)
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Updates particulate tracer concentrations
  !  </DESCRIPTION>
  !  <TEMPLATE>
  !   call generic_ERGOM_vmove
  !  </TEMPLATE>
  ! </SUBROUTINE>

  subroutine generic_ERGOM_vmove(move, field, flux_field, numerator, denominator, dzt, dt, ilb, iub, jlb, jub)
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: move        ! vertical speed [m/s]
    real, dimension(ilb:,jlb:,:) ,intent(inout) :: field       ! tracer array to be moved
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: flux_field  ! (old) tracer array from which the fluxes result
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: numerator   ! vertical transports are multiplied by this scaling factor
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: denominator ! vertical transports are divided by this scaling factor
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: dzt         ! cell heights [m]
    real,                   intent(in)          :: dt          ! timestep [s]
    integer,                intent(in)          :: ilb, iub, jlb, jub

    real, dimension(:,:,:)  , pointer  :: tmask
    integer, dimension(:,:) , pointer  :: kmt
    real, dimension(ilb:iub,jlb:jub)   :: ft1, ft2
    real                               :: velocity, wpos, wneg
    integer                            :: k, i, j, kp1
    integer                            :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau
    character(len=fm_string_len), parameter :: sub_name = 'generic_ERGOM_vmove'

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,&
         grid_tmask=tmask, grid_kmt=kmt)
    ft1  = 0.0
    do k = 1, nk-1
      kp1 = k+1
      do j = jsc, jec ; do i = isc, iec
        velocity = 0.5*move(i,j,k)
        wpos     = velocity + abs(velocity)                    ! velocity if upward, else 0.0   [m/s]
        wneg     = velocity - abs(velocity)                    ! velocity if downward, else 0.0 [m/s]
        if (numerator(i,j,k) .eq. denominator(i,j,k)) then
           ft2(i,j) = (wneg*max(flux_field(i,j,k),0.0)+wpos*max(flux_field(i,j,kp1),0.0) ) &
                   *tmask(i,j,k)*tmask(i,j,kp1)                ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        else
           ft2(i,j) = (wneg*max(flux_field(i,j,k),0.0)*max(numerator(i,j,k),0.0)/max(denominator(i,j,k),1e-20) &
                   +wpos*max(flux_field(i,j,kp1),0.0)*max(numerator(i,j,kp1),0.0)/max(denominator(i,j,kp1),1e-20) ) &
                   *tmask(i,j,k)*tmask(i,j,kp1)                ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        endif
        field(i,j,k) = field(i,j,k) - dt*tmask(i,j,k)*(ft1(i,j)-ft2(i,j))/dzt(i,j,k)  ! change in the t-cell due to transports through lower and upper boundary [mol/kg/s]
        ft1(i,j) = ft2(i,j)
      enddo; enddo
    enddo
    k = nk
    do j = jsc, jec ; do i = isc, iec
      field(i,j,k) = field(i,j,k) - dt*tmask(i,j,k)*ft1(i,j)/dzt(i,j,k)
    enddo; enddo
  end subroutine generic_ERGOM_vmove

    ! <SUBROUTINE NAME="generic_ERGOM_sedmove">
  !  <OVERVIEW>
  !   Performs vertical movement (up or down)
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Updates particulate tracer concentrations
  !  </DESCRIPTION>
  !  <TEMPLATE>
  !   call generic_ERGOM_vmove
  !  </TEMPLATE>
  ! </SUBROUTINE>

  subroutine generic_ERGOM_sedmove(move, field, flux_field, numerator, denominator, dzt, dztop, dt, ilb, iub, jlb, jub, kmax)
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: move        ! vertical speed [m/s]
    real, dimension(ilb:,jlb:,:) ,intent(inout) :: field       ! tracer array to be moved
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: flux_field  ! (old) tracer array from which the fluxes result
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: numerator   ! vertical transports are multiplied by this scaling factor
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: denominator ! vertical transports are divided by this scaling factor
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: dzt         ! cell heights [m]
    real, dimension(ilb:,jlb:)   ,intent(in)    :: dztop       ! cell height of fluffy layer [m]
    real,                   intent(in)          :: dt          ! timestep [s]
    integer,                intent(in)          :: ilb, iub, jlb, jub, kmax

    real, dimension(:,:,:)  , pointer  :: tmask
    integer, dimension(:,:) , pointer  :: kmt
    real, dimension(ilb:iub,jlb:jub)   :: ft1, ft2
    real                               :: velocity, wpos, wneg
    integer                            :: k, i, j, kp1
    integer                            :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau
    character(len=fm_string_len), parameter :: sub_name = 'generic_ERGOM_vmove'

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,&
         grid_tmask=tmask, grid_kmt=kmt)
    ft1  = 0.0
    do k = 1, 1
      kp1 = k+1
      do j = jsc, jec ; do i = isc, iec
        if (0.5e-9 .gt. dzt(i,j,2)) cycle
        velocity = 0.5*move(i,j,k)
        wpos     = velocity + abs(velocity)                    ! velocity if upward, else 0.0   [m/s]
        wneg     = velocity - abs(velocity)                    ! velocity if downward, else 0.0 [m/s]
        if (numerator(i,j,k) .eq. denominator(i,j,k)) then
           ft2(i,j) = (wneg*max(flux_field(i,j,k),0.0)+wpos*max(flux_field(i,j,kp1),0.0) ) ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        else
           ft2(i,j) = (wneg*max(flux_field(i,j,k),0.0)*max(numerator(i,j,k),0.0)/max(denominator(i,j,k),1e-20) &
                   +wpos*max(flux_field(i,j,kp1),0.0)*max(numerator(i,j,kp1),0.0)/max(denominator(i,j,kp1),1e-20) ) ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        endif
        field(i,j,k) = field(i,j,k) - dt*(ft1(i,j)-ft2(i,j))/dztop(i,j)  ! change in the t-cell due to transports through lower and upper boundary [mol/kg/s]
        ft1(i,j) = ft2(i,j)
      enddo; enddo
    enddo
    do k = 2, kmax-1
      kp1 = k+1
      do j = jsc, jec ; do i = isc, iec
        if (0.5e-9 .gt. dzt(i,j,2)) cycle
        velocity = 0.5*move(i,j,k)
        wpos     = velocity + abs(velocity)                    ! velocity if upward, else 0.0   [m/s]
        wneg     = velocity - abs(velocity)                    ! velocity if downward, else 0.0 [m/s]
        if (numerator(i,j,k) .eq. denominator(i,j,k)) then
           ft2(i,j) = (wneg*max(flux_field(i,j,k),0.0)+wpos*max(flux_field(i,j,kp1),0.0) ) ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        else
           ft2(i,j) = (wneg*max(flux_field(i,j,k),0.0)*max(numerator(i,j,k),0.0)/max(denominator(i,j,k),1e-20) &
                   +wpos*max(flux_field(i,j,kp1),0.0)*max(numerator(i,j,kp1),0.0)/max(denominator(i,j,kp1),1e-20) ) ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        endif
        field(i,j,k) = field(i,j,k) - dt*(ft1(i,j)-ft2(i,j))/dzt(i,j,k)  ! change in the t-cell due to transports through lower and upper boundary [mol/kg/s]
        ft1(i,j) = ft2(i,j)
      enddo; enddo
    enddo
    do k = kmax, kmax
      do j = jsc, jec ; do i = isc, iec
        if (0.5e-9 .gt. dzt(i,j,2)) cycle
        velocity = 0.5*move(i,j,k)
        wpos     = velocity + abs(velocity)                    ! velocity if upward, else 0.0   [m/s]
        wneg     = velocity - abs(velocity)                    ! velocity if downward, else 0.0 [m/s]
        if (numerator(i,j,k) .eq. denominator(i,j,k)) then
           ft2(i,j) = (wneg*max(flux_field(i,j,k),0.0))    ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        else
           ft2(i,j) = (wneg*max(flux_field(i,j,k),0.0)*max(numerator(i,j,k),0.0)/max(denominator(i,j,k),1e-20)) ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        endif
        field(i,j,k) = field(i,j,k) - dt*(ft1(i,j)-ft2(i,j))/dzt(i,j,k)  ! change in the t-cell due to transports through lower and upper boundary [mol/kg/s]
        ft1(i,j) = ft2(i,j)
      enddo; enddo
    enddo
  end subroutine generic_ERGOM_sedmove

  ! <SUBROUTINE NAME="generic_ERGOM_vdiff">
  !  <OVERVIEW>
  !   Performs additional vertical diffusion for some tracers
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Updates particulate tracer concentrations
  !  </DESCRIPTION>
  !  <TEMPLATE>
  !   call generic_ERGOM_vmove
  !  </TEMPLATE>
  ! </SUBROUTINE>

  subroutine generic_ERGOM_vdiff(diff, field, flux_field, numerator, denominator, dzt, dt, ilb, iub, jlb, jub)
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: diff        ! vertical diffusivity [m2/s]
    real, dimension(ilb:,jlb:,:) ,intent(inout) :: field       ! tracer array to be diffused
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: flux_field  ! (old) tracer array from which the fluxes are calculated
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: numerator   ! vertical transports are multiplied by this scaling factor
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: denominator ! vertical transports are divided by this scaling factor
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: dzt         ! cell heights [m]
    real,                   intent(in)          :: dt          ! timestep [s]
    integer,                intent(in)          :: ilb, iub, jlb, jub

    real, dimension(:,:,:)  , pointer  :: tmask
    integer, dimension(:,:) , pointer  :: kmt
    real, dimension(ilb:iub,jlb:jub)   :: ft1, ft2
    real                               :: diffusivity
    integer                            :: k, i, j, kp1
    integer                            :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau
    character(len=fm_string_len), parameter :: sub_name = 'generic_ERGOM_vdiff'

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,&
         grid_tmask=tmask, grid_kmt=kmt)
    ft1  = 0.0
    do k = 1, nk-1
      kp1 = k+1
      do j = jsc, jec ; do i = isc, iec
        diffusivity = 0.5*diff(i,j,k)+0.5*diff(i,j,kp1)
        if (numerator(i,j,k) .eq. denominator(i,j,k)) then
           ft2(i,j) = (max(flux_field(i,j,kp1),0.0)-max(flux_field(i,j,k),0.0))*diffusivity/(0.5*(dzt(i,j,k)+dzt(i,j,kp1))) &
                   *tmask(i,j,k)*tmask(i,j,kp1)                ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        else
           ft2(i,j) = (max(flux_field(i,j,kp1),0.0)*max(numerator(i,j,kp1),0.0)/max(denominator(i,j,kp1),1e-20) &
                      -max(flux_field(i,j,k),0.0)*max(numerator(i,j,k),0.0)/max(denominator(i,j,k),1e-20)) &
                      *diffusivity/(0.5*(dzt(i,j,k)+dzt(i,j,kp1))) &
                      *tmask(i,j,k)*tmask(i,j,kp1)                ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        endif
        field(i,j,k) = field(i,j,k) - dt*tmask(i,j,k)*(ft1(i,j)-ft2(i,j))/dzt(i,j,k)  ! change in the t-cell due to transports through lower and upper boundary [mol/kg/s]
        ft1(i,j) = ft2(i,j)
      enddo; enddo
    enddo
    k = nk
    do j = jsc, jec ; do i = isc, iec
      field(i,j,k) = field(i,j,k) - dt*tmask(i,j,k)*ft1(i,j)/dzt(i,j,k)
    enddo; enddo
  end subroutine generic_ERGOM_vdiff

  ! <SUBROUTINE NAME="generic_ERGOM_seddiff">
  !  <OVERVIEW>
  !   Performs vertical diffusion of tracers in the sediment
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Updates pore water concentrations
  !   A separate routine is needed for two reasons:
  !   Firstly, there may be less sediment layers than water layers.
  !   Secondly, the gradient between bottom water and uppermost porewater cell is higher than the
  !   distance between its cell centers suggests
  !  </DESCRIPTION>
  !  <TEMPLATE>
  !   call generic_ERGOM_vdiff
  !  </TEMPLATE>
  ! </SUBROUTINE>

  subroutine generic_ERGOM_seddiff(diff, field, flux_field, numerator, denominator, dzt, dztop, dt, ilb, iub, jlb, jub, kmax)
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: diff        ! vertical diffusivity [m2/s]
    real, dimension(ilb:,jlb:,:) ,intent(inout) :: field       ! tracer array to be diffused
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: flux_field  ! (old) tracer array from which the fluxes are calculated
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: numerator   ! vertical transports are multiplied by this scaling factor
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: denominator ! vertical transports are divided by this scaling factor
    real, dimension(ilb:,jlb:,:) ,intent(in)    :: dzt         ! cell heights [m]
    real, dimension(ilb:,jlb:)   ,intent(in)    :: dztop       ! distance across which the concentration gradient extends in the top layer [m]
    real,                   intent(in)          :: dt          ! timestep [s]
    integer,                intent(in)          :: ilb, iub, jlb, jub

    real, dimension(:,:,:)  , pointer  :: tmask
    integer, dimension(:,:) , pointer  :: kmt
    real, dimension(ilb:iub,jlb:jub)   :: ft1, ft2
    real                               :: diffusivity
    integer                            :: k, i, j, kp1
    integer                            :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,kmax
    character(len=fm_string_len), parameter :: sub_name = 'generic_ERGOM_seddiff'

    call g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,&
         grid_tmask=tmask, grid_kmt=kmt)
    ft1  = 0.0
    do k = 1, 1
      kp1 = k+1
      do j = jsc, jec ; do i = isc, iec
        if (0.5e-9 .gt. dzt(i,j,2)) cycle
        diffusivity = diff(i,j,kp1)
        if (numerator(i,j,k) .eq. denominator(i,j,k)) then
           ft2(i,j) = (max(flux_field(i,j,kp1),0.0)-max(flux_field(i,j,k),0.0))*diffusivity/(0.5*dzt(i,j,kp1)+dztop(i,j))  ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        else
           ft2(i,j) = (max(flux_field(i,j,kp1),0.0)*max(numerator(i,j,kp1),0.0)/max(denominator(i,j,kp1),1e-20) &
                      -max(flux_field(i,j,k),0.0)*max(numerator(i,j,k),0.0)/max(denominator(i,j,k),1e-20)) &
                      *diffusivity/(0.5*dzt(i,j,kp1)+dztop(i,j))                                                           ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        endif
        field(i,j,k) = field(i,j,k) - dt*(ft1(i,j)-ft2(i,j))/dzt(i,j,k)  ! change in the t-cell due to transports through lower and upper boundary [mol/kg/s]
        ft1(i,j) = ft2(i,j)
      enddo; enddo
    enddo
    do k = 2, kmax-1
      kp1 = k+1
      do j = jsc, jec ; do i = isc, iec
        if (0.5e-9 .gt. dzt(i,j,2)) cycle
        diffusivity = 0.5*diff(i,j,k)+0.5*diff(i,j,kp1)
        if (numerator(i,j,k) .eq. denominator(i,j,k)) then
           ft2(i,j) = (max(flux_field(i,j,kp1),0.0)-max(flux_field(i,j,k),0.0))*diffusivity/(0.5*(dzt(i,j,k)+dzt(i,j,kp1))) ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        else
           ft2(i,j) = (max(flux_field(i,j,kp1),0.0)*max(numerator(i,j,kp1),0.0)/max(denominator(i,j,kp1),1e-20) &
                      -max(flux_field(i,j,k),0.0)*max(numerator(i,j,k),0.0)/max(denominator(i,j,k),1e-20)) &
                      *diffusivity/(0.5*(dzt(i,j,k)+dzt(i,j,kp1)))  ! upward transport through lower boundary of the t cell [mol*m/kg/s]
        endif
        field(i,j,k) = field(i,j,k) - dt*(ft1(i,j)-ft2(i,j))/dzt(i,j,k)  ! change in the t-cell due to transports through lower and upper boundary [mol/kg/s]
        ft1(i,j) = ft2(i,j)
      enddo; enddo
    enddo
    k = kmax
    do j = jsc, jec ; do i = isc, iec
      if (0.5e-9 .gt. dzt(i,j,2)) cycle
      field(i,j,k) = field(i,j,k) - dt*ft1(i,j)/dzt(i,j,k)
    enddo; enddo
  end subroutine generic_ERGOM_seddiff

  ! <SUBROUTINE NAME="generic_ERGOM_end">
  !  <OVERVIEW>
  !   End the module.
  !  </OVERVIEW>
  !  <DESCRIPTION>
  !   Deallocate all work arrays
  !  </DESCRIPTION>
  !  <TEMPLATE>
  !   call generic_ERGOM_end
  !  </TEMPLATE>
  ! </SUBROUTINE>

  subroutine generic_ERGOM_end
    character(len=fm_string_len), parameter :: sub_name = 'generic_ERGOM_end'
    call user_deallocate_arrays
  end subroutine generic_ERGOM_end

  subroutine user_deallocate_arrays
  integer :: n

  <constants dependsOn=xyz>
    deallocate(ergom%<name>)
  </constants>
  <constants dependsOn=xyzt>
    deallocate(ergom%<name>)
  </constants>
  <auxiliaries vertLoc=WAT; isUsedElsewhere=1>
    deallocate(ergom%<name>)
  </auxiliaries>
  <auxiliaries vertLoc=WAT; isUsedElsewhere=0>
    if (ergom%id_<name> .gt. 0) deallocate(ergom%<name>)
  </auxiliaries>
  <auxiliaries vertLoc/=WAT>
    if (ergom%id_<name> .gt. 0) deallocate(ergom%<name>)
  </auxiliaries>
  <auxiliaries isZIntegral=1>
    deallocate(ergom%zintegralarray_<name>)
  </auxiliaries>
  <tracers vertLoc=FIS>
    deallocate(ergom%cumulated_change_of_<name>)
  </tracers>
  <celements isAging/=0>
    deallocate(ergom%old_<aged>)
  </celements>

  <processes>
    if (ergom%id_<name> .gt. 0) deallocate(ergom%<name>)
  </processes>

    if (implicit_surfacefluxes) then
  <tracers solubility/=0; childOf=none>
       deallocate(ergom%<gasName>_alpha)
       deallocate(ergom%<gasName>_sc_no)
  </tracers>
  <tracers solubility/=0>
       deallocate(ergom%<gasName>_<trimName>_csurf)
  </tracers>
    endif

  <tracers vertLoc/=WAT>
    deallocate(ergom%<name>)
</tracers>

    if (implicit_bottomfluxes) then
<tracers vertLoc=WAT>
      deallocate(ergom%btf_<name>)
</tracers>
    endif

    deallocate(tracers_2d)
    if (NUM_SEDIMENT_LAYERS .gt. 1) then
       deallocate(ergom%temp3darray)
       deallocate(ergom%vertspeedarray)
<tracers vertLoc=WAT; isInPorewater=1>
       deallocate(ergom%sed_<name>)
</tracers>
<tracers vertLoc=WAT; isInPorewater=1; molDiff/=0>
       deallocate(ergom%moldiff_<name>)
</tracers>
<tracers vertLoc=SED>
       deallocate(ergom%sed_<name>)
</tracers>
<celements>
       deallocate(ergom%sed_<total>)
</celements>
       deallocate(ergom%sed_cellheights          )
       deallocate(ergom%sed_diffusivity_porewater)
       deallocate(ergom%sed_diffusivity_solids   )
       deallocate(ergom%sed_2d_params            )
       deallocate(ergom%sed_inert_ratio          )
    endif

    if (RUNGE_KUTTA_DEPTH .eq. 2) then
<tracers>
       deallocate(ergom%tracer_array_intermediate_<name>)
       deallocate(ergom%patankar_modification_<name>)
</tracers>
       if (NUM_SEDIMENT_LAYERS .gt. 1) then
<tracers vertLoc=WAT; isInPorewater=1>
          deallocate(ergom%tracer_array_intermediate_sed_<name>)
          deallocate(ergom%patankar_modification_sed_<name>)
</tracers>
<tracers vertLoc=SED>
          deallocate(ergom%tracer_array_intermediate_sed_<name>)
          deallocate(ergom%patankar_modification_sed_<name>)
</tracers>
       endif
<processes>
       deallocate(ergom%saved_rate_<name>)
</processes>
    endif
  end subroutine user_deallocate_arrays

end module generic_ERGOM
