module biobenthic

  use constants, only : zero, one, half

  implicit none

  !---------------------------------------------------------
  ! This module contains the Benthic biogeochemical submodel
  ! part of the ecological model used in CMOD at DMI
  !---------------------------------------------------------

  private

  public :: bio_benthic  

contains

  function theta(x) 
    real(8)             :: theta
    real(8), intent(in) :: x
      theta = 0.0_8
    if (x .gt. 0.0_8) then
      theta = 1.0_8
    end if
  end function theta

  subroutine bio_benthic(eco,bottemp,ben,deriv,benderiv,rfr,hbot,kb)
    use bioparam, only :                                              &
                         dn_sed,dn_sed_anox,q10_recs, NOR,            &
                         ONamup, ONnitr, ksdenit, NNdenit,            &
                         fnitr, pfrac, Kspburial,                     &
                         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,                         &
                         ntbio,                                       &
                         bidx_nitr           ,                        &
                         bidx_sili           ,                        &
                         bidx_nitr_with_totatmos_N,                        &
                         bidx_nitr_with_agri_N,                        &
                         ntben 
    implicit none

    real(8),    intent(in)    :: eco(:,0:), ben(0:), rfr, hbot
    real(8),    intent(in)    :: bottemp
    real(8),    intent(inout) :: deriv(0:), benderiv(0:)
    integer,    intent(in)    :: kb

    !-----------------
    ! local variables
    !-----------------
    ! fix parameters
    real(8), parameter :: max_exp     = 0.4_8
    real(8), parameter :: burialcoeff = 5.0e-6_8

    ! quantities, independent of cgt
    real(8)    :: exparg
    integer(4) :: iben
    real(8)    :: burial(0:ntben)


    ! auxiliary parameters
    real(8) :: Kspburial2       !  
    real(8) :: oxy2             !  
    real(8) :: recs             ! auxiliary for calculating mineralization and nitrification 
    real(8) :: denitscal        ! using nitrate (and suphalte) to oxidize benthic organic material 
    real(8) :: ldn_N_sed        ! Under oxic conditions benthic OM is oxidated by oxygen. Under anoxic conditions benthic OM is oxidated by nitrate, if available. Otherwise by sulfate as negative o2 
    real(8) :: ldn_O_sed        ! Under oxic conditions benthic OM is oxidated by oxygen. Under anoxic conditions benthic OM is oxidated by nitrate, if available. Otherwise by sulfate as negative o2 
    real(8) :: one_pburial      ! Phosphate is buried under oxic conditions but no P burial under anoxic conditions 

    ! process rates
    real(8) :: nitrification_amm_sed  ! nitrification
    real(8) :: mineralization_nitr_sed  ! mineralization of bentic nitrogen, mmol N *m^-2*d^-1
    real(8) :: mineralization_sili_sed  ! mineralization of bentic silicate, mmol N *m^-2*d^-1
    real(8) :: nitrification_amm_sed_totatmos_N  ! nitrification; sub-process for totatmos nitrogen
    real(8) :: mineralization_nitr_sed_totatmos_N  ! mineralization of bentic nitrogen, mmol N *m^-2*d^-1; sub-process for totatmos nitrogen
    real(8) :: nitrification_amm_sed_agri_N  ! nitrification; sub-process for agri nitrogen
    real(8) :: mineralization_nitr_sed_agri_N  ! mineralization of bentic nitrogen, mmol N *m^-2*d^-1; sub-process for agri nitrogen

    ! tracers
    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) :: nitr             ! nitrogen in the sediment
    real(8) :: sili             ! silicate in the sediment
    real(8) :: nitr_with_totatmos_N  ! nitrogen in the sediment; containing 
    real(8) :: nitr_with_agri_N  ! nitrogen in the sediment; containing 


    !--------------
    ! calculations
    !--------------
    ! avoid unnecessary array access:
    amm              = eco(kb,idx_amm)
    nit              = eco(kb,idx_nit)
    phos             = eco(kb,idx_phos)
    sil              = eco(kb,idx_sil)
    dia              = eco(kb,idx_dia)
    flag             = eco(kb,idx_flag)
    cyano            = eco(kb,idx_cyano)
    mez              = eco(kb,idx_mez)
    miz              = eco(kb,idx_miz)
    det              = eco(kb,idx_det)
    dets             = eco(kb,idx_dets)
    ldon             = eco(kb,idx_ldon)
    oxy              = eco(kb,idx_oxy)
    amm_with_totatmos_N  = eco(kb,idx_amm_with_totatmos_N)
    nit_with_totatmos_N  = eco(kb,idx_nit_with_totatmos_N)
    dia_with_totatmos_N  = eco(kb,idx_dia_with_totatmos_N)
    flag_with_totatmos_N  = eco(kb,idx_flag_with_totatmos_N)
    cyano_with_totatmos_N  = eco(kb,idx_cyano_with_totatmos_N)
    mez_with_totatmos_N  = eco(kb,idx_mez_with_totatmos_N)
    miz_with_totatmos_N  = eco(kb,idx_miz_with_totatmos_N)
    det_with_totatmos_N  = eco(kb,idx_det_with_totatmos_N)
    ldon_with_totatmos_N  = eco(kb,idx_ldon_with_totatmos_N)
    amm_with_agri_N  = eco(kb,idx_amm_with_agri_N)
    nit_with_agri_N  = eco(kb,idx_nit_with_agri_N)
    dia_with_agri_N  = eco(kb,idx_dia_with_agri_N)
    flag_with_agri_N  = eco(kb,idx_flag_with_agri_N)
    cyano_with_agri_N  = eco(kb,idx_cyano_with_agri_N)
    mez_with_agri_N  = eco(kb,idx_mez_with_agri_N)
    miz_with_agri_N  = eco(kb,idx_miz_with_agri_N)
    det_with_agri_N  = eco(kb,idx_det_with_agri_N)
    ldon_with_agri_N  = eco(kb,idx_ldon_with_agri_N)
    nitr             = ben(bidx_nitr)
    sili             = ben(bidx_sili)
    nitr_with_totatmos_N  = ben(bidx_nitr_with_totatmos_N)
    nitr_with_agri_N  = ben(bidx_nitr_with_agri_N)

    ! Mineralized detritus becomes ammonium and nitrified under oxic conditions
    ! The produced nitrate is assumed to be used in coupled denitrification
    ! Thus no nitrate is produced but N2 which is lost from the system
    ! 
    ! Under oxic conditions benthic OM is oxidated by oxygen.
    ! Under anoxic conditions benthic OM is oxidated by nitrate, if available,
    ! otherwise by sulfate as negative o2
    ! set up some numbers
    Kspburial2      = Kspburial*Kspburial           
    oxy2            = oxy*oxy                       
    recs            = exp(q10_recs*bottemp) * (theta(oxy)*dn_sed + (1.0_8-theta(oxy))*dn_sed_anox)
    denitscal       = (1.0_8-theta(oxy)) * nit/(ksdenit+nit)
    ldn_N_sed       = NNdenit * denitscal           
    ldn_O_sed       = ONamup * (one - denitscal)    
    one_pburial     = 1.0_8 - theta(oxy) * pfrac * oxy2 * (oxy2 + Kspburial2)

    !- permanent burial of benthos when the concentration is very large
    do iben=1,ntben
      exparg = ben(iben)*ben(iben)*burialcoeff
      if (exparg > max_exp) then
        burial(iben) = half*ben(iben)
      else
        burial(iben) = (exp(exparg)-one)*ben(iben)
      endif
    enddo


    !fixme: denitrification in oxic conditions not same as Neumann
    nitrification_amm_sed  = fnitr*recs*nitr*theta(oxy)     ! nitrification
    mineralization_nitr_sed  = recs * nitr                    ! mineralization of bentic nitrogen, mmol N *m^-2*d^-1
    mineralization_sili_sed  = recs * sili                    ! mineralization of bentic silicate, mmol N *m^-2*d^-1
    nitrification_amm_sed_totatmos_N  = nitrification_amm_sed * ((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)) ! nitrification; sub-process for totatmos nitrogen
    mineralization_nitr_sed_totatmos_N  = mineralization_nitr_sed * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,nitr_with_totatmos_N/max(0.00000000001_8,nitr)))+(ldn_N_sed)*(1)*max(0.0_8,min(1.0_8,nit_with_totatmos_N/max(0.00000000001_8,nit)))) / ((1.0_8)*(1)+(ldn_N_sed)*(1)) ! mineralization of bentic nitrogen, mmol N *m^-2*d^-1; sub-process for totatmos nitrogen
    nitrification_amm_sed_agri_N  = nitrification_amm_sed * ((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)) ! nitrification; sub-process for agri nitrogen
    mineralization_nitr_sed_agri_N  = mineralization_nitr_sed * ((1.0_8)*(1)*max(0.0_8,min(1.0_8,nitr_with_agri_N/max(0.00000000001_8,nitr)))+(ldn_N_sed)*(1)*max(0.0_8,min(1.0_8,nit_with_agri_N/max(0.00000000001_8,nit)))) / ((1.0_8)*(1)+(ldn_N_sed)*(1)) ! mineralization of bentic nitrogen, mmol N *m^-2*d^-1; sub-process for agri nitrogen

    !------------
    ! Net fluxes
    !------------
    deriv(idx_amm) = deriv(idx_amm) + ( zero &
                          + mineralization_nitr_sed &
                          - nitrification_amm_sed &
                          )/hbot
    deriv(idx_nit) = deriv(idx_nit) + ( zero &
                          - (mineralization_nitr_sed)*(ldn_N_sed) &
                          )/hbot
    deriv(idx_phos) = deriv(idx_phos) + ( zero &
                          + (mineralization_nitr_sed)*(rfr*one_pburial) &
                          )/hbot
    deriv(idx_sil) = deriv(idx_sil) + ( zero &
                          )/hbot
    deriv(idx_dia) = deriv(idx_dia) + ( zero &
                          )/hbot
    deriv(idx_flag) = deriv(idx_flag) + ( zero &
                          )/hbot
    deriv(idx_cyano) = deriv(idx_cyano) + ( zero &
                          )/hbot
    deriv(idx_mez) = deriv(idx_mez) + ( zero &
                          )/hbot
    deriv(idx_miz) = deriv(idx_miz) + ( zero &
                          )/hbot
    deriv(idx_det) = deriv(idx_det) + ( zero &
                          )/hbot
    deriv(idx_dets) = deriv(idx_dets) + ( zero &
                          )/hbot
    deriv(idx_ldon) = deriv(idx_ldon) + ( zero &
                          )/hbot
    deriv(idx_oxy) = deriv(idx_oxy) + ( zero &
                          - (nitrification_amm_sed)*(ONnitr*NOR) &
                          - (mineralization_nitr_sed)*(NOR*ldn_O_sed) &
                          )/hbot
    deriv(idx_amm_with_totatmos_N) = deriv(idx_amm_with_totatmos_N) + ( zero &
                          + mineralization_nitr_sed_totatmos_N &
                          - nitrification_amm_sed*max(0.0_8,min(1.0_8,amm_with_totatmos_N/max(0.00000000001_8,amm))) &
                          )/hbot
    deriv(idx_nit_with_totatmos_N) = deriv(idx_nit_with_totatmos_N) + ( zero &
                          - (mineralization_nitr_sed)*(ldn_N_sed)*max(0.0_8,min(1.0_8,nit_with_totatmos_N/max(0.00000000001_8,nit))) &
                          )/hbot
    deriv(idx_dia_with_totatmos_N) = deriv(idx_dia_with_totatmos_N) + ( zero &
                          )/hbot
    deriv(idx_flag_with_totatmos_N) = deriv(idx_flag_with_totatmos_N) + ( zero &
                          )/hbot
    deriv(idx_cyano_with_totatmos_N) = deriv(idx_cyano_with_totatmos_N) + ( zero &
                          )/hbot
    deriv(idx_mez_with_totatmos_N) = deriv(idx_mez_with_totatmos_N) + ( zero &
                          )/hbot
    deriv(idx_miz_with_totatmos_N) = deriv(idx_miz_with_totatmos_N) + ( zero &
                          )/hbot
    deriv(idx_det_with_totatmos_N) = deriv(idx_det_with_totatmos_N) + ( zero &
                          )/hbot
    deriv(idx_ldon_with_totatmos_N) = deriv(idx_ldon_with_totatmos_N) + ( zero &
                          )/hbot
    deriv(idx_amm_with_agri_N) = deriv(idx_amm_with_agri_N) + ( zero &
                          + mineralization_nitr_sed_agri_N &
                          - nitrification_amm_sed*max(0.0_8,min(1.0_8,amm_with_agri_N/max(0.00000000001_8,amm))) &
                          )/hbot
    deriv(idx_nit_with_agri_N) = deriv(idx_nit_with_agri_N) + ( zero &
                          - (mineralization_nitr_sed)*(ldn_N_sed)*max(0.0_8,min(1.0_8,nit_with_agri_N/max(0.00000000001_8,nit))) &
                          )/hbot
    deriv(idx_dia_with_agri_N) = deriv(idx_dia_with_agri_N) + ( zero &
                          )/hbot
    deriv(idx_flag_with_agri_N) = deriv(idx_flag_with_agri_N) + ( zero &
                          )/hbot
    deriv(idx_cyano_with_agri_N) = deriv(idx_cyano_with_agri_N) + ( zero &
                          )/hbot
    deriv(idx_mez_with_agri_N) = deriv(idx_mez_with_agri_N) + ( zero &
                          )/hbot
    deriv(idx_miz_with_agri_N) = deriv(idx_miz_with_agri_N) + ( zero &
                          )/hbot
    deriv(idx_det_with_agri_N) = deriv(idx_det_with_agri_N) + ( zero &
                          )/hbot
    deriv(idx_ldon_with_agri_N) = deriv(idx_ldon_with_agri_N) + ( zero &
                          )/hbot

    benderiv(bidx_nitr) = &
                                - mineralization_nitr_sed      &
                                - burial(bidx_nitr)
    benderiv(bidx_sili) = &
                                - mineralization_sili_sed      &
                                - burial(bidx_sili)
    benderiv(bidx_nitr_with_totatmos_N) = &
                                - mineralization_nitr_sed*max(0.0_8,min(1.0_8,nitr_with_totatmos_N/max(0.00000000001_8,nitr))) &
                                - burial(bidx_nitr_with_totatmos_N)
    benderiv(bidx_nitr_with_agri_N) = &
                                - mineralization_nitr_sed*max(0.0_8,min(1.0_8,nitr_with_agri_N/max(0.00000000001_8,nitr))) &
                                - burial(bidx_nitr_with_agri_N)
  end subroutine bio_benthic

end module biobenthic
