Commit 84ed20a57205cc73a065326abfffceeb1ab1f8b8

Authored by Dany Dumont
1 parent 8d24f9ca
Exists in master and in 1 other branch snow

Construction d'un nouveau modele biogeochimique bio_gsj, base sur bio_ismer, com…

…prenant une classe hydrocarbure et une deuxieme classe de bacteries pouvant degrader les hydrocarbures. Le traceur hydrocarbon (hcb) est calque sur le traceur nitrate auquel on a ajoute la possibilite de specifier, via la namelist bio_gsj, une vitesse vertical (settling or rising). Trois fichers namelist sont affectes par ce changement (bio, obs et bio_gsj).
nml/bio.nml
... ... @@ -10,6 +10,7 @@
10 10 ! 4: Fasham et al. 1990 (7 variables)
11 11 ! 5: IOW-ERGOM MaBenE version (9 variables)
12 12 ! 6: ISMER model (9 variables)
  13 +! 7: GSJ model (11 variables)
13 14 !
14 15 ! bio_eulerian -> state variables are Eulerian (.true./.false.)
15 16 !
... ...
nml/bio_gsj.nml 0 → 100644
... ... @@ -0,0 +1,134 @@
  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 +! vb = maximum bacterial uptake rate [1/day]
  59 +! k4 = half saturation bacterial uptake [mmol n/m3]
  60 +! mu3 = bacteria excretion rate [1/day]
  61 +! eta = uptake ratio ammonium:LDON [-]
  62 +! mu4 = detritus breakdown rate [1/day]
  63 +! mu5 = nitrification rate [1/day]
  64 +! w_d = detritus settling velocity [m/day]
  65 +! kc = attenuation constant for the self shading effect [m**2/mmol N]
  66 +!-------------------------------------------------------------------------------
  67 + &bio_gsj_nml
  68 + numc = 11
  69 + p1_init = 0.012
  70 + p2_init = 0.012
  71 + z1_init = 0.012
  72 + z2_init = 0.012
  73 + b_init = 0.001
  74 + d_init = 0.01
  75 + l_init = 0.1
  76 + p0 = 0.0001
  77 + z0 = 0.0001
  78 + b0 = 0.0001
  79 + mte = .true.
  80 + ca1 = 3.61
  81 + ca2 = 14.58
  82 + ch1 = 3.265
  83 + ch2 = 24.923
  84 + amratio = 200
  85 + hmratio = 1000
  86 + vp1 = 0.02
  87 + vp2 = 0.8
  88 + alpha1 = 0.02
  89 + alpha2 = 0.04
  90 + inib1 = 0.0
  91 + inib2 = 0.006
  92 + kn1 = 1.0
  93 + ka1 = 0.8
  94 + kn2 = 1.0
  95 + ka2 = 0.8
  96 + mu11 = 0.05
  97 + mu12 = 0.05
  98 + k5 = 0.2
  99 + gamma = 0.05
  100 + w_p1 =-0.38
  101 + w_p2 =-0.00
  102 + theta = 0.0
  103 + w_p1min =-0.01
  104 + w_p1max =-0.10
  105 + w_p2min =-0.05
  106 + w_p2max =-0.38
  107 + g1max = 1.0
  108 + g2max = 1.0
  109 + k3 = 1.0
  110 + beta = 0.625
  111 + mu21 = 0.3
  112 + mu22 = 0.3
  113 + k6 = 0.2
  114 + delta = 0.1
  115 + epsi = 0.70
  116 + r11 = 0.55
  117 + r12 = 0.30
  118 + r13 = 0.05
  119 + r14 = 0.10
  120 + r21 = 0.50
  121 + r22 = 0.30
  122 + r23 = 0.05
  123 + r24 = 0.15
  124 + vb = 0.24
  125 + k4 = 0.5
  126 + k10 = 0.5
  127 + w_h = 100.0
  128 + mu3 = 0.03
  129 + eta = 0.0
  130 + mu4 = 0.02
  131 + mu5 = 0.00
  132 + w_d = -5.0
  133 + kc = 0.03
  134 + /
... ...
nml/obs.nml
... ... @@ -64,6 +64,68 @@
64 64 /
65 65  
66 66 !-------------------------------------------------------------------------------
  67 +! observed or prescribed potential temperature profiles
  68 +!
  69 +! t_prof_method -> method to create initial or observed temperature profiles
  70 +! 0: no initial values, T-equation is not solved
  71 +! 1: use analytically prescribed initial profile
  72 +! 2: read profiles at different dates from "t_prof_file"
  73 +! and interpolate to GOTM timestep
  74 +!
  75 +! t_analyt_method -> method to create analytically precribed inital profile
  76 +! 1: set profile to constant value s_1
  77 +! 2: set "two layer" stratification (see user's guide)
  78 +! 3: set profile with constant N^2 (see user's guide)
  79 +! This option can only be used toghether with
  80 +! s_analyt_method=1 (constant salinity).
  81 +!
  82 +! z_t1 -> upper layer thickness if t_analyt_method=2
  83 +!
  84 +! t_1 -> constant temperature if t_analyt_method=1
  85 +! upper layer temperature if t_analyt_method=2
  86 +! surface temperature if t_analyt_method=3
  87 +!
  88 +! z_t2 -> depth of top of lower layer if t_analyt_method=2
  89 +!
  90 +! t_2 -> lower layer temperature if t_analyt_method=2
  91 +!
  92 +! t_obs_NN -> constant value N^2 corresponding to temperature profile
  93 +! if t_analyt_method=3
  94 +
  95 +! t_prof_file -> filename of file with temperature profiles
  96 +! if t_prof_method=2
  97 +!
  98 +! The computed profiles can be relaxed towards observed or prescribed
  99 +! profiles with a certain relaxation time. If you do not want relaxation,
  100 +! set the relaxation times to 1.e15 (something large). It is possible to choose
  101 +! different relaxation times in a surface and bottom layer.
  102 +!
  103 +! TRelaxTauM -> relaxation time for bulk of the flow
  104 +! TRelaxTauB -> relaxation time for bottom layer
  105 +! TRelaxTauS -> relaxation time for surface layer
  106 +! TRelaxBott -> height of bottom relaxation layer
  107 +! (set to 0. if not used)
  108 +! TRelaxSurf -> height of surface relaxation layer
  109 +! (set to 0. if not used)
  110 +!
  111 +!-------------------------------------------------------------------------------
  112 + &tprofile
  113 + t_prof_method= 2
  114 + t_analyt_method= 2
  115 + z_t1= 25.
  116 + t_1= -1.0
  117 + z_t2= 35.
  118 + t_2= 0.0
  119 + t_obs_NN= 2.56e-4
  120 + t_prof_file= 'franklin_tprof_ctd.dat'
  121 + TRelaxTauM= 1209600.
  122 + TRelaxTauB= 1209600.
  123 + TRelaxTauS= 1209600.
  124 + TRelaxBott= 0.
  125 + TRelaxSurf= 0.
  126 + /
  127 +
  128 +!-------------------------------------------------------------------------------
67 129 ! observed or prescribed nitrate profiles
68 130 !
69 131 ! n_prof_method -> method to create initial or observed nitrate profiles
... ... @@ -85,7 +147,6 @@
85 147 !
86 148 ! n_2 -> lower layer nitrate if n_analyt_method=2
87 149 !
88   -
89 150 ! n_prof_file -> filename of file with nitrate profiles
90 151 ! if n_prof_method=2
91 152 !
... ... @@ -172,65 +233,57 @@
172 233 /
173 234  
174 235 !-------------------------------------------------------------------------------
175   -! observed or prescribed potential temperature profiles
  236 +! observed or prescribed hydrocarbon profiles
176 237 !
177   -! t_prof_method -> method to create initial or observed temperature profiles
178   -! 0: no initial values, T-equation is not solved
  238 +! hc_prof_method -> method to create initial or observed hydroarbon profiles
  239 +! 0: no initial values, N-equation is not solved
179 240 ! 1: use analytically prescribed initial profile
180   -! 2: read profiles at different dates from "t_prof_file"
  241 +! 2: read profiles at different dates from "hc_prof_file"
181 242 ! and interpolate to GOTM timestep
182 243 !
183   -! t_analyt_method -> method to create analytically precribed inital profile
184   -! 1: set profile to constant value s_1
  244 +! hc_analyt_method -> method to create analytically precribed inital profile
  245 +! 1: set profile to constant value hc_1
185 246 ! 2: set "two layer" stratification (see user's guide)
186   -! 3: set profile with constant N^2 (see user's guide)
187   -! This option can only be used toghether with
188   -! s_analyt_method=1 (constant salinity).
189 247 !
190   -! z_t1 -> upper layer thickness if t_analyt_method=2
  248 +! z_hc1 -> upper layer thickness if hc_analyt_method=2
191 249 !
192   -! t_1 -> constant temperature if t_analyt_method=1
193   -! upper layer temperature if t_analyt_method=2
194   -! surface temperature if t_analyt_method=3
  250 +! hc_1 -> constant nitrate if hc_analyt_method=1
  251 +! upper layer nitrate if hc_analyt_method=2
195 252 !
196   -! z_t2 -> depth of top of lower layer if t_analyt_method=2
  253 +! z_hc2 -> depth of top of lower layer if hc_analyt_method=2
197 254 !
198   -! t_2 -> lower layer temperature if t_analyt_method=2
  255 +! hc_2 -> lower layer nitrate if hc_analyt_method=2
199 256 !
200   -! t_obs_NN -> constant value N^2 corresponding to temperature profile
201   -! if t_analyt_method=3
202   -
203   -! t_prof_file -> filename of file with temperature profiles
204   -! if t_prof_method=2
  257 +! hc_prof_file -> filename of file with hydrocarbon profiles
  258 +! if hc_prof_method=2
205 259 !
206 260 ! The computed profiles can be relaxed towards observed or prescribed
207 261 ! profiles with a certain relaxation time. If you do not want relaxation,
208 262 ! set the relaxation times to 1.e15 (something large). It is possible to choose
209 263 ! different relaxation times in a surface and bottom layer.
210 264 !
211   -! TRelaxTauM -> relaxation time for bulk of the flow
212   -! TRelaxTauB -> relaxation time for bottom layer
213   -! TRelaxTauS -> relaxation time for surface layer
214   -! TRelaxBott -> height of bottom relaxation layer
  265 +! HCRelaxTauM -> relaxation time for bulk of the flow
  266 +! HCRelaxTauB -> relaxation time for bottom layer
  267 +! HCRelaxTauS -> relaxation time for surface layer
  268 +! HCRelaxBott -> height of bottom relaxation layer
215 269 ! (set to 0. if not used)
216   -! TRelaxSurf -> height of surface relaxation layer
  270 +! HCRelaxSurf -> height of surface relaxation layer
217 271 ! (set to 0. if not used)
218 272 !
219 273 !-------------------------------------------------------------------------------
220   - &tprofile
221   - t_prof_method= 2
222   - t_analyt_method= 2
223   - z_t1= 25.
224   - t_1= -1.0
225   - z_t2= 35.
226   - t_2= 0.0
227   - t_obs_NN= 2.56e-4
228   - t_prof_file= 'franklin_tprof_ctd.dat'
229   - TRelaxTauM= 1209600.
230   - TRelaxTauB= 1209600.
231   - TRelaxTauS= 1209600.
232   - TRelaxBott= 0.
233   - TRelaxSurf= 0.
  274 + &hcprofile
  275 + hc_prof_method= 1
  276 + hc_analyt_method= 2
  277 + z_hc1= 35.
  278 + hc_1= 3.0
  279 + z_hc2= 70.
  280 + hc_2= 15.0
  281 + hc_prof_file= 'hcprof.dat'
  282 + HCRelaxTauM= 1.e15
  283 + HCRelaxTauB= 1.e15
  284 + HCRelaxTauS= 1.e15
  285 + HCRelaxBott= 0.
  286 + HCRelaxSurf= 0.
234 287 /
235 288  
236 289 !-------------------------------------------------------------------------------
... ...
src/Makefile
1   -#$Id: Makefile,v 1.9 2005-08-16 14:38:37 kbk Exp $
  1 +#Id: Makefile,v 1.9 2005-08-16 14:38:37 kbk Exp $
2 2 #
3 3 # Master Makefile for making the GOTM executable.
4 4 #
... ...
src/extras/bio/Makefile
... ... @@ -15,9 +15,11 @@ bio_iow.F90 \
15 15 bio_sed.F90 \
16 16 bio_fasham.F90 \
17 17 bio_ismer.F90 \
  18 +bio_gsj.F90 \
18 19 bio_save.F90 \
19 20 nitrate.F90 \
20   -ammonium.F90
  21 +ammonium.F90 \
  22 +hydrocarbon.F90
21 23  
22 24 OBJ = \
23 25 ${LIB}(bio_var.o) \
... ... @@ -27,11 +29,13 @@ ${LIB}(bio_npzd.o) \
27 29 ${LIB}(bio_iow.o) \
28 30 ${LIB}(bio_mab.o) \
29 31 ${LIB}(bio_ismer.o) \
  32 +${LIB}(bio_gsj.o) \
30 33 ${LIB}(bio_fasham.o) \
31 34 ${LIB}(bio_sed.o) \
32 35 ${LIB}(bio_save.o) \
33 36 ${LIB}(nitrate.o) \
34 37 ${LIB}(ammonium.o) \
  38 +${LIB}(hydrocarbon.o) \
35 39 ${LIB}(bio.o)
36 40  
37 41 all: ${OBJ}
... ...
src/extras/bio/ammonium.F90
... ... @@ -172,6 +172,10 @@
172 172 do i=1,nlev
173 173 amm(i) = cc(9,i)
174 174 end do
  175 + else if (bio_model.eq.7) then
  176 + do i=1,nlev
  177 + amm(i) = cc(9,i)
  178 + end do
175 179 end if
176 180 #endif
177 181  
... ...
src/extras/bio/bio.F90
... ... @@ -30,16 +30,19 @@
30 30  
31 31 use bio_fasham, only : init_bio_fasham,init_var_fasham,var_info_fasham
32 32 use bio_fasham, only : light_fasham,do_bio_fasham
33   -!DD
  33 +
34 34 use bio_ismer, only : init_bio_ismer,init_var_ismer,var_info_ismer
35 35 use bio_ismer, only : light_ismer,do_bio_ismer
36 36  
37   - use bio_sed, only : init_bio_sed,init_var_sed,var_info_sed
  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 +
  40 + use bio_sed, only : init_bio_sed,init_var_sed,var_info_sed
38 41  
39   - use bio_mab, only : init_bio_mab,init_var_mab,var_info_mab
40   - use bio_mab, only : light_mab,surface_fluxes_mab,do_bio_mab
  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
41 44  
42   - use output, only: out_fmt,write_results,ts
  45 + use output, only: out_fmt,write_results,ts
43 46  
44 47 use util
45 48 !
... ... @@ -314,6 +317,16 @@
314 317  
315 318 call var_info_ismer()
316 319  
  320 + case (7) ! The GSJ model, modified from ISMER
  321 +
  322 + call init_bio_gsj(namlst,'bio_gsj.nml',unit)
  323 +
  324 + call allocate_memory(nlev)
  325 +
  326 + call init_var_gsj(nlev)
  327 +
  328 + call var_info_gsj()
  329 +
317 330 case default
318 331 stop "bio: no valid biomodel specified in bio.nml !"
319 332 end select
... ... @@ -409,8 +422,8 @@
409 422 ! modules
410 423 !
411 424 ! !INTERFACE:
412   - subroutine set_env_bio(nlev,h_,t_,s_,nit_,amm_,rho_,nuh_,rad_,wind_,I_0_, & !CHG3-5
413   - w_,w_adv_ctr_)
  425 + subroutine set_env_bio(nlev,h_,t_,s_,nit_,amm_,hcb_,rho_,nuh_,rad_,wind_, &
  426 + I_0_,w_,w_adv_ctr_)
414 427 !
415 428 ! !DESCRIPTION:
416 429 !
... ... @@ -423,8 +436,9 @@
423 436 REALTYPE, intent(in) :: nuh_(0:nlev)
424 437 REALTYPE, intent(in) :: t_(0:nlev)
425 438 REALTYPE, intent(in) :: s_(0:nlev)
426   - REALTYPE, intent(in) :: nit_(0:nlev) !CHG3
427   - REALTYPE, intent(in) :: amm_(0:nlev) !CHG5
  439 + REALTYPE, intent(in) :: nit_(0:nlev)
  440 + REALTYPE, intent(in) :: amm_(0:nlev)
  441 + REALTYPE, intent(in) :: hcb_(0:nlev)
428 442 REALTYPE, intent(in) :: rho_(0:nlev)
429 443 REALTYPE, intent(in) :: rad_(0:nlev)
430 444 REALTYPE, intent(in) :: wind_
... ... @@ -443,9 +457,10 @@
443 457 h = h_
444 458 t = t_
445 459 s = s_
446   - nit = nit_ !CHG3
447   - amm = amm_ !CHG5
  460 + nit = nit_
  461 + amm = amm_
448 462 rho = rho_
  463 + hcb = hcb_
449 464 nuh = nuh_
450 465 rad = rad_
451 466 wind = wind_
... ... @@ -565,7 +580,19 @@
565 580 if (bio_eulerian) then
566 581 do j=1,numcc
567 582  
568   - if (j==1 .and. bio_model==6) then
  583 + if (j==1 .and. bio_model==7) then
  584 + do i=1,nlev
  585 + cc(j,i) = nit(i)
  586 + end do
  587 + else if (j==9 .and. bio_model==7) then
  588 + do i=1,nlev
  589 + cc(j,i) = amm(i)
  590 + end do
  591 + else if (j==10 .and. bio_model==7) then
  592 + do i=1,nlev
  593 + cc(j,i) = hcb(i)
  594 + end do
  595 + else if (j==1 .and. bio_model==6) then
569 596 do i=1,nlev
570 597 cc(j,i) = nit(i)
571 598 end do
... ... @@ -679,6 +706,9 @@
679 706 case (6)
680 707 call light_ismer(nlev,bioshade_feedback)
681 708 call ode_solver(ode_method,numc,nlev,dt_eff,cc,do_bio_ismer)
  709 + case (7)
  710 + call light_gsj(nlev,bioshade_feedback)
  711 + call ode_solver(ode_method,numc,nlev,dt_eff,cc,do_bio_gsj)
682 712 end select
683 713  
684 714 end do
... ... @@ -772,6 +802,7 @@
772 802 if (allocated(s)) deallocate(s)
773 803 if (allocated(nit)) deallocate(nit) !CHG3
774 804 if (allocated(amm)) deallocate(amm) !CHG5
  805 + if (allocated(hcb)) deallocate(hcb)
775 806 if (allocated(rho)) deallocate(rho)
776 807 if (allocated(rad)) deallocate(rad)
777 808 if (allocated(w)) deallocate(w)
... ... @@ -892,6 +923,12 @@
892 923 allocate(aprof(0:nlev),stat=rc) !CHG5
893 924 if (rc /= 0) stop 'init_bio(): Error allocating (aprof)'
894 925  
  926 + allocate(hcb(0:nlev),stat=rc)
  927 + if (rc /= 0) stop 'init_bio(): Error allocating (hcb)'
  928 +
  929 + allocate(hcprof(0:nlev),stat=rc)
  930 + if (rc /= 0) stop 'init_bio(): Error allocating (hcprof)'
  931 +
895 932 allocate(rho(0:nlev),stat=rc)
896 933 if (rc /= 0) stop 'init_bio(): Error allocating (rho)'
897 934  
... ...
src/extras/bio/bio_gsj.F90 0 → 100644
... ... @@ -0,0 +1,772 @@
  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 :: vb = 1.2
  113 + REALTYPE :: k4 = 0.5
  114 + REALTYPE :: k10 = 0.5
  115 + REALTYPE :: w_h = 0.0
  116 + REALTYPE :: mu3 = 0.15
  117 + REALTYPE :: eta = 0.0
  118 + REALTYPE :: mu4 = 0.02
  119 + REALTYPE :: mu5 = 0.00
  120 + REALTYPE :: w_d = -2.0
  121 + REALTYPE, public :: kc = 0.03
  122 + integer :: out_unit
  123 + integer, parameter :: n=1,p1=2,p2=3,z1=4,z2=5,d=6,l=7,b1=8,a=9,hc=10,b2=11
  124 +!EOP
  125 +!-----------------------------------------------------------------------
  126 +
  127 + contains
  128 +
  129 +!-----------------------------------------------------------------------
  130 +!BOP
  131 +!
  132 +! !IROUTINE: Initialise the bio module
  133 +!
  134 +! !INTERFACE:
  135 + subroutine init_bio_gsj(namlst,fname,unit)
  136 +!
  137 +! !DESCRIPTION:
  138 +! Here, the bio namelist {\tt bio\_fasham.nml} is read and
  139 +! various variables (rates and settling velocities)
  140 +! are transformed into SI units.
  141 +!
  142 +! !USES:
  143 + IMPLICIT NONE
  144 +!
  145 +! !INPUT PARAMETERS:
  146 + integer, intent(in) :: namlst
  147 + character(len=*), intent(in) :: fname
  148 + character(len=20) :: pfile
  149 + integer, intent(in) :: unit
  150 +!
  151 +! !REVISION HISTORY:
  152 +! Original author(s): Hans Burchard & Karsten Bolding
  153 +!
  154 +! !LOCAL VARIABLES:
  155 + namelist /bio_gsj_nml/ numc, &
  156 + p1_init,p2_init,z1_init,z2_init, &
  157 + b_init,d_init,l_init, &
  158 + p0,z0,b0,vp1,alpha1,inib1,vp2,alpha2,inib2, &
  159 + kn1,ka1,kn2,ka2,mu11,mu12,k5,gamma,w_p1,w_p2, &
  160 + g1max,g2max,k3,beta,mu21,mu22,k6,delta,epsi, &
  161 + r11,r12,r13,r14,r21,r22,r23,r24, &
  162 + vb,k4,k10,w_h,mu3,eta,mu4,w_d,kc,mu5, &
  163 + theta,w_p1max,w_p1min,w_p2min,w_p2max, &
  164 + mte,ca1,ca2,ch1,ch2,amratio,hmratio
  165 +
  166 +!EOP
  167 +!-----------------------------------------------------------------------
  168 +!BOC
  169 + LEVEL2 'init_bio_gsj'
  170 +
  171 + open(namlst,file=fname,action='read',status='old',err=98)
  172 + read(namlst,nml=bio_gsj_nml,err=99)
  173 + close(namlst)
  174 +
  175 + numcc=numc
  176 +
  177 +! Print some parameter values in standard output
  178 +! and save them in a separate file [out_fn]_ismer.par
  179 + pfile = trim(out_fn) // '_gsj.par'
  180 + open(10,status='unknown',action='write',file=pfile)
  181 + LEVEL3 'GSJ parameters saved in ', pfile
  182 + write(10,901) vp1
  183 + write(10,901) vp2
  184 + write(10,901) alpha1
  185 + write(10,901) alpha2
  186 + write(10,901) inib1
  187 + write(10,901) inib2
  188 + write(10,901) kn1
  189 + write(10,901) ka1
  190 + write(10,901) kn2
  191 + write(10,901) ka2
  192 + write(10,901) mu11
  193 + write(10,901) mu12
  194 + write(10,901) k5
  195 + write(10,901) gamma
  196 + write(10,901) w_p1
  197 + write(10,901) w_p2
  198 + write(10,901) g1max
  199 + write(10,901) g2max
  200 + write(10,901) k3
  201 + write(10,901) beta
  202 + write(10,901) mu21
  203 + write(10,901) mu22
  204 + write(10,901) k6
  205 + write(10,901) delta
  206 + write(10,901) epsi
  207 + write(10,901) r11
  208 + write(10,901) r12
  209 + write(10,901) r13
  210 + write(10,901) r14
  211 + write(10,901) r21
  212 + write(10,901) r22
  213 + write(10,901) r23
  214 + write(10,901) r24
  215 + write(10,901) vb
  216 + write(10,901) k4
  217 + write(10,901) mu3
  218 + write(10,901) eta
  219 + write(10,901) mu4
  220 + write(10,901) mu5
  221 + write(10,901) w_d
  222 + write(10,901) kc
  223 + if (mte) then
  224 + write(10,901) ca1
  225 + write(10,901) ca2
  226 + write(10,901) ch1
  227 + write(10,901) ch2
  228 + write(10,901) amratio
  229 + write(10,901) hmratio
  230 + endif
  231 +
  232 +
  233 +900 format (a,f8.5)
  234 +901 format (f8.5)
  235 +
  236 +! Conversion from day to second
  237 + vp1 = vp1 /secs_pr_day
  238 + vp2 = vp2 /secs_pr_day
  239 + vb = vb /secs_pr_day
  240 + mu11 = mu11 /secs_pr_day
  241 + mu12 = mu12 /secs_pr_day
  242 + mu21 = mu21 /secs_pr_day
  243 + mu22 = mu22 /secs_pr_day
  244 + mu3 = mu3 /secs_pr_day
  245 + mu4 = mu4 /secs_pr_day
  246 + mu5 = mu5 /secs_pr_day
  247 + g1max = g1max /secs_pr_day
  248 + g2max = g2max /secs_pr_day
  249 + w_p1 = w_p1 /secs_pr_day
  250 + w_p2 = w_p2 /secs_pr_day
  251 + w_p1min = w_p1min /secs_pr_day
  252 + w_p1max = w_p1max /secs_pr_day
  253 + w_p2min = w_p2min /secs_pr_day
  254 + w_p2max = w_p2max /secs_pr_day
  255 + theta = theta /secs_pr_day
  256 + w_d = w_d /secs_pr_day
  257 + alpha1 = alpha1 /secs_pr_day
  258 + inib1 = inib1 /secs_pr_day
  259 + alpha2 = alpha2 /secs_pr_day
  260 + inib2 = inib2 /secs_pr_day
  261 +
  262 + out_unit=unit
  263 +
  264 + LEVEL3 'GSJ bio module initialised ...'
  265 +
  266 + return
  267 +
  268 +98 LEVEL2 'I could not open bio_gsj.nml'
  269 + LEVEL2 'If thats not what you want you have to supply bio_gsj.nml'
  270 + LEVEL2 'See the bio example on www.gotm.net for a working bio_gsj.nml'
  271 + return
  272 +99 FATAL 'I could not read bio_gsj.nml'
  273 + stop 'init_bio_gsj'
  274 + end subroutine init_bio_gsj
  275 +!EOC
  276 +
  277 +!-----------------------------------------------------------------------
  278 +!BOP
  279 +!
  280 +! !IROUTINE: Initialise the concentration variables
  281 +!
  282 +! !INTERFACE:
  283 + subroutine init_var_gsj(nlev)
  284 +!
  285 +! !DESCRIPTION:
  286 +! Here, the the initial conditions are set and the settling velocities are
  287 +! transferred to all vertical levels. All concentrations are declared
  288 +! as non-negative variables, and it is defined which variables would be
  289 +! taken up by benthic filter feeders.
  290 +!
  291 +! !USES:
  292 + use observations, only: nprof,aprof,hcprof
  293 + use meanflow, only: nit,amm,hcb,T,S
  294 +
  295 + IMPLICIT NONE
  296 +
  297 +!
  298 +! !INPUT PARAMETERS:
  299 + integer, intent(in) :: nlev
  300 +!
  301 +! !REVISION HISTORY:
  302 +! Original author(s): Hans Burchard & Karsten Bolding
  303 +
  304 +! !LOCAL VARIABLES:
  305 + integer :: i
  306 +!EOP
  307 +!-----------------------------------------------------------------------
  308 +!BOC
  309 + do i=1,nlev
  310 + cc(n,i) = nprof(i) !CHG3
  311 + cc(p1,i)= p1_init
  312 + cc(p2,i)= p2_init
  313 + cc(z1,i)= z1_init
  314 + cc(z2,i)= z2_init
  315 + cc(d,i) = d_init
  316 + cc(l,i) = l_init
  317 + cc(b1,i)= b_init
  318 + cc(a,i) = aprof(i) !CHG5
  319 + cc(hc,i)= hcprof(i)
  320 + cc(b2,i)= b_init
  321 + end do
  322 +
  323 + do i=0,nlev
  324 + ws(n,i) = _ZERO_
  325 + ws(p1,i) = w_p1
  326 + ws(p2,i) = w_p2
  327 + ws(z1,i) = _ZERO_
  328 + ws(z2,i) = _ZERO_
  329 + ws(d,i) = w_d
  330 + ws(l,i) = _ZERO_
  331 + ws(b1,i) = _ZERO_
  332 + ws(a,i) = _ZERO_
  333 + ws(hc,i) = _ZERO_
  334 + ws(b2,i) = _ZERO_
  335 + end do
  336 +
  337 + posconc(n) = 1
  338 + posconc(p1) = 1
  339 + posconc(p2) = 1
  340 + posconc(z1) = 1
  341 + posconc(z2) = 1
  342 + posconc(d) = 1
  343 + posconc(l) = 1
  344 + posconc(b1) = 1
  345 + posconc(a) = 1
  346 + posconc(hc) = 1
  347 + posconc(b2) = 1
  348 +
  349 + LEVEL3 'GSJ variables initialised ...'
  350 +
  351 + return
  352 +
  353 + end subroutine init_var_gsj
  354 +!EOC
  355 +
  356 +!-----------------------------------------------------------------------
  357 +!BOP
  358 +!
  359 +! !IROUTINE: Providing info on variables
  360 +!
  361 +! !INTERFACE:
  362 + subroutine var_info_gsj()
  363 +!
  364 +! !DESCRIPTION:
  365 +! This subroutine provides information about the variable names as they
  366 +! will be used when storing data in NetCDF files.
  367 +!
  368 +! !USES:
  369 + IMPLICIT NONE
  370 +!
  371 +! !REVISION HISTORY:
  372 +! Original author(s): Hans Burchard & Karsten Bolding
  373 +!
  374 +! !LOCAL VARIABLES:
  375 +!EOP
  376 +!-----------------------------------------------------------------------
  377 +!BOC
  378 + var_names(1) = 'nit'
  379 + var_units(1) = 'mmol/m**3'
  380 + var_long(1) = 'nitrate'
  381 +
  382 + var_names(2) = 'fla'
  383 + var_units(2) = 'mmol/m**3'
  384 + var_long(2) = 'flagellates'
  385 +
  386 + var_names(3) = 'dia'
  387 + var_units(3) = 'mmol/m**3'
  388 + var_long(3) = 'diatoms'
  389 +
  390 + var_names(4) = 'mcz'
  391 + var_units(4) = 'mmol/m**3'
  392 + var_long(4) = 'micro-zooplankton'
  393 +
  394 + var_names(5) = 'msz'
  395 + var_units(5) = 'mmol/m**3'
  396 + var_long(5) = 'meso-zooplankton'
  397 +
  398 + var_names(6) = 'det'
  399 + var_units(6) = 'mmol/m**3'
  400 + var_long(6) = 'detritus'
  401 +
  402 + var_names(7) = 'ldn'
  403 + var_units(7) = 'mmol/m**3'
  404 + var_long(7) = 'labile_dissolved_organic_nitrogen'
  405 +
  406 + var_names(8) = 'bac1'
  407 + var_units(8) = 'mmol/m**3'
  408 + var_long(8) = 'bacteria'
  409 +
  410 + var_names(9) = 'amm'
  411 + var_units(9) = 'mmol/m**3'
  412 + var_long(9) = 'ammonium'
  413 +
  414 + var_names(10)= 'hcb'
  415 + var_units(10)= 'mmol/m**3'
  416 + var_long(10) = 'hydrocarbon'
  417 +
  418 + var_names(11)= 'bac2'
  419 + var_units(11)= 'mmol/m**3'
  420 + var_long(11) = 'bacteria'
  421 +
  422 + return
  423 + end subroutine var_info_gsj
  424 +!EOC
  425 +
  426 +!-----------------------------------------------------------------------
  427 +!BOP
  428 +!
  429 +! !IROUTINE: Light properties for the ISMER model
  430 +!
  431 +! !INTERFACE:
  432 + subroutine light_gsj(nlev,bioshade_feedback)
  433 +!
  434 +! !DESCRIPTION:
  435 +! Here, the photosynthetically available radiation is calculated
  436 +! by simply assuming that the short wave part of the total
  437 +! radiation is available for photosynthesis.
  438 +! The photosynthetically
  439 +! available radiation, $I_{PAR}$, follows from (\ref{light}).
  440 +! The user should make
  441 +! sure that this is consistent with the light class given in the
  442 +! {\tt extinct} namelist of the {\tt obs.nml} file.
  443 +! The self-shading effect is also calculated in this subroutine,
  444 +! which may be used to consider the effect of bio-turbidity also
  445 +! in the temperature equation (if {\tt bioshade\_feedback} is set
  446 +! to true in {\tt bio.nml}).
  447 +! For details, see section \ref{sec:do-bio}.
  448 +!
  449 +! !USES:
  450 + IMPLICIT NONE
  451 +!
  452 +! !INPUT PARAMETERS:
  453 + integer, intent(in) :: nlev
  454 + logical, intent(in) :: bioshade_feedback
  455 +!
  456 +! !REVISION HISTORY:
  457 +! Original author(s): Hans Burchard, Karsten Bolding
  458 +!
  459 +! !LOCAL VARIABLES:
  460 + integer :: i
  461 + REALTYPE :: zz,add
  462 +!EOP
  463 +!-----------------------------------------------------------------------
  464 +!BOC
  465 + zz = _ZERO_
  466 + add = _ZERO_
  467 + do i=nlev,1,-1
  468 + add=add+0.5*h(i)*(cc(p1,i)+cc(p2,i)+2*p0)
  469 + zz=zz+0.5*h(i)
  470 + par(i)=rad(nlev)*(1.-aa)*exp(-zz/g2)*exp(-kc*add)
  471 + add=add+0.5*h(i)*(cc(p1,i)+cc(p2,i)+2*p0)
  472 + zz=zz+0.5*h(i)
  473 + if (bioshade_feedback) bioshade_(i)=exp(-kc*add)
  474 + end do
  475 +
  476 +
  477 + return
  478 + end subroutine light_gsj
  479 +!EOC
  480 +
  481 +!-----------------------------------------------------------------------
  482 +!BOP
  483 +!
  484 +! !IROUTINE: Right hand sides of geobiochemical model \label{sec:bio-fasham-rhs}
  485 +!
  486 +! !INTERFACE:
  487 + subroutine do_bio_gsj(first,numc,nlev,cc,pp,dd)
  488 +!
  489 +! !DESCRIPTION:
  490 +!
  491 +! The \cite{Fashametal1990} model consisting of the $I=7$
  492 +! state variables phytoplankton, bacteria, detritus, zooplankton,
  493 +! nitrate, ammonium and dissolved organic nitrogen is described here
  494 +! in detail.
  495 +!
  496 +! Phytoplankton mortality and zooplankton grazing loss of phytoplankton:
  497 +! \begin{equation}\label{d13}