Commit 9cc189c0 authored by Jorn Bruggeman's avatar Jorn Bruggeman
Browse files

model clean-up

parent 50317b5f
module bb_model_library
use fabm_types, only: type_base_model_factory,type_base_model
use fabm_types, only: type_base_model_factory, type_base_model
implicit none
private
type,extends(type_base_model_factory) :: type_factory
type, extends(type_base_model_factory) :: type_factory
contains
procedure :: create
end type
type (type_factory),save,target,public :: bb_model_factory
type (type_factory), save, target, public :: bb_model_factory
contains
subroutine create(self,name,model)
subroutine create(self, name, model)
use bb_filter_feeder
use bb_lorenz63
use bb_passive
! Add new BB models here
class (type_factory),intent(in) :: self
class (type_factory), intent(in) :: self
character(*), intent(in) :: name
class (type_base_model),pointer :: model
class (type_base_model), pointer :: model
select case (name)
case ('passive'); allocate(type_bb_passive::model)
......
module examples_model_library
use fabm_types, only: type_base_model_factory,type_base_model
use fabm_types, only: type_base_model_factory, type_base_model
use examples_benthic_predator
use examples_duplicator
......@@ -14,13 +14,13 @@ module examples_model_library
private
type,extends(type_base_model_factory) :: type_factory
type, extends(type_base_model_factory) :: type_factory
contains
procedure :: initialize
procedure :: create
end type
type (type_factory),save,target,public :: examples_model_factory
type (type_factory), save, target, public :: examples_model_factory
contains
......@@ -34,10 +34,10 @@ contains
call self%type_base_model_factory%initialize()
end subroutine initialize
subroutine create(self,name,model)
class (type_factory),intent(in) :: self
subroutine create(self, name, model)
class (type_factory), intent(in) :: self
character(*), intent(in) :: name
class (type_base_model),pointer :: model
class (type_base_model), pointer :: model
select case (name)
case ('benthic_predator'); allocate(type_examples_benthic_predator::model)
......@@ -48,7 +48,7 @@ contains
case ('light_cycle'); allocate(type_examples_light_cycle::model)
! Add individual example models here
case default
call self%type_base_model_factory%create(name,model)
call self%type_base_model_factory%create(name, model)
end select
end subroutine
......
......@@ -17,14 +17,14 @@ module examples_npzd_model_library
procedure :: create
end type
type (type_factory),save,target,public :: examples_npzd_model_factory
type (type_factory), save, target, public :: examples_npzd_model_factory
contains
subroutine create(self,name,model)
class (type_factory),intent(in) :: self
subroutine create(self, name, model)
class (type_factory), intent(in) :: self
character(*), intent(in) :: name
class (type_base_model),pointer :: model
class (type_base_model), pointer :: model
select case (name)
case ('nut'); allocate(type_examples_npzd_nut::model)
......
......@@ -10,15 +10,8 @@
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: gotm\_ergom --- GOTM implementation of the biogeochemical model ERGOM
!
! !INTERFACE:
module gotm_ergom
! GOTM implementation of the biogeochemical model ERGOM
!
! !DESCRIPTION:
! The biogeochemical model by
! \cite{Neumannetal2002} consists of $I=10$
! state variables. The nutrient state variables are dissolved
......@@ -63,19 +56,16 @@
! }\label{fig_neumann}
! \end{center}
! \end{figure}
!
! !USES:
module gotm_ergom
use fabm_types
implicit none
private
!
! !PUBLIC MEMBER FUNCTIONS:
public type_gotm_ergom
!
! !PUBLIC_DERIVED_TYPES:
type, extends(type_base_model) :: type_gotm_ergom
type, extends(type_base_model), public :: type_gotm_ergom
! Variable identifiers
type (type_state_variable_id) :: id_p1,id_p2,id_p3,id_zo,id_de,id_am,id_ni,id_po,id_o2
type (type_bottom_state_variable_id) :: id_fl
......@@ -97,39 +87,21 @@
procedure :: do_surface
end type
!EOP
!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the ergom model
!
! !INTERFACE:
subroutine initialize(self, configunit)
!
! !DESCRIPTION:
! Here, parameter values are read and the variables exported by the model are registered with FABM
!
! !INPUT PARAMETERS:
class (type_gotm_ergom), intent(inout), target :: self
integer, intent(in) :: configunit
!
! !LOCAL VARIABLES:
real(rk) :: w_p1
real(rk) :: w_p2
real(rk) :: w_p3
real(rk) :: w_de
real(rk), parameter :: secs_pr_day=86400.0_rk
!EOP
!-----------------------------------------------------------------------
!BOC
! Store parameter values in our own derived type
! NB! All rates must be provided in values per day,
! and are converted here to values per second
! NB! All rates must be provided in values per day in the configuration file,
! and are converted here to values per second by specifying scale_factor=1.0_rk/secs_pr_day
call self%get_parameter(self%sfl_po, 'sfl_po', 'mmol P/m2/d', 'surface phosphate flux', default=0.0015_rk, scale_factor=1.0_rk/secs_pr_day)
call self%get_parameter(self%sfl_ni, 'sfl_ni', 'mmol N/m2/d', 'surface nitrate flux', default=0.09_rk, scale_factor=1.0_rk/secs_pr_day)
call self%get_parameter(self%sfl_am, 'sfl_am', 'mmol N/m2/d', 'surface ammonium flux', default=0.07_rk, scale_factor=1.0_rk/secs_pr_day)
......@@ -229,44 +201,31 @@ contains
call self%add_to_aggregate_variable(standard_variables%attenuation_coefficient_of_photosynthetic_radiative_flux, self%id_p3, scale_factor=self%kc)
call self%add_to_aggregate_variable(standard_variables%attenuation_coefficient_of_photosynthetic_radiative_flux, self%id_de, scale_factor=self%kc)
call self%add_to_aggregate_variable(standard_variables%attenuation_coefficient_of_photosynthetic_radiative_flux, (self%p10 + self%p20 + self%p30) * self%kc)
end subroutine initialize
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Right hand sides of ergom model
!
! !INTERFACE:
subroutine do(self, _ARGUMENTS_DO_)
!
! !INPUT PARAMETERS:
class (type_gotm_ergom), intent(in) :: self
_DECLARE_ARGUMENTS_DO_
!
! !LOCAL VARIABLES:
real(rk) :: p1,p2,p3,zo,am,ni,po,de,o2,par,I_0
real(rk) :: iopt,ppi,temp,psum,llda,llan,lp,r1,r2,r3
real(rk) :: wo=30.0_rk,wn=0.1_rk
real(rk) :: thopnp,thomnp,thomnm,thsum
real(rk) :: secs_pr_day=86400.0_rk
!EOP
!-----------------------------------------------------------------------
!BOC
! Enter spatial_loops (if any)
_LOOP_BEGIN_
! Retrieve current (local) state variable values
_GET_(self%id_p1,p1) !diatoms
_GET_(self%id_p2,p2) !flagellates
_GET_(self%id_p3,p3) !cyanobacteria
_GET_(self%id_zo,zo) !zooplankton
_GET_(self%id_de,de) !detritus
_GET_(self%id_am,am) !ammonium
_GET_(self%id_ni,ni) !nitrate
_GET_(self%id_po,po) !phosphate
_GET_(self%id_o2,o2) !oxygen
_GET_(self%id_p1,p1) ! diatoms
_GET_(self%id_p2,p2) ! flagellates
_GET_(self%id_p3,p3) ! cyanobacteria
_GET_(self%id_zo,zo) ! zooplankton
_GET_(self%id_de,de) ! detritus
_GET_(self%id_am,am) ! ammonium
_GET_(self%id_ni,ni) ! nitrate
_GET_(self%id_po,po) ! phosphate
_GET_(self%id_o2,o2) ! oxygen
! Retrieve current environmental conditions
_GET_(self%id_par,par) ! local photosynthetically active radiation
......@@ -298,10 +257,10 @@ contains
r3=self%r3max*1.0_rk/(1.0_rk+exp(self%beta_bg*(self%tbg-temp))) &
*min(yy(self%sr*self%alpha3,po),ppi)*(p3+self%p30)
_SET_ODE_(self%id_p1,r1-fpz(self,self%g1max,temp,self%topt,psum)*p1/psum*(zo+self%zo0)-lp*p1)
_SET_ODE_(self%id_p2,r2-fpz(self,self%g2max,temp,self%topt,psum)*p2/psum*(zo+self%zo0)-lp*p2)
_SET_ODE_(self%id_p3,r3-fpz(self,self%g3max,temp,self%topt,psum)*p3/psum*(zo+self%zo0)-lp*p3)
_SET_ODE_(self%id_zo,(fpz(self,self%g1max,temp,self%topt,psum)*p1+fpz(self,self%g2max,temp,self%topt,psum)*p2+fpz(self,self%g3max,temp,self%topt,psum)*p3)*(zo+self%zo0)/psum-self%lza*zo*(zo+self%zo0)-self%lzd*zo*(zo+self%zo0))
_SET_ODE_(self%id_p1,r1-fpz(self%iv,self%g1max,temp,self%topt,psum)*p1/psum*(zo+self%zo0)-lp*p1)
_SET_ODE_(self%id_p2,r2-fpz(self%iv,self%g2max,temp,self%topt,psum)*p2/psum*(zo+self%zo0)-lp*p2)
_SET_ODE_(self%id_p3,r3-fpz(self%iv,self%g3max,temp,self%topt,psum)*p3/psum*(zo+self%zo0)-lp*p3)
_SET_ODE_(self%id_zo,(fpz(self%iv,self%g1max,temp,self%topt,psum)*p1+fpz(self%iv,self%g2max,temp,self%topt,psum)*p2+fpz(self%iv,self%g3max,temp,self%topt,psum)*p3)*(zo+self%zo0)/psum-self%lza*zo*(zo+self%zo0)-self%lzd*zo*(zo+self%zo0))
_SET_ODE_(self%id_de,self%lpd*p1+self%lpd*p2+self%lpd*p3+self%lzd*zo*(zo+self%zo0)-llda*de)
_SET_ODE_(self%id_am,llda*de-llan*am-am/(am+ni)*r1-am/(am+ni)*r2+self%lpa*(p1+p2+p3)+self%lza*zo*(zo+self%zo0))
_SET_ODE_(self%id_ni,llan*am-ni/(am+ni)*r1-ni/(am+ni)*r2-self%s1*llda*de*thomnp)
......@@ -317,29 +276,16 @@ contains
! Leave spatial loops (if any)
_LOOP_END_
end subroutine do
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE:
!
! !INTERFACE:
subroutine do_bottom(self,_ARGUMENTS_DO_BOTTOM_)
!
! !INPUT PARAMETERS:
subroutine do_bottom(self, _ARGUMENTS_DO_BOTTOM_)
class (type_gotm_ergom), intent(in) :: self
_DECLARE_ARGUMENTS_DO_BOTTOM_
!
! !LOCAL VARIABLES:
real(rk) :: fl,amb,nib,pob,deb,oxb,taub,temp
real(rk) :: llds,llsd,llsa,wo=30.0_rk,wn=0.1_rk,dot2=0.2_rk
real(rk) :: thopnp,thomnp,thomnm,thsum
!EOP
!-----------------------------------------------------------------------
!BOC
if (.not. self%fluff) return
! Enter spatial loops over the horizontal domain (if any).
......@@ -389,53 +335,40 @@ contains
! Leave spatial loops over the horizontal domain (if any).
_HORIZONTAL_LOOP_END_
end subroutine do_bottom
!EOC
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: Weiss formula for the saturation oxygen (osat)
!
! !INTERFACE:
real(rk) function osat_weiss(t,s)
!
! !DESCRIPTION:
! Weiss formula for the saturation oxygen (osat) \cite{Weiss1970}:
!
! \begin{equation}\label{osat_weiss}
! 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}
!
! !INPUT PARAMETERS:
! Weiss formula for the saturation oxygen (osat) \cite{Weiss1970}:
!
! \begin{equation}\label{osat_weiss}
! 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}
elemental real(rk) function osat_weiss(t,s)
real(rk), intent(in) :: t, s
!
! !LOCAL VARIABLES:
real(rk) :: tk
real(rk), parameter :: aa1=-173.4292_rk
real(rk), parameter :: aa2=249.6339_rk
......@@ -446,36 +379,22 @@ contains
real(rk), parameter :: b3=-0.001700_rk
real(rk), parameter :: kelvin=273.16_rk
real(rk), parameter :: mol_per_liter=44.661_rk
!EOP
!-----------------------------------------------------------------------
!BOC
tk=(t+kelvin)*0.01_rk
osat_weiss=exp(aa1+aa2/tk+a3*log(tk)+a4*tk &
+s*(b1+(b2+b3*tk)*tk))*mol_per_liter
end function osat_weiss
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Surface fluxes for the ergom model
!
! !INTERFACE:
subroutine do_surface(self,_ARGUMENTS_DO_SURFACE_)
class (type_gotm_ergom),intent(in) :: self
subroutine do_surface(self, _ARGUMENTS_DO_SURFACE_)
class (type_gotm_ergom), intent(in) :: self
_DECLARE_ARGUMENTS_DO_SURFACE_
!
! !LOCAL VARIABLES:
real(rk) :: temp,wnd,salt,o2,ni,am,po
real(rk), parameter :: secs_pr_day=86400.0_rk
real(rk) :: p_vel,sc,flo2
integer, parameter :: newflux=1
!EOP
!-----------------------------------------------------------------------
!BOC
_HORIZONTAL_LOOP_BEGIN_
_HORIZONTAL_LOOP_BEGIN_
_GET_(self%id_temp,temp)
_GET_(self%id_salt,salt)
_GET_SURFACE_(self%id_wind,wnd)
......@@ -512,33 +431,20 @@ contains
_SET_SURFACE_EXCHANGE_(self%id_po,self%sfl_po)
_HORIZONTAL_LOOP_END_
end subroutine do_surface
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Step function
!
! !INTERFACE:
pure real(rk) 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}
!
! !INPUT PARAMETERS:
! Step function
!
! 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}
elemental real(rk) function th(x,w,min,max)
real(rk), intent(in) :: x,w,min,max
!
!EOP
!-----------------------------------------------------------------------
!BOC
if (w .gt. 1.e-10_rk) then
th=min+(max-min)*0.5_rk*(1.0_rk+tanh(x/w))
else
......@@ -548,70 +454,36 @@ contains
th=0.0_rk
end if
end if
end function th
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Michaelis-Menten formulation for nutrient uptake
!
! !INTERFACE:
pure real(rk) 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}
!
! !INPUT PARAMETERS:
! Functional response for nutrient uptake
!
! 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}
elemental real(rk) function yy(a,x)
real(rk), intent(in) :: a,x
!
! !REVISION HISTORY:
!
!EOP
!-----------------------------------------------------------------------
!BOC
yy=x**2/(a**2+x**2)
yy=x**2/(a**2+x**2)
end function yy
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: Ivlev formulation for zooplankton grazing on phytoplankton
!
! !INTERFACE:
pure real(rk) function fpz(self,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}
!
! !INPUT PARAMETERS:
class (type_gotm_ergom), intent(in) :: self
real(rk), intent(in) :: g,t,topt,psum
!
!EOP
!-----------------------------------------------------------------------
!BOC
fpz=g*(1.0_rk+t**2/topt**2*exp(1.0_rk-2.0_rk*t/topt))* &
(1.0_rk-exp(-self%iv**2*psum**2))
! Ivlev formulation for zooplankton grazing on phytoplankton
!
! 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}
pure real(rk) function fpz(iv,g,t,topt,psum)
real(rk), intent(in) :: iv,g,t,topt,psum
fpz = g * (1.0_rk+t**2/topt**2*exp(1.0_rk-2.0_rk*t/topt)) * (1.0_rk-exp(-iv**2*psum**2))
end function fpz
!EOC
!-----------------------------------------------------------------------
end module gotm_ergom
......
#include "fabm_driver.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: fabm_fasham --- Fasham et al. (1990, DSR I) biogeochemical model,
! Fasham et al. (1990, DSR I) biogeochemical model,
! with slight modifications by Kuehn \& Radach (1997, JRM),
! and minimum concentration for phytoplankton, zooplankton and bacteria
! as decribed in Burchard et al. (2006, JMS)
! taken from GOTM and adapted for FABM by Jorn Bruggeman
!
! !INTERFACE:
module gotm_fasham
!
! !DESCRIPTION:
! The model developed by \cite{Fashametal1990}
! uses nitrogen as 'currency' according to the evidence that in
! most cases nitrogen is the limiting macronutrient. It consists of
......@@ -39,22 +32,94 @@
! slight modifications by \cite{KuehnRadach1997} and has been
! included into GOTM by \cite{Burchardetal05}.
!
! !USES:
! The \cite{Fashametal1990} model consisting of the $I=7$
! state variables phytoplankton, bacteria, detritus, zooplankton,
! nitrate, ammonium and dissolved organic nitrogen is described here
! in detail.
!
! Phytoplankton mortality and zooplankton grazing loss of phytoplankton:
! \begin{equation}\label{d13}
! d_{1,3} = \mu_1 \frac{c_1+c_{1}^{\min}}{K_5+c_1+c_{1}^{\min}}c_1+
! (1-\beta)\frac{g\rho_1 c_1^2}{K_3 \sum_{j=1}^3 \rho_jc_j
! + \sum_{j=1}^3 \rho_jc_j^2} (c_4+c_{4}^{\min}).
! \end{equation}
! Phytoplankton loss to LDON (labile dissolved organic nitrogen):
! \begin{equation}\label{d17}
! d_{1,7} = \gamma
! F(I_{PAR})\frac{\frac{c_5}{K_1}
! +\frac{c_6}{K_2}}{1+\frac{c_5}{K_1}+\frac{c_6}{K_2}}c_1,
! \end{equation}
! with
! \begin{equation}\label{FI}
! F(I_{PAR}) = \frac{V_p\alpha I_{PAR}(z)}{\left(V_p^2+\alpha^2(I_{PAR}(z))^2
! \right)^{1/2}}.
! \end{equation}
! With $I_{PAR}$ from (\ref{light}).