Compare View

switch
from
...
to
 
Commits (2)
nml/bio_gsj.nml 100644 → 100755
... ... @@ -64,6 +64,7 @@
64 64 ! mu4 = detritus breakdown rate [1/day]
65 65 ! mu5 = nitrification rate [1/day]
66 66 ! w_d = detritus settling velocity [m/day]
  67 +! w_h = hydrocarbon settling velocity [m/day]
67 68 ! kc = attenuation constant for the self shading effect [m**2/mmol N]
68 69 !-------------------------------------------------------------------------------
69 70 &bio_gsj_nml
... ... @@ -111,6 +112,8 @@
111 112 k3 = 1.0
112 113 beta = 0.625
113 114 mu21 = 0.3
  115 + mu21max = 1.5
  116 + bhc = 11.1
114 117 mu22 = 0.3
115 118 k6 = 0.2
116 119 delta = 0.1
... ... @@ -130,7 +133,7 @@
130 133 w_h = 100.0
131 134 mu3 = 0.03
132 135 etaa = 0.0
133   - etah = 0.0
  136 + etah = 0.2
134 137 mu4 = 0.02
135 138 mu5 = 0.00
136 139 w_d = -5.0
... ...
nml/bio_gsj_sansHCsurzoo.nml 0 → 100755
... ... @@ -0,0 +1,139 @@
  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_init = initial flagellate concentration [mmol n/m3]
  8 +! p2_init = initial diatom concentration [mmol n/m3]
  9 +! z1_init = initial micro-zooplakton concentration [mmol n/m3]
  10 +! z2_init = initial meso-zooplakton concentration [mmol n/m3]
  11 +! b_init = initial bacteria concentration [mmol n/m3]
  12 +! d_init = initial detritus concentration [mmol n/m3]
  13 +! l_init = initial LDON concentration [mmol n/m3]
  14 +! p0 = minimum phytoplankton concentration [mmol n/m3]
  15 +! z0 = minimum zooplakton concentration [mmol n/m3]
  16 +! b0 = minimum bacteria concentration [mmol n/m3]
  17 +! mte = if .true. use temperature-dependent metabolic rates
  18 +! ca1 = temp-dependence coeff for p1
  19 +! ca2 = temp-dependence coeff for p2
  20 +! ch1 = temp-dependence coeff for z1
  21 +! ch2 = temp-dependence coeff for z2
  22 +! amratio = Mass ratio between p2 and p1
  23 +! hmratio = Mass ratio between z2 and z1
  24 +! vp1 = maximum flagellate uptake rate by flagellates [1/day]
  25 +! vp2 = maximum diatom uptake rate by diatoms [1/day]
  26 +! alpha1 = slope of the flagellate PI-curve [m2/(W day)]
  27 +! alpha2 = slope of the diatom PI-curve [m2/(W day)]
  28 +! inib1 = inhibition slope of the flagellate PI-curve (pos.) [m2/(W day)]
  29 +! inib2 = inhibition slope of the PI-curve (pos.) [m2/(W day)]
  30 +! kn1 = half sat. constant nitrate uptake by p1 [mmol n/m3]
  31 +! ka1 = half sat. constant ammonium uptake by p1 [mmol n/m3]
  32 +! kn2 = half sat. constant nitrate uptake by p2 [mmol n/m3]
  33 +! ka2 = half sat. constant ammonium uptake by p2 [mmol n/m3]
  34 +! mu11 = mortality rate for p1 [1/day]
  35 +! mu12 = mortality rate for p2 [1/day]
  36 +! k5 = half sat. constant phy. mortality [mmol n/m3]
  37 +! gamma = exudation fraction [-]
  38 +! w_p1 = settling velocity for p1 [m/day]
  39 +! w_p2 = settling velocity for p2 [m/day]
  40 +! theta = phytoplancton buoyancy parameter [m3 day/(mmol N)]
  41 +! g1max = maximum ingestion rate for z1 [1/day]
  42 +! g2max = maximum ingestion rate for z2 [1/day]
  43 +! k3 = half saturation constant ingestion [mmol n/m3]
  44 +! beta = grazing efficiency [-]
  45 +! k6 = half saturation zooplankton loss (z1 & z2) [mmol n/m3]
  46 +! mu21 = maximum loss rate for z1 [1/day]
  47 +! mu22 = maximum loss rate for z2 [1/day]
  48 +! delta = fractional zooplankton loss to LDON (z1 & z2) [-]
  49 +! epsi = fractional zooplankton loss to ammonium (z1 & z2) [-]
  50 +! r11 = z1 grazing preference on p1 [-]
  51 +! r12 = z1 grazing preference on p2 [-]
  52 +! r13 = z1 grazing preference on bacteria [-]
  53 +! r14 = z1 grazing preference on detritus [-]
  54 +! r21 = z2 grazing preference on p1 [-]
  55 +! r22 = z2 grazing preference on p2 [-]
  56 +! r23 = z2 grazing preference on detritus [-]
  57 +! r24 = z2 grazing preference on z1 [-]
  58 +! vb1 = maximum bac1 uptake rate [1/day]
  59 +! vb2 = maximum bac2 uptake rate [1/day]
  60 +! k4 = half saturation bacterial uptake [mmol n/m3]
  61 +! mu3 = bacteria excretion rate [1/day]
  62 +! etaa = uptake ratio ammonium:LDON [-]
  63 +! etah = uptake ratio hydrocarbon:LDON [-]
  64 +! mu4 = detritus breakdown rate [1/day]
  65 +! mu5 = nitrification rate [1/day]
  66 +! w_d = detritus settling velocity [m/day]
  67 +! w_h = hydrocarbon settling velocity [m/day]
  68 +! kc = attenuation constant for the self shading effect [m**2/mmol N]
  69 +!-------------------------------------------------------------------------------
  70 + &bio_gsj_nml
  71 + numc = 11
  72 + p1_init = 0.012
  73 + p2_init = 0.012
  74 + z1_init = 0.012
  75 + z2_init = 0.012
  76 + b_init = 0.001
  77 + d_init = 0.01
  78 + l_init = 0.1
  79 + p0 = 0.0001
  80 + z0 = 0.0001
  81 + b0 = 0.0001
  82 + mte = .true.
  83 + ca1 = 3.61
  84 + ca2 = 14.58
  85 + ch1 = 3.265
  86 + ch2 = 24.923
  87 + amratio = 200
  88 + hmratio = 1000
  89 + vp1 = 0.02
  90 + vp2 = 0.8
  91 + alpha1 = 0.02
  92 + alpha2 = 0.04
  93 + inib1 = 0.0
  94 + inib2 = 0.006
  95 + kn1 = 1.0
  96 + ka1 = 0.8
  97 + kn2 = 1.0
  98 + ka2 = 0.8
  99 + mu11 = 0.05
  100 + mu12 = 0.05
  101 + k5 = 0.2
  102 + gamma = 0.05
  103 + w_p1 =-0.38
  104 + w_p2 =-0.00
  105 + theta = 0.0
  106 + w_p1min =-0.01
  107 + w_p1max =-0.10
  108 + w_p2min =-0.05
  109 + w_p2max =-0.38
  110 + g1max = 1.0
  111 + g2max = 1.0
  112 + k3 = 1.0
  113 + beta = 0.625
  114 + mu21 = 0.3
  115 + mu22 = 0.3
  116 + k6 = 0.2
  117 + delta = 0.1
  118 + epsi = 0.70
  119 + r11 = 0.55
  120 + r12 = 0.30
  121 + r13 = 0.05
  122 + r14 = 0.10
  123 + r21 = 0.50
  124 + r22 = 0.30
  125 + r23 = 0.05
  126 + r24 = 0.15
  127 + vb1 = 0.24
  128 + vb2 = 0.24
  129 + k4 = 0.5
  130 + k10 = 0.5
  131 + w_h = 100.0
  132 + mu3 = 0.03
  133 + etaa = 0.0
  134 + etah = 0.2
  135 + mu4 = 0.02
  136 + mu5 = 0.00
  137 + w_d = -5.0
  138 + kc = 0.03
  139 + /
... ...
src/Rules.make 0 → 100755
... ... @@ -0,0 +1,157 @@
  1 +#$iD: Rules.make,v 1.17 2006-10-26 13:12:46 kbk Exp $
  2 +
  3 +SHELL = /bin/sh
  4 +
  5 +# The compilation mode is obtained from $COMPILATION_MODE -
  6 +# default production - else debug or profiling
  7 +ifndef COMPILATION_MODE
  8 +compilation=production
  9 +else
  10 +compilation=$(COMPILATION_MODE)
  11 +endif
  12 +
  13 +# Force this here. Could be done in bashrc.
  14 +#FORTRAN_COMPILER=PGF90
  15 +#FORTRAN_COMPILER=GFORTRAN
  16 +FORTRAN_COMPILER=IFORT
  17 +
  18 +DEFINES=-DNUDGE_VEL
  19 +DEFINES=-D$(FORTRAN_COMPILER)
  20 +
  21 +# What do we include in this compilation
  22 +NetCDF = true
  23 +SEDIMENT = false
  24 +SEAGRASS = false
  25 +BIO = true
  26 +
  27 +FEATURES =
  28 +FEATURE_LIBS =
  29 +EXTRA_LIBS =
  30 +INCDIRS =
  31 +LDFLAGS =
  32 +
  33 +# If we want NetCDF - where are the include files and the library
  34 +NETCDFINC = /usr/local/include
  35 +NETCDFLIBDIR = /usr/local/lib
  36 +
  37 +ifdef NETCDFINC
  38 +INCDIRS += -I$(NETCDFINC)
  39 +endif
  40 +ifdef NETCDFLIBNAME
  41 +NETCDFLIB = $(NETCDFLIBNAME)
  42 +else
  43 +NETCDFLIB = -lnetcdf
  44 +endif
  45 +ifdef NETCDFLIBDIR
  46 +LDFLAGS += -L$(NETCDFLIBDIR)
  47 +endif
  48 +
  49 +# phony targets
  50 +.PHONY: clean realclean distclean dummy
  51 +
  52 +# Top of this version of GOTM.
  53 +ifndef GOTMDIR
  54 +GOTMDIR := $(HOME)/git/gotm_ismer
  55 +endif
  56 +
  57 +CPP = /usr/bin/cpp
  58 +
  59 +# Here you can put defines for the [c|f]pp - some will also be set depending
  60 +# on compilation mode.
  61 +ifeq ($(NetCDF),true)
  62 +DEFINES += -DNETCDF_FMT
  63 +EXTRA_LIBS += $(NETCDFLIB)
  64 +endif
  65 +ifeq ($(SEDIMENT),true)
  66 +DEFINES += -DSEDIMENT
  67 +FEATURES += extras/sediment
  68 +FEATURE_LIBS += -lsediment$(buildtype)
  69 +endif
  70 +ifeq ($(SEAGRASS),true)
  71 +DEFINES += -DSEAGRASS
  72 +FEATURES += extras/seagrass
  73 +FEATURE_LIBS += -lseagrass$(buildtype)
  74 +endif
  75 +ifeq ($(BIO),true)
  76 +DEFINES += -DBIO
  77 +FEATURES += extras/bio
  78 +FEATURE_LIBS += -lbio$(buildtype)
  79 +endif
  80 +
  81 +# Directory related settings.
  82 +
  83 +ifndef BINDIR
  84 +BINDIR = $(GOTMDIR)/bin
  85 +endif
  86 +
  87 +ifndef LIBDIR
  88 +LIBDIR += $(GOTMDIR)/lib/$(FORTRAN_COMPILER)
  89 +endif
  90 +
  91 +ifndef MODDIR
  92 +#MODDIR = $(GOTMDIR)/modules
  93 +MODDIR = $(GOTMDIR)/modules/$(FORTRAN_COMPILER)
  94 +endif
  95 +INCDIRS += -I/usr/local/include -I$(GOTMDIR)/include -I$(MODDIR)
  96 +
  97 +# Normaly this should not be changed - unless you want something very specific.
  98 +
  99 +# The Fortran compiler is determined from the EV FORTRAN_COMPILER - options
  100 +# sofar NAG(linux), FUJITSU(Linux), DECF90 (OSF1 and likely Linux on alpha),
  101 +# SunOS, PGF90 - Portland Group Fortran Compiler (on Intel Linux).
  102 +
  103 +# Sets options for debug compilation
  104 +ifeq ($(compilation),debug)
  105 +buildtype = _debug
  106 +DEFINES += -DDEBUG $(STATIC)
  107 +FLAGS = $(DEBUG_FLAGS)
  108 +endif
  109 +
  110 +# Sets options for profiling compilation
  111 +ifeq ($(compilation),profiling)
  112 +buildtype = _prof
  113 +DEFINES += -DPROFILING $(STATIC)
  114 +FLAGS = $(PROF_FLAGS)
  115 +endif
  116 +
  117 +# Sets options for production compilation
  118 +ifeq ($(compilation),production)
  119 +buildtype = _prod
  120 +DEFINES += -DPRODUCTION $(STATIC)
  121 +FLAGS = $(PROD_FLAGS)
  122 +endif
  123 +
  124 +include $(GOTMDIR)/compilers/compiler.$(FORTRAN_COMPILER)
  125 +
  126 +#DEFINES += -DREAL_4B=$(REAL_4B)
  127 +#ifeq ($(FORTRAN_COMPILER),XLF)
  128 +#DEFINES:=-WF,"$(DEFINES)"
  129 +#DEFINES=-WF,"-DXLF -DNETCDF_FMT -DSEAGRASS -DFORTRAN95 -DPRODUCTION -DREAL_4B=real\(4\)"
  130 +#endif
  131 +
  132 +# For making the source code documentation.
  133 +PROTEX = protex -b -n -s
  134 +
  135 +.SUFFIXES:
  136 +.SUFFIXES: .F90
  137 +
  138 +LINKDIR = -L$(LIBDIR)
  139 +
  140 +CPPFLAGS = $(DEFINES) $(INCDIRS)
  141 +FFLAGS = $(DEFINES) $(FLAGS) $(MODULES) $(INCDIRS) $(EXTRAS)
  142 +F90FLAGS = $(FFLAGS)
  143 +LDFLAGS += $(FFLAGS) $(LINKDIR)
  144 +
  145 +#
  146 +# Common rules
  147 +#
  148 +ifeq ($(can_do_F90),true)
  149 +%.o: %.F90
  150 + $(FC) $(F90FLAGS) $(EXTRA_FFLAGS) -c $< -o $@
  151 +else
  152 +%.f90: %.F90
  153 +# $(CPP) $(CPPFLAGS) $< -o $@
  154 + $(F90_to_f90)
  155 +%.o: %.f90
  156 + $(FC) $(F90FLAGS) $(EXTRA_FFLAGS) -c $< -o $@
  157 +endif
... ...
src/extras/bio/bio_gsj.F90 100644 → 100755
... ... @@ -97,8 +97,13 @@
97 97 REALTYPE :: k3 = 1.0
98 98 REALTYPE :: beta = 0.625
99 99 REALTYPE :: mu21 = 0.3
100   - REALTYPE :: mu22 = 0.3
  100 + REALTYPE :: mu21max = 1.5
101 101 REALTYPE :: k6 = 0.2
  102 + REALTYPE :: k6hc = 31.4
  103 + REALTYPE :: bhc = 11.1
  104 + REALTYPE :: mhc = 5
  105 + REALTYPE :: mu22 = 0.3
  106 + REALTYPE :: mu22max = 1.5
102 107 REALTYPE :: delta = 0.1
103 108 REALTYPE :: epsi = 0.70
104 109 REALTYPE :: r11 = 0.55
... ... @@ -124,6 +129,79 @@
124 129 integer :: out_unit
125 130 integer, parameter :: n=1,p1=2,p2=3,z1=4,z2=5,d=6,l=7,b1=8,a=9,hc=10,b2=11
126 131 !EOP
  132 +!! from a namelist
  133 +! p1_init : Small phytoplankton initial concentration
  134 +! p2_init : Large phytoplankton initial concentration
  135 +! z1_init : Small zooplankton initial concentration
  136 +! z2_init : Large zooplankton initial concentration
  137 +! b_init : Bacteria (Bac1) and OHC Bacteria (Bac2) initial concentration
  138 +! d_init : Detritus initial concentration
  139 +! l_init : Lablie Dissolved Organic Nitrogen (LDON) initial concentration
  140 +! p0 : Phytoplankton minimal concentration
  141 +! z0 : Zooplankton minimal concentration
  142 +! b0 : Bacteria minimal concentration
  143 +! mte : Temperature effect on phytoplankton and zooplankton loop
  144 +! hcz : Hydrocarbon mortality rate for large and small zooplankton
  145 +! ca1 :
  146 +! ca2 :
  147 +! ch1 :
  148 +! ch2 :
  149 +! amratio :
  150 +! hmratio :
  151 +! vp1 : Maximum uptake rate of small phytoplankton
  152 +! alpha1 : Slope of the small phytoplankton PI-curve
  153 +! inib1 : Inhibition slope of the small phytoplankton
  154 +! vp2 : Maximum uptake rate of large phytoplankton
  155 +! alpha2 : Slope of the large phytoplankton PI-curve
  156 +! inib2 : Inhibition slope of the large phytoplankton
  157 +! theta : Phytoplankton buoyancy parameter
  158 +! w_p1min : Maximum small phytoplankton settling velocity
  159 +! w_p1max : Minmimum small phytoplankton settling velocity
  160 +! w_p2min : Maximum large phytoplankton settling velocity
  161 +! w_p2max : Minmimum large phytoplankton settling velocity
  162 +! kn1 : Half saturation constant of nitrate uptake by small phytoplankton
  163 +! ka1 : Half saturation constant of ammonium uptake by small phytoplankton
  164 +! kn2 : Half saturation constant of nitrate uptake by large phytoplankton
  165 +! ka2 : Half saturation constant of ammonium uptake by large phytoplankton
  166 +! mu11 : Small phytoplankton mortality
  167 +! mu12 : Large phytoplankton mortality
  168 +! k5 : Half saturation of small and large phytoplankton mortality
  169 +! gamma : Excsudation fraction
  170 +! w_p1 : Small phytoplankton settling velocity
  171 +! w_p2 : Large phytoplankton settling velocity
  172 +! g1max : Maximum small zooplankton ingestion
  173 +! g2max : Maximum large zooplankton ingestion
  174 +! k3 : Half saturation constant of ingestion
  175 +! beta : Grazing efficiency
  176 +! mu21 : Maximum small zooplankton loss rate (variable function of hydrocarbon concentration)
  177 +! mu22 : Maximum large zooplankton loss rate
  178 +! k6 : Half saturation constant of zooplankton loss (small and large)
  179 +!! k7 : Half saturation constant of zooplankton loss by hydrocarbon toxicity (small and large)
  180 +! bhc : Slope factor of zooplankton mortality for sigmoïd exponential type curve
  181 +! mhc : Factor of zooplankton mortality for Michaelis-Menten exponential type curve
  182 +! hc0 : Hydrocarbon sill concentration for zooplankton mortality by hydrocarbons
  183 +! delta : Fractional zooplankton loss of LDON
  184 +! epsi : Fractional zooplankton loss of ammonuim
  185 +! r11 : Small zooplankton preference on small phytoplantkon
  186 +! r12 : Small zooplankton preference on small phytoplantkon
  187 +! r13 : Small zooplankton preference on small phytoplantkon
  188 +! r14 : Small zooplankton preference on small phytoplantkon
  189 +! r21 : Small zooplankton preference on small phytoplantkon
  190 +! r22 : Small zooplankton preference on small phytoplantkon
  191 +! r23 : Small zooplankton preference on small phytoplantkon
  192 +! r24 : Small zooplankton preference on small phytoplantkon
  193 +! vb1 : Maximum bacteria (b1) uptake rate
  194 +! vb2 : Maximum hydrocarbon degrading bacteria (b2) uptake rate
  195 +! k4 : Half saturation of bacteria (b1) uptake
  196 +! k10 : Half saturation of hydrocarbon degrading bacteria (b2) uptake
  197 +! w_h : Hydrocarbon settling velocity
  198 +! mu3 : Bacterial excretion rate
  199 +! etaa : Uptake ratio ammonium : LDON
  200 +! etah : Uptake ratio hydrocarbon : LDON
  201 +! mu4 : Detritus breakdown rate
  202 +! mu5 : Nitrification rate
  203 +! w_d : Detritus settling velocity
  204 +! kc :
127 205 !-----------------------------------------------------------------------
128 206  
129 207 contains
... ... @@ -159,7 +237,7 @@
159 237 b_init,d_init,l_init, &
160 238 p0,z0,b0,vp1,alpha1,inib1,vp2,alpha2,inib2, &
161 239 kn1,ka1,kn2,ka2,mu11,mu12,k5,gamma,w_p1,w_p2, &
162   - g1max,g2max,k3,beta,mu21,mu22,k6,delta,epsi, &
  240 + g1max,g2max,k3,beta,mu21,mu22,k6,bhc,delta,epsi, &
163 241 r11,r12,r13,r14,r21,r22,r23,r24, &
164 242 vb1,vb2,k4,k10,w_h,mu3,etaa,etah,mu4,w_d,kc,mu5, &
165 243 theta,w_p1max,w_p1min,w_p2min,w_p2max, &
... ... @@ -201,9 +279,8 @@
201 279 write(10,901) g2max
202 280 write(10,901) k3
203 281 write(10,901) beta
204   - write(10,901) mu21
205 282 write(10,901) mu22
206   - write(10,901) k6
  283 + write(10,901) bhc
207 284 write(10,901) delta
208 285 write(10,901) epsi
209 286 write(10,901) r11
... ... @@ -588,6 +665,7 @@
588 665 REALTYPE :: minal,minhl,qn1,qa1,qn2,qa2
589 666 REALTYPE :: ps1,ps2,ff1,ff2
590 667 REALTYPE :: Ea,Eh,kBeV,T0
  668 + REALTYPE :: k6hc,mu21hc,mu22hc
591 669 integer :: i,j,ci
592 670 !EOP
593 671 !-----------------------------------------------------------------------
... ... @@ -637,6 +715,10 @@
637 715 hmr2 = 1.0
638 716 endif
639 717  
  718 +!-----------------------------------------------------------------------
  719 +
  720 + ! Light limitation factor (PI curve)
  721 +
640 722 ! Smith (1936) - saturation (default)
641 723 ! ff= vp*alpha*par(ci)/sqrt(vp**2+alpha**2*par(ci)**2)
642 724 ! Blackman (1919)
... ... @@ -651,7 +733,6 @@
651 733 ! ff= vp*((par(ci)/I_opt)*exp(1-(par(ci)/I_opt)))**2
652 734  
653 735 ! Platt et al. (1980) - inhibition
654   - ! Light limitation factor (PI curve)
655 736 ps1= vp1/((alpha1/(alpha1+inib1))*(alpha1/(alpha1+inib1))**(inib1/alpha1))
656 737 ff1= ps1*(1.-exp(-1.*alpha1*par(ci)/ps1))*exp(-1.*inib1*par(ci)/ps1)
657 738  
... ... @@ -675,7 +756,7 @@
675 756 +r23*cc(d,ci)**2+r24*cc(z1,ci)**2)
676 757  
677 758 minal=min(cc(a,ci),etaa*cc(l,ci))
678   - minhl=min(cc(hc,ci),etah*cc(l,ci))
  759 + minhl=min(cc(hc,ci)/(k10+cc(hc,ci)),cc(l,ci)/(k4+cc(l,ci)))
679 760  
680 761 ! Light and nutrient limitation factors
681 762 lumlim1(ci) =amr1*ff1
... ... @@ -684,6 +765,11 @@
684 765 lumlim2(ci) =amr2*ff2
685 766 nitlim2(ci) =qn2
686 767 ammlim2(ci) =qa2
  768 +
  769 + ! Hydrocarbon inhibition factors
  770 + !mu21hc = mu21+(mu21max-mu21)/(1d0 + exp(-1d0*(cc(hc,ci)-k6hc)/bhc))
  771 + mu21hc = mu21+((mu21max-mu21)*(cc(hc,ci))**mhc)/((cc(hc,ci))**mhc + k6hc**mhc)
  772 + mu22hc = mu21hc
687 773  
688 774 ! Nutrient uptake by pico- and nano-phytoplankton
689 775 dd(n,p1,ci) =amr1*ff1*qn1*(cc(p1,ci)+p0)
... ... @@ -701,6 +787,8 @@
701 787  
702 788 dd(b1,d,ci) =(1.-beta)*g1max*r13*cc(b1,ci)**2*fac1
703 789  
  790 + ! Grazing and uptake by zooplankton
  791 +
704 792 dd(p1,z1,ci)=hmr1*beta*g1max*r11*cc(p1,ci)**2*fac1
705 793 dd(p2,z1,ci)=hmr1*beta*g1max*r12*cc(p2,ci)**2*fac1
706 794 dd(b1,z1,ci)=hmr1*beta*g1max*0.5*r13*cc(b1,ci)**2*fac1
... ... @@ -711,25 +799,32 @@
711 799 dd(p2,z2,ci)=hmr2*beta*g2max*r22*cc(p2,ci)**2*fac2
712 800 dd(d,z2,ci) =hmr2*beta*g2max*r23*cc(d,ci)**2*fac2
713 801 dd(z1,z2,ci)=hmr2*beta*g2max*r24*cc(z1,ci)**2*fac2
  802 +
  803 + ! Other processes
714 804  
715 805 dd(b1,a,ci) =mu3*cc(b1,ci)
716 806 dd(b2,a,ci) =mu3*cc(b2,ci)
717 807 dd(d,l,ci) =mu4*cc(d,ci)
718 808 dd(a,n,ci) =mu5*cc(a,ci)
719 809  
  810 + ! Zooplankton mortality and exsudation
  811 +
720 812 dd(z1,d,ci) =(1.-beta)*hmr2*g2max*r24*cc(z1,ci)**2*fac2 &
721   - +hmr1*(1.-epsi-delta)*mu21*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
722   - dd(z1,a,ci) =hmr1*epsi*mu21*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
723   - dd(z1,l,ci) =hmr1*delta*mu21*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  813 + +hmr1*(1.-epsi-delta)*mu21hc*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  814 + dd(z1,a,ci) =hmr1*epsi*mu21hc*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  815 + dd(z1,l,ci) =hmr1*delta*mu21hc*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
724 816  
725   - dd(z2,d,ci) =hmr2*(1.-epsi-delta)*mu22*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
726   - dd(z2,a,ci) =hmr2*epsi*mu22*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
727   - dd(z2,l,ci) =hmr2*delta*mu22*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  817 + dd(z2,d,ci) =hmr2*(1.-epsi-delta)*mu22hc*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  818 + dd(z2,a,ci) =hmr2*epsi*mu22hc*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  819 + dd(z2,l,ci) =hmr2*delta*mu22hc*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  820 +
  821 + ! Nutrient and Hydrocarbon uptake by bacteria and OHC bacteria
728 822  
729 823 dd(a,b1,ci) =vb1*minal/(k4+minal+cc(l,ci))*(cc(b1,ci)+b0)
730 824 dd(l,b1,ci) =vb1*cc(l,ci)/(k4+minal+cc(l,ci))*(cc(b1,ci)+b0)
731   - dd(hc,b2,ci)=vb2*minhl/(k10+minhl+cc(l,ci))*(cc(b2,ci)+b0)
732   - dd(l,b2,ci) =vb2*cc(l,ci)/(k10+minhl+cc(l,ci))*(cc(b2,ci)+b0)
  825 + dd(hc,b2,ci)=vb2*minhl*(cc(b2,ci)+b0)
  826 + dd(l,b2,ci) =vb2*etah*minhl*(cc(b2,ci)+b0)
  827 + !dd(l,b2,ci) =vb2*cc(l,ci)/(k10+minhl+cc(l,ci))*(cc(b2,ci)+b0)
733 828 !dd(hc,b2,ci)=vb*cc(hc,ci)/(k10+cc(hc,ci))*(cc(b2,ci)+b0)
734 829  
735 830  
... ...
src/extras/bio/bio_gsj_sansHCsurzoo.F90 0 → 100644
... ... @@ -0,0 +1,779 @@
  1 +!$Id: bio_ismer.F90,v 1.11 2007-01-06 11:49:15 kbk Exp $
  2 +#include"cppdefs.h"
  3 +!-----------------------------------------------------------------------
  4 +!BOP
  5 +!
  6 +! !MODULE: bio_gsj --- Modified from Fasham et al. biological model \label{sec:bio-fasham}
  7 +!
  8 +! !INTERFACE:
  9 + module bio_gsj
  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_gsj, init_var_gsj, var_info_gsj, &
  47 + light_gsj, do_bio_gsj, end_bio_gsj
  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_init = 0.05
  58 + REALTYPE :: p2_init = 0.05
  59 + REALTYPE :: z1_init = 0.05
  60 + REALTYPE :: z2_init = 0.05
  61 + REALTYPE :: b_init = 0.001
  62 + REALTYPE :: d_init = 0.4
  63 + REALTYPE :: l_init = 0.14
  64 + REALTYPE :: p0 = 0.0
  65 + REALTYPE :: z0 = 0.0
  66 + REALTYPE :: b0 = 0.0
  67 + LOGICAL :: mte = .true.
  68 + REALTYPE :: ca1 = 3.61
  69 + REALTYPE :: ca2 = 14.58
  70 + REALTYPE :: ch1 = 3.265
  71 + REALTYPE :: ch2 = 24.923
  72 + REALTYPE :: amratio = 200.0
  73 + REALTYPE :: hmratio = 1000.0
  74 + REALTYPE :: vp1 = 1.5
  75 + REALTYPE :: alpha1 = 0.065
  76 + REALTYPE :: inib1 = 0.05
  77 + REALTYPE :: vp2 = 1.5
  78 + REALTYPE :: alpha2 = 0.065
  79 + REALTYPE :: inib2 = 0.05
  80 + REALTYPE :: theta = 0.0
  81 + REALTYPE :: w_p1min = -0.06
  82 + REALTYPE :: w_p1max = -0.38
  83 + REALTYPE :: w_p2min = -0.06
  84 + REALTYPE :: w_p2max = -0.38
  85 + REALTYPE :: kn1 = 0.2
  86 + REALTYPE :: ka1 = 0.8
  87 + REALTYPE :: kn2 = 0.2
  88 + REALTYPE :: ka2 = 0.8
  89 + REALTYPE :: mu11 = 0.05
  90 + REALTYPE :: mu12 = 0.05
  91 + REALTYPE :: k5 = 0.2
  92 + REALTYPE :: gamma = 0.05
  93 + REALTYPE :: w_p1 = -0.5
  94 + REALTYPE :: w_p2 = -0.5
  95 + REALTYPE :: g1max = 1.0
  96 + REALTYPE :: g2max = 1.0
  97 + REALTYPE :: k3 = 1.0
  98 + REALTYPE :: beta = 0.625
  99 + REALTYPE :: mu21 = 0.3
  100 + REALTYPE :: mu22 = 0.3
  101 + REALTYPE :: k6 = 0.2
  102 + REALTYPE :: delta = 0.1
  103 + REALTYPE :: epsi = 0.70
  104 + REALTYPE :: r11 = 0.55
  105 + REALTYPE :: r12 = 0.30
  106 + REALTYPE :: r13 = 0.05
  107 + REALTYPE :: r14 = 0.10
  108 + REALTYPE :: r21 = 0.50
  109 + REALTYPE :: r22 = 0.30
  110 + REALTYPE :: r23 = 0.05
  111 + REALTYPE :: r24 = 0.15
  112 + REALTYPE :: vb1 = 1.2
  113 + REALTYPE :: vb2 = 1.2
  114 + REALTYPE :: k4 = 0.5
  115 + REALTYPE :: k10 = 0.15
  116 + REALTYPE :: w_h = 0.0
  117 + REALTYPE :: mu3 = 0.15
  118 + REALTYPE :: etaa = 0.0
  119 + REALTYPE :: etah = 0.0
  120 + REALTYPE :: mu4 = 0.02
  121 + REALTYPE :: mu5 = 0.00
  122 + REALTYPE :: w_d = -2.0
  123 + REALTYPE, public :: kc = 0.03
  124 + integer :: out_unit
  125 + integer, parameter :: n=1,p1=2,p2=3,z1=4,z2=5,d=6,l=7,b1=8,a=9,hc=10,b2=11
  126 +!EOP
  127 +!-----------------------------------------------------------------------
  128 +
  129 + contains
  130 +
  131 +!-----------------------------------------------------------------------
  132 +!BOP
  133 +!
  134 +! !IROUTINE: Initialise the bio module
  135 +!
  136 +! !INTERFACE:
  137 + subroutine init_bio_gsj(namlst,fname,unit)
  138 +!
  139 +! !DESCRIPTION:
  140 +! Here, the bio namelist {\tt bio\_fasham.nml} is read and
  141 +! various variables (rates and settling velocities)
  142 +! are transformed into SI units.
  143 +!
  144 +! !USES:
  145 + IMPLICIT NONE
  146 +!
  147 +! !INPUT PARAMETERS:
  148 + integer, intent(in) :: namlst
  149 + character(len=*), intent(in) :: fname
  150 + character(len=20) :: pfile
  151 + integer, intent(in) :: unit
  152 +!
  153 +! !REVISION HISTORY:
  154 +! Original author(s): Hans Burchard & Karsten Bolding
  155 +!
  156 +! !LOCAL VARIABLES:
  157 + namelist /bio_gsj_nml/ numc, &
  158 + p1_init,p2_init,z1_init,z2_init, &
  159 + b_init,d_init,l_init, &
  160 + p0,z0,b0,vp1,alpha1,inib1,vp2,alpha2,inib2, &
  161 + kn1,ka1,kn2,ka2,mu11,mu12,k5,gamma,w_p1,w_p2, &
  162 + g1max,g2max,k3,beta,mu21,mu22,k6,delta,epsi, &
  163 + r11,r12,r13,r14,r21,r22,r23,r24, &
  164 + vb1,vb2,k4,k10,w_h,mu3,etaa,etah,mu4,w_d,kc,mu5, &
  165 + theta,w_p1max,w_p1min,w_p2min,w_p2max, &
  166 + mte,ca1,ca2,ch1,ch2,amratio,hmratio
  167 +
  168 +!EOP
  169 +!-----------------------------------------------------------------------
  170 +!BOC
  171 + LEVEL2 'init_bio_gsj'
  172 +
  173 + open(namlst,file=fname,action='read',status='old',err=98)
  174 + read(namlst,nml=bio_gsj_nml,err=99)
  175 + close(namlst)
  176 +
  177 + numcc=numc
  178 +
  179 +! Print some parameter values in standard output
  180 +! and save them in a separate file [out_fn]_ismer.par
  181 + pfile = trim(out_fn) // '_gsj.par'
  182 + open(10,status='unknown',action='write',file=pfile)
  183 + LEVEL3 'GSJ parameters saved in ', pfile
  184 + write(10,901) vp1
  185 + write(10,901) vp2
  186 + write(10,901) alpha1
  187 + write(10,901) alpha2
  188 + write(10,901) inib1
  189 + write(10,901) inib2
  190 + write(10,901) kn1
  191 + write(10,901) ka1
  192 + write(10,901) kn2
  193 + write(10,901) ka2
  194 + write(10,901) mu11
  195 + write(10,901) mu12
  196 + write(10,901) k5
  197 + write(10,901) gamma
  198 + write(10,901) w_p1
  199 + write(10,901) w_p2
  200 + write(10,901) g1max
  201 + write(10,901) g2max
  202 + write(10,901) k3
  203 + write(10,901) beta
  204 + write(10,901) mu21
  205 + write(10,901) mu22
  206 + write(10,901) k6
  207 + write(10,901) delta
  208 + write(10,901) epsi
  209 + write(10,901) r11
  210 + write(10,901) r12
  211 + write(10,901) r13
  212 + write(10,901) r14
  213 + write(10,901) r21
  214 + write(10,901) r22
  215 + write(10,901) r23
  216 + write(10,901) r24
  217 + write(10,901) vb1
  218 + write(10,901) vb2
  219 + write(10,901) k4
  220 + write(10,901) mu3
  221 + write(10,901) etaa
  222 + write(10,901) etah
  223 + write(10,901) mu4
  224 + write(10,901) mu5
  225 + write(10,901) w_d
  226 + write(10,901) kc
  227 + if (mte) then
  228 + write(10,901) ca1
  229 + write(10,901) ca2
  230 + write(10,901) ch1
  231 + write(10,901) ch2
  232 + write(10,901) amratio
  233 + write(10,901) hmratio
  234 + endif
  235 +
  236 +
  237 +901 format (f8.5)
  238 +
  239 +! Conversion from day to second
  240 + vp1 = vp1 /secs_pr_day
  241 + vp2 = vp2 /secs_pr_day
  242 + vb1 = vb1 /secs_pr_day
  243 + vb2 = vb2 /secs_pr_day
  244 + mu11 = mu11 /secs_pr_day
  245 + mu12 = mu12 /secs_pr_day
  246 + mu21 = mu21 /secs_pr_day
  247 + mu22 = mu22 /secs_pr_day
  248 + mu3 = mu3 /secs_pr_day
  249 + mu4 = mu4 /secs_pr_day
  250 + mu5 = mu5 /secs_pr_day
  251 + g1max = g1max /secs_pr_day
  252 + g2max = g2max /secs_pr_day
  253 + w_p1 = w_p1 /secs_pr_day
  254 + w_p2 = w_p2 /secs_pr_day
  255 + w_p1min = w_p1min /secs_pr_day
  256 + w_p1max = w_p1max /secs_pr_day
  257 + w_p2min = w_p2min /secs_pr_day
  258 + w_p2max = w_p2max /secs_pr_day
  259 + theta = theta /secs_pr_day
  260 + w_d = w_d /secs_pr_day
  261 + w_h = w_h /secs_pr_day
  262 + alpha1 = alpha1 /secs_pr_day
  263 + inib1 = inib1 /secs_pr_day
  264 + alpha2 = alpha2 /secs_pr_day
  265 + inib2 = inib2 /secs_pr_day
  266 +
  267 + out_unit=unit
  268 +
  269 + LEVEL3 'GSJ bio module initialised ...'
  270 +
  271 + return
  272 +
  273 +98 LEVEL2 'I could not open bio_gsj.nml'
  274 + LEVEL2 'If thats not what you want you have to supply bio_gsj.nml'
  275 + LEVEL2 'See the bio example on www.gotm.net for a working bio_gsj.nml'
  276 + return
  277 +99 FATAL 'I could not read bio_gsj.nml'
  278 + stop 'init_bio_gsj'
  279 + end subroutine init_bio_gsj
  280 +!EOC
  281 +
  282 +!-----------------------------------------------------------------------
  283 +!BOP
  284 +!
  285 +! !IROUTINE: Initialise the concentration variables
  286 +!
  287 +! !INTERFACE:
  288 + subroutine init_var_gsj(nlev)
  289 +!
  290 +! !DESCRIPTION:
  291 +! Here, the the initial conditions are set and the settling velocities are
  292 +! transferred to all vertical levels. All concentrations are declared
  293 +! as non-negative variables, and it is defined which variables would be
  294 +! taken up by benthic filter feeders.
  295 +!
  296 +! !USES:
  297 + use observations, only: nprof,aprof,hcprof
  298 + use meanflow, only: nit,amm,hcb,T,S
  299 +
  300 + IMPLICIT NONE
  301 +
  302 +!
  303 +! !INPUT PARAMETERS:
  304 + integer, intent(in) :: nlev
  305 +!
  306 +! !REVISION HISTORY:
  307 +! Original author(s): Hans Burchard & Karsten Bolding
  308 +
  309 +! !LOCAL VARIABLES:
  310 + integer :: i
  311 +!EOP
  312 +!-----------------------------------------------------------------------
  313 +!BOC
  314 + do i=1,nlev
  315 + cc(n,i) = nprof(i) !CHG3
  316 + cc(p1,i)= p1_init
  317 + cc(p2,i)= p2_init
  318 + cc(z1,i)= z1_init
  319 + cc(z2,i)= z2_init
  320 + cc(d,i) = d_init
  321 + cc(l,i) = l_init
  322 + cc(b1,i)= b_init
  323 + cc(a,i) = aprof(i) !CHG5
  324 + cc(hc,i)= hcprof(i)
  325 + cc(b2,i)= b_init
  326 + end do
  327 +
  328 + do i=0,nlev
  329 + ws(n,i) = _ZERO_
  330 + ws(p1,i) = w_p1
  331 + ws(p2,i) = w_p2
  332 + ws(z1,i) = _ZERO_
  333 + ws(z2,i) = _ZERO_
  334 + ws(d,i) = w_d
  335 + ws(l,i) = _ZERO_
  336 + ws(b1,i) = _ZERO_
  337 + ws(a,i) = _ZERO_
  338 + ws(hc,i) = w_h
  339 + ws(b2,i) = _ZERO_
  340 + end do
  341 +
  342 + posconc(n) = 1
  343 + posconc(p1) = 1
  344 + posconc(p2) = 1
  345 + posconc(z1) = 1
  346 + posconc(z2) = 1
  347 + posconc(d) = 1
  348 + posconc(l) = 1
  349 + posconc(b1) = 1
  350 + posconc(a) = 1
  351 + posconc(hc) = 1
  352 + posconc(b2) = 1
  353 +
  354 + LEVEL3 'GSJ variables initialised ...'
  355 +
  356 + return
  357 +
  358 + end subroutine init_var_gsj
  359 +!EOC
  360 +
  361 +!-----------------------------------------------------------------------
  362 +!BOP
  363 +!
  364 +! !IROUTINE: Providing info on variables
  365 +!
  366 +! !INTERFACE:
  367 + subroutine var_info_gsj()
  368 +!
  369 +! !DESCRIPTION:
  370 +! This subroutine provides information about the variable names as they
  371 +! will be used when storing data in NetCDF files.
  372 +!
  373 +! !USES:
  374 + IMPLICIT NONE
  375 +!
  376 +! !REVISION HISTORY:
  377 +! Original author(s): Hans Burchard & Karsten Bolding
  378 +!
  379 +! !LOCAL VARIABLES:
  380 +!EOP
  381 +!-----------------------------------------------------------------------
  382 +!BOC
  383 + var_names(1) = 'nit'
  384 + var_units(1) = 'mmol/m**3'
  385 + var_long(1) = 'nitrate'
  386 +
  387 + var_names(2) = 'fla'
  388 + var_units(2) = 'mmol/m**3'
  389 + var_long(2) = 'flagellates'
  390 +
  391 + var_names(3) = 'dia'
  392 + var_units(3) = 'mmol/m**3'
  393 + var_long(3) = 'diatoms'
  394 +
  395 + var_names(4) = 'mcz'
  396 + var_units(4) = 'mmol/m**3'
  397 + var_long(4) = 'micro-zooplankton'
  398 +
  399 + var_names(5) = 'msz'
  400 + var_units(5) = 'mmol/m**3'
  401 + var_long(5) = 'meso-zooplankton'
  402 +
  403 + var_names(6) = 'det'
  404 + var_units(6) = 'mmol/m**3'
  405 + var_long(6) = 'detritus'
  406 +
  407 + var_names(7) = 'ldn'
  408 + var_units(7) = 'mmol/m**3'
  409 + var_long(7) = 'labile_dissolved_organic_nitrogen'
  410 +
  411 + var_names(8) = 'bac1'
  412 + var_units(8) = 'mmol/m**3'
  413 + var_long(8) = 'bacteria'
  414 +
  415 + var_names(9) = 'amm'
  416 + var_units(9) = 'mmol/m**3'
  417 + var_long(9) = 'ammonium'
  418 +
  419 + var_names(10)= 'hcb'
  420 + var_units(10)= 'mmol/m**3'
  421 + var_long(10) = 'hydrocarbon'
  422 +
  423 + var_names(11)= 'bac2'
  424 + var_units(11)= 'mmol/m**3'
  425 + var_long(11) = 'bacteria'
  426 +
  427 + return
  428 + end subroutine var_info_gsj
  429 +!EOC
  430 +
  431 +!-----------------------------------------------------------------------
  432 +!BOP
  433 +!
  434 +! !IROUTINE: Light properties for the ISMER model
  435 +!
  436 +! !INTERFACE:
  437 + subroutine light_gsj(nlev,bioshade_feedback)
  438 +!
  439 +! !DESCRIPTION:
  440 +! Here, the photosynthetically available radiation is calculated
  441 +! by simply assuming that the short wave part of the total
  442 +! radiation is available for photosynthesis.
  443 +! The photosynthetically
  444 +! available radiation, $I_{PAR}$, follows from (\ref{light}).
  445 +! The user should make
  446 +! sure that this is consistent with the light class given in the
  447 +! {\tt extinct} namelist of the {\tt obs.nml} file.
  448 +! The self-shading effect is also calculated in this subroutine,
  449 +! which may be used to consider the effect of bio-turbidity also
  450 +! in the temperature equation (if {\tt bioshade\_feedback} is set
  451 +! to true in {\tt bio.nml}).
  452 +! For details, see section \ref{sec:do-bio}.
  453 +!
  454 +! !USES:
  455 + IMPLICIT NONE
  456 +!
  457 +! !INPUT PARAMETERS:
  458 + integer, intent(in) :: nlev
  459 + logical, intent(in) :: bioshade_feedback
  460 +!
  461 +! !REVISION HISTORY:
  462 +! Original author(s): Hans Burchard, Karsten Bolding
  463 +!
  464 +! !LOCAL VARIABLES:
  465 + integer :: i
  466 + REALTYPE :: zz,add
  467 +!EOP
  468 +!-----------------------------------------------------------------------
  469 +!BOC
  470 + zz = _ZERO_
  471 + add = _ZERO_
  472 + do i=nlev,1,-1
  473 + add=add+0.5*h(i)*(cc(p1,i)+cc(p2,i)+2*p0)
  474 + zz=zz+0.5*h(i)
  475 + par(i)=rad(nlev)*(1.-aa)*exp(-zz/g2)*exp(-kc*add)
  476 + add=add+0.5*h(i)*(cc(p1,i)+cc(p2,i)+2*p0)
  477 + zz=zz+0.5*h(i)
  478 + if (bioshade_feedback) bioshade_(i)=exp(-kc*add)
  479 + end do
  480 +
  481 +
  482 + return
  483 + end subroutine light_gsj
  484 +!EOC
  485 +
  486 +!-----------------------------------------------------------------------
  487 +!BOP
  488 +!
  489 +! !IROUTINE: Right hand sides of geobiochemical model \label{sec:bio-fasham-rhs}
  490 +!
  491 +! !INTERFACE:
  492 + subroutine do_bio_gsj(first,numc,nlev,cc,pp,dd)
  493 +!
  494 +! !DESCRIPTION:
  495 +!
  496 +! The \cite{Fashametal1990} model consisting of the $I=7$
  497 +! state variables phytoplankton, bacteria, detritus, zooplankton,
  498 +! nitrate, ammonium and dissolved organic nitrogen is described here
  499 +! in detail.
  500 +!
  501 +! Phytoplankton mortality and zooplankton grazing loss of phytoplankton:
  502 +! \begin{equation}\label{d13}
  503 +! d_{1,3} = \mu_1 \frac{c_1+c_{1}^{\min}}{K_5+c_1+c_{1}^{\min}}c_1+
  504 +! (1-\beta)\frac{g\rho_1 c_1^2}{K_3 \sum_{j=1}^3 \rho_jc_j
  505 +! + \sum_{j=1}^3 \rho_jc_j^2} (c_4+c_{4}^{\min}).
  506 +! \end{equation}
  507 +! Phytoplankton loss to LDON (labile dissolved organic nitrogen):
  508 +! \begin{equation}\label{d17}
  509 +! d_{1,7} = \gamma
  510 +! F(I_{PAR})\frac{\frac{c_5}{K_1}
  511 +! +\frac{c_6}{K_2}}{1+\frac{c_5}{K_1}+\frac{c_6}{K_2}}c_1,
  512 +! \end{equation}
  513 +! with
  514 +! \begin{equation}\label{FI}
  515 +! F(I_{PAR}) = \frac{V_p\alpha I_{PAR}(z)}{\left(V_p^2+\alpha^2(I_{PAR}(z))^2
  516 +! \right)^{1/2}}.
  517 +! \end{equation}
  518 +! With $I_{PAR}$ from (\ref{light}).
  519 +!
  520 +! Zooplankton grazing loss:
  521 +! \begin{equation}\label{di3}
  522 +! d_{2,3} = (1-\beta)\frac{g\rho_2 c_2^2}{K_3 \sum_{j=1}^3 \rho_jc_j
  523 +! + \sum_{j=1}^3 \rho_jc_j^2} (c_4+c_{4}^{\min}).
  524 +! \end{equation}
  525 +! Zooplankton grazing:
  526 +! \begin{equation}\label{di4}
  527 +! d_{i,4} = \beta\frac{g\rho_i c_i^2}{K_3 \sum_{j=1}^3 \rho_jc_j
  528 +! + \sum_{j=1}^3 \rho_jc_j^2} (c_4+c_{4}^{\min}), \quad i=1,\dots,3.
  529 +! \end{equation}
  530 +! Bacteria excretion rate:
  531 +! \begin{equation}\label{d26}
  532 +! d_{2,6} = \mu_3 c_2.
  533 +! \end{equation}
  534 +! Detritus breakdown rate:
  535 +! \begin{equation}\label{d37}
  536 +! d_{3,7} = \mu_4 c_3.
  537 +! \end{equation}
  538 +! Zooplankton losses to detritus, ammonium and LDON:
  539 +! \begin{equation}\label{d43}
  540 +! d_{4,3} = (1-\epsilon-\delta)\mu_2
  541 +! \frac{c_4+c_{4}^{\min}}{K_6+c_4+c_{4}^{\min}}c_4.
  542 +! \end{equation}
  543 +! \begin{equation}\label{d46}
  544 +! d_{4,6} = \epsilon\mu_2 \frac{c_4+c_{4}^{\min}}{K_6+c_4+c_{4}^{\min}}c_4.
  545 +! \end{equation}
  546 +! \begin{equation}\label{d47}
  547 +! d_{4,7} = \delta\mu_2 \frac{c_4+c_{4}^{\min}}{K_6+c_4+c_{4}^{\min}}c_4.
  548 +! \end{equation}
  549 +! Nitrate uptake by phytoplankton:
  550 +! \begin{equation}\label{d51}
  551 +! d_{5,1} = F(I_{PAR})\frac{\frac{c_5}{K_1}}{1+\frac{c_5}{K_1}
  552 +! +\frac{c_6}{K_2}}(c_1+c_{1}^{\min}).
  553 +! \end{equation}
  554 +! Ammonium uptake by phytoplankton:
  555 +! \begin{equation}\label{d61}
  556 +! d_{6,1} = F(I_{PAR})\frac{\frac{c_6}{K_2}}{1+\frac{c_5}{K_1}
  557 +! +\frac{c_6}{K_2}}(c_1+c_{1}^{\min}).
  558 +! \end{equation}
  559 +! Ammonium uptake by bacteria:
  560 +! \begin{equation}\label{d62}
  561 +! d_{6,2} = V_b \frac{\min(c_6,\eta c_7)}{K_4+\min(c_6,\eta c_7)+c_7}
  562 +! (c_2+c_{2}^{\min}).
  563 +! \end{equation}
  564 +! LDON uptake by bacteria:
  565 +! \begin{equation}\label{d72}
  566 +! d_{7,2} = V_b \frac{c_7}{K_4+\min(c_6,\eta c_7)+c_7} (c_2+c_{2}^{\min}).
  567 +! \end{equation}
  568 +!
  569 +! !USES:
  570 + IMPLICIT NONE
  571 +!
  572 +! !INPUT PARAMETERS:
  573 + logical, intent(in) :: first
  574 + integer, intent(in) :: numc,nlev
  575 + REALTYPE, intent(in) :: cc(1:numc,0:nlev)
  576 +!
  577 +! !OUTPUT PARAMETERS:
  578 + REALTYPE, intent(out) :: pp(1:numc,1:numc,0:nlev)
  579 + REALTYPE, intent(out) :: dd(1:numc,1:numc,0:nlev)
  580 +!
  581 +! !REVISION HISTORY:
  582 +! Original author(s): Hans Burchard, Karsten Bolding
  583 +!
  584 +! !LOCAL VARIABLES:
  585 + REALTYPE :: amr1,amr2
  586 + REALTYPE :: hmr1,hmr2
  587 + REALTYPE :: fac1,fac2,fac3,fac4
  588 + REALTYPE :: minal,minhl,qn1,qa1,qn2,qa2
  589 + REALTYPE :: ps1,ps2,ff1,ff2
  590 + REALTYPE :: Ea,Eh,kBeV,T0
  591 + integer :: i,j,ci
  592 +!EOP
  593 +!-----------------------------------------------------------------------
  594 +!BOC
  595 +!KBK - is it necessary to initialise every time - expensive in a 3D model
  596 + pp = _ZERO_
  597 + dd = _ZERO_
  598 +
  599 + do ci=1,nlev
  600 +
  601 + if (mte) then
  602 + ! Temperature-dependent metabolic rates
  603 + ! Gillooly et al. (2002)
  604 + Eh = 0.65 ! Activation energy for heterotrophs (eV)
  605 + Ea = 0.32 ! Activation energy for autotrophs (eV)
  606 + kBeV = 8.62e-5 ! Boltzmann constant (eV K-1)
  607 + T0 = 273.15-1.9 ! Temperature at which metabolism stops
  608 +
  609 + ! metabolic rate coefficient (a(T0) in Gillooly et al. 2002 in kg^1/4 s^-1)
  610 + ! they are tuned to
  611 + ca1 = 3.61 ! pico-phytoplankton
  612 + ca2 = 14.58 ! micro-phytoplankton
  613 + ch1 = 3.265 ! micro-zooplankton
  614 + ch2 = 24.923 ! meso-zooplankton
  615 +
  616 + ! The mass ratio between pico-phytoplankton (p1) and
  617 + ! nano- and micro-phytoplankton (p2) is amratio,
  618 + ! which translates into a ratio of fac3=amratio**0.25 between
  619 + ! their respective metabolic rates, according to the MTE.
  620 + ! The same applies to micro-zooplankton (z1) versus
  621 + ! meso-zooplankton (z2) through fac4=hmratio**0.25.
  622 +
  623 + fac3 = amratio**0.25
  624 + fac4 = hmratio**0.25
  625 +
  626 + ! Autotroph metabolic rate
  627 + amr1 = max(0d0,ca1*0.25*exp(Ea/(kBeV*T0**2)*(T(ci)/(1+T(ci)/T0))) )
  628 + amr2 = max(0d0,ca2*0.25*exp(Ea/(kBeV*T0**2)*(T(ci)/(1+T(ci)/T0)))/fac3)
  629 +
  630 + ! Heterotroph metabolic rate
  631 + hmr1 = max(0d0,ch1*0.25*exp(Eh/(kBeV*T0**2)*(T(ci)/(1+T(ci)/T0))) )
  632 + hmr2 = max(0d0,ch2*0.25*exp(Eh/(kBeV*T0**2)*(T(ci)/(1+T(ci)/T0)))/fac4)
  633 + else
  634 + amr1 = 1.0
  635 + amr2 = 1.0
  636 + hmr1 = 1.0
  637 + hmr2 = 1.0
  638 + endif
  639 +
  640 + ! Smith (1936) - saturation (default)
  641 + ! ff= vp*alpha*par(ci)/sqrt(vp**2+alpha**2*par(ci)**2)
  642 + ! Blackman (1919)
  643 + ! if (par(ci) .lt. vp/alpha) then
  644 + ! ff=alpha*par(ci)
  645 + ! else
  646 + ! ff=vp
  647 + ! endif
  648 + ! Steele (1962) - inhibition
  649 + ! ff= vp*((par(ci)/I_opt)*exp(1-(par(ci)/I_opt)))
  650 + ! Parker (1974) - inhibition
  651 + ! ff= vp*((par(ci)/I_opt)*exp(1-(par(ci)/I_opt)))**2
  652 +
  653 + ! Platt et al. (1980) - inhibition
  654 + ! Light limitation factor (PI curve)
  655 + ps1= vp1/((alpha1/(alpha1+inib1))*(alpha1/(alpha1+inib1))**(inib1/alpha1))
  656 + ff1= ps1*(1.-exp(-1.*alpha1*par(ci)/ps1))*exp(-1.*inib1*par(ci)/ps1)
  657 +
  658 + ps2= vp2/((alpha2/(alpha2+inib2))*(alpha2/(alpha2+inib2))**(inib2/alpha2))
  659 + ff2= ps2*(1.-exp(-1.*alpha2*par(ci)/ps2))*exp(-1.*inib2*par(ci)/ps2)
  660 +
  661 + ! Nutrient limitation factors
  662 + qn1=(cc(n,ci)/kn1)/(1.+cc(n,ci)/kn1+cc(a,ci)/ka1)
  663 + qa1=(cc(a,ci)/ka1)/(1.+cc(n,ci)/kn1+cc(a,ci)/ka1)
  664 +
  665 + qn2=(cc(n,ci)/kn2)/(1.+cc(n,ci)/kn2+cc(a,ci)/ka2)
  666 + qa2=(cc(a,ci)/ka2)/(1.+cc(n,ci)/kn2+cc(a,ci)/ka2)
  667 +
  668 + ! Grazing preference normalization factors
  669 + fac1=(cc(z1,ci)+z0)/(k3*(r11*cc(p1,ci)+r12*cc(p2,ci)+0.5*r13*(cc(b1,ci)+cc(b2,ci))+r14*cc(d,ci)) &
  670 + +r11*cc(p1,ci)**2+r12*cc(p2,ci)**2+0.5*r13*(cc(b1,ci)+cc(b2,ci))**2+r14*cc(d,ci)**2)
  671 +
  672 + fac2=(cc(z2,ci)+z0)/(k3*(r21*cc(p1,ci)+r22*cc(p2,ci) &
  673 + +r23*cc(d,ci)+r24*cc(z1,ci)) &
  674 + +r21*cc(p1,ci)**2+r22*cc(p2,ci)**2 &
  675 + +r23*cc(d,ci)**2+r24*cc(z1,ci)**2)
  676 +
  677 + minal=min(cc(a,ci),etaa*cc(l,ci))
  678 + minhl=min(cc(hc,ci)/(k10+cc(hc,ci)),cc(l,ci)/(k4+cc(l,ci)))
  679 +
  680 + ! Light and nutrient limitation factors
  681 + lumlim1(ci) =amr1*ff1
  682 + nitlim1(ci) =qn1
  683 + ammlim1(ci) =qa1
  684 + lumlim2(ci) =amr2*ff2
  685 + nitlim2(ci) =qn2
  686 + ammlim2(ci) =qa2
  687 +
  688 + ! Nutrient uptake by pico- and nano-phytoplankton
  689 + dd(n,p1,ci) =amr1*ff1*qn1*(cc(p1,ci)+p0)
  690 + dd(a,p1,ci) =amr1*ff1*qa1*(cc(p1,ci)+p0)
  691 + dd(n,p2,ci) =amr2*ff2*qn2*(cc(p2,ci)+p0)
  692 + dd(a,p2,ci) =amr2*ff2*qa2*(cc(p2,ci)+p0)
  693 +
  694 + dd(p1,l,ci) =amr1*gamma*ff1*(qn1+qa1)*cc(p1,ci)
  695 + dd(p2,l,ci) =amr2*gamma*ff2*(qn2+qa2)*cc(p2,ci)
  696 +
  697 + dd(p1,d,ci) =amr1*mu11*(cc(p1,ci)+p0)/(k5+cc(p1,ci)+p0)*cc(p1,ci) &
  698 + +(1.-beta)*cc(p1,ci)**2*(hmr1*g1max*r11*fac1+hmr2*g2max*r21*fac2)
  699 + dd(p2,d,ci) =amr2*mu12*(cc(p2,ci)+p0)/(k5+cc(p2,ci)+p0)*cc(p2,ci) &
  700 + +(1.-beta)*cc(p2,ci)**2*(hmr1*g1max*r12*fac1+hmr2*g2max*r22*fac2)
  701 +
  702 + dd(b1,d,ci) =(1.-beta)*g1max*r13*cc(b1,ci)**2*fac1
  703 +
  704 + dd(p1,z1,ci)=hmr1*beta*g1max*r11*cc(p1,ci)**2*fac1
  705 + dd(p2,z1,ci)=hmr1*beta*g1max*r12*cc(p2,ci)**2*fac1
  706 + dd(b1,z1,ci)=hmr1*beta*g1max*0.5*r13*cc(b1,ci)**2*fac1
  707 + dd(b2,z1,ci)=hmr1*beta*g1max*0.5*r13*cc(b2,ci)**2*fac1
  708 + dd(d,z1,ci) =hmr1*beta*g1max*r14*cc(d,ci)**2*fac1
  709 +
  710 + dd(p1,z2,ci)=hmr2*beta*g2max*r21*cc(p1,ci)**2*fac2
  711 + dd(p2,z2,ci)=hmr2*beta*g2max*r22*cc(p2,ci)**2*fac2
  712 + dd(d,z2,ci) =hmr2*beta*g2max*r23*cc(d,ci)**2*fac2
  713 + dd(z1,z2,ci)=hmr2*beta*g2max*r24*cc(z1,ci)**2*fac2
  714 +
  715 + dd(b1,a,ci) =mu3*cc(b1,ci)
  716 + dd(b2,a,ci) =mu3*cc(b2,ci)
  717 + dd(d,l,ci) =mu4*cc(d,ci)
  718 + dd(a,n,ci) =mu5*cc(a,ci)
  719 +
  720 + dd(z1,d,ci) =(1.-beta)*hmr2*g2max*r24*cc(z1,ci)**2*fac2 &
  721 + +hmr1*(1.-epsi-delta)*mu21*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  722 + dd(z1,a,ci) =hmr1*epsi*mu21*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  723 + dd(z1,l,ci) =hmr1*delta*mu21*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  724 +
  725 + dd(z2,d,ci) =hmr2*(1.-epsi-delta)*mu22*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  726 + dd(z2,a,ci) =hmr2*epsi*mu22*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  727 + dd(z2,l,ci) =hmr2*delta*mu22*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  728 +
  729 + dd(a,b1,ci) =vb1*minal/(k4+minal+cc(l,ci))*(cc(b1,ci)+b0)
  730 + dd(l,b1,ci) =vb1*cc(l,ci)/(k4+minal+cc(l,ci))*(cc(b1,ci)+b0)
  731 + dd(hc,b2,ci)=vb2*minhl*(cc(b2,ci)+b0)
  732 + dd(l,b2,ci) =vb2*etah*minhl*(cc(b2,ci)+b0)
  733 + !dd(l,b2,ci) =vb2*cc(l,ci)/(k10+minhl+cc(l,ci))*(cc(b2,ci)+b0)
  734 + !dd(hc,b2,ci)=vb*cc(hc,ci)/(k10+cc(hc,ci))*(cc(b2,ci)+b0)
  735 +
  736 +
  737 + do i=1,numc
  738 + do j=1,numc
  739 + pp(i,j,ci)=dd(j,i,ci)
  740 + end do
  741 + end do
  742 + end do
  743 +
  744 + return
  745 + end subroutine do_bio_gsj
  746 +!EOC
  747 +
  748 +!-----------------------------------------------------------------------
  749 +!BOP
  750 +!
  751 +! !IROUTINE: Finish the bio calculations
  752 +!
  753 +! !INTERFACE:
  754 + subroutine end_bio_gsj
  755 +!
  756 +! !DESCRIPTION:
  757 +! Nothing done yet --- supplied for completeness.
  758 +!
  759 +! !USES:
  760 + IMPLICIT NONE
  761 +!
  762 +! !REVISION HISTORY:
  763 +! Original author(s): Hans Burchard & Karsten Bolding
  764 +!
  765 +!EOP
  766 +!-----------------------------------------------------------------------
  767 +!BOC
  768 +
  769 + return
  770 + end subroutine end_bio_gsj
  771 +!EOC
  772 +
  773 +!-----------------------------------------------------------------------
  774 +
  775 + end module bio_gsj
  776 +
  777 +!-----------------------------------------------------------------------
  778 +! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
  779 +!-----------------------------------------------------------------------
... ...
src/extras/bio/bio_gsj_zoo.F90 0 → 100755
... ... @@ -0,0 +1,922 @@
  1 +!$Id: bio_ismer.F90,v 1.11 2007-01-06 11:49:15 kbk Exp $
  2 +#include"cppdefs.h"
  3 +!-----------------------------------------------------------------------
  4 +!BOP
  5 +!
  6 +! !MODULE: bio_gsj --- Modified from Fasham et al. biological model \label{sec:bio-fasham}
  7 +!
  8 +! !INTERFACE:
  9 + module bio_gsj
  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_gsj, init_var_gsj, var_info_gsj, &
  47 + light_gsj, do_bio_gsj, end_bio_gsj
  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_init = 0.05
  58 + REALTYPE :: p2_init = 0.05
  59 + REALTYPE :: z1_init = 0.05
  60 + REALTYPE :: z2_init = 0.05
  61 + REALTYPE :: b_init = 0.001
  62 + REALTYPE :: d_init = 0.4
  63 + REALTYPE :: l_init = 0.14
  64 + REALTYPE :: p0 = 0.0
  65 + REALTYPE :: z0 = 0.0
  66 + REALTYPE :: b0 = 0.0
  67 + LOGICAL :: mte = .true.
  68 + REALTYPE :: ca1 = 3.61
  69 + REALTYPE :: ca2 = 14.58
  70 + REALTYPE :: ch1 = 3.265
  71 + REALTYPE :: ch2 = 24.923
  72 + REALTYPE :: amratio = 200.0
  73 + REALTYPE :: hmratio = 1000.0
  74 + REALTYPE :: vp1 = 1.5
  75 + REALTYPE :: alpha1 = 0.065
  76 + REALTYPE :: inib1 = 0.05
  77 + REALTYPE :: vp2 = 1.5
  78 + REALTYPE :: alpha2 = 0.065
  79 + REALTYPE :: inib2 = 0.05
  80 + REALTYPE :: theta = 0.0
  81 + REALTYPE :: w_p1min = -0.06
  82 + REALTYPE :: w_p1max = -0.38
  83 + REALTYPE :: w_p2min = -0.06
  84 + REALTYPE :: w_p2max = -0.38
  85 + REALTYPE :: kn1 = 0.2
  86 + REALTYPE :: ka1 = 0.8
  87 + REALTYPE :: kn2 = 0.2
  88 + REALTYPE :: ka2 = 0.8
  89 + REALTYPE :: mu11 = 0.05
  90 + REALTYPE :: mu12 = 0.05
  91 + REALTYPE :: k5 = 0.2
  92 + REALTYPE :: gamma = 0.05
  93 + REALTYPE :: w_p1 = -0.5
  94 + REALTYPE :: w_p2 = -0.5
  95 + REALTYPE :: g1max = 1.0
  96 + REALTYPE :: g2max = 1.0
  97 + REALTYPE :: k3 = 1.0
  98 + REALTYPE :: beta = 0.625
  99 + REALTYPE :: mu22 = 0.3
  100 +! REALTYPE :: k6 = 0.2
  101 +! REALTYPE :: lc50 =
  102 +! REALYTPE :: bhc = 11.1
  103 + REALTYPE :: delta = 0.1
  104 + REALTYPE :: epsi = 0.70
  105 + REALTYPE :: r11 = 0.55
  106 + REALTYPE :: r12 = 0.30
  107 + REALTYPE :: r13 = 0.05
  108 + REALTYPE :: r14 = 0.10
  109 + REALTYPE :: r21 = 0.50
  110 + REALTYPE :: r22 = 0.30
  111 + REALTYPE :: r23 = 0.05
  112 + REALTYPE :: r24 = 0.15
  113 + REALTYPE :: vb1 = 1.2
  114 + REALTYPE :: vb2 = 1.2
  115 + REALTYPE :: k4 = 0.5
  116 + REALTYPE :: k10 = 0.15
  117 + REALTYPE :: w_h = 0.0
  118 + REALTYPE :: mu3 = 0.15
  119 + REALTYPE :: etaa = 0.0
  120 + REALTYPE :: etah = 0.0
  121 + REALTYPE :: mu4 = 0.02
  122 + REALTYPE :: mu5 = 0.00
  123 + REALTYPE :: w_d = -2.0
  124 + REALTYPE, public :: kc = 0.03
  125 + integer :: out_unit
  126 + integer, parameter :: n=1,p1=2,p2=3,z1=4,z2=5,d=6,l=7,b1=8,a=9,hc=10,b2=11
  127 +!EOP
  128 +!! from a namelist
  129 +! p1_init : Small phytoplankton initial concentration
  130 +! p2_init : Large phytoplankton initial concentration
  131 +! z1_init : Small zooplankton initial concentration
  132 +! z2_init : Large zooplankton initial concentration
  133 +! b_init : Bacteria (Bac1) and OHC Bacteria (Bac2) initial concentration
  134 +! d_init : Detritus initial concentration
  135 +! l_init : Lablie Dissolved Organic Nitrogen (LDON) initial concentration
  136 +! p0 : Phytoplankton minimal concentration
  137 +! z0 : Zooplankton minimal concentration
  138 +! b0 : Bacteria minimal concentration
  139 +! mte : Temperature effect on phytoplankton and zooplankton loop
  140 +! hcz : Hydrocarbon mortality rate for large and small zooplankton
  141 +! ca1 :
  142 +! ca2 :
  143 +! ch1 :
  144 +! ch2 :
  145 +! amratio :
  146 +! hmratio :
  147 +! vp1 : Maximum uptake rate of small phytoplankton
  148 +! alpha1 : Slope of the small phytoplankton PI-curve
  149 +! inib1 : Inhibition slope of the small phytoplankton
  150 +! vp2 : Maximum uptake rate of large phytoplankton
  151 +! alpha2 : Slope of the large phytoplankton PI-curve
  152 +! inib2 : Inhibition slope of the large phytoplankton
  153 +! theta : Phytoplankton buoyancy parameter
  154 +! w_p1min : Maximum small phytoplankton settling velocity
  155 +! w_p1max : Minmimum small phytoplankton settling velocity
  156 +! w_p2min : Maximum large phytoplankton settling velocity
  157 +! w_p2max : Minmimum large phytoplankton settling velocity
  158 +! kn1 : Half saturation constant of nitrate uptake by small phytoplankton
  159 +! ka1 : Half saturation constant of ammonium uptake by small phytoplankton
  160 +! kn2 : Half saturation constant of nitrate uptake by large phytoplankton
  161 +! ka2 : Half saturation constant of ammonium uptake by large phytoplankton
  162 +! mu11 : Small phytoplankton mortality
  163 +! mu12 : Large phytoplankton mortality
  164 +! k5 : Half saturation of small and large phytoplankton mortality
  165 +! gamma : Excsudation fraction
  166 +! w_p1 : Small phytoplankton settling velocity
  167 +! w_p2 : Large phytoplankton settling velocity
  168 +! g1max : Maximum small zooplankton ingestion
  169 +! g2max : Maximum large zooplankton ingestion
  170 +! k3 : Half saturation constant of ingestion
  171 +! beta : Grazing efficiency
  172 +! mu21 : Maximum small zooplankton loss rate (variable function of hydrocarbon concentration)
  173 +! mu22 : Maximum large zooplankton loss rate
  174 +! k6 : Half saturation constant of zooplankton loss (small and large)
  175 +! k7 : Half saturation constant of zooplankton loss by hydrocarbon toxicity (small and large)
  176 +! lc50 : Half saturation concentration of hydrocarbon toxicity on zooplankton
  177 +! bhc : Slope factor of zooplankton mortality
  178 +! hc0 : Hydrocarbon sill concentration for zooplankton mortality by hydrocarbons
  179 +! delta : Fractional zooplankton loss of LDON
  180 +! epsi : Fractional zooplankton loss of ammonuim
  181 +! r11 : Small zooplankton preference on small phytoplantkon
  182 +! r12 : Small zooplankton preference on small phytoplantkon
  183 +! r13 : Small zooplankton preference on small phytoplantkon
  184 +! r14 : Small zooplankton preference on small phytoplantkon
  185 +! r21 : Small zooplankton preference on small phytoplantkon
  186 +! r22 : Small zooplankton preference on small phytoplantkon
  187 +! r23 : Small zooplankton preference on small phytoplantkon
  188 +! r24 : Small zooplankton preference on small phytoplantkon
  189 +! vb1 : Maximum bacteria (b1) uptake rate
  190 +! vb2 : Maximum hydrocarbon degrading bacteria (b2) uptake rate
  191 +! k4 : Half saturation of bacteria (b1) uptake
  192 +! k10 : Half saturation of hydrocarbon degrading bacteria (b2) uptake
  193 +! w_h : Hydrocarbon settling velocity
  194 +! mu3 : Bacterial excretion rate
  195 +! etaa : Uptake ratio ammonium : LDON
  196 +! etah : Uptake ratio hydrocarbon : LDON
  197 +! mu4 : Detritus breakdown rate
  198 +! mu5 : Nitrification rate
  199 +! w_d : Detritus settling velocity
  200 +! kc :
  201 +!-----------------------------------------------------------------------
  202 +
  203 + contains
  204 +
  205 +!-----------------------------------------------------------------------
  206 +!BOP
  207 +!
  208 +! !IROUTINE: Initialise the bio module
  209 +!
  210 +! !INTERFACE:
  211 + subroutine init_bio_gsj(namlst,fname,unit)
  212 +!
  213 +! !DESCRIPTION:
  214 +! Here, the bio namelist {\tt bio\_fasham.nml} is read and
  215 +! various variables (rates and settling velocities)
  216 +! are transformed into SI units.
  217 +!
  218 +! !USES:
  219 + IMPLICIT NONE
  220 +!
  221 +! !INPUT PARAMETERS:
  222 + integer, intent(in) :: namlst
  223 + character(len=*), intent(in) :: fname
  224 + character(len=20) :: pfile
  225 + integer, intent(in) :: unit
  226 +!
  227 +! !REVISION HISTORY:
  228 +! Original author(s): Hans Burchard & Karsten Bolding
  229 +!
  230 +! !LOCAL VARIABLES:
  231 + namelist /bio_gsj_nml/ numc, &
  232 + p1_init,p2_init,z1_init,z2_init, &
  233 + b_init,d_init,l_init, &
  234 + p0,z0,b0,vp1,alpha1,inib1,vp2,alpha2,inib2, &
  235 + kn1,ka1,kn2,ka2,mu11,mu12,k5,gamma,w_p1,w_p2, &
  236 + g1max,g2max,k3,beta,mu21,mu22,k6,delta,epsi, &
  237 + r11,r12,r13,r14,r21,r22,r23,r24, &
  238 + vb1,vb2,k4,k10,w_h,mu3,etaa,etah,mu4,w_d,kc,mu5, &
  239 + theta,w_p1max,w_p1min,w_p2min,w_p2max, &
  240 + mte,ca1,ca2,ch1,ch2,amratio,hmratio
  241 +
  242 +!EOP
  243 +!-----------------------------------------------------------------------
  244 +!BOC
  245 + LEVEL2 'init_bio_gsj'
  246 +
  247 + open(namlst,file=fname,action='read',status='old',err=98)
  248 + read(namlst,nml=bio_gsj_nml,err=99)
  249 + close(namlst)
  250 +
  251 + numcc=numc
  252 +
  253 +! Print some parameter values in standard output
  254 +! and save them in a separate file [out_fn]_ismer.par
  255 + pfile = trim(out_fn) // '_gsj.par'
  256 + open(10,status='unknown',action='write',file=pfile)
  257 + LEVEL3 'GSJ parameters saved in ', pfile
  258 + write(10,901) vp1
  259 + write(10,901) vp2
  260 + write(10,901) alpha1
  261 + write(10,901) alpha2
  262 + write(10,901) inib1
  263 + write(10,901) inib2
  264 + write(10,901) kn1
  265 + write(10,901) ka1
  266 + write(10,901) kn2
  267 + write(10,901) ka2
  268 + write(10,901) mu11
  269 + write(10,901) mu12
  270 + write(10,901) k5
  271 + write(10,901) gamma
  272 + write(10,901) w_p1
  273 + write(10,901) w_p2
  274 + write(10,901) g1max
  275 + write(10,901) g2max
  276 + write(10,901) k3
  277 + write(10,901) beta
  278 + write(10,901) mu21
  279 + write(10,901) mu22
  280 + write(10,901) k6
  281 + write(10,901) k6z
  282 + write(10,901) lc50
  283 + write(10,901) bhc
  284 + write(10,901) hc0
  285 + write(10,901) delta
  286 + write(10,901) epsi
  287 + write(10,901) r11
  288 + write(10,901) r12
  289 + write(10,901) r13
  290 + write(10,901) r14
  291 + write(10,901) r21
  292 + write(10,901) r22
  293 + write(10,901) r23
  294 + write(10,901) r24
  295 + write(10,901) vb1
  296 + write(10,901) vb2
  297 + write(10,901) k4
  298 + write(10,901) mu3
  299 + write(10,901) etaa
  300 + write(10,901) etah
  301 + write(10,901) mu4
  302 + write(10,901) mu5
  303 + write(10,901) w_d
  304 + write(10,901) kc
  305 + if (mte) then
  306 + write(10,901) ca1
  307 + write(10,901) ca2
  308 + write(10,901) ch1
  309 + write(10,901) ch2
  310 + write(10,901) amratio
  311 + write(10,901) hmratio
  312 + endif
  313 +
  314 +
  315 +901 format (f8.5)
  316 +
  317 +! Conversion from day to second
  318 + vp1 = vp1 /secs_pr_day
  319 + vp2 = vp2 /secs_pr_day
  320 + vb1 = vb1 /secs_pr_day
  321 + vb2 = vb2 /secs_pr_day
  322 + mu11 = mu11 /secs_pr_day
  323 + mu12 = mu12 /secs_pr_day
  324 +! mu21 = mu21 /secs_pr_day
  325 + mu22 = mu22 /secs_pr_day
  326 + mu3 = mu3 /secs_pr_day
  327 + mu4 = mu4 /secs_pr_day
  328 + mu5 = mu5 /secs_pr_day
  329 + g1max = g1max /secs_pr_day
  330 + g2max = g2max /secs_pr_day
  331 + w_p1 = w_p1 /secs_pr_day
  332 + w_p2 = w_p2 /secs_pr_day
  333 + w_p1min = w_p1min /secs_pr_day
  334 + w_p1max = w_p1max /secs_pr_day
  335 + w_p2min = w_p2min /secs_pr_day
  336 + w_p2max = w_p2max /secs_pr_day
  337 + theta = theta /secs_pr_day
  338 + w_d = w_d /secs_pr_day
  339 + w_h = w_h /secs_pr_day
  340 + alpha1 = alpha1 /secs_pr_day
  341 + inib1 = inib1 /secs_pr_day
  342 + alpha2 = alpha2 /secs_pr_day
  343 + inib2 = inib2 /secs_pr_day
  344 +
  345 + out_unit=unit
  346 +
  347 + LEVEL3 'GSJ bio module initialised ...'
  348 +
  349 + return
  350 +
  351 +98 LEVEL2 'I could not open bio_gsj.nml'
  352 + LEVEL2 'If thats not what you want you have to supply bio_gsj.nml'
  353 + LEVEL2 'See the bio example on www.gotm.net for a working bio_gsj.nml'
  354 + return
  355 +99 FATAL 'I could not read bio_gsj.nml'
  356 + stop 'init_bio_gsj'
  357 + end subroutine init_bio_gsj
  358 +!EOC
  359 +
  360 +!-----------------------------------------------------------------------
  361 +!BOP
  362 +!
  363 +! !IROUTINE: Initialise the concentration variables
  364 +!
  365 +! !INTERFACE:
  366 + subroutine init_var_gsj(nlev)
  367 +!
  368 +! !DESCRIPTION:
  369 +! Here, the the initial conditions are set and the settling velocities are
  370 +! transferred to all vertical levels. All concentrations are declared
  371 +! as non-negative variables, and it is defined which variables would be
  372 +! taken up by benthic filter feeders.
  373 +!
  374 +! !USES:
  375 + use observations, only: nprof,aprof,hcprof
  376 + use meanflow, only: nit,amm,hcb,T,S
  377 +
  378 + IMPLICIT NONE
  379 +
  380 +!
  381 +! !INPUT PARAMETERS:
  382 + integer, intent(in) :: nlev
  383 +!
  384 +! !REVISION HISTORY:
  385 +! Original author(s): Hans Burchard & Karsten Bolding
  386 +
  387 +! !LOCAL VARIABLES:
  388 + integer :: i
  389 +!EOP
  390 +!-----------------------------------------------------------------------
  391 +!BOC
  392 + do i=1,nlev
  393 + cc(n,i) = nprof(i) !CHG3
  394 + cc(p1,i)= p1_init
  395 + cc(p2,i)= p2_init
  396 + cc(z1,i)= z1_init
  397 + cc(z2,i)= z2_init
  398 + cc(d,i) = d_init
  399 + cc(l,i) = l_init
  400 + cc(b1,i)= b_init
  401 + cc(a,i) = aprof(i) !CHG5
  402 + cc(hc,i)= hcprof(i)
  403 + cc(b2,i)= b_init
  404 + end do
  405 +
  406 + do i=0,nlev
  407 + ws(n,i) = _ZERO_
  408 + ws(p1,i) = w_p1
  409 + ws(p2,i) = w_p2
  410 + ws(z1,i) = _ZERO_
  411 + ws(z2,i) = _ZERO_
  412 + ws(d,i) = w_d
  413 + ws(l,i) = _ZERO_
  414 + ws(b1,i) = _ZERO_
  415 + ws(a,i) = _ZERO_
  416 + ws(hc,i) = w_h
  417 + ws(b2,i) = _ZERO_
  418 + end do
  419 +
  420 + posconc(n) = 1
  421 + posconc(p1) = 1
  422 + posconc(p2) = 1
  423 + posconc(z1) = 1
  424 + posconc(z2) = 1
  425 + posconc(d) = 1
  426 + posconc(l) = 1
  427 + posconc(b1) = 1
  428 + posconc(a) = 1
  429 + posconc(hc) = 1
  430 + posconc(b2) = 1