Commit e57d1d1d9ad69d293ac0b42a58d399d92bb954cb

Authored by Dany Dumont
2 parents 0f87998d 5ab2e8dc
Exists in master and in 1 other branch snow

Merge branch 'master' of gitlasso.uqar.ca:dumoda01/gotm_ismer

Conflicts:
	src/Rules.make
compilers/compiler.GFORTRAN
... ... @@ -8,7 +8,7 @@ can_do_F90 = true
8 8 F90_to_f90 = $(FC) -E $(F90FLAGS) $(EXTRA_FFLAGS) $< > $@
9 9 #MODULES = -M $(MODDIR)
10 10 MODULES = -J $(MODDIR)
11   -EXTRAS = -ffree-form -Wall -std=f95 -pedantic
  11 +EXTRAS = -ffree-form -std=f95 -pedantic
12 12 DEBUG_FLAGS = -g
13 13 PROF_FLAGS = -p
14 14 PROD_FLAGS = -O0
... ...
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
... ... @@ -12,8 +12,8 @@ endif
12 12  
13 13 # Force this here. Could be done in bashrc.
14 14 #FORTRAN_COMPILER=PGF90
15   -#FORTRAN_COMPILER=GFORTRAN
16   -FORTRAN_COMPILER=IFORT
  15 +FORTRAN_COMPILER=GFORTRAN
  16 +#FORTRAN_COMPILER=IFORT
17 17  
18 18 DEFINES=-DNUDGE_VEL
19 19 DEFINES=-D$(FORTRAN_COMPILER)
... ...
src/extras/bio/Makefile
... ... @@ -16,6 +16,7 @@ bio_sed.F90 \
16 16 bio_fasham.F90 \
17 17 bio_ismer.F90 \
18 18 bio_gsj.F90 \
  19 +bio_nocera.F90 \
19 20 bio_save.F90 \
20 21 nitrate.F90 \
21 22 ammonium.F90 \
... ... @@ -31,6 +32,7 @@ ${LIB}(bio_mab.o) \
31 32 ${LIB}(bio_ismer.o) \
32 33 ${LIB}(bio_gsj.o) \
33 34 ${LIB}(bio_fasham.o) \
  35 +${LIB}(bio_nocera.o) \
34 36 ${LIB}(bio_sed.o) \
35 37 ${LIB}(bio_save.o) \
36 38 ${LIB}(nitrate.o) \
... ...
src/extras/bio/ammonium.F90
... ... @@ -176,6 +176,10 @@
176 176 do i=1,nlev
177 177 amm(i) = cc(9,i)
178 178 end do
  179 + else if (bio_model.eq.8) then
  180 + do i=1,nlev
  181 + amm(i) = cc(6,i)
  182 + end do
179 183 end if
180 184 #endif
181 185  
... ...
src/extras/bio/bio.F90
... ... @@ -22,27 +22,30 @@
22 22 use bio_template, only : init_bio_template,init_var_template
23 23 use bio_template, only : var_info_template,light_template
24 24  
25   - use bio_npzd, only : init_bio_npzd,init_var_npzd,var_info_npzd
26   - use bio_npzd, only : light_npzd, do_bio_npzd
  25 + use bio_npzd, only : init_bio_npzd,init_var_npzd,var_info_npzd
  26 + use bio_npzd, only : light_npzd, do_bio_npzd
27 27  
28   - use bio_iow, only : init_bio_iow,init_var_iow,var_info_iow
29   - use bio_iow, only : light_iow,surface_fluxes_iow,do_bio_iow
  28 + use bio_iow, only : init_bio_iow,init_var_iow,var_info_iow
  29 + use bio_iow, only : light_iow,surface_fluxes_iow,do_bio_iow
30 30  
31   - use bio_fasham, only : init_bio_fasham,init_var_fasham,var_info_fasham
32   - use bio_fasham, only : light_fasham,do_bio_fasham
  31 + use bio_fasham, only : init_bio_fasham,init_var_fasham,var_info_fasham
  32 + use bio_fasham, only : light_fasham,do_bio_fasham
33 33  
34   - use bio_ismer, only : init_bio_ismer,init_var_ismer,var_info_ismer
35   - use bio_ismer, only : light_ismer,do_bio_ismer
  34 + use bio_ismer, only : init_bio_ismer,init_var_ismer,var_info_ismer
  35 + use bio_ismer, only : light_ismer,do_bio_ismer
36 36  
37   - use bio_gsj, only : init_bio_gsj,init_var_gsj,var_info_gsj
38   - use bio_gsj, only : light_gsj,do_bio_gsj
  37 + use bio_gsj, only : init_bio_gsj,init_var_gsj,var_info_gsj
  38 + use bio_gsj, only : light_gsj,do_bio_gsj
39 39  
40   - use bio_sed, only : init_bio_sed,init_var_sed,var_info_sed
  40 + use bio_nocera, only : init_bio_nocera,init_var_nocera,var_info_nocera
  41 + use bio_nocera, only : light_nocera,do_bio_nocera
41 42  
42   - use bio_mab, only : init_bio_mab,init_var_mab,var_info_mab
43   - use bio_mab, only : light_mab,surface_fluxes_mab,do_bio_mab
  43 + use bio_sed, only : init_bio_sed,init_var_sed,var_info_sed
44 44  
45   - use output, only: out_fmt,write_results,ts
  45 + use bio_mab, only : init_bio_mab,init_var_mab,var_info_mab
  46 + use bio_mab, only : light_mab,surface_fluxes_mab,do_bio_mab
  47 +
  48 + use output, only : out_fmt,write_results,ts
46 49  
47 50 use util
48 51 !
... ... @@ -327,6 +330,16 @@
327 330  
328 331 call var_info_gsj()
329 332  
  333 + case (8) ! The model of Ariadna Celina Nocera, modified from Fasham
  334 +
  335 + call init_bio_nocera(namlst,'bio_nocera.nml',unit)
  336 +
  337 + call allocate_memory(nlev)
  338 +
  339 + call init_var_nocera(nlev)
  340 +
  341 + call var_info_nocera()
  342 +
330 343 case default
331 344 stop "bio: no valid biomodel specified in bio.nml !"
332 345 end select
... ... @@ -710,6 +723,9 @@
710 723 case (7)
711 724 call light_gsj(nlev,bioshade_feedback)
712 725 call ode_solver(ode_method,numc,nlev,dt_eff,cc,do_bio_gsj)
  726 + case (8)
  727 + call light_nocera(nlev,bioshade_feedback)
  728 + call ode_solver(ode_method,numc,nlev,dt_eff,cc,do_bio_nocera)
713 729 end select
714 730  
715 731 end do
... ...
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: