#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: bio_mab --- IOW model ERGOM, MaBenE version \label{sec:bio-mab}
!
! !INTERFACE:
   module bio_mab
!
! !DESCRIPTION:
! Special version of the biogeochemical model by \cite{Neumannetal2002}
! for the MaBenE project. 
! 
!
! !USES:
!  default: all is private.
   use bio_var
   use observations, only : aa,g2
   private
!
! !PUBLIC MEMBER FUNCTIONS:
   public init_bio_mab, init_var_mab, var_info_mab, &
          surface_fluxes_mab,light_mab, do_bio_mab, end_bio_mab
!
! !PRIVATE DATA MEMBERS:
!
! !REVISION HISTORY:!
!  Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
!  from a namelist
   REALTYPE                  :: p1_initial=4.5
   REALTYPE                  :: p2_initial=4.5
   REALTYPE                  :: p3_initial=4.5
   REALTYPE                  :: zo_initial=4.5
   REALTYPE                  :: de_initial=4.5
   REALTYPE                  :: am_initial=4.5
   REALTYPE                  :: ni_initial=4.5
   REALTYPE                  :: po_initial=4.5
   REALTYPE                  :: o2_initial=4.5
   REALTYPE                  :: se_initial=4.5
   REALTYPE                  :: sfl_po=0.0015
   REALTYPE                  :: sfl_am=0.07
   REALTYPE                  :: sfl_ni=0.09
   logical                   :: fluff=.false.
   REALTYPE                  :: fl_initial=0.0
   REALTYPE, public          :: p10=0.0225
   REALTYPE, public          :: p20=0.0225
   REALTYPE, public          :: p30=0.0225
   REALTYPE                  :: zo0=0.0225
   REALTYPE                  :: w_p1=-1.157407e-05
   REALTYPE                  :: w_p2=-5.787037e-05
   REALTYPE                  :: w_p3=-5.787037e-05
   REALTYPE                  :: w_de=-3.
   REALTYPE                  :: w_se=-3.
   REALTYPE, public          :: kc=0.03
   REALTYPE                  :: i_min=25.
   REALTYPE                  :: r1max=1.
   REALTYPE                  :: r2max=1.
   REALTYPE                  :: r3max=1.
   REALTYPE                  :: alpha1=0.3
   REALTYPE                  :: alpha2=0.15
   REALTYPE                  :: alpha3=0.5
   REALTYPE                  :: lpa=0.01
   REALTYPE                  :: lpd=0.02
   REALTYPE                  :: Tf=10.
   REALTYPE                  :: Tbg=16.
   REALTYPE                  :: beta_bg=1.
   REALTYPE                  :: g1max=0.5
   REALTYPE                  :: g2max=0.5
   REALTYPE                  :: g3max=0.25
   REALTYPE                  :: lza=0.3
   REALTYPE                  :: lzd=0.6
   REALTYPE, public          :: iv=1.2
   REALTYPE                  :: topt=20.
   REALTYPE                  :: lan=0.1
   REALTYPE                  :: oan=0.01
   REALTYPE                  :: beta_an=0.11
   REALTYPE                  :: lda=0.003
   REALTYPE                  :: Tda=13.
   REALTYPE                  :: beta_da=20.
   REALTYPE                  :: lds=4.05e-5
   REALTYPE                  :: lsa=1.16e-8
   REALTYPE                  :: bsa=0.15
   REALTYPE                  :: ph1=0.15
   REALTYPE                  :: ph2=0.1
   REALTYPE                  :: pvel=5.
   REALTYPE                  :: sr=0.0625
   REALTYPE                  :: s1=5.3
   REALTYPE                  :: s2=6.625
   REALTYPE                  :: s3=8.125
   REALTYPE                  :: s4=0.666666666
   REALTYPE                  :: a0=31.25
   REALTYPE                  :: a1=14.603
   REALTYPE                  :: a2=0.4025
   integer                   :: out_unit
   integer, parameter        :: p1=1,p2=2,p3=3,zo=4,de=5,     &
                                am=6,ni=7,po=8,o2=9,se=10,fl=11
   REALTYPE, allocatable     :: ppi(:)
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the bio module
!
! !INTERFACE:
   subroutine init_bio_mab(namlst,fname,unit)
!
! !DESCRIPTION:
!  Here, the bio namelist {\tt bio\_mab.nml} is read and
!  various variables (rates and settling velocities)
!  are transformed into SI units.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer,          intent(in)   :: namlst
   character(len=*), intent(in)   :: fname
   integer,          intent(in)   :: unit
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
   namelist /bio_mab_nml/ numc,p1_initial,p2_initial,p3_initial,zo_initial,  &
                      de_initial,am_initial,ni_initial,po_initial,           &
                      o2_initial,se_initial,sfl_po,sfl_am,sfl_ni,            &
                      surface_flux_method,   &
                      fluff,fl_initial,p10,p20,p30,zo0,                      &
                      w_p1,w_p2,w_p3,                                        &
                      w_de,w_se,kc,i_min,r1max,r2max,r3max,alpha1,alpha2,         &
                      alpha3,lpa,lpd,tf,tbg,beta_bg,g1max,g2max,             &
                      g3max,lza,lzd,iv,topt,lan,oan,beta_an,lda,             &
                      tda,beta_da,lds,lsa,bsa,ph1,ph2,pvel,sr,               &
                      s1,s2,s3,s4,a0,a1,a2
!EOP
!-----------------------------------------------------------------------
!BOC
   LEVEL2 'init_bio_mab'

   numc=10
   open(namlst,file=fname,action='read',status='old',err=98)
   read(namlst,nml=bio_mab_nml,err=99)
   close(namlst)

   n_surface_fluxes=3

   numcc=numc
   if (fluff) numc=numc+1

!  Conversion from day to second
   w_p1   = w_p1    /secs_pr_day
   w_p2   = w_p2    /secs_pr_day
   w_p3   = w_p3    /secs_pr_day
   w_de   = w_de    /secs_pr_day
   w_se   = w_se    /secs_pr_day
   r1max  = r1max   /secs_pr_day
   r2max  = r2max   /secs_pr_day
   r3max  = r3max   /secs_pr_day
   lpa    = lpa     /secs_pr_day
   lpd    = lpd     /secs_pr_day
   g1max  = g1max   /secs_pr_day
   g2max  = g2max   /secs_pr_day
   g3max  = g3max   /secs_pr_day
   lza    = lza     /secs_pr_day
   lzd    = lzd     /secs_pr_day
   lan    = lan     /secs_pr_day
   lda    = lda     /secs_pr_day
   lds    = lds     /secs_pr_day
   lsa    = lsa     /secs_pr_day
   pvel   = pvel    /secs_pr_day

   out_unit=unit

   LEVEL3 'IOW bio module initialised ...'

   return

98 LEVEL2 'I could not open bio_mab.nml'
   LEVEL2 'If thats not what you want you have to supply bio_mab.nml'
   LEVEL2 'See the bio example on www.gotm.net for a working bio_mab.nml'
   return
99 FATAL 'I could not read bio_mab.nml'
   stop 'init_bio_mab'
   end subroutine init_bio_mab
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the concentration variables
!
! !INTERFACE:
   subroutine init_var_mab(nlev)
!
! !DESCRIPTION:
!  Here, the the initial conditions are set and the settling velocities are
!  transferred to all vertical levels. All concentrations except oxygen 
!  are declared
!  as non-negative variables, and it is defined which variables would be
!  taken up by benthic filter feeders.
!  Furthermore, the primary production {\tt ppi} is allocated.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer, intent(in)                 :: nlev
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding

! !LOCAL VARIABLES:
  integer                    :: i,rc
!EOP
!-----------------------------------------------------------------------
!BOC
   do i=1,nlev
      cc(p1,i)=p1_initial
      cc(p2,i)=p2_initial
      cc(p3,i)=p3_initial
      cc(zo,i)=zo_initial
      cc(de,i)=de_initial
      cc(am,i)=am_initial
      cc(ni,i)=ni_initial
      cc(po,i)=po_initial
      cc(o2,i)=o2_initial
      cc(se,i)=se_initial
      if (fluff) then 
         if (i .eq. 1) then
            cc(fl,i)=fl_initial+1.e-10
         else
            cc(fl,i)=1.e-10
         end if
      end if
   end do

   do i=0,nlev
      ws(p1,i) = w_p1
      ws(p2,i) = w_p2
      ws(p3,i) = w_p3
      ws(zo,i) = _ZERO_
      ws(de,i) = w_de
      ws(am,i) = _ZERO_
      ws(ni,i) = _ZERO_
      ws(po,i) = _ZERO_
      ws(o2,i) = _ZERO_
      ws(se,i) = w_se
   end do

   sfl = _ZERO_

   posconc(p1) = 1
   posconc(p2) = 1
   posconc(p3) = 1
   posconc(zo) = 1
   posconc(de) = 1
   posconc(am) = 1
   posconc(ni) = 1
   posconc(po) = 1
   posconc(o2) = 0

#if 0
   mussels_inhale(p1) = .true.
   mussels_inhale(p2) = .true.
   mussels_inhale(p3) = .true.
   mussels_inhale(zo) = .true.
   mussels_inhale(de) = .true.
   mussels_inhale(am) = .false.
   mussels_inhale(ni) = .false.
   mussels_inhale(po) = .false.
   mussels_inhale(o2) = .false.
#endif

   allocate(ppi(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_var_mab(): Error allocating ppi)'

   select case (surface_flux_method)
      case (-1)! absolutely nothing
      case (0) ! constant

         sfl(po)=-sfl_po /secs_pr_day
         sfl(am)=-sfl_am /secs_pr_day
         sfl(ni)=-sfl_ni /secs_pr_day

      case (2) ! from file via sfl_read

      case default
   end select

   LEVEL3 'IOW variables initialised ...'

   return
   end subroutine init_var_mab
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Providing info on variables
!
! !INTERFACE:
   subroutine var_info_mab()
!
! !DESCRIPTION:
!  This subroutine provides information about the variable names as they
!  will be used when storing data in NetCDF files.
!
! !USES:
   IMPLICIT NONE
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
!EOP
!-----------------------------------------------------------------------
!BOC
   var_names(1) = 'dia'
   var_units(1) = 'mmol n/m**3'
   var_long(1)  = 'diatoms'

   var_names(2) = 'fla'
   var_units(2) = 'mmol n/m**3'
   var_long(2)  = 'flagellates'

   var_names(3) = 'cya'
   var_units(3) = 'mmol n/m**3'
   var_long(3)  = 'cyanobacteria'

   var_names(4) = 'zoo'
   var_units(4) = 'mmol n/m**3'
   var_long(4)  = 'zooplankton'

   var_names(5) = 'det'
   var_units(5) = 'mmol n/m**3'
   var_long(5)  = 'detritus'

   var_names(6) = 'amm'
   var_units(6) = 'mmol n/m**3'
   var_long(6)  = 'ammonium'

   var_names(7) = 'nit'
   var_units(7) = 'mmol n/m**3'
   var_long(7)  = 'nitrate'

   var_names(8) = 'pho'
   var_units(8) = 'mmol p/m**3'
   var_long(8)  = 'phosphate'

   var_names(9) = 'oxy'
   var_units(9) = 'mmol o2/m**3'
   var_long(9)  = 'oxygen'   

   var_names(10) = 'slt'
   var_units(10) = 'g/m**3'
   var_long(10)  = 'silt'   

   if (fluff) then
      var_names(11) = 'flf'
      var_units(11) = 'mmol n/m**2'
      var_long(11)  = 'fluff'
   end if

   return
   end subroutine var_info_mab
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Step function
!
! !INTERFACE:
   REALTYPE function th(x,w,min,max)
!
! !DESCRIPTION:
! Instead of the
! heavyside switches used by \cite{Neumannetal2002}, we apply here a smoothed
! {\it tangens hyperbolicus} transition with prescribed width $x_w$:
! \begin{equation}\label{theta}
! \theta (x,x_w,y_{\min},y_{\max})= y_{\min}+(y_{\max}-y_{\min})
! \frac12\left(1-\tanh \left(\frac{x}{x_w}   \right)      \right).
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: x,w,min,max
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
!EOP
!-----------------------------------------------------------------------
!BOC
   if (w .gt. 1.e-10) then
      th=min+(max-min)*0.5*(1.+tanh(x/w))
   else
      if (x .gt. _ZERO_) then
         th=_ONE_
      else
         th=_ZERO_
      end if    
   end if
   return
   end function th
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Saturation function squared
!
! !INTERFACE:
   REALTYPE function yy(a,x)
!
! !DESCRIPTION:
! This is a squared Michaelis-Menten type of limiter:
! \begin{equation}\label{Y}
! Y(x_w,x) = \frac{x^2}{x_w^2+x^2}.
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: a,x
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
!EOP
!-----------------------------------------------------------------------
!BOC
   yy=x**2/(a**2+x**2)
   return
   end function yy
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Ivlev formulation for zooplankton grazing on phytoplankton
!
! !INTERFACE:
   REALTYPE function fpz(g,t,topt,psum)
!
! !DESCRIPTION:
! The Ivlev formulation for zooplankton grazing on the three phytoplankton 
! classes $c_1$, $c_2$, and $c_3$ is given here as a function:
! \begin{equation}\label{neu_di4}
! d_{i,4}=g_i^{\max}\left(1+\frac{T^2}{T_{opt}^2}\exp
! \left(1-\frac{2T}{T_{opt}} \right)\right)
! \left( 1-\exp\left(-I_v^2 \left( \sum_{
! j=1}^3c_j \right)^2\right)  \right)
! \frac{c_i}{\sum_{j=1}^3c_j}\left( c_4+c_4^{\min} \right)
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   REALTYPE, intent(in)                :: g,t,topt,psum
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
!EOP
!-----------------------------------------------------------------------
!BOC
   fpz=g*(1.+t**2/topt**2*exp(1.-2.*t/topt))*               &
        (1.-exp(-iv**2*psum**2))
   return
   end function fpz
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Surface fluxes for the IOW model
!
! !INTERFACE:
   subroutine surface_fluxes_mab(nlev,t,s)
!
! !DESCRIPTION:
! Here, those surface fluxes which have been read from a file are transformed
! to SI units, and the surface oxygen flux is calculated by means of the 
! following formula:
! \begin{equation}\label{o2flux}
! F^s_9 = p_{vel} \left(O_{sat}-c_9 \right)
! \end{equation}
! with the \cite{Weiss1970} formula
! \begin{equation}\label{osat}
! O_{sat}= \exp\left[a_1 +a_2\frac{100}{T}+a_3\ln\left(\frac{T}{100}\right)
! +a_4\frac{T}{100}+S \left\{b_1+b_2\frac{T}{100}
! +b_3\left(\frac{T}{100}\right)^2   \right\}\right],
! \end{equation}
!
! where $T$ is the temperature in Kelvin and the empirical constants are
! given in table \ref{table_weiss}.
! 
! \begin{table}[h]
! \begin{center}
! \begin{tabular}{|l|l|l|l|l|l|l|}
! \hline
! $a_1$ & $a_2$ & $a_3$ & $a_4$ & $b_1$ & $b_2$ & $b_3$  \\ \hline
! -173.4292 &
! 249.6339 &
! 143.3483 &
! -21.8492 &
! -0.033096 &
! 0.014259 &
! -0.001700 \\ \hline
! \end{tabular}
! \caption{Constants for the oxygen saturation formula by \cite{Weiss1970},
! see equation (\ref{osat_weiss}).}
! \label{table_weiss}
! \end{center}
! \end{table}
! 
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
  integer                              :: nlev
  REALTYPE, intent(in)                 :: t,s

!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  REALTYPE                 :: osat_weiss,tk
  REALTYPE                 :: a1=-173.4292
  REALTYPE                 :: a2=249.6339
  REALTYPE                 :: a3=143.3483
  REALTYPE                 :: a4=-21.8492
  REALTYPE                 :: b1=-0.033096
  REALTYPE                 :: b2=0.014259
  REALTYPE                 :: b3=-0.001700
  REALTYPE                 :: kelvin=273.16
  REALTYPE                 :: mol_per_liter=44.661

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

   select case (surface_flux_method)
      case (-1)! absolutely nothing
      case (0) ! constant
      case (2) ! from file via sfl_read
         sfl(ni) =-sfl_read(1)/secs_pr_day
         sfl(am) =-sfl_read(2)/secs_pr_day
         sfl(po) =-sfl_read(3)/secs_pr_day
      case (3) ! sfl array filled externally - for 3D models
      case default
   end select

! surface oxygen flux
      tk=t+kelvin
      osat_weiss=a1+a2*100./tk+a3*log(0.01*tk)+a4*0.01*tk
      osat_weiss=osat_weiss+s*(b1+b2*0.01*tk+b3*(0.01*tk)**2)
      osat_weiss=exp(osat_weiss)*mol_per_liter

   sfl(o2) = pvel*(osat_weiss-cc(o2,nlev))
   return
   end subroutine surface_fluxes_mab
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Light properties for the IOW model
!
! !INTERFACE:
   subroutine light_mab(nlev,bioshade_feedback)
!
! !DESCRIPTION:
! Here, the photosynthetically available radiation is calculated
! by simply assuming that the short wave part of the total
! radiation is available for photosynthesis. The user should make
! sure that this is consistent with the light class given in the
! {\tt extinct} namelist of the {\tt obs.nml} file.
! The self-shading effect is also calculated in this subroutine,
! which may be used to consider the effect of bio-turbidity also
! in the temperature equation (if {\tt bioshade\_feedback} is set
! to true in {\tt bio.nml}).
! For details, see section \ref{sec:do-bio}.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer, intent(in)                 :: nlev
   logical, intent(in)                 :: bioshade_feedback
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
   integer                   :: i
   REALTYPE                  :: zz,add
!EOP
!-----------------------------------------------------------------------
!BOC
   zz = _ZERO_
   add = _ZERO_
   do i=nlev,1,-1
      add=add+0.5*h(i)*(cc(de,i)+cc(p1,i)+cc(p2,i)+cc(p3,i)+p10+p20+p30)
      zz=zz+0.5*h(i)
      par(i)=rad(nlev)*(1.-aa)*exp(-zz/g2)*exp(-kc*add)
      add=add+0.5*h(i)*(cc(de,i)+cc(p1,i)+cc(p2,i)+cc(p3,i)+p10+p20+p30)
      zz=zz+0.5*h(i)
      if (bioshade_feedback) bioshade_(i)=exp(-kc*add)
   end do

   return
   end subroutine light_mab
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Right hand sides of the IOW geobiochemical model\label{sec:bio-mab-details}
!
! !INTERFACE:
   subroutine do_bio_mab(first,numc,nlev,cc,pp,dd)
!
! !DESCRIPTION:
! The right hand sides of the \cite{Neumannetal2002} biogeochemical model are 
! coded in this soubroutine.
! First of all, based on (\ref{theta}) and (\ref{Y}), 
! we construct limiters for chemical
! reactions which depend on the availability of oxygen ($c_9$) and nitrate
! ($c_7$) and have to add up to unity:
! \begin{equation}\label{limits}
! \begin{array}{rcl}
! l^+_+ &=& \theta(c_9,c_9^t,0,1)Y(c_7^t,c_7), \\ \\
! l^-_+ &=& \theta(-c_9,c_9^t,0,1)Y(c_7^t,c_7), \\ \\
! l^-_- &=& \theta(-c_9,c_9^t,0,1)(1-Y(c_7^t,c_7)), \\ \\
! L^+_+ &=& \frac{l^+_+}{l^+_+ + l^-_+ + l^-_-}, \\ \\
! L^-_+ &=& \frac{l^-_+}{l^+_+ + l^-_+ + l^-_-}, \\ \\
! L^-_- &=& \frac{l^-_-}{l^+_+ + l^-_+ + l^-_-}. \\ \\
! \end{array}
! \end{equation}
! 
! Mortality of the three phytoplankton classes $c_i$, $i=1,\dots,3$:
! \begin{equation}\label{neu_di5}
! d_{i,5}=l_{PD} c_i
! \end{equation}
! 
! Respiration of the three phytoplankton classes $c_i$, $i=1,\dots,3$
! into ammonium:
! \begin{equation}\label{neu_di6}
! d_{i,6}=l_{PA} c_i
! \end{equation}
! 
! Zooplankton mortality:
! \begin{equation}\label{neu_d45}
! d_{4,5}=l_{ZD}(c_4+c_4^{\min})c_4
! \end{equation}
! 
! Zooplankton exudation into ammonium:
! \begin{equation}\label{neu_d46}
! d_{4,6}=l_{ZA}(c_4+c_4^{\min})c_4
! \end{equation}
! 
! Detritus mineralisation:
! \begin{equation}\label{neu_d56}
! d_{5,6}=L_{DA}c_5
! \end{equation}
! with
! \begin{equation}\label{LDA}
! L_{DA} = l_{DA} \left(1+\beta_{DA}Y(T_{DA},T)\right).
! \end{equation}
! 
! Ammonium uptake by phytoplankta $c_i$, $i=1,2$:
! \begin{equation}\label{neu_d6i}
! d_{6,i}=R_i\frac{c_6}{c_6+c_7}\left(c_i+c_i^{\min} \right)
! \end{equation}
! with the growth rate for diatoms,
! \begin{equation}\label{r1}
! R_1=r_1^{\max} \min\left\{ 
! Y(\alpha_1,c_6+c_7), Y(s_R\alpha_1,c_8), PPI  
! \right\}
! \end{equation}
! and the growth rate for flagellates,
! \begin{equation}\label{r2}
! R_2=r_2^{\max}\left(1+Y\left(T_f,T \right)\right)\min\left\{
! Y(\alpha_2,c_6+c_7),Y(s_R\alpha_2,c_8),PPI 
! \right\}.
! \end{equation}
! Here, 
! \begin{equation}\label{ppi}
! PPI=\frac{I_{PAR}}{I_{opt}}\exp\left(1-\frac{I_{PAR}}{I_{opt}}  \right)
! \end{equation}
! with
! \begin{equation}\label{iopt}
! I_{opt}=\max\left\{\frac{I_0}{4},I_{\min}   \right\}
! \end{equation}
! and $I_{PAR}$ from (\ref{light}).
! 
! Nitrification of ammonium to nitrate:
! \begin{equation}\label{neu_d67}
! d_{6,7}=L_{AN}c_6
! \end{equation}
! with
! \begin{equation}\label{LAN}
! L_{AN}=l_{AN}\theta(c_9,0,0,1)\frac{c_9}{O_{AN}+c_9}\exp\left(\beta_{AN}T\right).
! \end{equation}
! 
! Nitrate uptake by phytoplankta $c_i$, $i=1,2$:
! \begin{equation}\label{neu_d7i}
! d_{7,i}=R_i\frac{c_7}{c_6+c_7}\left(c_i+c_i^{\min} \right).
! \end{equation}
! 
! Settling of detritus into sediment:
! \begin{equation}\label{neu_d510}
! d_{5,10}=l_{DS} \frac{c_5}{h_1}\delta_{k,1}
! \end{equation}
! 
! Mineralisation of sediment into ammonium:
! \begin{equation}\label{neu_d106}
! d_{10,6}=L_{SA} c_{10}
! \end{equation} 
! with
! \begin{equation}\label{LSA}
! L_{SA}=l_{SA} \exp\left(\beta_{SA}T \right) \theta(c_9,c_9^t,0.2,1)
! \end{equation}
! 
! From the above sink terms, respective source terms are calculated by means of (\ref{eq:am:symmetry}),
! except for settling of detritus into sediment and mineralisation of sediment into
! ammonium, for which we have:
! \begin{equation}\label{neu_p105}
! p_{10,5}=h_1 d_{5,10}, \quad p_{6,10}=\frac{d_{10,6}}{h_1}.
! \end{equation}
! 
! Denitrification in water column:
! \begin{equation}\label{neu_d77}
! d_{7,7}=s_1 \left(L_{DA} c_5 +L_{SA} \frac{c_{10}}{h_1}\delta_{k,1} \right)L^-_+.
! \end{equation}
! 
! Denitrification in sediment:
! \begin{equation}\label{neu_d1010}
! d_{10,10}=\theta(c_9,c_9^t,0,1) L_{SA} c_{10}
! \end{equation}
! 
! Phosphorus uptake by the three phytoplankton classes $c_i$, $i=1,\dots,3$:
! \begin{equation}\label{neu_d88}
! d_{8,8}=s_R \left(\sum_{j=1}^3 R_j \left(c_j+c_j^{\min}   \right)     \right). 
! \end{equation}
! 
! Nitrogen fixation:
! \begin{equation}\label{neu_p33}
! p_{3,3}=R_3\left(c_3+c_3^{\min}\right)
! \end{equation}
! with
! \begin{equation}\label{r3}
! R_3=r_3^{\max}\frac{1}{1+\exp\left(\beta_{bg}\left(T_{bg}-T  \right)   \right)}\min\left\{Y\left(s_R\alpha_3,c_8\right),PPI   \right\}
! \end{equation}
! 
! Respiration of the three phytoplankton classes $c_i$, $i=1,\dots,3$
! into phosphorus:
! \begin{equation}\label{neu_p8i}
! p_{8,i}=s_R d_{i,6}.
! \end{equation}
! 
! Zooplankton exudation into phosphorus:
! \begin{equation}\label{neu_p84}
! p_{8,4}=s_R d_{4,6}.
! \end{equation}
! 
! Oxygen production due to ammonium uptake by phytoplankton classes $c_i$, $i=1,2$and nitrification of ammonium into nitrate:
! \begin{equation}\label{neu_p96}
! p_{9,6}= s_2 \left(d_{6,1}+d_{6,2} \right)-s_4 d_{6,7}.
! \end{equation}
! 
! Oxygen production due to nitrate uptake by phytoplankton classes $c_i$, $i=1,2$:
! \begin{equation}\label{neu_p97}
! p_{9,7}= s_3 \left(d_{7,1}+d_{7,2} \right).
! \end{equation}
! 
! Oxygen production due to nitrogen fixation by blue-greens:
! \begin{equation}\label{neu_p99}
! p_{9,9}=s_2 p{3,3}
! \end{equation}
! 
! Oxygen demand due to
! respiration of the three phytoplankton classes $c_i$, $i=1,\dots,3$:
! \begin{equation}\label{neu_p9i}
! p_{9,i}=-s_2 d_{i,6}.
! \end{equation}
! 
! Oxygen demand of zooplankton exudation:
! \begin{equation}\label{neu_p94}
! p_{9,4}=-s_2 d_{4,6}.
! \end{equation}
! 
! Oxygen demand of mineralisation of detritus into ammonium:
! \begin{equation}\label{neu_p95}
! p_{9,5}=-s_2\left(L_+^++L_-^- \right) d_{5,6}.
! \end{equation}
! 
! Oxygen demand of mineralisation of sediment into ammonium:
! \begin{equation}\label{neu_p910}
! p_{9,10}=-\left(s_4+s_2\left(L_+^++L_-^- \right)\right) \frac{d_{10,6}}{h_1}\delta_{k,1}.
! \end{equation}
! 
! Phosphate release due to mineralisation of sediment into ammonium:
! \begin{equation}\label{neu_p88}
! p_{8,8}=(1-p_1\theta\left(c_9,c_9^t,0,1\right)Y(p_2,c_9))\frac{d_{10,6}}{h_1}\delta_{k,1}.
! \end{equation}
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   logical, intent(in)                 :: first
   integer, intent(in)                 :: numc,nlev
   REALTYPE, intent(in)                :: cc(1:numc,0:nlev)
!
! !OUTPUT PARAMETERS:
   REALTYPE, intent(out)               :: pp(1:numc,1:numc,0:nlev)
   REALTYPE, intent(out)               :: dd(1:numc,1:numc,0:nlev)
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard, Karsten Bolding
!
! !LOCAL VARIABLES:
  REALTYPE, save             :: iopt
  REALTYPE                   :: rat(0:nlev,0:nlev)
  REALTYPE                   :: psum,llda,llan,llsa,r1,r2,r3
  REALTYPE                   :: wo=30.,wn=0.1,dot2=0.2
  REALTYPE                   :: thopnp,thomnp,thomnm,thsum
  integer                    :: i,j,ci
!EOP
!-----------------------------------------------------------------------
!BOC
!KBK - is it necessary to initialise every time - expensive in a 3D model
   pp = _ZERO_
   dd = _ZERO_

   if (first) then
      iopt=max(0.25*I_0,I_min)
      do ci=1,nlev
         ppi(ci)=par(ci)/iopt*exp(1.-par(ci)/iopt)
      end do
   end if

   rat=1.         ! fixed (in time  space) ratio between sink and source
   rat(de,fl)=h(1)
   rat(fl,am)=1./h(1)

   do ci=1,nlev

      thopnp=th( cc(o2,ci),wo,_ZERO_,_ONE_)*yy(wn,cc(ni,ci))
      thomnp=th(-cc(o2,ci),wo,_ZERO_,_ONE_)*yy(wn,cc(ni,ci))
      thomnm=th(-cc(o2,ci),wo,_ZERO_,_ONE_)*(1.-yy(wn,cc(ni,ci)))
      thsum=thopnp+thomnp+thomnm
      thopnp=thopnp/thsum
      thomnp=thomnp/thsum
      thomnm=thomnm/thsum

      psum=cc(p1,ci)+cc(p2,ci)+cc(p3,ci)+p10+p20+p30 
      llda=lda*(1.+beta_da*yy(tda,t(ci)))
      llan=th(cc(o2,ci),_ZERO_,_ZERO_,_ONE_)*cc(o2,ci)/(oan+cc(o2,ci))      &
              *lan*exp(beta_an*t(ci))
      if ((fluff).and.(ci.eq.1)) then
         llsa=lsa*exp(bsa*t(ci))*(th(cc(o2,ci),wo,dot2,_ONE_))
      end if
      r1=r1max*min(yy(alpha1,cc(am,ci)+cc(ni,ci)),yy(sr*alpha1,cc(po,ci)),   &
                   ppi(ci))
      r2=r2max*(1.+yy(tf,t(ci)))*                                            &
               min(yy(alpha2,cc(am,ci)+cc(ni,ci)),yy(sr*alpha2,cc(po,ci)),   &
                   ppi(ci))
      r3=r3max*1./(1.+exp(beta_bg*(tbg-t(ci))))                              &
                    *min(yy(sr*alpha3,cc(po,ci)),ppi(ci))

!  Sink terms for non-negative compartments, which appear exactly
!  as or proportional to source terms for other compartments:
      dd(p1,zo,ci)=fpz(g1max,t(ci),topt,psum)*cc(p1,ci)/psum*(cc(zo,ci)+zo0)
      dd(p1,de,ci)=lpd*cc(p1,ci)
      dd(p1,am,ci)=lpa*cc(p1,ci)
      dd(p2,zo,ci)=fpz(g2max,t(ci),topt,psum)*cc(p2,ci)/psum*(cc(zo,ci)+zo0)
      dd(p2,de,ci)=lpd*cc(p2,ci)
      dd(p2,am,ci)=lpa*cc(p2,ci)
      dd(p3,zo,ci)=fpz(g3max,t(ci),topt,psum)*cc(p3,ci)/psum*(cc(zo,ci)+zo0)
      dd(p3,de,ci)=lpd*cc(p3,ci)
      dd(p3,am,ci)=lpa*cc(p3,ci)
      dd(zo,de,ci)=lzd*(cc(zo,ci)+zo0)*cc(zo,ci)
      dd(zo,am,ci)=lza*(cc(zo,ci)+zo0)*cc(zo,ci)
      dd(de,am,ci)=llda*cc(de,ci)
      dd(am,p1,ci)=cc(am,ci)/(cc(am,ci)+cc(ni,ci))*r1*(cc(p1,ci)+p10)
      dd(am,p2,ci)=cc(am,ci)/(cc(am,ci)+cc(ni,ci))*r2*(cc(p2,ci)+p20)
      dd(am,ni,ci)=llan*cc(am,ci)
      dd(ni,p1,ci)=cc(ni,ci)/(cc(ni,ci)+cc(am,ci))*r1*(cc(p1,ci)+p10)
      dd(ni,p2,ci)=cc(ni,ci)/(cc(ni,ci)+cc(am,ci))*r2*(cc(p2,ci)+p20)
      if ((fluff).and.(ci.eq.1)) then
         dd(de,fl,ci)=lds*cc(de,ci)/h(ci)
         dd(fl,am,ci)=llsa*cc(fl,ci)
      end if

!  Sink terms for positive compartments, which do not appear 
!  as source terms for other compartments:
      dd(ni,ni,ci)=s1*llda*cc(de,ci)*thomnp    ! denitrification
      dd(po,po,ci)=sr*( r1*(cc(p1,ci)+p10)+r2*(cc(p2,ci)+p20)                &
                       +r3*(cc(p3,ci)+p30)) 
      
      if ((fluff).and.(ci.eq.1)) then
         dd(fl,fl,ci)=th(cc(o2,ci),wo,_ZERO_,_ONE_)*dd(fl,am,ci)
         dd(ni,ni,ci)=dd(ni,ni,ci)+s1*thomnp*dd(fl,am,ci)/h(ci)
      end if

!  Source terms which are exactly sinks terms of other compartments or
!   proportional to them:
      do i=1,numc
         do j=1,numc
            if (i.ne.j) pp(i,j,ci)=rat(j,i)*dd(j,i,ci)
         end do
      end do

!   Non-conservative source terms or source and sink terms which are 
!   stoichiometrically related to other source terms:
      pp(p3,p3,ci)=r3*(cc(p3,ci)+p30)     ! nitrogen fixation
      pp(po,p1,ci)=sr*dd(p1,am,ci)   
      pp(po,p2,ci)=sr*dd(p2,am,ci)  
      pp(po,p3,ci)=sr*dd(p3,am,ci)  
      pp(po,de,ci)=sr*dd(de,am,ci)   
      pp(po,zo,ci)=sr*dd(zo,am,ci)    
      pp(o2,am,ci)=s2*(dd(am,p1,ci)+dd(am,p2,ci))-s4*dd(am,ni,ci)    
      pp(o2,ni,ci)=s3*(dd(ni,p1,ci)+dd(ni,p2,ci))  
      pp(o2,o2,ci)=s2*pp(p3,p3,ci)        ! nitrogen fixation
      pp(o2,p1,ci)=-s2*dd(p1,am,ci)  
      pp(o2,p2,ci)=-s2*dd(p2,am,ci)  
      pp(o2,p3,ci)=-s2*dd(p3,am,ci)  
      pp(o2,zo,ci)=-s2*dd(zo,am,ci)  
!      pp(o2,am,ci)=-s4*dd(am,ni,ci) 
      pp(o2,de,ci)=-s2*(thopnp+thomnm)*dd(de,am,ci)
      if ((fluff).and.(ci.eq.1)) then
         pp(o2,fl,ci)=-(s4+s2*(thopnp+thomnm))*dd(fl,am,ci)/h(ci)
         pp(po,po,ci)=(1.-ph1*th(cc(o2,ci),wo,_ZERO_,_ONE_)* &
                      yy(ph2,cc(o2,ci)))*dd(fl,am,ci)/h(i)
      end if
   end do

   return
   end subroutine do_bio_mab
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Finish the bio calculations
!
! !INTERFACE:
   subroutine end_bio_mab
!
! !DESCRIPTION:
!  Nothing done yet --- supplied for completeness.
!
! !USES:
   IMPLICIT NONE
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
!EOP
!-----------------------------------------------------------------------
!BOC

   return
   end subroutine end_bio_mab
!EOC

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

   end module bio_mab

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