Commit e74c0f0a955296af2392036e37d48959401205cc

Authored by dumoda01
1 parent 55e3f794
Exists in master and in 1 other branch snow

Ajout d'une routine et d'un fichier namelist pour le nouveau modele biologique b…

…io_ismer. Celui-ci est une version modifiee de Fasham auquel on a ajoute une classe de phyto et une classe de zoo supplementaires. Un terme de nitrification a ete ajoute. Les parametres par defaut devront etre ajustes.
Showing 2 changed files with 763 additions and 0 deletions   Show diff stats
nml/bio_ismer.nml 0 → 100644
... ... @@ -0,0 +1,112 @@
  1 +#$Id$
  2 +!-------------------------------------------------------------------------------
  3 +! Fasham et al. biological model with modifications by Kuehn and Radach
  4 +!
  5 +! numc= number of compartments for geobiochemical model
  6 +!
  7 +! p1_initial= initial flagellate concentration [mmol n/m3]
  8 +! p2_initial= initial diatom concentration [mmol n/m3]
  9 +! z1_initial= initial micro-zooplakton concentration [mmol n/m3]
  10 +! z2_initial= initial meso-zooplakton concentration [mmol n/m3]
  11 +! b_initial= initial bacteria concentration [mmol n/m3]
  12 +! d_initial= initial detritus concentration [mmol n/m3]
  13 +! n_initial= *** see obs.nml *** [mmol n/m3]
  14 +! a_initial= *** see obs.nml *** [mmol n/m3]
  15 +! l_initial= initial LDON concentration [mmol n/m3]
  16 +! p0 = minimum phytoplankton concentration [mmol n/m3]
  17 +! z0 = minimum zooplakton concentration [mmol n/m3]
  18 +! b0 = minimum bacteria concentration [mmol n/m3]
  19 +! theta = phytoplancton buoyancy parameter [m3 day/(mmol N)] !CHG2
  20 +! vp1 = maximum flagellate uptake rate by flagellates [1/day]
  21 +! vp2 = maximum diatom uptake rate by diatoms [1/day]
  22 +! alpha1 = slope of the flagellate PI-curve [m2/(W day)]
  23 +! alpha2 = slope of the diatom PI-curve [m2/(W day)]
  24 +! inib1 = inhibition slope of the flagellate PI-curve (pos.) [m2/(W day)] !CHG1
  25 +! inib2 = inhibition slope of the PI-curve (pos.) [m2/(W day)] !CHG1
  26 +! kn1 = half sat. constant nitrate uptake by fla [mmol n/m3]
  27 +! ka1 = half sat. constant ammonium uptake by fla [mmol n/m3]
  28 +! kn2 = half sat. constant nitrate uptake by diatoms [mmol n/m3]
  29 +! ka2 = half sat. constant ammonium uptake by diatoms [mmol n/m3]
  30 +! mu1 = phytoplankton (fla & dia) mortality rate [1/day]
  31 +! k5 = half sat. constant phy. mortality (fla & dia) [mmol n/m3]
  32 +! gamma = exudation fraction [-]
  33 +! w_p1 = flagellate settling velocity [m/day]
  34 +! w_p2 = diatom settling velocity [m/day]
  35 +! g1max = maximum microzooplankton ingestion rate [1/day]
  36 +! g2max = maximum mesozooplankton ingestion rate [1/day]
  37 +! k3 = half saturation constant ingestion [mmol n/m3]
  38 +! beta = grazing efficiency [-]
  39 +! mu2 = maximum zooplankton loss rate (mcz & msz) [1/day]
  40 +! k6 = half saturation zooplankton loss (mcz & msz) [mmol n/m3]
  41 +! delta = fractional zooplankton loss to LDON (mcz & msz) [-]
  42 +! epsi = fractional zooplankton loss to ammonium (mcz & msz) [-]
  43 +! r11 = mcz grazing preference on flagellates [-]
  44 +! r12 = mcz grazing preference on bacteria [-]
  45 +! r13 = mcz grazing preference on detritus [-]
  46 +! r21 = msz grazing preference on flagellates [-]
  47 +! r22 = msz grazing preference on diatoms [-]
  48 +! r23 = msz grazing preference on detritus [-]
  49 +! vb = maximum bacterial uptake rate [1/day]
  50 +! k4 = half saturation bacterial uptake [mmol n/m3]
  51 +! mu3 = bacteria excretion rate [1/day]
  52 +! eta = uptake ratio ammonium:LDON [-]
  53 +! mu4 = detritus breakdown rate [1/day]
  54 +! mu5 = nitrification rate [1/day]
  55 +! w_d = detritus settling velocity [m/day]
  56 +! kc = attenuation constant for the self shading effect [m**2/mmol N]
  57 +!-------------------------------------------------------------------------------
  58 + &bio_ismer_nml
  59 + numc= 9
  60 + p1_initial= 0.012
  61 + p2_initial= 0.012
  62 + z1_initial= 0.012
  63 + z2_initial= 0.012
  64 + b_initial= 0.001
  65 + d_initial= 0.01
  66 + l_initial= 0.1
  67 + p0= 0.0001
  68 + z0= 0.0001
  69 + b0= 0.0001
  70 + vp1= 0.3
  71 + vp2= 0.3
  72 + alpha1= 0.04
  73 + alpha2= 0.04
  74 + inib1= 0.06
  75 + inib2= 0.06
  76 + kn1= 1.0
  77 + ka1= 0.8
  78 + kn2= 1.0
  79 + ka2= 0.8
  80 + mu1= 0.05
  81 + k5= 0.2
  82 + gamma= 0.05
  83 + w_p1= -0.38
  84 + w_p2= -0.10
  85 + theta= 0.0
  86 + w_p1min= -0.06
  87 + w_p1max= -0.38
  88 + w_p2min= -0.06
  89 + w_p2max= -0.38
  90 + g1max= 1.0
  91 + g2max= 1.0
  92 + k3= 1.0
  93 + beta= 0.625
  94 + mu2= 0.3
  95 + k6= 0.2
  96 + delta= 0.1
  97 + epsi= 0.70
  98 + r11= 0.55
  99 + r12= 0.4
  100 + r13= 0.05
  101 + r21= 0.55
  102 + r22= 0.4
  103 + r23= 0.05
  104 + vb= 0.24
  105 + k4= 0.5
  106 + mu3= 0.03
  107 + eta= 0.0
  108 + mu4= 0.02
  109 + mu5= 0.00
  110 + w_d= -5.0
  111 + kc= 0.03
  112 + /
... ...
src/extras/bio/bio_ismer.F90 0 → 100644
... ... @@ -0,0 +1,651 @@
  1 +!$Id: bio_fasham.F90,v 1.11 2007-01-06 11:49:15 kbk Exp $
  2 +#include"cppdefs.h"
  3 +!-----------------------------------------------------------------------
  4 +!BOP
  5 +!
  6 +! !MODULE: bio_ismer --- Modified from Fasham et al. biological model \label{sec:bio-fasham}
  7 +!
  8 +! !INTERFACE:
  9 + module bio_ismer
  10 +!
  11 +! !DESCRIPTION:
  12 +! The model developed by \cite{Fashametal1990}
  13 +! uses nitrogen as 'currency' according to the evidence that in
  14 +! most cases nitrogen is the limiting macronutrient. It consists of
  15 +! seven state variables: phytoplankton, zooplankton, bacteria,
  16 +! particulate organic matter (detritus), dissolved organic matter
  17 +! and the nutrients nitrate and ammonium.
  18 +! The structure of the \cite{Fashametal1990} biogeochemical model
  19 +! is given in figure \ref{fig_fasham}.
  20 +! \begin{figure}
  21 +! \begin{center}
  22 +! \scalebox{0.5}{\includegraphics{figures/fasham_structure.eps}}
  23 +! \caption{Structure of the \cite{Fashametal1990} model with bacteria (bac),
  24 +! phytoplankton (phy), detritus (det), zooplankton (zoo), labile dissolved
  25 +! organic nitrogen (don), ammonium (amm) and nitrate (nit) as the seven
  26 +! state variables.
  27 +! The concentrations are in mmol N\,m$^{-3}$,
  28 +! all fluxes (green arrows) are conservative.
  29 +! }\label{fig_fasham}
  30 +! \end{center}
  31 +! \end{figure}
  32 +! A detailed mathematical description of all
  33 +! processes is given in section \ref{sec:bio-fasham-rhs}.
  34 +! The version of the \cite{Fashametal1990} model which is implemented includes
  35 +! slight modifications by \cite{KuehnRadach1997} and has been
  36 +! included into GOTM by \cite{Burchardetal05}.
  37 +
  38 +! !USES:
  39 +! default: all is private.
  40 + use bio_var
  41 + use output
  42 + use observations, only : aa,g2
  43 + private
  44 +!
  45 +! !PUBLIC MEMBER FUNCTIONS:
  46 + public init_bio_ismer, init_var_ismer, var_info_ismer, &
  47 + light_ismer, do_bio_ismer, end_bio_ismer
  48 +!
  49 +! !PRIVATE DATA MEMBERS:
  50 +!
  51 +! !REVISION HISTORY:!
  52 +! Original author(s): Hans Burchard & Karsten Bolding
  53 +!
  54 +!
  55 +! !LOCAL VARIABLES:
  56 +! from a namelist
  57 + REALTYPE :: p1_initial= 0.05
  58 + REALTYPE :: p2_initial= 0.05
  59 + REALTYPE :: z1_initial= 0.05
  60 + REALTYPE :: z2_initial= 0.05
  61 + REALTYPE :: b_initial= 0.001
  62 + REALTYPE :: d_initial= 0.416666666
  63 +! Nitrate and ammonium are initialized within the GOTM observation module
  64 +! REALTYPE :: n_initial= 8.3
  65 +! REALTYPE :: a_initial= 0.22
  66 + REALTYPE :: l_initial= 0.14
  67 + REALTYPE :: p0 = 0.0
  68 + REALTYPE :: z0 = 0.0
  69 + REALTYPE :: b0 = 0.0
  70 + REALTYPE :: vp1 = 1.5
  71 + REALTYPE :: alpha1 = 0.065
  72 + REALTYPE :: inib1 = 0.05
  73 + REALTYPE :: vp2 = 1.5
  74 + REALTYPE :: alpha2 = 0.065
  75 + REALTYPE :: inib2 = 0.05
  76 + REALTYPE :: theta = 0.0
  77 + REALTYPE :: w_p1min = -0.06
  78 + REALTYPE :: w_p1max = -0.38
  79 + REALTYPE :: w_p2min = -0.06
  80 + REALTYPE :: w_p2max = -0.38
  81 + REALTYPE :: kn1 = 0.2
  82 + REALTYPE :: ka1 = 0.8
  83 + REALTYPE :: kn2 = 0.2
  84 + REALTYPE :: ka2 = 0.8
  85 + REALTYPE :: mu1 = 0.05
  86 + REALTYPE :: k5 = 0.2
  87 + REALTYPE :: gamma = 0.05
  88 + REALTYPE :: w_p1 = -0.5
  89 + REALTYPE :: w_p2 = -0.5
  90 + REALTYPE :: g1max = 1.0
  91 + REALTYPE :: g2max = 1.0
  92 + REALTYPE :: k3 = 1.0
  93 + REALTYPE :: beta = 0.625
  94 + REALTYPE :: mu2 = 0.3
  95 + REALTYPE :: k6 = 0.2
  96 + REALTYPE :: delta = 0.1
  97 + REALTYPE :: epsi = 0.70
  98 + REALTYPE :: r11 = 0.55
  99 + REALTYPE :: r12 = 0.4
  100 + REALTYPE :: r13 = 0.05
  101 + REALTYPE :: r21 = 0.55
  102 + REALTYPE :: r22 = 0.4
  103 + REALTYPE :: r23 = 0.05
  104 + REALTYPE :: vb = 1.2
  105 + REALTYPE :: k4 = 0.5
  106 + REALTYPE :: mu3 = 0.15
  107 + REALTYPE :: eta = 0.0
  108 + REALTYPE :: mu4 = 0.02
  109 + REALTYPE :: mu5 = 0.00
  110 + REALTYPE :: w_d = -2.0
  111 + REALTYPE, public :: kc = 0.03
  112 + integer :: out_unit
  113 + integer, parameter :: n=1,p1=2,p2=3,z1=4,z2=5,d=6,l=7,b=8,a=9
  114 +!EOP
  115 +!-----------------------------------------------------------------------
  116 +
  117 + contains
  118 +
  119 +!-----------------------------------------------------------------------
  120 +!BOP
  121 +!
  122 +! !IROUTINE: Initialise the bio module
  123 +!
  124 +! !INTERFACE:
  125 + subroutine init_bio_ismer(namlst,fname,unit)
  126 +!
  127 +! !DESCRIPTION:
  128 +! Here, the bio namelist {\tt bio\_fasham.nml} is read and
  129 +! various variables (rates and settling velocities)
  130 +! are transformed into SI units.
  131 +!
  132 +! !USES:
  133 + IMPLICIT NONE
  134 +!
  135 +! !INPUT PARAMETERS:
  136 + integer, intent(in) :: namlst
  137 + character(len=*), intent(in) :: fname
  138 + character(len=20) :: pfile
  139 + integer, intent(in) :: unit
  140 +!
  141 +! !REVISION HISTORY:
  142 +! Original author(s): Hans Burchard & Karsten Bolding
  143 +!
  144 +! !LOCAL VARIABLES:
  145 + namelist /bio_ismer_nml/ numc, &
  146 + p1_initial,p2_initial,z1_initial,z2_initial, &
  147 + b_initial,d_initial,l_initial, &
  148 + p0,z0,b0,vp1,alpha1,inib1,vp2,alpha2,inib2, &
  149 + kn1,ka1,kn2,ka2,mu1,k5,gamma,w_p1,w_p2, &
  150 + g1max,g2max,k3,beta,mu2,k6,delta,epsi, &
  151 + r11,r12,r13,r21,r22,r23, &
  152 + vb,k4,mu3,eta,mu4,w_d,kc,mu5, &
  153 + theta,w_p1max,w_p1min,w_p2min,w_p2max !CHG2
  154 +!EOP
  155 +!-----------------------------------------------------------------------
  156 +!BOC
  157 + LEVEL2 'init_bio_ismer'
  158 +
  159 + open(namlst,file=fname,action='read',status='old',err=98)
  160 + read(namlst,nml=bio_ismer_nml,err=99)
  161 + close(namlst)
  162 +
  163 + numcc=numc
  164 +
  165 +! Print some parameter values in standard output
  166 +! and save them in a separate file [out_fn]_ismer.par
  167 + pfile = trim(out_fn) // '_ismer.par'
  168 + open(10,status='unknown',action='write',file=pfile)
  169 + LEVEL3 'ISMER parameters saved in ', pfile
  170 +! write(*,900) ' vp = ',vp
  171 +! write(10,901) vp
  172 +! write(*,900) ' alpha = ',alpha
  173 +! write(10,901) alpha
  174 +! write(*,900) ' inib = ',inib
  175 +! write(10,901) inib
  176 +
  177 +900 format (a,f8.5)
  178 +901 format (f8.5)
  179 +
  180 +! Conversion from day to second
  181 + vp1 = vp1 /secs_pr_day
  182 + vp2 = vp2 /secs_pr_day
  183 + vb = vb /secs_pr_day
  184 + mu1 = mu1 /secs_pr_day
  185 + mu2 = mu2 /secs_pr_day
  186 + mu3 = mu3 /secs_pr_day
  187 + mu4 = mu4 /secs_pr_day
  188 + mu5 = mu5 /secs_pr_day !DD
  189 + g1max = g1max /secs_pr_day
  190 + g2max = g2max /secs_pr_day
  191 + w_p1 = w_p1 /secs_pr_day
  192 + w_p2 = w_p2 /secs_pr_day
  193 + w_p1min = w_p1min /secs_pr_day !DD
  194 + w_p1max = w_p1max /secs_pr_day !DD
  195 + w_p2min = w_p2min /secs_pr_day !DD
  196 + w_p2max = w_p2max /secs_pr_day !DD
  197 + theta = theta /secs_pr_day !DD
  198 + w_d = w_d /secs_pr_day
  199 + alpha1 = alpha1/secs_pr_day
  200 + inib1 = inib1 /secs_pr_day !CHG1
  201 + alpha2 = alpha2/secs_pr_day
  202 + inib2 = inib2 /secs_pr_day !CHG1
  203 +
  204 + out_unit=unit
  205 +
  206 + LEVEL3 'ISMER bio module initialised ...'
  207 +
  208 + return
  209 +
  210 +98 LEVEL2 'I could not open bio_ismer.nml'
  211 + LEVEL2 'If thats not what you want you have to supply bio_ismer.nml'
  212 + LEVEL2 'See the bio example on www.gotm.net for a working bio_ismer.nml'
  213 + return
  214 +99 FATAL 'I could not read bio_ismer.nml'
  215 + stop 'init_bio_ismer'
  216 + end subroutine init_bio_ismer
  217 +!EOC
  218 +
  219 +!-----------------------------------------------------------------------
  220 +!BOP
  221 +!
  222 +! !IROUTINE: Initialise the concentration variables
  223 +!
  224 +! !INTERFACE:
  225 + subroutine init_var_ismer(nlev)
  226 +!
  227 +! !DESCRIPTION:
  228 +! Here, the the initial conditions are set and the settling velocities are
  229 +! transferred to all vertical levels. All concentrations are declared
  230 +! as non-negative variables, and it is defined which variables would be
  231 +! taken up by benthic filter feeders.
  232 +!
  233 +! !USES:
  234 + use observations, only: nprof,aprof !CHG3-5
  235 + use meanflow, only: nit,amm,T,S !CHG3-5
  236 +
  237 + IMPLICIT NONE
  238 +
  239 +!
  240 +! !INPUT PARAMETERS:
  241 + integer, intent(in) :: nlev
  242 +!
  243 +! !REVISION HISTORY:
  244 +! Original author(s): Hans Burchard & Karsten Bolding
  245 +
  246 +! !LOCAL VARIABLES:
  247 + integer :: i
  248 +!EOP
  249 +!-----------------------------------------------------------------------
  250 +!BOC
  251 + do i=1,nlev
  252 + cc(n,i) = nprof(i) !CHG3
  253 + cc(p1,i)= p1_initial
  254 + cc(p2,i)= p2_initial
  255 + cc(z1,i)= z1_initial
  256 + cc(z2,i)= z2_initial
  257 + cc(d,i) = d_initial
  258 + cc(l,i) = l_initial
  259 + cc(b,i) = b_initial
  260 + cc(a,i) = aprof(i) !CHG5
  261 + end do
  262 +
  263 + do i=0,nlev
  264 + ws(n,i) = _ZERO_
  265 + ws(p1,i) = w_p1
  266 + ws(p2,i) = w_p2
  267 + ws(z1,i) = _ZERO_
  268 + ws(z2,i) = _ZERO_
  269 + ws(d,i) = w_d
  270 + ws(l,i) = _ZERO_
  271 + ws(b,i) = _ZERO_
  272 + ws(a,i) = _ZERO_
  273 + end do
  274 +
  275 + posconc(n) = 1
  276 + posconc(p1) = 1
  277 + posconc(p2) = 1
  278 + posconc(z1) = 1
  279 + posconc(z2) = 1
  280 + posconc(d) = 1
  281 + posconc(l) = 1
  282 + posconc(b) = 1
  283 + posconc(a) = 1
  284 +
  285 + LEVEL3 'ISMER variables initialised ...'
  286 +
  287 + return
  288 +
  289 + end subroutine init_var_ismer
  290 +!EOC
  291 +
  292 +!-----------------------------------------------------------------------
  293 +!BOP
  294 +!
  295 +! !IROUTINE: Providing info on variables
  296 +!
  297 +! !INTERFACE:
  298 + subroutine var_info_ismer()
  299 +!
  300 +! !DESCRIPTION:
  301 +! This subroutine provides information about the variable names as they
  302 +! will be used when storing data in NetCDF files.
  303 +!
  304 +! !USES:
  305 + IMPLICIT NONE
  306 +!
  307 +! !REVISION HISTORY:
  308 +! Original author(s): Hans Burchard & Karsten Bolding
  309 +!
  310 +! !LOCAL VARIABLES:
  311 +!EOP
  312 +!-----------------------------------------------------------------------
  313 +!BOC
  314 + var_names(1) = 'nit'
  315 + var_units(1) = 'mmol/m**3'
  316 + var_long(1) = 'nitrate'
  317 +
  318 + var_names(2) = 'fla'
  319 + var_units(2) = 'mmol/m**3'
  320 + var_long(2) = 'flagellates'
  321 +
  322 + var_names(3) = 'dia'
  323 + var_units(3) = 'mmol/m**3'
  324 + var_long(3) = 'diatoms'
  325 +
  326 + var_names(4) = 'mcz'
  327 + var_units(4) = 'mmol/m**3'
  328 + var_long(4) = 'micro-zooplankton'
  329 +
  330 + var_names(5) = 'msz'
  331 + var_units(5) = 'mmol/m**3'
  332 + var_long(5) = 'meso-zooplankton'
  333 +
  334 + var_names(6) = 'det'
  335 + var_units(6) = 'mmol/m**3'
  336 + var_long(6) = 'detritus'
  337 +
  338 + var_names(7) = 'ldn'
  339 + var_units(7) = 'mmol/m**3'
  340 + var_long(7) = 'labile_dissolved_organic_nitrogen'
  341 +
  342 + var_names(8) = 'bac'
  343 + var_units(8) = 'mmol/m**3'
  344 + var_long(8) = 'bacteria'
  345 +
  346 + var_names(9) = 'amm'
  347 + var_units(9) = 'mmol/m**3'
  348 + var_long(9) = 'ammonium'
  349 +
  350 + return
  351 + end subroutine var_info_ismer
  352 +!EOC
  353 +
  354 +!-----------------------------------------------------------------------
  355 +!BOP
  356 +!
  357 +! !IROUTINE: Light properties for the ISMER model
  358 +!
  359 +! !INTERFACE:
  360 + subroutine light_ismer(nlev,bioshade_feedback)
  361 +!
  362 +! !DESCRIPTION:
  363 +! Here, the photosynthetically available radiation is calculated
  364 +! by simply assuming that the short wave part of the total
  365 +! radiation is available for photosynthesis.
  366 +! The photosynthetically
  367 +! available radiation, $I_{PAR}$, follows from (\ref{light}).
  368 +! The user should make
  369 +! sure that this is consistent with the light class given in the
  370 +! {\tt extinct} namelist of the {\tt obs.nml} file.
  371 +! The self-shading effect is also calculated in this subroutine,
  372 +! which may be used to consider the effect of bio-turbidity also
  373 +! in the temperature equation (if {\tt bioshade\_feedback} is set
  374 +! to true in {\tt bio.nml}).
  375 +! For details, see section \ref{sec:do-bio}.
  376 +!
  377 +! !USES:
  378 + IMPLICIT NONE
  379 +!
  380 +! !INPUT PARAMETERS:
  381 + integer, intent(in) :: nlev
  382 + logical, intent(in) :: bioshade_feedback
  383 +!
  384 +! !REVISION HISTORY:
  385 +! Original author(s): Hans Burchard, Karsten Bolding
  386 +!
  387 +! !LOCAL VARIABLES:
  388 + integer :: i
  389 + REALTYPE :: zz,add
  390 +!EOP
  391 +!-----------------------------------------------------------------------
  392 +!BOC
  393 + zz = _ZERO_
  394 + add = _ZERO_
  395 + do i=nlev,1,-1
  396 + add=add+0.5*h(i)*(cc(p1,i)+cc(p2,i)+p0)
  397 + zz=zz+0.5*h(i)
  398 + par(i)=rad(nlev)*(1.-aa)*exp(-zz/g2)*exp(-kc*add)
  399 + add=add+0.5*h(i)*(cc(p1,i)+cc(p2,i)+p0)
  400 + zz=zz+0.5*h(i)
  401 + if (bioshade_feedback) bioshade_(i)=exp(-kc*add)
  402 + end do
  403 +
  404 +
  405 + return
  406 + end subroutine light_ismer
  407 +!EOC
  408 +
  409 +!-----------------------------------------------------------------------
  410 +!BOP
  411 +!
  412 +! !IROUTINE: Right hand sides of geobiochemical model \label{sec:bio-fasham-rhs}
  413 +!
  414 +! !INTERFACE:
  415 + subroutine do_bio_ismer(first,numc,nlev,cc,pp,dd)
  416 +!
  417 +! !DESCRIPTION:
  418 +!
  419 +! The \cite{Fashametal1990} model consisting of the $I=7$
  420 +! state variables phytoplankton, bacteria, detritus, zooplankton,
  421 +! nitrate, ammonium and dissolved organic nitrogen is described here
  422 +! in detail.
  423 +!
  424 +! Phytoplankton mortality and zooplankton grazing loss of phytoplankton:
  425 +! \begin{equation}\label{d13}
  426 +! d_{1,3} = \mu_1 \frac{c_1+c_{1}^{\min}}{K_5+c_1+c_{1}^{\min}}c_1+
  427 +! (1-\beta)\frac{g\rho_1 c_1^2}{K_3 \sum_{j=1}^3 \rho_jc_j
  428 +! + \sum_{j=1}^3 \rho_jc_j^2} (c_4+c_{4}^{\min}).
  429 +! \end{equation}
  430 +! Phytoplankton loss to LDON (labile dissolved organic nitrogen):
  431 +! \begin{equation}\label{d17}
  432 +! d_{1,7} = \gamma
  433 +! F(I_{PAR})\frac{\frac{c_5}{K_1}
  434 +! +\frac{c_6}{K_2}}{1+\frac{c_5}{K_1}+\frac{c_6}{K_2}}c_1,
  435 +! \end{equation}
  436 +! with
  437 +! \begin{equation}\label{FI}
  438 +! F(I_{PAR}) = \frac{V_p\alpha I_{PAR}(z)}{\left(V_p^2+\alpha^2(I_{PAR}(z))^2
  439 +! \right)^{1/2}}.
  440 +! \end{equation}
  441 +! With $I_{PAR}$ from (\ref{light}).
  442 +!
  443 +! Zooplankton grazing loss:
  444 +! \begin{equation}\label{di3}
  445 +! d_{2,3} = (1-\beta)\frac{g\rho_2 c_2^2}{K_3 \sum_{j=1}^3 \rho_jc_j
  446 +! + \sum_{j=1}^3 \rho_jc_j^2} (c_4+c_{4}^{\min}).
  447 +! \end{equation}
  448 +! Zooplankton grazing:
  449 +! \begin{equation}\label{di4}
  450 +! d_{i,4} = \beta\frac{g\rho_i c_i^2}{K_3 \sum_{j=1}^3 \rho_jc_j
  451 +! + \sum_{j=1}^3 \rho_jc_j^2} (c_4+c_{4}^{\min}), \quad i=1,\dots,3.
  452 +! \end{equation}
  453 +! Bacteria excretion rate:
  454 +! \begin{equation}\label{d26}
  455 +! d_{2,6} = \mu_3 c_2.
  456 +! \end{equation}
  457 +! Detritus breakdown rate:
  458 +! \begin{equation}\label{d37}
  459 +! d_{3,7} = \mu_4 c_3.
  460 +! \end{equation}
  461 +! Zooplankton losses to detritus, ammonium and LDON:
  462 +! \begin{equation}\label{d43}
  463 +! d_{4,3} = (1-\epsilon-\delta)\mu_2
  464 +! \frac{c_4+c_{4}^{\min}}{K_6+c_4+c_{4}^{\min}}c_4.
  465 +! \end{equation}
  466 +! \begin{equation}\label{d46}
  467 +! d_{4,6} = \epsilon\mu_2 \frac{c_4+c_{4}^{\min}}{K_6+c_4+c_{4}^{\min}}c_4.
  468 +! \end{equation}
  469 +! \begin{equation}\label{d47}
  470 +! d_{4,7} = \delta\mu_2 \frac{c_4+c_{4}^{\min}}{K_6+c_4+c_{4}^{\min}}c_4.
  471 +! \end{equation}
  472 +! Nitrate uptake by phytoplankton:
  473 +! \begin{equation}\label{d51}
  474 +! d_{5,1} = F(I_{PAR})\frac{\frac{c_5}{K_1}}{1+\frac{c_5}{K_1}
  475 +! +\frac{c_6}{K_2}}(c_1+c_{1}^{\min}).
  476 +! \end{equation}
  477 +! Ammonium uptake by phytoplankton:
  478 +! \begin{equation}\label{d61}
  479 +! d_{6,1} = F(I_{PAR})\frac{\frac{c_6}{K_2}}{1+\frac{c_5}{K_1}
  480 +! +\frac{c_6}{K_2}}(c_1+c_{1}^{\min}).
  481 +! \end{equation}
  482 +! Ammonium uptake by bacteria:
  483 +! \begin{equation}\label{d62}
  484 +! d_{6,2} = V_b \frac{\min(c_6,\eta c_7)}{K_4+\min(c_6,\eta c_7)+c_7}
  485 +! (c_2+c_{2}^{\min}).
  486 +! \end{equation}
  487 +! LDON uptake by bacteria:
  488 +! \begin{equation}\label{d72}
  489 +! d_{7,2} = V_b \frac{c_7}{K_4+\min(c_6,\eta c_7)+c_7} (c_2+c_{2}^{\min}).
  490 +! \end{equation}
  491 +!
  492 +! !USES:
  493 + IMPLICIT NONE
  494 +!
  495 +! !INPUT PARAMETERS:
  496 + logical, intent(in) :: first
  497 + integer, intent(in) :: numc,nlev
  498 + REALTYPE, intent(in) :: cc(1:numc,0:nlev)
  499 +!
  500 +! !OUTPUT PARAMETERS:
  501 + REALTYPE, intent(out) :: pp(1:numc,1:numc,0:nlev)
  502 + REALTYPE, intent(out) :: dd(1:numc,1:numc,0:nlev)
  503 +!
  504 +! !REVISION HISTORY:
  505 +! Original author(s): Hans Burchard, Karsten Bolding
  506 +!
  507 +! !LOCAL VARIABLES:
  508 + REALTYPE :: fac1,fac2,fac3,minal,qn1,qa1,qn2,qa2
  509 + REALTYPE :: Ps1,Ps2,ff1,ff2 !CHG1
  510 + integer :: i,j,ci
  511 +!EOP
  512 +!-----------------------------------------------------------------------
  513 +!BOC
  514 +!KBK - is it necessary to initialise every time - expensive in a 3D model
  515 + pp = _ZERO_
  516 + dd = _ZERO_
  517 +
  518 + do ci=1,nlev
  519 +
  520 +!CHG1
  521 +! Smith (1936) - saturation (default)
  522 +! ff= vp*alpha*par(ci)/sqrt(vp**2+alpha**2*par(ci)**2)
  523 +! Blackman (1919)
  524 +! if (par(ci) .lt. vp/alpha) then
  525 +! ff=alpha*par(ci)
  526 +! else
  527 +! ff=vp
  528 +! endif
  529 +! Steele (1962) - inhibition
  530 +! ff= vp*((par(ci)/I_opt)*exp(1-(par(ci)/I_opt)))
  531 +! Parker (1974) - inhibition
  532 +! ff= vp*((par(ci)/I_opt)*exp(1-(par(ci)/I_opt)))**2
  533 +! Platt et al. (1980) - inhibition
  534 + Ps1= vp1/((alpha1/(alpha1+inib1))*(alpha1/(alpha1+inib1))**(inib1/alpha1))
  535 + ff1= Ps1*(1.-exp(-1.*alpha1*par(ci)/Ps1))*exp(-1.*inib1*par(ci)/Ps1)
  536 + Ps2= vp2/((alpha2/(alpha2+inib2))*(alpha2/(alpha2+inib2))**(inib2/alpha2))
  537 + ff2= Ps2*(1.-exp(-1.*alpha2*par(ci)/Ps2))*exp(-1.*inib2*par(ci)/Ps2)
  538 +! --------------------------------------------------------------------
  539 +
  540 +
  541 + qn1=(cc(n,ci)/kn1)/(1.+cc(n,ci)/kn1+cc(a,ci)/ka1)
  542 + qa1=(cc(a,ci)/ka1)/(1.+cc(n,ci)/kn1+cc(a,ci)/ka1)
  543 +
  544 + qn2=(cc(n,ci)/kn2)/(1.+cc(n,ci)/kn2+cc(a,ci)/ka2)
  545 + qa2=(cc(a,ci)/ka2)/(1.+cc(n,ci)/kn2+cc(a,ci)/ka2)
  546 +
  547 + fac1=(cc(z1,ci)+z0)/(k3*(r11*cc(p1,ci)+r12*cc(b,ci)+r13*cc(d,ci))+ &
  548 + r11*cc(p1,ci)**2+r12*cc(b,ci)**2+r13*cc(d,ci)**2)
  549 +
  550 + fac2=(cc(z2,ci)+z0)/(k3*(r21*cc(p1,ci)+r22*cc(p2,ci)+r23*cc(d,ci))+ &
  551 + r21*cc(p1,ci)**2+r22*cc(p2,ci)**2+r23*cc(d,ci)**2)
  552 +
  553 + minal=min(cc(a,ci),eta*cc(l,ci))
  554 +
  555 + ! Temperature-dependent function for zooplankton grazing (for test purposes)
  556 + ! This function should come from thermodynamical arguments
  557 + !fac3 = max(0.0,tanh(0.5*(T(ci)+1.9)))
  558 + fac3 = 1.0
  559 +
  560 + ! Light and nutrient limitation factors
  561 + lumlim(ci) =ff1
  562 + nitlim(ci) =qn1
  563 + ammlim(ci) =qa1
  564 +
  565 + dd(p1,d,ci)=mu1*(cc(p1,ci)+p0)/(k5+cc(p1,ci)+p0)*cc(p1,ci) &
  566 + +fac3*(1.-beta)*cc(p1,ci)**2*(g1max*r11*fac1 + g2max*r21*fac2)
  567 + dd(p2,d,ci)=mu1*(cc(p2,ci)+p0)/(k5+cc(p2,ci)+p0)*cc(p2,ci) &
  568 + +fac3*(1.-beta)*g2max*r12*cc(p2,ci)**2*fac2
  569 +
  570 + dd(p1,l,ci)=gamma*ff1*(qn1+qa1)*cc(p1,ci)
  571 + dd(p2,l,ci)=gamma*ff2*(qn2+qa2)*cc(p2,ci)
  572 +
  573 + dd(b,d,ci) =fac3*(1.-beta)*g1max*r12*cc(b,ci)**2*fac1
  574 +
  575 + dd(p1,z1,ci)=fac3*beta*g1max*r11*cc(p1,ci)**2*fac1
  576 + dd(b,z1,ci)=fac3*beta*g1max*r12*cc(b,ci)**2*fac1
  577 + dd(d,z1,ci)=fac3*beta*g1max*r13*cc(d,ci)**2*fac1
  578 +
  579 + dd(p1,z2,ci)=fac3*beta*g2max*r21*cc(p1,ci)**2*fac2
  580 + dd(p2,z2,ci)=fac3*beta*g2max*r22*cc(p2,ci)**2*fac2
  581 + dd(d,z2,ci)=fac3*beta*g2max*r23*cc(d,ci)**2*fac2
  582 +
  583 + dd(b,a,ci) =mu3*cc(b,ci)
  584 + dd(d,l,ci) =mu4*cc(d,ci)
  585 + dd(a,n,ci) =mu5*cc(a,ci)
  586 +
  587 + dd(z1,d,ci)=(1.-epsi-delta)*mu2*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  588 + dd(z1,a,ci)=epsi*mu2*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  589 + dd(z1,l,ci)=delta*mu2*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  590 +
  591 + dd(z2,d,ci)=(1.-epsi-delta)*mu2*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  592 + dd(z2,a,ci)=epsi*mu2*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  593 + dd(z2,l,ci)=delta*mu2*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  594 +
  595 + dd(n,p1,ci) =ff1*qn1*(cc(p1,ci)+p0)
  596 + dd(a,p1,ci) =ff1*qa1*(cc(p1,ci)+p0)
  597 + dd(n,p2,ci) =ff2*qn2*(cc(p2,ci)+p0)
  598 + dd(a,p2,ci) =ff2*qa2*(cc(p2,ci)+p0)
  599 +
  600 + dd(a,b,ci) =vb*minal/(k4+minal+cc(l,ci))*(cc(b,ci)+b0)
  601 + dd(l,b,ci) =vb*cc(l,ci)/(k4+minal+cc(l,ci))*(cc(b,ci)+b0)
  602 +
  603 +! Taux de chute du phytoplancton en fonction de la croissance
  604 +! ws(p,ci)=w_pmin+(w_pmax-w_pmin)*exp(-ff/theta)
  605 +! ws(p,ci)=w_pmax-(w_pmax-w_pmin)*(ff/theta/(1.+ff/theta))
  606 + ws(p1,ci)=w_p1
  607 + ws(p2,ci)=w_p2
  608 +
  609 + do i=1,numc
  610 + do j=1,numc
  611 + pp(i,j,ci)=dd(j,i,ci)
  612 + end do
  613 + end do
  614 + end do
  615 +
  616 + return
  617 + end subroutine do_bio_ismer
  618 +!EOC
  619 +
  620 +!-----------------------------------------------------------------------
  621 +!BOP
  622 +!
  623 +! !IROUTINE: Finish the bio calculations
  624 +!
  625 +! !INTERFACE:
  626 + subroutine end_bio_ismer
  627 +!
  628 +! !DESCRIPTION:
  629 +! Nothing done yet --- supplied for completeness.
  630 +!
  631 +! !USES:
  632 + IMPLICIT NONE
  633 +!
  634 +! !REVISION HISTORY:
  635 +! Original author(s): Hans Burchard & Karsten Bolding
  636 +!
  637 +!EOP
  638 +!-----------------------------------------------------------------------
  639 +!BOC
  640 +
  641 + return
  642 + end subroutine end_bio_ismer
  643 +!EOC
  644 +
  645 +!-----------------------------------------------------------------------
  646 +
  647 + end module bio_ismer
  648 +
  649 +!-----------------------------------------------------------------------
  650 +! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
  651 +!-----------------------------------------------------------------------
... ...