module biomod

  !============================================================================!
  ! A module containing core subroutines for the eco model. The model is based 
  ! on Neumann 2000, Neumann 2002 and DMI developments.
  !============================================================================!

  implicit none
  private
  
  public :: bio_dynamics

contains

  subroutine bio_dynamics(iw2,ke,msrf,mcol,ind,kh,h,u,v,wu,wv,t,s,humid,cloud, &
                          ts, benthos,xkoor,ykoor, dispvc,casus,npr, dts, idx, &
                          irc,rivmap,diag_light, diag_kpar                     &
                          , srfammflx                                          &
                          , srfnitflx                                          &
                          , srfphosflx                                          &
                          , srfamm_with_totatmos_Nflx                                          &
                          , srfnit_with_totatmos_Nflx                                          &
                          , srfamm_with_agri_Nflx                                          &
                          , srfnit_with_agri_Nflx                                          &
                          )

    !==========================================================================!
    !
    ! Biogeochemical dynamics
    !
    ! Biogeochemical-active tracers (physically passive) are solved with an 
    ! operator splitting scheme. tflow solves the tracer advection while the 
    ! change in concentration due to biogeochemical intercations and vertical 
    ! movement are calculated in the subroutine bio_dynamics. 
    !
    ! bio_dynamics:
    !
    ! The subroutine is called from the coupler at every #NTBIOMOD timestep 
    !
    !   i m                   zonal dimension (n-->s)
    !   i n                   merid dimension (w-->e)
    !
    ! input:
    !
    !   i iw2                 2d (surface) wet points
    !   i ke                  number of layers
    !   i msrf(0:m+1,0:n+1), mcol(0:iw2) padded land/sea mask
    !   i ind(2,iw2)          index i,j of each surface point
    !   i kh(0:iw2)           No. of active layers at each wet surface point
    !   r h(0:iw3)            layer thicknesses
    !   r wu(0:iw2)           meridional wind speed 
    !   r wv(0:iw2)           zonal wind speed 
    !   r ts(1:2+nrt,0:iw3)   temperature, salinity and ecotracers
    !   r benthos(0:iw2,1:ntben)      benthos nitrogen
    !   r swrad        shortwave radiation under the air-sea interface
    !   r dts                 time step  (seconds)
    !
    !   parameter:
    !
    !   i ntbio=9             number of ecological variables
    !
    !   returns:
    !
    !   r ts                  with updated ecological tracer concentrations
    !   r benthos             with uptaded concentrations
    !
    !   ecological variables, optional tracers in ():
    !
    !   0     dummy values
    !   1     ammonium
    !   2     nitrate
    !   3     phosphate
    !   4     diatoms
    !   5     flagellates
    !   6     cyano bacteria
    !   7     mesozooplankton (or bulk zooplankton)
    !   (8)   Optional microzooplankton
    !   8(9)  detritus
    !   9(10) oxygen
    !
    ! TODO:
    !   * Sediment erosion
    !   * Implicit scheme
    !   * Should nutrients be accumulated on the ice?
    !   * Should we really use the NP array in the surface loop if it is const.?
    !
    !==========================================================================!

    !--------------------------------------------------------------------------!
    ! Biomod params: constants, water type specific
    !--------------------------------------------------------------------------!
    use cmod_timer,     only : timer
    use dmi_omp,        only : domp_get_domain
    use dmi_mpi_global, only : iu06
    use bioparam,       only : &
                               zero, one, half, four, ten, tenth,             &
                               idx_amm            ,                          &
                               idx_nit            ,                          &
                               idx_phos           ,                          &
                               idx_sil            ,                          &
                               idx_dia            ,                          &
                               idx_flag           ,                          &
                               idx_cyano          ,                          &
                               idx_mez            ,                          &
                               idx_miz            ,                          &
                               idx_det            ,                          &
                               idx_dets           ,                          &
                               idx_ldon           ,                          &
                               idx_oxy            ,                          &
                               idx_amm_with_totatmos_N,                          &
                               idx_nit_with_totatmos_N,                          &
                               idx_dia_with_totatmos_N,                          &
                               idx_flag_with_totatmos_N,                          &
                               idx_cyano_with_totatmos_N,                          &
                               idx_mez_with_totatmos_N,                          &
                               idx_miz_with_totatmos_N,                          &
                               idx_det_with_totatmos_N,                          &
                               idx_ldon_with_totatmos_N,                          &
                               idx_amm_with_agri_N,                          &
                               idx_nit_with_agri_N,                          &
                               idx_dia_with_agri_N,                          &
                               idx_flag_with_agri_N,                          &
                               idx_cyano_with_agri_N,                          &
                               idx_mez_with_agri_N,                          &
                               idx_miz_with_agri_N,                          &
                               idx_det_with_agri_N,                          &
                               idx_ldon_with_agri_N,                          &
                               bidx_nitr           ,                         &
                               bidx_sili           ,                         &
                               bidx_nitr_with_totatmos_N,                         &
                               bidx_nitr_with_agri_N,                         &
                               alphab, alphaf, alphap, ntbio, deltao,         &
                               cyanosul, b0, anit, cyanotll, dn,ds, f0, flagtll,&
                               ksdenit, ksnit, nb,  nor, onnitr, ONamup, ntben, &
                               onniup, Onorm_inv, p0, q10_rec, rb0, rf0,rp0,  &
                               NNdenit, rnit,                                 &
                               wdiaz,                                        &
                               wcyanoz,                                        &
                               wdetz,                                        &
                               wdetsz,                                        &
                               Nnorm_inv, min_odial, min_oflagl, min_ocyanol, &
                               ppispmmod, rfs, flagtsc, ldon_frac,            &
                               ldon_tor, do_diag_secchi, par_calc_don,        &
                               miztk, mizgraz,mizgrazf,mezgraz, mezgrazf,     &
                               mezprefmiz,                                    &
                               mezcl, ldmez, lnmez,ldmiz, lnmiz,              &
                               sigma_b,                                       &
                               mizprefdia, mizprefflag, mizprefcyano,         &
                               mezprefdia, mezprefflag, mezprefcyano,         &
                               mez_tscale, meztk, warn_bad_solver
    use biogas,         only : oxygensat, pistonvel
    use biolight,       only : bio_light
    use biovmove,       only : vmove, sedimentation_resuspension, drying_case
    use biobenthic,     only : bio_benthic
    use biosolver,      only : solvebio
    use biosun,         only : shortwaverad

    !--------------------------------------------------------------------------!
    ! Cmod params: dimension, depth, lsm
    !--------------------------------------------------------------------------!
    integer(4), intent(in)    :: iw2, ke
    integer(4), intent(in)    :: msrf(0:,0:), mcol(0:)
    integer(4), intent(in)    :: ind(:,:)         ! (2,iw2)
    integer(4), intent(inout) :: idx(1:,1:)       ! (2,iw2) thread decompo
    integer(4), intent(in)    :: kh(0:)           ! (0:iw2)
    integer(4), intent(inout) :: irc              ! error status
    real(8),    intent(in)    :: xkoor(:),ykoor(:)! (n),(m)
    real(8),    intent(in)    :: h(0:)            ! (0:iw3)
    real(8),    intent(in)    :: u(0:), v(0:)     ! (0:iw3)
    real(8),    intent(in)    :: wu(0:), wv(0:)   ! (0:iw2)
    real(8),    intent(in)    :: dispvc(0:)       ! (0:iw3) eddyvdiff for salt
    real(8),    intent(in)    :: t(0:), s(0:)     ! (0:iw3) T&S
    real(8),    intent(in)    :: humid(0:), cloud(0:) ! (0:iw3) humidity&cloud.
    real(8),    intent(inout) :: ts(0:,1:)        ! (0:iw3,1:ntbio)  eco
    logical,    intent(in)    :: casus(0:)        ! (0:iw2) ice cover
    real(8),    intent(inout) :: benthos(0:,1:)   ! sediment
    real(8),    intent(in)    :: npr(0:)          ! N/P ratio
    integer(4), intent(in)    :: rivmap(0:)       !index list of deactivated riv
    real(8),    intent(inout) :: diag_light(0:)   !(0:iw3)
    real(8),    intent(out)   :: diag_kpar(0:)    !(0:iw3)
    real(8),    intent(in)    :: srfammflx(0:)     !(0:iw2)
    real(8),    intent(in)    :: srfnitflx(0:)     !(0:iw2)
    real(8),    intent(in)    :: srfphosflx(0:)     !(0:iw2)
    real(8),    intent(in)    :: srfamm_with_totatmos_Nflx(0:)     !(0:iw2)
    real(8),    intent(in)    :: srfnit_with_totatmos_Nflx(0:)     !(0:iw2)
    real(8),    intent(in)    :: srfamm_with_agri_Nflx(0:)     !(0:iw2)
    real(8),    intent(in)    :: srfnit_with_agri_Nflx(0:)     !(0:iw2)

    !--------------------------------------------------------------------------!
    ! Model time step (in seconds)
    !--------------------------------------------------------------------------!
    real(8), intent(in) :: dts

    !--------------------------------------------------------------------------!
    ! local arrays for the whole water column 
    !--------------------------------------------------------------------------!
    real(8) :: ecoa(1:ke,0:ntbio)   ! eco variables (inactive vars at 0)
    real(8) :: ben(0:ntben)         ! benthos variable
    real(8) :: deriv(0:ntbio)       ! Derivatves of eco variables
    real(8) :: benderiv(0:ntben)    ! Deriviate of benthos
    real(8) :: light(1:ke)          ! PAR in water column
    real(8) :: kpar(1:ke)           ! Kpar in water column
    real(8) :: light_o(1:ke)        ! PAR in water column as input for light3
    real(8) :: hi(1:ke)             ! grid depth
    real(8) :: t_k(1:ke), s_k(1:ke) ! salt and temp
    real(8) :: eddyvdiff(1:ke)      ! vertical diffusivity
   
    !--------------------------------------------------------------------------!
    ! local variables
    !--------------------------------------------------------------------------!
    real(8) :: dt             ! time step in days
    real(8) :: rfr, rfr2      ! N/P ratio, N/P ratio squared
    real(8) :: amm            ! ammonium
    real(8) :: nit            ! nitrate
    real(8) :: phos           ! phosphate
    real(8) :: sil            ! silicate
    real(8) :: dia            ! diatoms
    real(8) :: flag           ! flagellates
    real(8) :: cyano          ! cyano bacteria
    real(8) :: mez            ! meso zooplankton
    real(8) :: miz            ! micro zooplakton
    real(8) :: det            ! Detritus
    real(8) :: dets           ! S-Detritus
    real(8) :: ldon           ! ldon is the labile DON fraction
    real(8) :: oxy            ! Oxygen
    real(8) :: amm_with_totatmos_N! ammonium; containing 
    real(8) :: nit_with_totatmos_N! nitrate; containing 
    real(8) :: dia_with_totatmos_N! diatoms; containing 
    real(8) :: flag_with_totatmos_N! flagellates; containing 
    real(8) :: cyano_with_totatmos_N! cyano bacteria; containing 
    real(8) :: mez_with_totatmos_N! meso zooplankton; containing 
    real(8) :: miz_with_totatmos_N! micro zooplakton; containing 
    real(8) :: det_with_totatmos_N! Detritus; containing 
    real(8) :: ldon_with_totatmos_N! ldon is the labile DON fraction; containing 
    real(8) :: amm_with_agri_N! ammonium; containing 
    real(8) :: nit_with_agri_N! nitrate; containing 
    real(8) :: dia_with_agri_N! diatoms; containing 
    real(8) :: flag_with_agri_N! flagellates; containing 
    real(8) :: cyano_with_agri_N! cyano bacteria; containing 
    real(8) :: mez_with_agri_N! meso zooplankton; containing 
    real(8) :: miz_with_agri_N! micro zooplakton; containing 
    real(8) :: det_with_agri_N! Detritus; containing 
    real(8) :: ldon_with_agri_N! ldon is the labile DON fraction; containing 
    real(8) :: temp           ! temperature
    real(8) :: lightk         ! light
    real(8) :: sali           ! salinity
    real(8) :: osat           ! oxygen saturation 
    real(8) :: hsrf, hbot     ! surfacei/bottom layer height
    real(8) :: swrad          ! short wave radiation
    real(8) :: ksdenit2       ! half saturation oxygen dependence of nitrif.
    real(8) :: wui, wvi, wind ! wind components and windspeed
                              ! dia and flag nutrent half saturation squared
    real(8) :: alphap2, alphaf2 
                              ! nutrent half saturation times redfield squared
    real(8) :: albrfr2, alprfs2, alphab2
    real(8) :: ust,vst,velkb
    
    logical :: icefree        ! ice cover
   

    !--------------------------------------------------------------------------!
    ! Rates
    !--------------------------------------------------------------------------!
    
    REAL(8) :: nutlimc          ! nutrient limitation phosphate
    REAL(8) :: ppikb            ! ratio of optimal light for bluegreen
    REAL(8) :: rb               ! uptake rate of cyano
    REAL(8) :: nutlimf          ! nutrient limitation phosphate or nitrogen
    REAL(8) :: ppikf            ! ratio of optimal light for flagellats
    REAL(8) :: rf               ! uptake rate of flag
    REAL(8) :: nutlimd          ! nutrient limitation phosphate, nitrogen or silicate
    REAL(8) :: ppikp            ! ratio of optimal light for diatoms
    REAL(8) :: rp               ! uptake rate of dia
    REAL(8) :: nf               ! nitrification rate
    REAL(8) :: foodmez          ! More help variables
    REAL(8) :: temp2            ! help variable
    REAL(8) :: mezgt            ! mez grazing temperature dependece
    REAL(8) :: mezgscal         ! food dependent grazing; these have been divided by food; mezfgrazf2 low compared to paper
    REAL(8) :: meztotgraz       ! total mez grazin
    REAL(8) :: foodmiz          ! More help variables
    REAL(8) :: mizgt            ! miz grazing temperature dependece
    REAL(8) :: mizgscal         ! food dependent grazing; these have been divided by food
    REAL(8) :: miztotgraz       ! total miz grazing
    REAL(8) :: lpn              ! phytoplankton respiration loss rate (to inorg. compounds)
    REAL(8) :: denitnscal       ! using nitrate (and suphalte) to oxidize organic material
    REAL(8) :: ldn_N            ! respiration of nitrate under hypoxic conditions; else = 0.0_8
    REAL(8) :: ldn_O            ! less hydrogen sulphite formation under hypoxic conditions; else = ONamup
    REAL(8) :: lpd              ! phytoplankton respiration loss rate (to detritus)
    REAL(8) :: tldmez           ! meso zoo plankton loss rates (mortality)
    REAL(8) :: tldmiz           ! micro zoo plankton loss rates (mortality)
    REAL(8) :: fac              ! factor for detritus recycling rates (lds and ldn)
    REAL(8) :: ldn              ! detritus recycling rate (nitrogen, D => Amm)
    REAL(8) :: lds              ! detritus recycling rate (silicate, Ds => S)
    REAL(8) :: invdin_eps       ! total inorganic N inversed and eps
    real(8) :: o2flx          ! O2 surface flux
    real(8) :: pistonvel0     ! piston velocity
    real(8) :: srfammflx_col  ! Value for single column deposition rate
    real(8) :: srfnitflx_col  ! Value for single column deposition rate
    real(8) :: srfphosflx_col  ! Value for single column deposition rate
    real(8) :: srfamm_with_totatmos_Nflx_col  ! Value for single column deposition rate
    real(8) :: srfnit_with_totatmos_Nflx_col  ! Value for single column deposition rate
    real(8) :: srfamm_with_agri_Nflx_col  ! Value for single column deposition rate
    real(8) :: srfnit_with_agri_Nflx_col  ! Value for single column deposition rate
    
    !--------------------------------------------------------------------------!
    ! processes
    !--------------------------------------------------------------------------!
    real(8) :: grazing_miz_on_dia     ! description
    real(8) :: grazing_miz_on_flag     ! description
    real(8) :: grazing_miz_on_cyano     ! description
    real(8) :: respiration_miz     ! description
    real(8) :: mortality_miz       ! description
    real(8) :: grazing_mez_on_dia     ! description
    real(8) :: grazing_mez_on_flag     ! description
    real(8) :: grazing_mez_on_cyano     ! description
    real(8) :: grazing_mez_on_miz     ! description
    real(8) :: respiration_mez     ! description
    real(8) :: mortality_mez       ! description
    real(8) :: uptake_amm_by_dia     ! description
    real(8) :: uptake_nit_by_dia     ! description
    real(8) :: respiration_dia     ! description
    real(8) :: mortality_dia       ! description
    real(8) :: uptake_amm_by_flag     ! description
    real(8) :: uptake_nit_by_flat     ! description
    real(8) :: respiration_flag     ! description
    real(8) :: mortality_flag      ! description
    real(8) :: uptake_phos_by_cyano     ! description
    real(8) :: respiration_cyano     ! description
    real(8) :: mortality_cyano     ! description
    real(8) :: recycling_detritus_to_n     ! description
    real(8) :: recycling_detritus_to_sil     ! description
    real(8) :: nitrification       ! description
    real(8) :: degradation_ldon     ! description
    real(8) :: grazing_miz_on_dia_totatmos_N     ! description
    real(8) :: grazing_miz_on_flag_totatmos_N     ! description
    real(8) :: grazing_miz_on_cyano_totatmos_N     ! description
    real(8) :: respiration_miz_totatmos_N     ! description
    real(8) :: mortality_miz_totatmos_N     ! description
    real(8) :: grazing_mez_on_dia_totatmos_N     ! description
    real(8) :: grazing_mez_on_flag_totatmos_N     ! description
    real(8) :: grazing_mez_on_cyano_totatmos_N     ! description
    real(8) :: grazing_mez_on_miz_totatmos_N     ! description
    real(8) :: respiration_mez_totatmos_N     ! description
    real(8) :: mortality_mez_totatmos_N     ! description
    real(8) :: uptake_amm_by_dia_totatmos_N     ! description
    real(8) :: uptake_nit_by_dia_totatmos_N     ! description
    real(8) :: respiration_dia_totatmos_N     ! description
    real(8) :: mortality_dia_totatmos_N     ! description
    real(8) :: uptake_amm_by_flag_totatmos_N     ! description
    real(8) :: uptake_nit_by_flat_totatmos_N     ! description
    real(8) :: respiration_flag_totatmos_N     ! description
    real(8) :: mortality_flag_totatmos_N     ! description
    real(8) :: respiration_cyano_totatmos_N     ! description
    real(8) :: mortality_cyano_totatmos_N     ! description
    real(8) :: recycling_detritus_to_n_totatmos_N     ! description
    real(8) :: nitrification_totatmos_N     ! description
    real(8) :: degradation_ldon_totatmos_N     ! description
    real(8) :: grazing_miz_on_dia_agri_N     ! description
    real(8) :: grazing_miz_on_flag_agri_N     ! description
    real(8) :: grazing_miz_on_cyano_agri_N     ! description
    real(8) :: respiration_miz_agri_N     ! description
    real(8) :: mortality_miz_agri_N     ! description
    real(8) :: grazing_mez_on_dia_agri_N     ! description
    real(8) :: grazing_mez_on_flag_agri_N     ! description
    real(8) :: grazing_mez_on_cyano_agri_N     ! description
    real(8) :: grazing_mez_on_miz_agri_N     ! description
    real(8) :: respiration_mez_agri_N     ! description
    real(8) :: mortality_mez_agri_N     ! description
    real(8) :: uptake_amm_by_dia_agri_N     ! description
    real(8) :: uptake_nit_by_dia_agri_N     ! description
    real(8) :: respiration_dia_agri_N     ! description
    real(8) :: mortality_dia_agri_N     ! description
    real(8) :: uptake_amm_by_flag_agri_N     ! description
    real(8) :: uptake_nit_by_flat_agri_N     ! description
    real(8) :: respiration_flag_agri_N     ! description
    real(8) :: mortality_flag_agri_N     ! description
    real(8) :: respiration_cyano_agri_N     ! description
    real(8) :: mortality_cyano_agri_N     ! description
    real(8) :: recycling_detritus_to_n_agri_N     ! description
    real(8) :: nitrification_agri_N     ! description
    real(8) :: degradation_ldon_agri_N     ! description
 
    !--------------------------------------------------------------------------!
    ! Loop counters
    !--------------------------------------------------------------------------!
    integer(4) :: i,k,j,nsrf,mi0,mi,mww,mnn,n2dl,n2du,kb,itra,mil,miu

    !--------------------------------------------------------------------------!
    ! Local parameters
    !--------------------------------------------------------------------------!
    real(8), parameter :: eps   = 0.00000000001_8            ! no division by 0
    real(8), parameter :: pd2ps = 1.0_8/(24.0_8*60.0_8*60.0_8)
    real(8), parameter :: pi    = 3.14159265358979323846_8 
    real(8), parameter :: hnull = 0.1_8

    !--------------------------------------------------------------------------!
    ! Bio model params are in units of days, so we scale the coefficients by
    ! a simple scale of the time step size and hope that time-params only
    ! enter the model as coefficients to dt:
    !--------------------------------------------------------------------------!
    dt = dts*pd2ps

    !--------------------------------------------------------------------------!
    ! some local parameters to ease computation
    !--------------------------------------------------------------------------!

    ksdenit2 = ksdenit*ksdenit
    alphap2  = alphap*alphap
    alphaf2  = alphaf*alphaf
    alphab2  = alphab*alphab
    alprfs2  = alphap2*rfs*rfs !silicate redfield ratio

    call domp_get_domain(kh, 1, iw2, n2dl, n2du, idx)

!$OMP MASTER
    !--------------------------------------------------------------------------!
    ! Storing input tracers in old array ecoa
    !--------------------------------------------------------------------------!
    ! init dummy vars for unused parameters
    ben(0)        = zero 
    benderiv(0)   = zero 
    deriv(0)      = zero
    ecoa(:,0)     = zero
    diag_light(0) = zero 
    diag_kpar(0)  = zero
!$OMP END MASTER

    light(:) = zero
  
    ! if (ntben == 2) then
      !- Handle drying cases, only valid for ntben==2
      ! FIXME: drying_case must be re-defined when ntben>2
      ! NEUMANND: solved
      call drying_case(n2dl, n2du, kh, h, hnull, ts, benthos, diag_kpar,      &
                       diag_light)
    ! endif

    do nsrf=n2dl,n2du
      kb = kh(nsrf)
     
      !- Skip drying situation (done above)
      if (kb <= 1 .and. h(nsrf) < hnull) cycle

      if (kb > 1) then
        mi0 = mcol(nsrf)
        mil = mi0
        miu = mi0 + kb - 2
      endif
      i = ind(1,nsrf)
      j = ind(2,nsrf)    
 
      ! Obtain solar insolation
      call shortwaverad(i, j, xkoor, ykoor, cloud(nsrf), humid(nsrf),          &
                        casus(nsrf), swrad )
      
      if (rivmap(nsrf) == 0) then
        !keeping the old array without going through the calculations for 
        !specified river areas
        ecoa(1,idx_amm            ) = max(ts(1,idx_amm),zero)
        ecoa(1,idx_nit            ) = max(ts(1,idx_nit),zero)
        ecoa(1,idx_phos           ) = max(ts(1,idx_phos),zero)
        ecoa(1,idx_sil            ) = max(ts(1,idx_sil),zero)
        ecoa(1,idx_dia            ) = max(ts(1,idx_dia),zero)
        ecoa(1,idx_flag           ) = max(ts(1,idx_flag),zero)
        ecoa(1,idx_cyano          ) = max(ts(1,idx_cyano),zero)
        ecoa(1,idx_mez            ) = max(ts(1,idx_mez),zero)
        ecoa(1,idx_miz            ) = max(ts(1,idx_miz),zero)
        ecoa(1,idx_det            ) = max(ts(1,idx_det),zero)
        ecoa(1,idx_dets           ) = max(ts(1,idx_dets),zero)
        ecoa(1,idx_ldon           ) = max(ts(1,idx_ldon),zero)
        !oxygen can be negative since it represent hydrogen sulphite
        ! colored tracer should be negative due to some numeric stability
        ecoa(1,idx_oxy            ) = ts(1,idx_oxy)
        ecoa(1,idx_amm_with_totatmos_N) = ts(1,idx_amm_with_totatmos_N)
        ecoa(1,idx_nit_with_totatmos_N) = ts(1,idx_nit_with_totatmos_N)
        ecoa(1,idx_dia_with_totatmos_N) = ts(1,idx_dia_with_totatmos_N)
        ecoa(1,idx_flag_with_totatmos_N) = ts(1,idx_flag_with_totatmos_N)
        ecoa(1,idx_cyano_with_totatmos_N) = ts(1,idx_cyano_with_totatmos_N)
        ecoa(1,idx_mez_with_totatmos_N) = ts(1,idx_mez_with_totatmos_N)
        ecoa(1,idx_miz_with_totatmos_N) = ts(1,idx_miz_with_totatmos_N)
        ecoa(1,idx_det_with_totatmos_N) = ts(1,idx_det_with_totatmos_N)
        ecoa(1,idx_ldon_with_totatmos_N) = ts(1,idx_ldon_with_totatmos_N)
        ecoa(1,idx_amm_with_agri_N) = ts(1,idx_amm_with_agri_N)
        ecoa(1,idx_nit_with_agri_N) = ts(1,idx_nit_with_agri_N)
        ecoa(1,idx_dia_with_agri_N) = ts(1,idx_dia_with_agri_N)
        ecoa(1,idx_flag_with_agri_N) = ts(1,idx_flag_with_agri_N)
        ecoa(1,idx_cyano_with_agri_N) = ts(1,idx_cyano_with_agri_N)
        ecoa(1,idx_mez_with_agri_N) = ts(1,idx_mez_with_agri_N)
        ecoa(1,idx_miz_with_agri_N) = ts(1,idx_miz_with_agri_N)
        ecoa(1,idx_det_with_agri_N) = ts(1,idx_det_with_agri_N)
        ecoa(1,idx_ldon_with_agri_N) = ts(1,idx_ldon_with_agri_N)
        s_k(1)        = s(nsrf)
        light_o(1)    = diag_light(nsrf)
        hi(1)         = h(nsrf)
        eddyvdiff(1)  = dispvc(nsrf)

        ! REPLACED DO LOOP ...
        !~ do k=2,kb
          !~ mi = mi0+k-2
          !~ 
          !~ ecoa(k,idx_amm            ) = max(ts(mi,idx_amm),zero)
          !~ ecoa(k,idx_nit            ) = max(ts(mi,idx_nit),zero)
          !~ ecoa(k,idx_phos           ) = max(ts(mi,idx_phos),zero)
          !~ ecoa(k,idx_sil            ) = max(ts(mi,idx_sil),zero)
          !~ ecoa(k,idx_dia            ) = max(ts(mi,idx_dia),zero)
          !~ ecoa(k,idx_flag           ) = max(ts(mi,idx_flag),zero)
          !~ ecoa(k,idx_cyano          ) = max(ts(mi,idx_cyano),zero)
          !~ ecoa(k,idx_mez            ) = max(ts(mi,idx_mez),zero)
          !~ ecoa(k,idx_miz            ) = max(ts(mi,idx_miz),zero)
          !~ ecoa(k,idx_det            ) = max(ts(mi,idx_det),zero)
          !~ ecoa(k,idx_dets           ) = max(ts(mi,idx_dets),zero)
          !~ ecoa(k,idx_ldon           ) = max(ts(mi,idx_ldon),zero)
          !~ !oxygen can be negative since it represent hydrogen sulphite
          !~ ! colored tracer should be negative due to some numeric stability
          !~ ecoa(k,idx_oxy            ) = ts(mi,idx_oxy)
          !~ ecoa(k,idx_amm_with_totatmos_N) = ts(mi,idx_amm_with_totatmos_N)
          !~ ecoa(k,idx_nit_with_totatmos_N) = ts(mi,idx_nit_with_totatmos_N)
          !~ ecoa(k,idx_dia_with_totatmos_N) = ts(mi,idx_dia_with_totatmos_N)
          !~ ecoa(k,idx_flag_with_totatmos_N) = ts(mi,idx_flag_with_totatmos_N)
          !~ ecoa(k,idx_cyano_with_totatmos_N) = ts(mi,idx_cyano_with_totatmos_N)
          !~ ecoa(k,idx_mez_with_totatmos_N) = ts(mi,idx_mez_with_totatmos_N)
          !~ ecoa(k,idx_miz_with_totatmos_N) = ts(mi,idx_miz_with_totatmos_N)
          !~ ecoa(k,idx_det_with_totatmos_N) = ts(mi,idx_det_with_totatmos_N)
          !~ ecoa(k,idx_ldon_with_totatmos_N) = ts(mi,idx_ldon_with_totatmos_N)
          !~ ecoa(k,idx_amm_with_agri_N) = ts(mi,idx_amm_with_agri_N)
          !~ ecoa(k,idx_nit_with_agri_N) = ts(mi,idx_nit_with_agri_N)
          !~ ecoa(k,idx_dia_with_agri_N) = ts(mi,idx_dia_with_agri_N)
          !~ ecoa(k,idx_flag_with_agri_N) = ts(mi,idx_flag_with_agri_N)
          !~ ecoa(k,idx_cyano_with_agri_N) = ts(mi,idx_cyano_with_agri_N)
          !~ ecoa(k,idx_mez_with_agri_N) = ts(mi,idx_mez_with_agri_N)
          !~ ecoa(k,idx_miz_with_agri_N) = ts(mi,idx_miz_with_agri_N)
          !~ ecoa(k,idx_det_with_agri_N) = ts(mi,idx_det_with_agri_N)
          !~ ecoa(k,idx_ldon_with_agri_N) = ts(mi,idx_ldon_with_agri_N)
          !~ s_k(k)        = s(mi)
          !~ light_o(k)    = diag_light(mi)
          !~ hi(k)         = h(mi)
          !~ eddyvdiff(k)  = dispvc(mi)
        !~ enddo
        ! ... BY ...
        if (kb > 1) then
          ecoa(2:kb,idx_amm            ) = max(ts(mil:miu,idx_amm),zero)
          ecoa(2:kb,idx_nit            ) = max(ts(mil:miu,idx_nit),zero)
          ecoa(2:kb,idx_phos           ) = max(ts(mil:miu,idx_phos),zero)
          ecoa(2:kb,idx_sil            ) = max(ts(mil:miu,idx_sil),zero)
          ecoa(2:kb,idx_dia            ) = max(ts(mil:miu,idx_dia),zero)
          ecoa(2:kb,idx_flag           ) = max(ts(mil:miu,idx_flag),zero)
          ecoa(2:kb,idx_cyano          ) = max(ts(mil:miu,idx_cyano),zero)
          ecoa(2:kb,idx_mez            ) = max(ts(mil:miu,idx_mez),zero)
          ecoa(2:kb,idx_miz            ) = max(ts(mil:miu,idx_miz),zero)
          ecoa(2:kb,idx_det            ) = max(ts(mil:miu,idx_det),zero)
          ecoa(2:kb,idx_dets           ) = max(ts(mil:miu,idx_dets),zero)
          ecoa(2:kb,idx_ldon           ) = max(ts(mil:miu,idx_ldon),zero)
          !oxygen can be negative since it represent hydrogen sulphite
          ! colored tracer should be negative due to some numeric stability
          ecoa(2:kb,idx_oxy            ) = ts(mil:miu,idx_oxy)
          ecoa(2:kb,idx_amm_with_totatmos_N) = ts(mil:miu,idx_amm_with_totatmos_N)
          ecoa(2:kb,idx_nit_with_totatmos_N) = ts(mil:miu,idx_nit_with_totatmos_N)
          ecoa(2:kb,idx_dia_with_totatmos_N) = ts(mil:miu,idx_dia_with_totatmos_N)
          ecoa(2:kb,idx_flag_with_totatmos_N) = ts(mil:miu,idx_flag_with_totatmos_N)
          ecoa(2:kb,idx_cyano_with_totatmos_N) = ts(mil:miu,idx_cyano_with_totatmos_N)
          ecoa(2:kb,idx_mez_with_totatmos_N) = ts(mil:miu,idx_mez_with_totatmos_N)
          ecoa(2:kb,idx_miz_with_totatmos_N) = ts(mil:miu,idx_miz_with_totatmos_N)
          ecoa(2:kb,idx_det_with_totatmos_N) = ts(mil:miu,idx_det_with_totatmos_N)
          ecoa(2:kb,idx_ldon_with_totatmos_N) = ts(mil:miu,idx_ldon_with_totatmos_N)
          ecoa(2:kb,idx_amm_with_agri_N) = ts(mil:miu,idx_amm_with_agri_N)
          ecoa(2:kb,idx_nit_with_agri_N) = ts(mil:miu,idx_nit_with_agri_N)
          ecoa(2:kb,idx_dia_with_agri_N) = ts(mil:miu,idx_dia_with_agri_N)
          ecoa(2:kb,idx_flag_with_agri_N) = ts(mil:miu,idx_flag_with_agri_N)
          ecoa(2:kb,idx_cyano_with_agri_N) = ts(mil:miu,idx_cyano_with_agri_N)
          ecoa(2:kb,idx_mez_with_agri_N) = ts(mil:miu,idx_mez_with_agri_N)
          ecoa(2:kb,idx_miz_with_agri_N) = ts(mil:miu,idx_miz_with_agri_N)
          ecoa(2:kb,idx_det_with_agri_N) = ts(mil:miu,idx_det_with_agri_N)
          ecoa(2:kb,idx_ldon_with_agri_N) = ts(mil:miu,idx_ldon_with_agri_N)
          s_k(2:kb)        = s(mil:miu)
          light_o(2:kb)    = diag_light(mil:miu)
          hi(2:kb)         = h(mil:miu)
          eddyvdiff(2:kb)  = dispvc(mil:miu)
        endif
        
        hsrf = max(hi(1),hnull)
        
        !! CURRENTLY, WE DO NOT USE benthos BELOW!!!
        ! not all benthos species should be set to zero if negativ
        ! benthos(nsrf,1:ntben) = max(benthos(nsrf,1:ntben),zero)
        ! benthos(nsrf,1:ntben) = max(benthos(nsrf,1:ntben),zero)
        
        !  warning if par_calc_don = .true. and ppismmod = .true. should
        !  be implemented

        if (par_calc_don) then
          call bio_light(kb,hi,hsrf,swrad,ecoa,s_k,light_o,kpar,light)
        else if (ppispmmod) then
          call bio_light(kb,hi,hsrf,swrad,ecoa,eddyvdiff(:),kpar,light)
        else
          call bio_light(kb,hi,hsrf,swrad,ecoa,kpar,light)
        endif

        diag_light(nsrf) = light(1)
        diag_kpar(nsrf)  = kpar(1)
        if (kb > 1) then
          diag_light(mi0:mi0+kb-2) = light(2:kb)
          diag_kpar(mi0:mi0+kb-2)  = kpar(2:kb)
        endif
      else

        icefree      = casus(nsrf)
        ben(1:ntben) = benthos(nsrf,1:ntben)
        rfr          = npr(nsrf)
        rfr2         = rfr*rfr
        albrfr2      = alphab2*rfr2

        hi(1)           = h(nsrf)
        t_k(1)          = max(t(nsrf),zero) ! all bioeq's are based on temp>0
        s_k(1)          = s(nsrf)
        light_o(1)      = diag_light(nsrf)
        ecoa(1,1:ntbio) = ts(nsrf,1:ntbio)
        eddyvdiff(1)    = dispvc(nsrf)

        srfammflx_col   = srfammflx(nsrf)
        srfnitflx_col   = srfnitflx(nsrf)
        srfphosflx_col   = srfphosflx(nsrf)
        srfamm_with_totatmos_Nflx_col   = srfamm_with_totatmos_Nflx(nsrf)
        srfnit_with_totatmos_Nflx_col   = srfnit_with_totatmos_Nflx(nsrf)
        srfamm_with_agri_Nflx_col   = srfamm_with_agri_Nflx(nsrf)
        srfnit_with_agri_Nflx_col   = srfnit_with_agri_Nflx(nsrf)
        
        !~ do k=2,kb
          !~ mi = mi0+k-2
          !~ hi(k)           = h(mi)
          !~ t_k(k)          = max(t(mi),zero) ! all bioeq's are based on temp>0
          !~ s_k(k)          = s(mi)
          !~ light_o(k)      = diag_light(mi)
          !~ ecoa(k,1:ntbio) = ts(mi,1:ntbio)
          !~ eddyvdiff(k)    = dispvc(mi)
        !~ enddo
        ! REPLACED BY
        if (kb > 1) then
          hi(2:kb)           = h(mil:miu)
          t_k(2:kb)          = max(t(mil:miu),zero) ! all bioeq's are based on temp>0
          s_k(2:kb)          = s(mil:miu)
          light_o(2:kb)      = diag_light(mil:miu)
          ecoa(2:kb,1:ntbio) = ts(mil:miu,1:ntbio)
          eddyvdiff(2:kb)    = dispvc(mil:miu)
        endif

        ! hsrf/hbot and bottom velocity
        hsrf = max(hi(1),hnull)
        mww = 0
        mnn = 0
        if (kb == 1) then
          hbot = hsrf
          mi  = msrf(i,j)
          mww = msrf(i,j-1)
          mnn = msrf(i-1,j)
        else
          hbot = hi(kb)
          mi = mi0+kb-2
          if(kh(msrf(i,j-1)) >= kb) mww = mcol(msrf(i,j-1))+kb-2
          if(kh(msrf(i-1,j)) >= kb) mnn = mcol(msrf(i-1,j))+kb-2
        endif
        ust   = half*(u(mi)+u(mww))
        vst   = half*(v(mi)+v(mnn))
        velkb = sqrt(ust*ust+vst*vst)

        !----------------------------------------------------------------------!
        ! Sinking and rising of particles and algae
        !----------------------------------------------------------------------!
        if ( kb > 1 ) then
          ! biological vertical movement neglected for single layer column
          call vmove(kb,hi,ecoa(:,idx_dia            ),wdiaz,dt,irc)
          call vmove(kb,hi,ecoa(:,idx_dia_with_totatmos_N),wdiaz,dt,irc)
          call vmove(kb,hi,ecoa(:,idx_dia_with_agri_N),wdiaz,dt,irc)
          call vmove(kb,hi,ecoa(:,idx_cyano          ),wcyanoz,dt,irc)
          call vmove(kb,hi,ecoa(:,idx_cyano_with_totatmos_N),wcyanoz,dt,irc)
          call vmove(kb,hi,ecoa(:,idx_cyano_with_agri_N),wcyanoz,dt,irc)
          call vmove(kb,hi,ecoa(:,idx_det            ),wdetz,dt,irc)
          call vmove(kb,hi,ecoa(:,idx_det_with_totatmos_N),wdetz,dt,irc)
          call vmove(kb,hi,ecoa(:,idx_det_with_agri_N),wdetz,dt,irc)
          call vmove(kb,hi,ecoa(:,idx_dets           ),wdetsz,dt,irc)
          if (irc /= 0) return
        endif
        call sedimentation_resuspension(hbot, dt, velkb, ecoa, ben, kb)

        !----------------------------------------------------------------------!
        !  Light limitation considering self shading
        !----------------------------------------------------------------------!
        !  warning if par_calc_don = .true. and ppismmod = .true. should
        !  be implemented
        if (par_calc_don) then
          call bio_light(kb,hi,hsrf,swrad,ecoa,s_k,light_o,kpar,light)
        else if (ppispmmod) then
          call bio_light(kb,hi,hsrf,swrad,ecoa,eddyvdiff(:),kpar,light)
        else
          call bio_light(kb,hi,hsrf,swrad,ecoa,kpar,light)
        endif

        !----------------------------------------------------------------------!
        !
        !            surface fluxes
        !
        !----------------------------------------------------------------------!
        !fixme: srf flux should probably be a function of the ice concentration
        !       nutrient should also be accumulated on the snow
        if (icefree) then
          osat = oxygensat(t_k(1),s_k(1)) * Onorm_inv  !using biogas
               
          !--------------------------------------------------------------------!
          !    Gas transfer at surface using a fixed piston-velocity 
          !--------------------------------------------------------------------!
          oxy  = ecoa(1,idx_oxy) ! NEUMANND: ntbio -> idx_oxy
          wui  = wu(nsrf)
          wvi  = wv(nsrf)
          wind = sqrt(wui*wui + wvi*wvi)      ! wind velocity
          pistonvel0 = pistonvel(wind,t_k(1)) ! using a variable piston-velocity
          o2flx      = pistonvel0*(osat - oxy) ! o2_rel * m/day
          ecoa(1,idx_oxy) = oxy + dt*o2flx/hsrf

          !--------------------------------------------------------------------!
          !   Atmospheric deposistion of nit, amm and phos
          !--------------------------------------------------------------------!
          ecoa(1,idx_amm)  = ecoa(1,idx_amm)  + srfammflx_col/hsrf*Nnorm_inv*dt
          ecoa(1,idx_nit)  = ecoa(1,idx_nit)  + srfnitflx_col/hsrf*Nnorm_inv*dt
          ecoa(1,idx_phos)  = ecoa(1,idx_phos)  + srfphosflx_col/hsrf*Nnorm_inv*dt
          ecoa(1,idx_amm_with_totatmos_N)  = ecoa(1,idx_amm_with_totatmos_N)  + srfamm_with_totatmos_Nflx_col/hsrf*Nnorm_inv*dt
          ecoa(1,idx_nit_with_totatmos_N)  = ecoa(1,idx_nit_with_totatmos_N)  + srfnit_with_totatmos_Nflx_col/hsrf*Nnorm_inv*dt
          ecoa(1,idx_amm_with_agri_N)  = ecoa(1,idx_amm_with_agri_N)  + srfamm_with_agri_Nflx_col/hsrf*Nnorm_inv*dt
          ecoa(1,idx_nit_with_agri_N)  = ecoa(1,idx_nit_with_agri_N)  + srfnit_with_agri_Nflx_col/hsrf*Nnorm_inv*dt
        endif
      
        !- done with surface fluxes, now loop through the water column
        do k=1,kb
          temp   = t_k(k)
          sali   = s_k(k)
          lightk = light(k)

          amm             =  ecoa(k,idx_amm            )
          nit             =  ecoa(k,idx_nit            )
          phos            =  ecoa(k,idx_phos           )
          sil             =  ecoa(k,idx_sil            )
          dia             =  ecoa(k,idx_dia            )
          flag            =  ecoa(k,idx_flag           )
          cyano           =  ecoa(k,idx_cyano          )
          mez             =  ecoa(k,idx_mez            )
          miz             =  ecoa(k,idx_miz            )
          det             =  ecoa(k,idx_det            )
          dets            =  ecoa(k,idx_dets           )
          ldon            =  ecoa(k,idx_ldon           )
          oxy             =  ecoa(k,idx_oxy            )
          amm_with_totatmos_N =  ecoa(k,idx_amm_with_totatmos_N)
          nit_with_totatmos_N =  ecoa(k,idx_nit_with_totatmos_N)
          dia_with_totatmos_N =  ecoa(k,idx_dia_with_totatmos_N)
          flag_with_totatmos_N =  ecoa(k,idx_flag_with_totatmos_N)
          cyano_with_totatmos_N =  ecoa(k,idx_cyano_with_totatmos_N)
          mez_with_totatmos_N =  ecoa(k,idx_mez_with_totatmos_N)
          miz_with_totatmos_N =  ecoa(k,idx_miz_with_totatmos_N)
          det_with_totatmos_N =  ecoa(k,idx_det_with_totatmos_N)
          ldon_with_totatmos_N =  ecoa(k,idx_ldon_with_totatmos_N)
          amm_with_agri_N =  ecoa(k,idx_amm_with_agri_N)
          nit_with_agri_N =  ecoa(k,idx_nit_with_agri_N)
          dia_with_agri_N =  ecoa(k,idx_dia_with_agri_N)
          flag_with_agri_N =  ecoa(k,idx_flag_with_agri_N)
          cyano_with_agri_N =  ecoa(k,idx_cyano_with_agri_N)
          mez_with_agri_N =  ecoa(k,idx_mez_with_agri_N)
          miz_with_agri_N =  ecoa(k,idx_miz_with_agri_N)
          det_with_agri_N =  ecoa(k,idx_det_with_agri_N)
          ldon_with_agri_N =  ecoa(k,idx_ldon_with_agri_N)

          if (oxy <= zero) then
            !------------------------------------------------------------------!
            !
            !            Hypoxia
            !
            !------------------------------------------------------------------!
            nutlimc         = zero
            ppikb           = zero
            rb              = zero
            nutlimf         = zero
            ppikf           = zero
            rf              = zero
            nutlimd         = zero
            ppikp           = zero
            rp              = zero
            nf              = zero
            foodmez         = zero
            temp2           = zero
            mezgt           = zero
            mezgscal        = zero
            meztotgraz      = zero
            foodmiz         = zero
            mizgt           = zero
            mizgscal        = zero
            miztotgraz      = zero
            lpn             = zero
            denitnscal      = nit*nit/(ksdenit2+nit*nit)
            ldn_N           = NNdenit * denitnscal
            ldn_O           = ONamup * (one - denitnscal)
            lpd             = ten * deltao
            tldmez          = ten*sigma_b
            tldmiz          = ten*sigma_b
            fac             = exp(q10_rec*temp)/ten
          ELSE
            !------------------------------------------------------------------!
            !
            !            Oxic conditions
            !
            !------------------------------------------------------------------!
            nutlimc         = phos*phos/(albrfr2+phos*phos) 
            ppikb           = min(lightk/min_ocyanol,one)   
            rb              = rb0*min(nutlimc, ppikb*exp(one-ppikb))*four/(four + exp(cyanotll - temp))*(atan(-sali+cyanosul)/pi+half)
            nutlimf         = min(((amm+nit)*(amm+nit))/(alphaf2+(amm+nit)*(amm+nit)), nutlimc)
            ppikf           = min(lightk/min_oflagl,one)    
            rf              = rf0*min(nutlimf,ppikf*exp(one-ppikf))*flagtsc*exp(flagtll*temp)
            nutlimd         = min(nutlimf, sil*sil/(alprfs2+sil*sil))
            ppikp           = min(lightk/min_odial,one)     
            rp              = rp0*min(nutlimd, ppikp*exp(one-ppikp))
            nf              = rnit*oxy/(ksnit + oxy)*exp(anit*temp)
            foodmez         = mezprefdia*dia + mezprefflag*flag + mezprefcyano*cyano + mezprefmiz*miz
            temp2           = temp * temp                   
            mezgt           = half + mez_tscale*temp2/(meztk+temp2)
            mezgscal        = mezgraz*foodmez/(foodmez*foodmez+mezgrazf*mezgrazf)
            meztotgraz      = mezgscal*mezgt                
            foodmiz         = mizprefdia*dia + mizprefflag*flag + mizprefcyano*cyano
            mizgt           = one  + temp2/(miztk+temp2)    
            mizgscal        = mizgraz*foodmiz/(foodmiz*foodmiz+mizgrazf*mizgrazf)
            miztotgraz      = mizgscal*mizgt                
            lpn             = nb                            
            denitnscal      = zero                          
            ldn_N           = zero                          
            ldn_O           = ONamup                        
            lpd             = deltao                        
            tldmez          = ldmez + mezcl*mez             
            tldmiz          = ldmiz                         
            fac             = exp(q10_rec*temp)             
         
            ! NOTE!!The light description below is problematic since we get very
            ! low ppik for cyano and dia at the surface...consider change the 
            ! formulation to this one:
            !   Ifac  = light/min_odial
            !   ppikf = Ifac/sqrt(1+Ifac*Ifac)

            ! fixme: check formulation of T dependence of bluegreen algae
          endif

          ldn             = dn*fac                        
          lds             = ds*fac                        
          invdin_eps      = one / (nit + amm + eps)       


          !--------------------------------------------------------------------!
          !
          !            biochemical tracer dynamics
          !
          !--------------------------------------------------------------------!
          grazing_miz_on_dia = mizprefdia  * miztotgraz * miz * dia
          grazing_miz_on_flag = mizprefflag * miztotgraz * miz * flag
          grazing_miz_on_cyano = mizprefcyano * miztotgraz * miz * cyano
          respiration_miz = lnmiz * miztotgraz * foodmiz * miz
          mortality_miz   = tldmiz * miz                  
          grazing_mez_on_dia = mezprefdia  * meztotgraz * mez * dia
          grazing_mez_on_flag = mezprefflag * meztotgraz * mez * flag
          grazing_mez_on_cyano = mezprefcyano * meztotgraz * mez * cyano
          grazing_mez_on_miz = mezprefmiz  * meztotgraz * mez * miz
          respiration_mez = lnmez * meztotgraz * foodmez * mez
          mortality_mez   = tldmez * mez                  
          uptake_amm_by_dia = rp * (dia+p0) * amm * invdin_eps
          uptake_nit_by_dia = rp * (dia+p0) * nit * invdin_eps
          respiration_dia = lpn * dia                     
          mortality_dia   = lpd * dia                     
          uptake_amm_by_flag = rf * (flag+f0) * amm * invdin_eps
          uptake_nit_by_flat = rf * (flag+f0) * nit * invdin_eps
          respiration_flag = lpn * flag                    
          mortality_flag  = lpd * flag                    
          uptake_phos_by_cyano = rb * (cyano + b0)             
          respiration_cyano = lpn * cyano                   
          mortality_cyano = lpd * cyano                   
          recycling_detritus_to_n = ldn * det                     
          recycling_detritus_to_sil = lds * dets                    
          nitrification   = nf * amm                      
          degradation_ldon = ldon_tor * ldon               
          grazing_miz_on_dia_totatmos_N = grazing_miz_on_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,dia_with_totatmos_N/max(0.00000000001_8,dia)))) / ((1.0_8)*(1))
          grazing_miz_on_flag_totatmos_N = grazing_miz_on_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,flag_with_totatmos_N/max(0.00000000001_8,flag)))) / ((1.0_8)*(1))
          grazing_miz_on_cyano_totatmos_N = grazing_miz_on_cyano * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,cyano_with_totatmos_N/max(0.00000000001_8,cyano)))) / ((1.0_8)*(1))
          respiration_miz_totatmos_N = respiration_miz * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,miz_with_totatmos_N/max(0.00000000001_8,miz)))) / ((1.0_8)*(1))
          mortality_miz_totatmos_N = mortality_miz * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,miz_with_totatmos_N/max(0.00000000001_8,miz)))) / ((1.0_8)*(1))
          grazing_mez_on_dia_totatmos_N = grazing_mez_on_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,dia_with_totatmos_N/max(0.00000000001_8,dia)))) / ((1.0_8)*(1))
          grazing_mez_on_flag_totatmos_N = grazing_mez_on_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,flag_with_totatmos_N/max(0.00000000001_8,flag)))) / ((1.0_8)*(1))
          grazing_mez_on_cyano_totatmos_N = grazing_mez_on_cyano * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,cyano_with_totatmos_N/max(0.00000000001_8,cyano)))) / ((1.0_8)*(1))
          grazing_mez_on_miz_totatmos_N = grazing_mez_on_miz * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,miz_with_totatmos_N/max(0.00000000001_8,miz)))) / ((1.0_8)*(1))
          respiration_mez_totatmos_N = respiration_mez * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,mez_with_totatmos_N/max(0.00000000001_8,mez)))) / ((1.0_8)*(1))
          mortality_mez_totatmos_N = mortality_mez * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,mez_with_totatmos_N/max(0.00000000001_8,mez)))) / ((1.0_8)*(1))
          uptake_amm_by_dia_totatmos_N = uptake_amm_by_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,amm_with_totatmos_N/max(0.00000000001_8,amm)))) / ((1.0_8)*(1))
          uptake_nit_by_dia_totatmos_N = uptake_nit_by_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,nit_with_totatmos_N/max(0.00000000001_8,nit)))) / ((1.0_8)*(1))
          respiration_dia_totatmos_N = respiration_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,dia_with_totatmos_N/max(0.00000000001_8,dia)))) / ((1.0_8)*(1))
          mortality_dia_totatmos_N = mortality_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,dia_with_totatmos_N/max(0.00000000001_8,dia)))) / ((1.0_8)*(1))
          uptake_amm_by_flag_totatmos_N = uptake_amm_by_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,amm_with_totatmos_N/max(0.00000000001_8,amm)))) / ((1.0_8)*(1))
          uptake_nit_by_flat_totatmos_N = uptake_nit_by_flat * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,nit_with_totatmos_N/max(0.00000000001_8,nit)))) / ((1.0_8)*(1))
          respiration_flag_totatmos_N = respiration_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,flag_with_totatmos_N/max(0.00000000001_8,flag)))) / ((1.0_8)*(1))
          mortality_flag_totatmos_N = mortality_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,flag_with_totatmos_N/max(0.00000000001_8,flag)))) / ((1.0_8)*(1))
          respiration_cyano_totatmos_N = respiration_cyano * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,cyano_with_totatmos_N/max(0.00000000001_8,cyano)))) / ((1.0_8)*(1))
          mortality_cyano_totatmos_N = mortality_cyano * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,cyano_with_totatmos_N/max(0.00000000001_8,cyano)))) / ((1.0_8)*(1))
          recycling_detritus_to_n_totatmos_N = recycling_detritus_to_n * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,det_with_totatmos_N/max(0.00000000001_8,det)))+(ldn_N)*(1)*max(0.0_8,min(1.0_8,nit_with_totatmos_N/max(0.00000000001_8,nit)))) / ((1.0_8)*(1)+(ldn_N)*(1))
          nitrification_totatmos_N = nitrification * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,amm_with_totatmos_N/max(0.00000000001_8,amm)))) / ((1.0_8)*(1))
          degradation_ldon_totatmos_N = degradation_ldon * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,ldon_with_totatmos_N/max(0.00000000001_8,ldon)))) / ((1.0_8)*(1))
          grazing_miz_on_dia_agri_N = grazing_miz_on_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,dia_with_agri_N/max(0.00000000001_8,dia)))) / ((1.0_8)*(1))
          grazing_miz_on_flag_agri_N = grazing_miz_on_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,flag_with_agri_N/max(0.00000000001_8,flag)))) / ((1.0_8)*(1))
          grazing_miz_on_cyano_agri_N = grazing_miz_on_cyano * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,cyano_with_agri_N/max(0.00000000001_8,cyano)))) / ((1.0_8)*(1))
          respiration_miz_agri_N = respiration_miz * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,miz_with_agri_N/max(0.00000000001_8,miz)))) / ((1.0_8)*(1))
          mortality_miz_agri_N = mortality_miz * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,miz_with_agri_N/max(0.00000000001_8,miz)))) / ((1.0_8)*(1))
          grazing_mez_on_dia_agri_N = grazing_mez_on_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,dia_with_agri_N/max(0.00000000001_8,dia)))) / ((1.0_8)*(1))
          grazing_mez_on_flag_agri_N = grazing_mez_on_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,flag_with_agri_N/max(0.00000000001_8,flag)))) / ((1.0_8)*(1))
          grazing_mez_on_cyano_agri_N = grazing_mez_on_cyano * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,cyano_with_agri_N/max(0.00000000001_8,cyano)))) / ((1.0_8)*(1))
          grazing_mez_on_miz_agri_N = grazing_mez_on_miz * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,miz_with_agri_N/max(0.00000000001_8,miz)))) / ((1.0_8)*(1))
          respiration_mez_agri_N = respiration_mez * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,mez_with_agri_N/max(0.00000000001_8,mez)))) / ((1.0_8)*(1))
          mortality_mez_agri_N = mortality_mez * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,mez_with_agri_N/max(0.00000000001_8,mez)))) / ((1.0_8)*(1))
          uptake_amm_by_dia_agri_N = uptake_amm_by_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,amm_with_agri_N/max(0.00000000001_8,amm)))) / ((1.0_8)*(1))
          uptake_nit_by_dia_agri_N = uptake_nit_by_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,nit_with_agri_N/max(0.00000000001_8,nit)))) / ((1.0_8)*(1))
          respiration_dia_agri_N = respiration_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,dia_with_agri_N/max(0.00000000001_8,dia)))) / ((1.0_8)*(1))
          mortality_dia_agri_N = mortality_dia * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,dia_with_agri_N/max(0.00000000001_8,dia)))) / ((1.0_8)*(1))
          uptake_amm_by_flag_agri_N = uptake_amm_by_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,amm_with_agri_N/max(0.00000000001_8,amm)))) / ((1.0_8)*(1))
          uptake_nit_by_flat_agri_N = uptake_nit_by_flat * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,nit_with_agri_N/max(0.00000000001_8,nit)))) / ((1.0_8)*(1))
          respiration_flag_agri_N = respiration_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,flag_with_agri_N/max(0.00000000001_8,flag)))) / ((1.0_8)*(1))
          mortality_flag_agri_N = mortality_flag * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,flag_with_agri_N/max(0.00000000001_8,flag)))) / ((1.0_8)*(1))
          respiration_cyano_agri_N = respiration_cyano * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,cyano_with_agri_N/max(0.00000000001_8,cyano)))) / ((1.0_8)*(1))
          mortality_cyano_agri_N = mortality_cyano * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,cyano_with_agri_N/max(0.00000000001_8,cyano)))) / ((1.0_8)*(1))
          recycling_detritus_to_n_agri_N = recycling_detritus_to_n * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,det_with_agri_N/max(0.00000000001_8,det)))+(ldn_N)*(1)*max(0.0_8,min(1.0_8,nit_with_agri_N/max(0.00000000001_8,nit)))) / ((1.0_8)*(1)+(ldn_N)*(1))
          nitrification_agri_N = nitrification * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,amm_with_agri_N/max(0.00000000001_8,amm)))) / ((1.0_8)*(1))
          degradation_ldon_agri_N = degradation_ldon * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,ldon_with_agri_N/max(0.00000000001_8,ldon)))) / ((1.0_8)*(1))

          !-----------------
          ! ammonium
          !-----------------
          deriv(idx_amm            ) = (                               &
                                       + (respiration_miz)*((one-ldon_frac))  &
                                       + (respiration_mez)*((one-ldon_frac))  &
                                       + (respiration_dia)*((one-ldon_frac))  &
                                       + (respiration_flag)*((one-ldon_frac))  &
                                       + (respiration_cyano)*((one-ldon_frac))  &
                                       + recycling_detritus_to_n       &
                                       + degradation_ldon              &
                                       - uptake_amm_by_dia             &
                                       - uptake_amm_by_flag            &
                                       - nitrification                 &
                                       )
          !-----------------
          ! nitrate
          !-----------------
          deriv(idx_nit            ) = (                               &
                                       + nitrification                 &
                                       - uptake_nit_by_dia             &
                                       - uptake_nit_by_flat            &
                                       - (recycling_detritus_to_n)*(ldn_N)  &
                                       )
          !-----------------
          ! phosphate
          !-----------------
          deriv(idx_phos           ) = (                               &
                                       + (respiration_miz)*(rfr)       &
                                       + (respiration_mez)*(rfr)       &
                                       + (respiration_dia)*(rfr)       &
                                       + (respiration_flag)*(rfr)      &
                                       + (respiration_cyano)*(rfr)     &
                                       + (recycling_detritus_to_n)*(rfr)  &
                                       - (uptake_amm_by_dia)*(rfr)     &
                                       - (uptake_nit_by_dia)*(rfr)     &
                                       - (uptake_amm_by_flag)*(rfr)    &
                                       - (uptake_nit_by_flat)*(rfr)    &
                                       - (uptake_phos_by_cyano)*(rfr)  &
                                       )
          !-----------------
          ! silicate
          !-----------------
          deriv(idx_sil            ) = (                               &
                                       + (respiration_dia)*(rfs)       &
                                       + recycling_detritus_to_sil     &
                                       - (uptake_amm_by_dia)*(rfs)     &
                                       - (uptake_nit_by_dia)*(rfs)     &
                                       )
          !-----------------
          ! diatoms
          !-----------------
          deriv(idx_dia            ) = (                               &
                                       + uptake_amm_by_dia             &
                                       + uptake_nit_by_dia             &
                                       - grazing_miz_on_dia            &
                                       - grazing_mez_on_dia            &
                                       - respiration_dia               &
                                       - mortality_dia                 &
                                       )
          !-----------------
          ! flagellates
          !-----------------
          deriv(idx_flag           ) = (                               &
                                       + uptake_amm_by_flag            &
                                       + uptake_nit_by_flat            &
                                       - grazing_miz_on_flag           &
                                       - grazing_mez_on_flag           &
                                       - respiration_flag              &
                                       - mortality_flag                &
                                       )
          !-----------------
          ! cyano bacteria
          !-----------------
          deriv(idx_cyano          ) = (                               &
                                       + uptake_phos_by_cyano          &
                                       - grazing_miz_on_cyano          &
                                       - grazing_mez_on_cyano          &
                                       - respiration_cyano             &
                                       - mortality_cyano               &
                                       )
          !-----------------
          ! meso zooplankton
          !-----------------
          deriv(idx_mez            ) = (                               &
                                       + grazing_mez_on_dia            &
                                       + grazing_mez_on_flag           &
                                       + grazing_mez_on_cyano          &
                                       + grazing_mez_on_miz            &
                                       - respiration_mez               &
                                       - mortality_mez                 &
                                       )
          !-----------------
          ! micro zooplakton
          !-----------------
          deriv(idx_miz            ) = (                               &
                                       + grazing_miz_on_dia            &
                                       + grazing_miz_on_flag           &
                                       + grazing_miz_on_cyano          &
                                       - respiration_miz               &
                                       - mortality_miz                 &
                                       - grazing_mez_on_miz            &
                                       )
          !-----------------
          ! Detritus
          !-----------------
          deriv(idx_det            ) = (                               &
                                       + mortality_miz                 &
                                       + mortality_mez                 &
                                       + mortality_dia                 &
                                       + mortality_flag                &
                                       + mortality_cyano               &
                                       - recycling_detritus_to_n       &
                                       )
          !-----------------
          ! S-Detritus
          !-----------------
          deriv(idx_dets           ) = (                               &
                                       + (grazing_miz_on_dia)*(rfs)    &
                                       + (grazing_mez_on_dia)*(rfs)    &
                                       + (mortality_dia)*(rfs)         &
                                       - recycling_detritus_to_sil     &
                                       )
          !-----------------
          ! ldon is the labile DON fraction
          !-----------------
          deriv(idx_ldon           ) = (                               &
                                       + (respiration_miz)*(ldon_frac)  &
                                       + (respiration_mez)*(ldon_frac)  &
                                       + (respiration_dia)*(ldon_frac)  &
                                       + (respiration_flag)*(ldon_frac)  &
                                       + (respiration_cyano)*(ldon_frac)  &
                                       - degradation_ldon              &
                                       )
          !-----------------
          ! Oxygen
          !-----------------
          deriv(idx_oxy            ) = (                               &
                                       + (uptake_amm_by_dia)*(NOR*ONamup)  &
                                       + (uptake_nit_by_dia)*(NOR*ONniup)  &
                                       + (uptake_amm_by_flag)*(NOR*ONamup)  &
                                       + (uptake_nit_by_flat)*(NOR*ONniup)  &
                                       + (uptake_phos_by_cyano)*(NOR*ONamup)  &
                                       - (respiration_miz)*(NOR*ONamup)  &
                                       - (respiration_mez)*(NOR*ONamup)  &
                                       - (respiration_dia)*(NOR*ONamup)  &
                                       - (respiration_flag)*(NOR*ONamup)  &
                                       - (respiration_cyano)*(NOR*ONamup)  &
                                       - (recycling_detritus_to_n)*(NOR*ldn_O)  &
                                       - (nitrification)*(NOR*ONnitr)  &
                                       )
          !-----------------
          ! ammonium; containing 
          !-----------------
          deriv(idx_amm_with_totatmos_N) = (                               &
                                       + (respiration_miz_totatmos_N)*((one-ldon_frac))  &
                                       + (respiration_mez_totatmos_N)*((one-ldon_frac))  &
                                       + (respiration_dia_totatmos_N)*((one-ldon_frac))  &
                                       + (respiration_flag_totatmos_N)*((one-ldon_frac))  &
                                       + (respiration_cyano_totatmos_N)*((one-ldon_frac))  &
                                       + recycling_detritus_to_n_totatmos_N  &
                                       + degradation_ldon_totatmos_N   &
                                       - uptake_amm_by_dia*max(0.0_8,min(1.0_8,amm_with_totatmos_N/max(0.00000000001_8,amm)))  &
                                       - uptake_amm_by_flag*max(0.0_8,min(1.0_8,amm_with_totatmos_N/max(0.00000000001_8,amm)))  &
                                       - nitrification*max(0.0_8,min(1.0_8,amm_with_totatmos_N/max(0.00000000001_8,amm)))  &
                                       )
          !-----------------
          ! nitrate; containing 
          !-----------------
          deriv(idx_nit_with_totatmos_N) = (                               &
                                       + nitrification_totatmos_N      &
                                       - uptake_nit_by_dia*max(0.0_8,min(1.0_8,nit_with_totatmos_N/max(0.00000000001_8,nit)))  &
                                       - uptake_nit_by_flat*max(0.0_8,min(1.0_8,nit_with_totatmos_N/max(0.00000000001_8,nit)))  &
                                       - (recycling_detritus_to_n)*(ldn_N)*max(0.0_8,min(1.0_8,nit_with_totatmos_N/max(0.00000000001_8,nit)))  &
                                       )
          !-----------------
          ! diatoms; containing 
          !-----------------
          deriv(idx_dia_with_totatmos_N) = (                               &
                                       + uptake_amm_by_dia_totatmos_N  &
                                       + uptake_nit_by_dia_totatmos_N  &
                                       - grazing_miz_on_dia*max(0.0_8,min(1.0_8,dia_with_totatmos_N/max(0.00000000001_8,dia)))  &
                                       - grazing_mez_on_dia*max(0.0_8,min(1.0_8,dia_with_totatmos_N/max(0.00000000001_8,dia)))  &
                                       - respiration_dia*max(0.0_8,min(1.0_8,dia_with_totatmos_N/max(0.00000000001_8,dia)))  &
                                       - mortality_dia*max(0.0_8,min(1.0_8,dia_with_totatmos_N/max(0.00000000001_8,dia)))  &
                                       )
          !-----------------
          ! flagellates; containing 
          !-----------------
          deriv(idx_flag_with_totatmos_N) = (                               &
                                       + uptake_amm_by_flag_totatmos_N  &
                                       + uptake_nit_by_flat_totatmos_N  &
                                       - grazing_miz_on_flag*max(0.0_8,min(1.0_8,flag_with_totatmos_N/max(0.00000000001_8,flag)))  &
                                       - grazing_mez_on_flag*max(0.0_8,min(1.0_8,flag_with_totatmos_N/max(0.00000000001_8,flag)))  &
                                       - respiration_flag*max(0.0_8,min(1.0_8,flag_with_totatmos_N/max(0.00000000001_8,flag)))  &
                                       - mortality_flag*max(0.0_8,min(1.0_8,flag_with_totatmos_N/max(0.00000000001_8,flag)))  &
                                       )
          !-----------------
          ! cyano bacteria; containing 
          !-----------------
          deriv(idx_cyano_with_totatmos_N) = (                               &
                                       - grazing_miz_on_cyano*max(0.0_8,min(1.0_8,cyano_with_totatmos_N/max(0.00000000001_8,cyano)))  &
                                       - grazing_mez_on_cyano*max(0.0_8,min(1.0_8,cyano_with_totatmos_N/max(0.00000000001_8,cyano)))  &
                                       - respiration_cyano*max(0.0_8,min(1.0_8,cyano_with_totatmos_N/max(0.00000000001_8,cyano)))  &
                                       - mortality_cyano*max(0.0_8,min(1.0_8,cyano_with_totatmos_N/max(0.00000000001_8,cyano)))  &
                                       )
          !-----------------
          ! meso zooplankton; containing 
          !-----------------
          deriv(idx_mez_with_totatmos_N) = (                               &
                                       + grazing_mez_on_dia_totatmos_N  &
                                       + grazing_mez_on_flag_totatmos_N  &
                                       + grazing_mez_on_cyano_totatmos_N  &
                                       + grazing_mez_on_miz_totatmos_N  &
                                       - respiration_mez*max(0.0_8,min(1.0_8,mez_with_totatmos_N/max(0.00000000001_8,mez)))  &
                                       - mortality_mez*max(0.0_8,min(1.0_8,mez_with_totatmos_N/max(0.00000000001_8,mez)))  &
                                       )
          !-----------------
          ! micro zooplakton; containing 
          !-----------------
          deriv(idx_miz_with_totatmos_N) = (                               &
                                       + grazing_miz_on_dia_totatmos_N  &
                                       + grazing_miz_on_flag_totatmos_N  &
                                       + grazing_miz_on_cyano_totatmos_N  &
                                       - respiration_miz*max(0.0_8,min(1.0_8,miz_with_totatmos_N/max(0.00000000001_8,miz)))  &
                                       - mortality_miz*max(0.0_8,min(1.0_8,miz_with_totatmos_N/max(0.00000000001_8,miz)))  &
                                       - grazing_mez_on_miz*max(0.0_8,min(1.0_8,miz_with_totatmos_N/max(0.00000000001_8,miz)))  &
                                       )
          !-----------------
          ! Detritus; containing 
          !-----------------
          deriv(idx_det_with_totatmos_N) = (                               &
                                       + mortality_miz_totatmos_N      &
                                       + mortality_mez_totatmos_N      &
                                       + mortality_dia_totatmos_N      &
                                       + mortality_flag_totatmos_N     &
                                       + mortality_cyano_totatmos_N    &
                                       - recycling_detritus_to_n*max(0.0_8,min(1.0_8,det_with_totatmos_N/max(0.00000000001_8,det)))  &
                                       )
          !-----------------
          ! ldon is the labile DON fraction; containing 
          !-----------------
          deriv(idx_ldon_with_totatmos_N) = (                               &
                                       + (respiration_miz_totatmos_N)*(ldon_frac)  &
                                       + (respiration_mez_totatmos_N)*(ldon_frac)  &
                                       + (respiration_dia_totatmos_N)*(ldon_frac)  &
                                       + (respiration_flag_totatmos_N)*(ldon_frac)  &
                                       + (respiration_cyano_totatmos_N)*(ldon_frac)  &
                                       - degradation_ldon*max(0.0_8,min(1.0_8,ldon_with_totatmos_N/max(0.00000000001_8,ldon)))  &
                                       )
          !-----------------
          ! ammonium; containing 
          !-----------------
          deriv(idx_amm_with_agri_N) = (                               &
                                       + (respiration_miz_agri_N)*((one-ldon_frac))  &
                                       + (respiration_mez_agri_N)*((one-ldon_frac))  &
                                       + (respiration_dia_agri_N)*((one-ldon_frac))  &
                                       + (respiration_flag_agri_N)*((one-ldon_frac))  &
                                       + (respiration_cyano_agri_N)*((one-ldon_frac))  &
                                       + recycling_detritus_to_n_agri_N  &
                                       + degradation_ldon_agri_N       &
                                       - uptake_amm_by_dia*max(0.0_8,min(1.0_8,amm_with_agri_N/max(0.00000000001_8,amm)))  &
                                       - uptake_amm_by_flag*max(0.0_8,min(1.0_8,amm_with_agri_N/max(0.00000000001_8,amm)))  &
                                       - nitrification*max(0.0_8,min(1.0_8,amm_with_agri_N/max(0.00000000001_8,amm)))  &
                                       )
          !-----------------
          ! nitrate; containing 
          !-----------------
          deriv(idx_nit_with_agri_N) = (                               &
                                       + nitrification_agri_N          &
                                       - uptake_nit_by_dia*max(0.0_8,min(1.0_8,nit_with_agri_N/max(0.00000000001_8,nit)))  &
                                       - uptake_nit_by_flat*max(0.0_8,min(1.0_8,nit_with_agri_N/max(0.00000000001_8,nit)))  &
                                       - (recycling_detritus_to_n)*(ldn_N)*max(0.0_8,min(1.0_8,nit_with_agri_N/max(0.00000000001_8,nit)))  &
                                       )
          !-----------------
          ! diatoms; containing 
          !-----------------
          deriv(idx_dia_with_agri_N) = (                               &
                                       + uptake_amm_by_dia_agri_N      &
                                       + uptake_nit_by_dia_agri_N      &
                                       - grazing_miz_on_dia*max(0.0_8,min(1.0_8,dia_with_agri_N/max(0.00000000001_8,dia)))  &
                                       - grazing_mez_on_dia*max(0.0_8,min(1.0_8,dia_with_agri_N/max(0.00000000001_8,dia)))  &
                                       - respiration_dia*max(0.0_8,min(1.0_8,dia_with_agri_N/max(0.00000000001_8,dia)))  &
                                       - mortality_dia*max(0.0_8,min(1.0_8,dia_with_agri_N/max(0.00000000001_8,dia)))  &
                                       )
          !-----------------
          ! flagellates; containing 
          !-----------------
          deriv(idx_flag_with_agri_N) = (                               &
                                       + uptake_amm_by_flag_agri_N     &
                                       + uptake_nit_by_flat_agri_N     &
                                       - grazing_miz_on_flag*max(0.0_8,min(1.0_8,flag_with_agri_N/max(0.00000000001_8,flag)))  &
                                       - grazing_mez_on_flag*max(0.0_8,min(1.0_8,flag_with_agri_N/max(0.00000000001_8,flag)))  &
                                       - respiration_flag*max(0.0_8,min(1.0_8,flag_with_agri_N/max(0.00000000001_8,flag)))  &
                                       - mortality_flag*max(0.0_8,min(1.0_8,flag_with_agri_N/max(0.00000000001_8,flag)))  &
                                       )
          !-----------------
          ! cyano bacteria; containing 
          !-----------------
          deriv(idx_cyano_with_agri_N) = (                               &
                                       - grazing_miz_on_cyano*max(0.0_8,min(1.0_8,cyano_with_agri_N/max(0.00000000001_8,cyano)))  &
                                       - grazing_mez_on_cyano*max(0.0_8,min(1.0_8,cyano_with_agri_N/max(0.00000000001_8,cyano)))  &
                                       - respiration_cyano*max(0.0_8,min(1.0_8,cyano_with_agri_N/max(0.00000000001_8,cyano)))  &
                                       - mortality_cyano*max(0.0_8,min(1.0_8,cyano_with_agri_N/max(0.00000000001_8,cyano)))  &
                                       )
          !-----------------
          ! meso zooplankton; containing 
          !-----------------
          deriv(idx_mez_with_agri_N) = (                               &
                                       + grazing_mez_on_dia_agri_N     &
                                       + grazing_mez_on_flag_agri_N    &
                                       + grazing_mez_on_cyano_agri_N   &
                                       + grazing_mez_on_miz_agri_N     &
                                       - respiration_mez*max(0.0_8,min(1.0_8,mez_with_agri_N/max(0.00000000001_8,mez)))  &
                                       - mortality_mez*max(0.0_8,min(1.0_8,mez_with_agri_N/max(0.00000000001_8,mez)))  &
                                       )
          !-----------------
          ! micro zooplakton; containing 
          !-----------------
          deriv(idx_miz_with_agri_N) = (                               &
                                       + grazing_miz_on_dia_agri_N     &
                                       + grazing_miz_on_flag_agri_N    &
                                       + grazing_miz_on_cyano_agri_N   &
                                       - respiration_miz*max(0.0_8,min(1.0_8,miz_with_agri_N/max(0.00000000001_8,miz)))  &
                                       - mortality_miz*max(0.0_8,min(1.0_8,miz_with_agri_N/max(0.00000000001_8,miz)))  &
                                       - grazing_mez_on_miz*max(0.0_8,min(1.0_8,miz_with_agri_N/max(0.00000000001_8,miz)))  &
                                       )
          !-----------------
          ! Detritus; containing 
          !-----------------
          deriv(idx_det_with_agri_N) = (                               &
                                       + mortality_miz_agri_N          &
                                       + mortality_mez_agri_N          &
                                       + mortality_dia_agri_N          &
                                       + mortality_flag_agri_N         &
                                       + mortality_cyano_agri_N        &
                                       - recycling_detritus_to_n*max(0.0_8,min(1.0_8,det_with_agri_N/max(0.00000000001_8,det)))  &
                                       )
          !-----------------
          ! ldon is the labile DON fraction; containing 
          !-----------------
          deriv(idx_ldon_with_agri_N) = (                               &
                                       + (respiration_miz_agri_N)*(ldon_frac)  &
                                       + (respiration_mez_agri_N)*(ldon_frac)  &
                                       + (respiration_dia_agri_N)*(ldon_frac)  &
                                       + (respiration_flag_agri_N)*(ldon_frac)  &
                                       + (respiration_cyano_agri_N)*(ldon_frac)  &
                                       - degradation_ldon*max(0.0_8,min(1.0_8,ldon_with_agri_N/max(0.00000000001_8,ldon)))  &
                                       )
        
          !--------------------------------------------------------------------!
          ! Integration 
          !--------------------------------------------------------------------!
          if (k==kb) then
            !-----------------
            ! Add benthic dynamics
            !-----------------
            call timer(1,'biobenthic')
            call bio_benthic(ecoa,temp,ben,deriv,benderiv,rfr,hbot,kb)
            call timer(2,'biobenthic')
          
            call timer(1,'biosolver')
            call solvebio(ecoa,deriv,dt,irc,kb,ben,benderiv)
            call timer(2,'biosolver')
          else
            ! pelagic 
           
            call timer(1,'biosolver')
            call solvebio(ecoa,deriv,dt,irc,k)
            call timer(2,'biosolver')
          endif
        
          if (irc /= 0) return
        enddo  ! k 

        ! Transfer results to full arrays --------------------------------------
        
        ! save light in the diagnostics array for later use in other modules
        ! light always gets the first index:
        diag_light(nsrf) = light(1)
        diag_kpar(nsrf)  = kpar(1)
        
        ! First, the surface:
        benthos(nsrf,bidx_nitr           ) = max(ben(bidx_nitr),zero)
        benthos(nsrf,bidx_sili           ) = max(ben(bidx_sili),zero)
        benthos(nsrf,bidx_nitr_with_totatmos_N) = ben(bidx_nitr_with_totatmos_N)
        benthos(nsrf,bidx_nitr_with_agri_N) = ben(bidx_nitr_with_agri_N)
        
        ts(nsrf,idx_amm            ) = max(ecoa(1,idx_amm),zero)
        ts(nsrf,idx_nit            ) = max(ecoa(1,idx_nit),zero)
        ts(nsrf,idx_phos           ) = max(ecoa(1,idx_phos),zero)
        ts(nsrf,idx_sil            ) = max(ecoa(1,idx_sil),zero)
        ts(nsrf,idx_dia            ) = max(ecoa(1,idx_dia),zero)
        ts(nsrf,idx_flag           ) = max(ecoa(1,idx_flag),zero)
        ts(nsrf,idx_cyano          ) = max(ecoa(1,idx_cyano),zero)
        ts(nsrf,idx_mez            ) = max(ecoa(1,idx_mez),zero)
        ts(nsrf,idx_miz            ) = max(ecoa(1,idx_miz),zero)
        ts(nsrf,idx_det            ) = max(ecoa(1,idx_det),zero)
        ts(nsrf,idx_dets           ) = max(ecoa(1,idx_dets),zero)
        ts(nsrf,idx_ldon           ) = max(ecoa(1,idx_ldon),zero)
        ! oxygen can be negative since it represent hydrogen sulphite:
        ts(nsrf,idx_oxy            ) = ecoa(1,idx_oxy)
        ts(nsrf,idx_amm_with_totatmos_N) = ecoa(1,idx_amm_with_totatmos_N)
        ts(nsrf,idx_nit_with_totatmos_N) = ecoa(1,idx_nit_with_totatmos_N)
        ts(nsrf,idx_dia_with_totatmos_N) = ecoa(1,idx_dia_with_totatmos_N)
        ts(nsrf,idx_flag_with_totatmos_N) = ecoa(1,idx_flag_with_totatmos_N)
        ts(nsrf,idx_cyano_with_totatmos_N) = ecoa(1,idx_cyano_with_totatmos_N)
        ts(nsrf,idx_mez_with_totatmos_N) = ecoa(1,idx_mez_with_totatmos_N)
        ts(nsrf,idx_miz_with_totatmos_N) = ecoa(1,idx_miz_with_totatmos_N)
        ts(nsrf,idx_det_with_totatmos_N) = ecoa(1,idx_det_with_totatmos_N)
        ts(nsrf,idx_ldon_with_totatmos_N) = ecoa(1,idx_ldon_with_totatmos_N)
        ts(nsrf,idx_amm_with_agri_N) = ecoa(1,idx_amm_with_agri_N)
        ts(nsrf,idx_nit_with_agri_N) = ecoa(1,idx_nit_with_agri_N)
        ts(nsrf,idx_dia_with_agri_N) = ecoa(1,idx_dia_with_agri_N)
        ts(nsrf,idx_flag_with_agri_N) = ecoa(1,idx_flag_with_agri_N)
        ts(nsrf,idx_cyano_with_agri_N) = ecoa(1,idx_cyano_with_agri_N)
        ts(nsrf,idx_mez_with_agri_N) = ecoa(1,idx_mez_with_agri_N)
        ts(nsrf,idx_miz_with_agri_N) = ecoa(1,idx_miz_with_agri_N)
        ts(nsrf,idx_det_with_agri_N) = ecoa(1,idx_det_with_agri_N)
        ts(nsrf,idx_ldon_with_agri_N) = ecoa(1,idx_ldon_with_agri_N)

        ! Then, the sub-surface:
        if (kb > 1) then
          ts(mi0:mi0+kb-2,idx_amm            ) = max(ecoa(2:kb,idx_amm),zero)
          ts(mi0:mi0+kb-2,idx_nit            ) = max(ecoa(2:kb,idx_nit),zero)
          ts(mi0:mi0+kb-2,idx_phos           ) = max(ecoa(2:kb,idx_phos),zero)
          ts(mi0:mi0+kb-2,idx_sil            ) = max(ecoa(2:kb,idx_sil),zero)
          ts(mi0:mi0+kb-2,idx_dia            ) = max(ecoa(2:kb,idx_dia),zero)
          ts(mi0:mi0+kb-2,idx_flag           ) = max(ecoa(2:kb,idx_flag),zero)
          ts(mi0:mi0+kb-2,idx_cyano          ) = max(ecoa(2:kb,idx_cyano),zero)
          ts(mi0:mi0+kb-2,idx_mez            ) = max(ecoa(2:kb,idx_mez),zero)
          ts(mi0:mi0+kb-2,idx_miz            ) = max(ecoa(2:kb,idx_miz),zero)
          ts(mi0:mi0+kb-2,idx_det            ) = max(ecoa(2:kb,idx_det),zero)
          ts(mi0:mi0+kb-2,idx_dets           ) = max(ecoa(2:kb,idx_dets),zero)
          ts(mi0:mi0+kb-2,idx_ldon           ) = max(ecoa(2:kb,idx_ldon),zero)
          ts(mi0:mi0+kb-2,idx_oxy            ) = ecoa(2:kb,idx_oxy)
          ts(mi0:mi0+kb-2,idx_amm_with_totatmos_N) = ecoa(2:kb,idx_amm_with_totatmos_N)
          ts(mi0:mi0+kb-2,idx_nit_with_totatmos_N) = ecoa(2:kb,idx_nit_with_totatmos_N)
          ts(mi0:mi0+kb-2,idx_dia_with_totatmos_N) = ecoa(2:kb,idx_dia_with_totatmos_N)
          ts(mi0:mi0+kb-2,idx_flag_with_totatmos_N) = ecoa(2:kb,idx_flag_with_totatmos_N)
          ts(mi0:mi0+kb-2,idx_cyano_with_totatmos_N) = ecoa(2:kb,idx_cyano_with_totatmos_N)
          ts(mi0:mi0+kb-2,idx_mez_with_totatmos_N) = ecoa(2:kb,idx_mez_with_totatmos_N)
          ts(mi0:mi0+kb-2,idx_miz_with_totatmos_N) = ecoa(2:kb,idx_miz_with_totatmos_N)
          ts(mi0:mi0+kb-2,idx_det_with_totatmos_N) = ecoa(2:kb,idx_det_with_totatmos_N)
          ts(mi0:mi0+kb-2,idx_ldon_with_totatmos_N) = ecoa(2:kb,idx_ldon_with_totatmos_N)
          ts(mi0:mi0+kb-2,idx_amm_with_agri_N) = ecoa(2:kb,idx_amm_with_agri_N)
          ts(mi0:mi0+kb-2,idx_nit_with_agri_N) = ecoa(2:kb,idx_nit_with_agri_N)
          ts(mi0:mi0+kb-2,idx_dia_with_agri_N) = ecoa(2:kb,idx_dia_with_agri_N)
          ts(mi0:mi0+kb-2,idx_flag_with_agri_N) = ecoa(2:kb,idx_flag_with_agri_N)
          ts(mi0:mi0+kb-2,idx_cyano_with_agri_N) = ecoa(2:kb,idx_cyano_with_agri_N)
          ts(mi0:mi0+kb-2,idx_mez_with_agri_N) = ecoa(2:kb,idx_mez_with_agri_N)
          ts(mi0:mi0+kb-2,idx_miz_with_agri_N) = ecoa(2:kb,idx_miz_with_agri_N)
          ts(mi0:mi0+kb-2,idx_det_with_agri_N) = ecoa(2:kb,idx_det_with_agri_N)
          ts(mi0:mi0+kb-2,idx_ldon_with_agri_N) = ecoa(2:kb,idx_ldon_with_agri_N)
          diag_light(mi0:mi0+kb-2)   = light(2:kb)
          diag_kpar(mi0:mi0+kb-2)    = kpar(2:kb)
        endif
        
        if (warn_bad_solver) then
          if ( ANY (ecoa(1:kb,idx_amm) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer amm although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_nit) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer nit although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_phos) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer phos although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_sil) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer sil although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_dia) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer dia although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_flag) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer flag although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_cyano) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer cyano although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_mez) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer mez although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_miz) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer miz although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_det) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer det although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_dets) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer dets although should be >=0._8'
          if ( ANY (ecoa(1:kb,idx_ldon) < 0) )     &
              write(iu06,*) 'Neg. val. in bio tracer ldon although should be >=0._8'
          if ( ANY (ecoa(1:kb,bidx_nitr) < 0) )     &
              write(iu06,*) 'Neg. val. in ben tracer nitr although should be >=0._8'
          if ( ANY (ecoa(1:kb,bidx_sili) < 0) )     &
              write(iu06,*) 'Neg. val. in ben tracer sili although should be >=0._8'
        endif
        
      endif ! rivmap
    enddo   ! nsrf

  end subroutine bio_dynamics

end module biomod
