Line data Source code
1 : !=======================================================================
2 : !
3 : ! Compute sea ice biogeochemistry (vertical or skeletal layer)
4 : !
5 : ! authors: Nicole Jeffery, LANL
6 : ! Scott Elliot, LANL
7 : ! Elizabeth C. Hunke, LANL
8 : !
9 : module icepack_algae
10 :
11 : use icepack_kinds
12 :
13 : use icepack_parameters, only: p05, p5, c0, c1, c2, c6, c10
14 : use icepack_parameters, only: pi, secday, puny
15 : use icepack_parameters, only: hs_ssl, sk_l
16 :
17 : use icepack_parameters, only: dEdd_algae, solve_zbgc
18 : use icepack_parameters, only: R_dFe2dust, dustFe_sol, algal_vel
19 : use icepack_parameters, only: bgc_flux_type
20 : use icepack_parameters, only: grid_o
21 : use icepack_parameters, only: T_max, fsal , fr_resp
22 : use icepack_parameters, only: op_dep_min , fr_graze_s
23 : use icepack_parameters, only: fr_graze_e , fr_mort2min
24 : use icepack_parameters, only: fr_dFe , k_nitrif
25 : use icepack_parameters, only: t_iron_conv , max_loss
26 : use icepack_parameters, only: max_dfe_doc1 , fr_resp_s
27 : use icepack_parameters, only: y_sk_DMS , t_sk_conv
28 : use icepack_parameters, only: t_sk_ox
29 :
30 : use icepack_tracers, only: ntrcr, bio_index
31 : use icepack_tracers, only: nt_bgc_N, nt_fbri, nt_zbgc_frac
32 : use icepack_tracers, only: tr_brine
33 : use icepack_tracers, only: tr_bgc_Nit, tr_bgc_Am, tr_bgc_Sil
34 : use icepack_tracers, only: tr_bgc_DMS, tr_bgc_PON
35 : use icepack_tracers, only: tr_bgc_N, tr_bgc_C, tr_bgc_chl
36 : use icepack_tracers, only: tr_bgc_DON, tr_bgc_Fe, tr_zaero
37 : use icepack_tracers, only: nlt_bgc_Nit, nlt_bgc_Am, nlt_bgc_Sil
38 : use icepack_tracers, only: nlt_bgc_DMS, nlt_bgc_PON
39 : use icepack_tracers, only: nlt_bgc_N, nlt_bgc_C, nlt_bgc_chl
40 : use icepack_tracers, only: nlt_bgc_DOC, nlt_bgc_DON, nlt_bgc_DIC
41 : use icepack_tracers, only: nlt_zaero , nlt_bgc_DMSPp,nlt_bgc_DMSPd
42 : use icepack_tracers, only: nlt_bgc_Fed, nlt_bgc_Fep
43 :
44 : use icepack_zbgc_shared, only: remap_zbgc, regrid_stationary
45 : use icepack_zbgc_shared, only: merge_bgc_fluxes
46 : use icepack_zbgc_shared, only: merge_bgc_fluxes_skl
47 : use icepack_zbgc_shared, only: phi_sk, bgc_tracer_type
48 : use icepack_zbgc_shared, only: zbgc_init_frac
49 : use icepack_zbgc_shared, only: zbgc_frac_init
50 : use icepack_zbgc_shared, only: tau_rel, tau_ret, thinS
51 : use icepack_zbgc_shared, only: r_Si2N, R_Fe2N, R_S2N, R_chl2N
52 : use icepack_zbgc_shared, only: chlabs, alpha2max_low, beta2max, mu_max
53 : use icepack_zbgc_shared, only: k_exude, K_Nit, K_Am, K_Sil, K_Fe
54 : use icepack_zbgc_shared, only: grow_Tdep, fr_graze, mort_pre, mort_Tdep
55 : use icepack_zbgc_shared, only: f_don, kn_bac, f_don_Am
56 : use icepack_zbgc_shared, only: f_doc, f_exude, k_bac, R_chl2N, R_C2N
57 :
58 : use icepack_warnings, only: warnstr, icepack_warnings_add
59 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
60 :
61 : use icepack_aerosol, only: update_snow_bgc
62 :
63 : implicit none
64 :
65 : private
66 : public :: zbio, sklbio
67 :
68 : real (kind=dbl_kind), parameter :: &
69 : exp_argmax = c10 ! maximum argument of exponential
70 :
71 : !=======================================================================
72 :
73 : contains
74 :
75 : !=======================================================================
76 :
77 0 : subroutine zbio (dt, nblyr, &
78 : nslyr, nilyr, & ! LCOV_EXCL_LINE
79 : meltt, melts, & ! LCOV_EXCL_LINE
80 : meltb, congel, & ! LCOV_EXCL_LINE
81 : snoice, nbtrcr, & ! LCOV_EXCL_LINE
82 : fsnow, ntrcr, & ! LCOV_EXCL_LINE
83 : trcrn, bio_index, & ! LCOV_EXCL_LINE
84 : aice_old, & ! LCOV_EXCL_LINE
85 : vice_old, vsno_old, & ! LCOV_EXCL_LINE
86 : vicen, vsnon, & ! LCOV_EXCL_LINE
87 : aicen, flux_bio_atm,& ! LCOV_EXCL_LINE
88 : n_cat, n_algae, & ! LCOV_EXCL_LINE
89 : n_doc, n_dic, & ! LCOV_EXCL_LINE
90 : n_don, & ! LCOV_EXCL_LINE
91 : n_fed, n_fep, & ! LCOV_EXCL_LINE
92 : n_zaero, first_ice, & ! LCOV_EXCL_LINE
93 : hice_old, ocean_bio, & ! LCOV_EXCL_LINE
94 : bphin, iphin, & ! LCOV_EXCL_LINE
95 : iDin, & ! LCOV_EXCL_LINE
96 : fswthrul, & ! LCOV_EXCL_LINE
97 : dh_top, dh_bot, & ! LCOV_EXCL_LINE
98 : zfswin, & ! LCOV_EXCL_LINE
99 : hbri, hbri_old, & ! LCOV_EXCL_LINE
100 : ! darcy_V, darcy_V_chl, & ! LCOV_EXCL_LINE
101 : darcy_V, & ! LCOV_EXCL_LINE
102 : bgrid, & ! LCOV_EXCL_LINE
103 : igrid, icgrid, & ! LCOV_EXCL_LINE
104 : bphi_min, & ! LCOV_EXCL_LINE
105 : iTin, & ! LCOV_EXCL_LINE
106 : Zoo, & ! LCOV_EXCL_LINE
107 : flux_bio, dh_direct, & ! LCOV_EXCL_LINE
108 : upNO, upNH, & ! LCOV_EXCL_LINE
109 : fbio_snoice, fbio_atmice, & ! LCOV_EXCL_LINE
110 : PP_net, ice_bio_net, & ! LCOV_EXCL_LINE
111 0 : snow_bio_net, grow_net )
112 :
113 : integer (kind=int_kind), intent(in) :: &
114 : nblyr, & ! number of bio layers ! LCOV_EXCL_LINE
115 : nslyr, & ! number of snow layers ! LCOV_EXCL_LINE
116 : nilyr, & ! number of ice layers ! LCOV_EXCL_LINE
117 : nbtrcr, & ! number of distinct bio tracers ! LCOV_EXCL_LINE
118 : n_cat, & ! category number ! LCOV_EXCL_LINE
119 : n_algae, & ! number of autotrophs ! LCOV_EXCL_LINE
120 : n_zaero, & ! number of aerosols ! LCOV_EXCL_LINE
121 : n_doc, n_dic, n_don, n_fed, n_fep, & ! LCOV_EXCL_LINE
122 : ntrcr ! number of tracers
123 :
124 : integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
125 : bio_index
126 :
127 : real (kind=dbl_kind), intent(in) :: &
128 : dt, & ! time step ! LCOV_EXCL_LINE
129 : hbri, & ! brine height (m) ! LCOV_EXCL_LINE
130 : bphi_min, & ! surface porosity ! LCOV_EXCL_LINE
131 : meltt, & ! thermodynamic melt/growth rates in dt (m) ! LCOV_EXCL_LINE
132 : melts, & ! LCOV_EXCL_LINE
133 : meltb, & ! LCOV_EXCL_LINE
134 : congel, & ! LCOV_EXCL_LINE
135 : snoice, & ! LCOV_EXCL_LINE
136 : fsnow, & ! snowfall rate (kg/m^2 s) ! LCOV_EXCL_LINE
137 : hice_old, & ! ice height (m) ! LCOV_EXCL_LINE
138 : vicen, & ! ice volume (m) ! LCOV_EXCL_LINE
139 : vsnon, & ! snow volume (m) ! LCOV_EXCL_LINE
140 : aicen, & ! ice area fraction ! LCOV_EXCL_LINE
141 : aice_old, & ! values prior to thermodynamic changes ! LCOV_EXCL_LINE
142 : vice_old, & ! LCOV_EXCL_LINE
143 : vsno_old, & ! LCOV_EXCL_LINE
144 : darcy_V, & ! darcy velocity ! LCOV_EXCL_LINE
145 : ! darcy_V_chl,& ! darcy velocity for algae ! LCOV_EXCL_LINE
146 : dh_bot, & ! change in brine bottom (m) ! LCOV_EXCL_LINE
147 : dh_top, & ! change in brine top (m) ! LCOV_EXCL_LINE
148 : dh_direct ! surface flooding or surface runoff (m)
149 :
150 : real (kind=dbl_kind), dimension (nbtrcr), intent(inout) :: &
151 : snow_bio_net,& ! net bio tracer in snow (mmol/m^2) ! LCOV_EXCL_LINE
152 : ice_bio_net, & ! net bio tracer in ice (mmol/m^2) ! LCOV_EXCL_LINE
153 : fbio_atmice, & ! bio flux from atm to ice (mmol/m^2/s) ! LCOV_EXCL_LINE
154 : fbio_snoice, & ! bio flux from snow to ice (mmol/m^2/s) ! LCOV_EXCL_LINE
155 : flux_bio ! total ocean tracer flux (mmol/m^2/s)
156 :
157 : real (kind=dbl_kind), intent(inout) :: &
158 : hbri_old ! brine height (m)
159 :
160 : real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
161 : bgrid ! biology nondimensional vertical grid points
162 :
163 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
164 : igrid , & ! biology vertical interface points ! LCOV_EXCL_LINE
165 : iTin , & ! salinity vertical interface points ! LCOV_EXCL_LINE
166 : iphin , & ! Porosity on the igrid ! LCOV_EXCL_LINE
167 : iDin ! Diffusivity/h on the igrid (1/s)
168 :
169 : real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
170 : icgrid , & ! CICE interface coordinate ! LCOV_EXCL_LINE
171 : fswthrul ! visible short wave radiation on icgrid (W/m^2)
172 :
173 : real (kind=dbl_kind), dimension(nbtrcr), intent(in) :: &
174 : flux_bio_atm ! aerosol/bgc deposition rate (mmol/m^2 s)
175 :
176 : real (kind=dbl_kind), dimension(ntrcr), intent(inout) :: &
177 : trcrn
178 :
179 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
180 : zfswin ! visible Short wave flux on igrid (W/m^2)
181 :
182 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
183 : Zoo ! N losses to the system from reaction terms
184 : ! (ie. zooplankton/bacteria) (mmol/m^3)
185 :
186 : real (kind=dbl_kind), dimension (nbtrcr), intent(in) :: &
187 : !change to inout when updating ocean fields
188 : ocean_bio ! ocean concentrations (mmol/m^3)
189 :
190 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
191 : bphin ! Porosity on the bgrid
192 :
193 : real (kind=dbl_kind), intent(inout):: &
194 : PP_net , & ! net PP (mg C/m^2/d) times aice ! LCOV_EXCL_LINE
195 : grow_net , & ! net specific growth (m/d) times vice ! LCOV_EXCL_LINE
196 : upNO , & ! tot nitrate uptake rate (mmol/m^2/d) times aice ! LCOV_EXCL_LINE
197 : upNH ! tot ammonium uptake rate (mmol/m^2/d) times aice
198 :
199 : logical (kind=log_kind), intent(in) :: &
200 : first_ice ! initialized values should be used
201 :
202 : ! local variables
203 :
204 : integer (kind=int_kind) :: &
205 : mm ! thickness category index
206 :
207 : real (kind=dbl_kind), dimension (nblyr+1,n_algae) :: &
208 : upNOn , & ! algal nitrate uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
209 : upNHn , & ! algal ammonium uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
210 0 : grow_alg ! algal growth rate (mmol/m^3/s)
211 :
212 : real (kind=dbl_kind), dimension (nbtrcr) :: &
213 0 : flux_bion !tracer flux to ocean
214 :
215 : real (kind=dbl_kind),dimension(nbtrcr) :: &
216 : zbgc_snown, & ! aerosol contribution from snow to ice ! LCOV_EXCL_LINE
217 0 : zbgc_atmn ! and atm to ice concentration * volume (mmol/m^3*m)
218 :
219 : real (kind=dbl_kind), dimension(nbtrcr) :: &
220 : Tot_BGC_i, & ! initial column sum, ice + snow, of BGC tracer (mmol/m^2) ! LCOV_EXCL_LINE
221 : Tot_BGC_f, & ! final column sum ! LCOV_EXCL_LINE
222 0 : flux_bio_sno !
223 :
224 : real (kind=dbl_kind) :: &
225 : hsnow_i, & ! initial snow thickness (m) ! LCOV_EXCL_LINE
226 0 : hsnow_f ! final snow thickness (m)
227 :
228 : logical (kind=log_kind) :: &
229 : write_flux_diag
230 :
231 : real (kind=dbl_kind) :: &
232 0 : a_ice
233 :
234 : character(len=*),parameter :: subname='(zbio)'
235 :
236 0 : zbgc_snown(:) = c0
237 0 : zbgc_atmn (:) = c0
238 0 : flux_bion (:) = c0
239 0 : flux_bio_sno(:) = c0
240 0 : Tot_BGC_i (:) = c0
241 0 : Tot_BGC_f (:) = c0
242 0 : hsnow_i = c0
243 0 : hsnow_f = c0
244 0 : write_flux_diag = .false.
245 :
246 0 : if (write_flux_diag) then
247 0 : if (aice_old > c0) then
248 0 : hsnow_i = vsno_old/aice_old
249 0 : do mm = 1,nbtrcr
250 : call bgc_column_sum (nblyr, nslyr, hsnow_i, hbri_old, &
251 : trcrn(bio_index(mm):bio_index(mm)+nblyr+2), & ! LCOV_EXCL_LINE
252 0 : Tot_BGC_i(mm))
253 0 : if (icepack_warnings_aborted(subname)) return
254 : enddo
255 : endif
256 : endif
257 :
258 : call update_snow_bgc (dt, nblyr, &
259 : nslyr, & ! LCOV_EXCL_LINE
260 : meltt, melts, & ! LCOV_EXCL_LINE
261 : meltb, congel, & ! LCOV_EXCL_LINE
262 : snoice, nbtrcr, & ! LCOV_EXCL_LINE
263 : fsnow, ntrcr, & ! LCOV_EXCL_LINE
264 : trcrn, bio_index, & ! LCOV_EXCL_LINE
265 : aice_old, zbgc_snown, & ! LCOV_EXCL_LINE
266 : vice_old, vsno_old, & ! LCOV_EXCL_LINE
267 : vicen, vsnon, & ! LCOV_EXCL_LINE
268 : aicen, flux_bio_atm, & ! LCOV_EXCL_LINE
269 0 : zbgc_atmn, flux_bio_sno)
270 0 : if (icepack_warnings_aborted(subname)) return
271 :
272 : call z_biogeochemistry (n_cat, dt, &
273 : nilyr, & ! LCOV_EXCL_LINE
274 : nblyr, nbtrcr, & ! LCOV_EXCL_LINE
275 : n_algae, n_doc, & ! LCOV_EXCL_LINE
276 : n_dic, n_don, & ! LCOV_EXCL_LINE
277 : n_fed, n_fep, & ! LCOV_EXCL_LINE
278 : n_zaero, first_ice, & ! LCOV_EXCL_LINE
279 : aicen, vicen, & ! LCOV_EXCL_LINE
280 : hice_old, ocean_bio, & ! LCOV_EXCL_LINE
281 : flux_bion, bphin, & ! LCOV_EXCL_LINE
282 : iphin, trcrn, & ! LCOV_EXCL_LINE
283 : iDin, & ! LCOV_EXCL_LINE
284 : fswthrul, grow_alg, & ! LCOV_EXCL_LINE
285 : upNOn, upNHn, & ! LCOV_EXCL_LINE
286 : dh_top, dh_bot, & ! LCOV_EXCL_LINE
287 : zfswin, hbri, & ! LCOV_EXCL_LINE
288 : hbri_old, darcy_V, & ! LCOV_EXCL_LINE
289 : ! darcy_V_chl, bgrid, & ! LCOV_EXCL_LINE
290 : bgrid, & ! LCOV_EXCL_LINE
291 : igrid, icgrid, & ! LCOV_EXCL_LINE
292 : bphi_min, zbgc_snown,& ! LCOV_EXCL_LINE
293 : zbgc_atmn, & ! LCOV_EXCL_LINE
294 : iTin, dh_direct, & ! LCOV_EXCL_LINE
295 : Zoo, meltb, & ! LCOV_EXCL_LINE
296 0 : congel )
297 0 : if (icepack_warnings_aborted(subname)) return
298 :
299 0 : do mm = 1,nbtrcr
300 0 : flux_bion(mm) = flux_bion(mm) + flux_bio_sno(mm)
301 : enddo
302 :
303 0 : if (write_flux_diag) then
304 0 : if (aicen > c0) then
305 0 : hsnow_f = vsnon/aicen
306 0 : do mm = 1,nbtrcr
307 : call bgc_column_sum (nblyr, nslyr, hsnow_f, hbri, &
308 : trcrn(bio_index(mm):bio_index(mm)+nblyr+2), & ! LCOV_EXCL_LINE
309 0 : Tot_BGC_f(mm))
310 0 : write(warnstr,*) subname, 'mm, Tot_BGC_i(mm), Tot_BGC_f(mm)'
311 0 : call icepack_warnings_add(warnstr)
312 0 : write(warnstr,*) subname, mm, Tot_BGC_i(mm), Tot_BGC_f(mm)
313 0 : call icepack_warnings_add(warnstr)
314 0 : write(warnstr,*) subname, 'flux_bion(mm), flux_bio_atm(mm)'
315 0 : call icepack_warnings_add(warnstr)
316 0 : write(warnstr,*) subname, flux_bion(mm), flux_bio_atm(mm)
317 0 : call icepack_warnings_add(warnstr)
318 0 : write(warnstr,*) subname, 'zbgc_snown(mm),zbgc_atmn(mm)'
319 0 : call icepack_warnings_add(warnstr)
320 0 : write(warnstr,*) subname, zbgc_snown(mm),zbgc_atmn(mm)
321 0 : call icepack_warnings_add(warnstr)
322 0 : write(warnstr,*) subname, 'Tot_BGC_i(mm) + flux_bio_atm(mm)*dt - flux_bion(mm)*dt'
323 0 : call icepack_warnings_add(warnstr)
324 0 : write(warnstr,*) subname, Tot_BGC_i(mm) + flux_bio_atm(mm)*dt - flux_bion(mm)*dt
325 0 : call icepack_warnings_add(warnstr)
326 : enddo
327 : endif
328 : endif
329 :
330 0 : if (icepack_warnings_aborted(subname)) return
331 :
332 : call merge_bgc_fluxes (dt, nblyr, &
333 : nslyr, & ! LCOV_EXCL_LINE
334 : bio_index, n_algae, & ! LCOV_EXCL_LINE
335 : nbtrcr, aicen, & ! LCOV_EXCL_LINE
336 : vicen, vsnon, & ! LCOV_EXCL_LINE
337 : iphin, & ! LCOV_EXCL_LINE
338 : trcrn, & ! LCOV_EXCL_LINE
339 : flux_bion, flux_bio, & ! LCOV_EXCL_LINE
340 : upNOn, upNHn, & ! LCOV_EXCL_LINE
341 : upNO, upNH, & ! LCOV_EXCL_LINE
342 : zbgc_snown, zbgc_atmn, & ! LCOV_EXCL_LINE
343 : fbio_snoice, fbio_atmice,& ! LCOV_EXCL_LINE
344 : PP_net, ice_bio_net,& ! LCOV_EXCL_LINE
345 : snow_bio_net, grow_alg, & ! LCOV_EXCL_LINE
346 0 : grow_net)
347 0 : if (icepack_warnings_aborted(subname)) return
348 :
349 0 : if (write_flux_diag) then
350 0 : if (aicen > c0) then
351 0 : if (n_cat .eq. 1) a_ice = c0
352 0 : a_ice = a_ice + aicen
353 0 : write(warnstr,*) subname, 'after merge_bgc_fluxes, n_cat:', n_cat
354 0 : call icepack_warnings_add(warnstr)
355 0 : do mm = 1,nbtrcr
356 0 : write(warnstr,*) subname, 'mm, flux_bio(mm):',mm,flux_bio(mm)
357 0 : call icepack_warnings_add(warnstr)
358 0 : write(warnstr,*) subname, 'fbio_snoice(mm)',fbio_snoice(mm)
359 0 : call icepack_warnings_add(warnstr)
360 0 : write(warnstr,*) subname, 'fbio_atmice(mm)',fbio_atmice(mm)
361 0 : call icepack_warnings_add(warnstr)
362 0 : write(warnstr,*) subname, 'flux_bio_atm(mm)', flux_bio_atm(mm)
363 0 : call icepack_warnings_add(warnstr)
364 0 : write(warnstr,*) subname, 'flux_bio_atm(mm)*a_ice', flux_bio_atm(mm)*a_ice
365 0 : call icepack_warnings_add(warnstr)
366 : enddo
367 : endif
368 : endif
369 :
370 : end subroutine zbio
371 :
372 : !=======================================================================
373 :
374 0 : subroutine sklbio (dt, ntrcr, &
375 : nbtrcr, n_algae, & ! LCOV_EXCL_LINE
376 : n_doc, & ! LCOV_EXCL_LINE
377 : n_dic, n_don, & ! LCOV_EXCL_LINE
378 : n_fed, n_fep, & ! LCOV_EXCL_LINE
379 : flux_bio, ocean_bio, & ! LCOV_EXCL_LINE
380 : aicen, & ! LCOV_EXCL_LINE
381 : meltb, congel, & ! LCOV_EXCL_LINE
382 : fswthru, first_ice, & ! LCOV_EXCL_LINE
383 : trcrn, & ! LCOV_EXCL_LINE
384 : PP_net, upNO, & ! LCOV_EXCL_LINE
385 : upNH, grow_net )
386 :
387 : integer (kind=int_kind), intent(in) :: &
388 : nbtrcr, & ! number of distinct bio tracers ! LCOV_EXCL_LINE
389 : n_algae, & ! number of autotrophs ! LCOV_EXCL_LINE
390 : n_doc, n_dic, & ! number of dissolved organic, inorganic carbon ! LCOV_EXCL_LINE
391 : n_don, & ! number of dissolved organic nitrogen ! LCOV_EXCL_LINE
392 : n_fed, n_fep, & ! number of iron ! LCOV_EXCL_LINE
393 : ntrcr ! number of tracers
394 :
395 : logical (kind=log_kind), intent(in) :: &
396 : first_ice ! initialized values should be used
397 :
398 : real (kind=dbl_kind), intent(in) :: &
399 : dt, & ! time step ! LCOV_EXCL_LINE
400 : ! hmix, & ! mixed layer depth (m) ! LCOV_EXCL_LINE
401 : aicen, & ! ice area fraction ! LCOV_EXCL_LINE
402 : meltb, & ! bottom melt (m) ! LCOV_EXCL_LINE
403 : congel, & ! bottom growth (m) ! LCOV_EXCL_LINE
404 : fswthru ! visible shortwave passing to ocean(W/m^2)
405 :
406 : real (kind=dbl_kind), dimension(ntrcr), intent(inout) :: &
407 : trcrn ! bulk concentration per m^3
408 :
409 : real (kind=dbl_kind), dimension (nbtrcr), intent(inout) :: &
410 : flux_bio ! ocean tracer flux (mmol/m^2/s) positive into ocean
411 :
412 : real (kind=dbl_kind), dimension (nbtrcr), intent(in) :: &
413 : ocean_bio ! ocean tracer concentration (mmol/m^3)
414 :
415 : ! history output
416 : real (kind=dbl_kind), intent(inout):: &
417 : PP_net , & ! Bulk net PP (mg C/m^2/s) ! LCOV_EXCL_LINE
418 : grow_net, & ! net specific growth (/s) ! LCOV_EXCL_LINE
419 : upNO , & ! tot nitrate uptake rate (mmol/m^2/s) ! LCOV_EXCL_LINE
420 : upNH ! tot ammonium uptake rate (mmol/m^2/s)
421 :
422 : ! local variables
423 :
424 : real (kind=dbl_kind), dimension (n_algae) :: &
425 : upNOn , & ! algal nitrate uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
426 : upNHn , & ! algal ammonium uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
427 0 : grow_alg ! algal growth rate (mmol/m^3/s)
428 :
429 : real (kind=dbl_kind), dimension (nbtrcr) :: &
430 0 : flux_bion !tracer flux to ocean
431 :
432 : character(len=*),parameter :: subname='(sklbio)'
433 :
434 0 : flux_bion (:) = c0
435 0 : upNOn (:) = c0
436 0 : upNHn (:) = c0
437 0 : grow_alg (:) = c0
438 :
439 : call skl_biogeochemistry (dt, &
440 : n_doc, & ! LCOV_EXCL_LINE
441 : n_dic, n_don, & ! LCOV_EXCL_LINE
442 : n_fed, n_fep, & ! LCOV_EXCL_LINE
443 : nbtrcr, n_algae, & ! LCOV_EXCL_LINE
444 : flux_bion, ocean_bio, & ! LCOV_EXCL_LINE
445 : ! hmix, aicen, & ! LCOV_EXCL_LINE
446 : meltb, congel, & ! LCOV_EXCL_LINE
447 : fswthru, first_ice, & ! LCOV_EXCL_LINE
448 : trcrn, upNOn, & ! LCOV_EXCL_LINE
449 0 : upNHn, grow_alg)
450 0 : if (icepack_warnings_aborted(subname)) return
451 :
452 : call merge_bgc_fluxes_skl (nbtrcr, n_algae, &
453 : aicen, trcrn, & ! LCOV_EXCL_LINE
454 : flux_bion, flux_bio, & ! LCOV_EXCL_LINE
455 : PP_net, upNOn, & ! LCOV_EXCL_LINE
456 : upNHn, upNO, & ! LCOV_EXCL_LINE
457 : upNH, grow_net, & ! LCOV_EXCL_LINE
458 0 : grow_alg)
459 0 : if (icepack_warnings_aborted(subname)) return
460 :
461 : end subroutine sklbio
462 :
463 : !=======================================================================
464 : !
465 : ! skeletal layer biochemistry
466 : !
467 0 : subroutine skl_biogeochemistry (dt, &
468 : n_doc, & ! LCOV_EXCL_LINE
469 : n_dic, n_don, & ! LCOV_EXCL_LINE
470 : n_fed, n_fep, & ! LCOV_EXCL_LINE
471 : nbtrcr, n_algae, & ! LCOV_EXCL_LINE
472 : flux_bio, ocean_bio, & ! LCOV_EXCL_LINE
473 : ! hmix, aicen, & ! LCOV_EXCL_LINE
474 : meltb, congel, & ! LCOV_EXCL_LINE
475 : fswthru, first_ice, & ! LCOV_EXCL_LINE
476 : trcrn, upNOn, & ! LCOV_EXCL_LINE
477 0 : upNHn, grow_alg_skl)
478 :
479 : integer (kind=int_kind), intent(in) :: &
480 : n_doc, n_dic, n_don, n_fed, n_fep, & ! LCOV_EXCL_LINE
481 : nbtrcr , n_algae ! number of bgc tracers and number algae
482 :
483 : real (kind=dbl_kind), intent(in) :: &
484 : dt , & ! time step ! LCOV_EXCL_LINE
485 : ! hmix , & ! mixed layer depth ! LCOV_EXCL_LINE
486 : ! aicen , & ! ice area ! LCOV_EXCL_LINE
487 : meltb , & ! bottom ice melt ! LCOV_EXCL_LINE
488 : congel , & ! bottom ice growth ! LCOV_EXCL_LINE
489 : fswthru ! shortwave passing through ice to ocean
490 :
491 : logical (kind=log_kind), intent(in) :: &
492 : first_ice ! initialized values should be used
493 :
494 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
495 : trcrn ! bulk concentration per m^3
496 :
497 : ! history variables
498 :
499 : real (kind=dbl_kind), dimension (:), intent(out) :: &
500 : flux_bio ! ocean tracer flux (mmol/m^2/s) positive into ocean
501 :
502 : real (kind=dbl_kind), dimension (:), intent(in) :: &
503 : ocean_bio ! ocean tracer concentration (mmol/m^3)
504 :
505 : real (kind=dbl_kind), dimension (:), intent(out) :: &
506 : grow_alg_skl, & ! tot algal growth rate (mmol/m^3/s) ! LCOV_EXCL_LINE
507 : upNOn , & ! algal NO uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
508 : upNHn ! algal NH uptake rate (mmol/m^3/s)
509 :
510 : ! local variables
511 :
512 : integer (kind=int_kind) :: nn
513 :
514 : real (kind=dbl_kind), dimension(nbtrcr):: &
515 : react , & ! biological sources and sinks (mmol/m^3) ! LCOV_EXCL_LINE
516 : cinit , & ! initial brine concentration*sk_l (mmol/m^2) ! LCOV_EXCL_LINE
517 : cinit_v , & ! initial brine concentration (mmol/m^3) ! LCOV_EXCL_LINE
518 : congel_alg , & ! congelation flux contribution to ice algae (mmol/m^2 s) ! LCOV_EXCL_LINE
519 : ! (used as initialization)
520 0 : f_meltn , & ! vertical melt fraction of skeletal layer in dt
521 : flux_bio_temp, & ! tracer flux to ocean (mmol/m^2 s) ! LCOV_EXCL_LINE
522 : PVflag , & ! 1 for tracers that flow with the brine, 0 otherwise ! LCOV_EXCL_LINE
523 0 : cling ! 1 for tracers that cling, 0 otherwise
524 :
525 : real (kind=dbl_kind) :: &
526 : Zoo_skl , & ! N losses from zooplankton/bacteria ... (mmol/m^3) ! LCOV_EXCL_LINE
527 : iTin , & ! LCOV_EXCL_LINE
528 : PVt , & ! type 'Jin2006' piston velocity (m/s) ! LCOV_EXCL_LINE
529 : ice_growth , & ! Jin2006 definition: either congel rate or bottom melt rate (m/s) ! LCOV_EXCL_LINE
530 : grow_val , & ! (m/x) ! LCOV_EXCL_LINE
531 : rphi_sk , & ! 1 / skeletal layer porosity ! LCOV_EXCL_LINE
532 : cinit_tmp , & ! temporary variable for concentration (mmol/m^2) ! LCOV_EXCL_LINE
533 0 : Nerror ! change in total nitrogen from reactions
534 :
535 : real (kind=dbl_kind), parameter :: &
536 : PVc = 1.e-6_dbl_kind , & ! type 'constant' piston velocity for interface (m/s) ! LCOV_EXCL_LINE
537 : PV_scale_growth = p5 , & ! scale factor in Jin code PV during ice growth ! LCOV_EXCL_LINE
538 : PV_scale_melt = p05 , & ! scale factor in Jin code PV during ice melt ! LCOV_EXCL_LINE
539 : growth_max = 1.85e-5_dbl_kind , & ! PVt function reaches maximum here. (m/s) ! LCOV_EXCL_LINE
540 : Tin_bot = -1.8_dbl_kind , & ! temperature of the ice bottom (oC) ! LCOV_EXCL_LINE
541 : MJ1 = 9.667e-9_dbl_kind , & ! (m/s) coefficients in Jin2008 ! LCOV_EXCL_LINE
542 : MJ2 = 38.8_dbl_kind , & ! (1) from:4.49e-4_dbl_kind*secday ! LCOV_EXCL_LINE
543 : MJ3 = 1.04e7_dbl_kind , & ! 1/(m/s) from: 1.39e-3_dbl_kind*secday^2 ! LCOV_EXCL_LINE
544 : PV_frac_max = 0.9_dbl_kind ! Maximum Piston velocity is 90% of skeletal layer/dt
545 :
546 : logical (kind=log_kind) :: conserve_N
547 :
548 : character(len=*),parameter :: subname='(skl_biogeochemistry)'
549 :
550 : !-----------------------------------------------------------------
551 : ! Initialize
552 : !-----------------------------------------------------------------
553 :
554 0 : conserve_N = .true.
555 0 : Zoo_skl = c0
556 0 : rphi_sk = c1/phi_sk
557 0 : PVt = c0
558 0 : iTin = Tin_bot
559 :
560 0 : do nn = 1, nbtrcr
561 0 : cinit (nn) = c0
562 0 : cinit_v (nn) = c0
563 0 : congel_alg(nn) = c0
564 0 : f_meltn (nn) = c0
565 0 : react (nn) = c0
566 0 : PVflag (nn) = c1
567 0 : cling (nn) = c0
568 :
569 : !-----------------------------------------------------------------
570 : ! only the dominant tracer_type affects behavior
571 : ! < 0 is purely mobile: > 0 stationary behavior
572 : ! NOTE: retention times are not used in skl model
573 : !-----------------------------------------------------------------
574 :
575 0 : if (bgc_tracer_type(nn) >= c0) then
576 0 : PVflag(nn) = c0
577 0 : cling (nn) = c1
578 : endif
579 :
580 0 : ice_growth = (congel-meltb)/dt
581 0 : if (first_ice) then
582 0 : trcrn(bio_index(nn)) = ocean_bio(nn) ! * sk_l*rphi_sk
583 : endif ! first_ice
584 0 : cinit (nn) = trcrn(bio_index(nn)) * sk_l * rphi_sk
585 0 : cinit_v(nn) = cinit(nn)/sk_l
586 0 : if (cinit(nn) < c0) then
587 0 : write(warnstr,*) subname,'initial sk_bgc < 0, nn,nbtrcr,cinit(nn)', &
588 0 : nn,nbtrcr,cinit(nn)
589 0 : call icepack_warnings_add(warnstr)
590 0 : call icepack_warnings_add(subname//' cinit < c0')
591 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
592 0 : return
593 : endif
594 : enddo ! nbtrcr
595 :
596 0 : if (icepack_warnings_aborted(subname)) return
597 :
598 0 : if (trim(bgc_flux_type) == 'Jin2006') then
599 :
600 : !-----------------------------------------------------------------
601 : ! 'Jin2006':
602 : ! 1. congel/melt dependent piston velocity (PV) for growth and melt
603 : ! 2. If congel > melt use 'congel'; if melt > congel use 'melt'
604 : ! 3. For algal N, PV for ice growth only provides a seeding concentration
605 : ! 4. Melt affects nutrients and algae in the same manner through PV(melt)
606 : !-----------------------------------------------------------------
607 :
608 0 : if (ice_growth > c0) then ! ice_growth = congel/dt
609 0 : grow_val = min(ice_growth,growth_max)
610 : PVt = -min(abs(PV_scale_growth*(MJ1 + MJ2*grow_val &
611 : - MJ3*grow_val**2)), & ! LCOV_EXCL_LINE
612 0 : PV_frac_max*sk_l/dt)
613 : else ! ice_growth = -meltb/dt
614 : PVt = min(abs(PV_scale_melt *( MJ2*ice_growth &
615 : - MJ3*ice_growth**2)), & ! LCOV_EXCL_LINE
616 0 : PV_frac_max*sk_l/dt)
617 : endif
618 0 : do nn = 1, nbtrcr
619 0 : if (bgc_tracer_type(nn) >= c0) then
620 0 : if (ice_growth < c0) then ! flux from ice to ocean
621 : ! Algae and clinging tracers melt like nutrients
622 0 : f_meltn(nn) = PVt*cinit_v(nn) ! for algae only
623 0 : elseif (ice_growth > c0 .AND. &
624 0 : cinit(nn) < ocean_bio(nn)*sk_l/phi_sk) then
625 : ! Growth only contributes to seeding from ocean
626 0 : congel_alg(nn) = (ocean_bio(nn)*sk_l/phi_sk - cinit(nn))/dt
627 : endif ! PVt > c0
628 : endif
629 : enddo
630 :
631 : else ! bgc_flux_type = 'constant'
632 :
633 : !-----------------------------------------------------------------
634 : ! 'constant':
635 : ! 1. Constant PV for congel > melt
636 : ! 2. For algae, PV for ice growth only provides a seeding concentration
637 : ! 3. Melt loss (f_meltn) affects algae only and is proportional to melt
638 : !-----------------------------------------------------------------
639 :
640 0 : if (ice_growth > c0) PVt = -PVc
641 0 : do nn = 1, nbtrcr
642 0 : if (bgc_tracer_type(nn) >= c0 ) then
643 0 : if (ice_growth >= c0 .AND. cinit_v(nn) < ocean_bio(nn)/phi_sk) then
644 0 : congel_alg(nn) = (ocean_bio(nn)*sk_l/phi_sk - cinit(nn))/dt
645 0 : elseif (ice_growth < c0) then
646 0 : f_meltn(nn) = min(c1, meltb/sk_l)*cinit(nn)/dt
647 : endif
648 : endif
649 : enddo ! nn
650 :
651 : endif ! bgc_flux_type
652 :
653 : !-----------------------------------------------------------------
654 : ! begin building biogeochemistry terms
655 : !-----------------------------------------------------------------
656 :
657 0 : react(:) = c0
658 0 : grow_alg_skl(:) = c0
659 :
660 : call algal_dyn (dt, &
661 : n_doc, n_dic, n_don, n_fed, n_fep, & ! LCOV_EXCL_LINE
662 : dEdd_algae, & ! LCOV_EXCL_LINE
663 : fswthru, react, & ! LCOV_EXCL_LINE
664 : cinit_v, & ! LCOV_EXCL_LINE
665 : grow_alg_skl(:), n_algae, & ! LCOV_EXCL_LINE
666 : iTin, & ! LCOV_EXCL_LINE
667 : upNOn(:), upNHn(:), & ! LCOV_EXCL_LINE
668 : Zoo_skl, & ! LCOV_EXCL_LINE
669 0 : Nerror, conserve_N)
670 0 : if (icepack_warnings_aborted(subname)) return
671 :
672 : !-----------------------------------------------------------------
673 : ! compute new tracer concencentrations
674 : !-----------------------------------------------------------------
675 :
676 0 : do nn = 1, nbtrcr
677 :
678 : !-----------------------------------------------------------------
679 : ! if PVt > 0, ie melt, then ocean_bio term drops out (MJ2006)
680 : ! Combine boundary fluxes
681 : !-----------------------------------------------------------------
682 :
683 0 : PVflag(nn) = SIGN(PVflag(nn),PVt)
684 0 : cinit_tmp = max(c0, cinit_v(nn) + react(nn))
685 0 : flux_bio_temp(nn) = (PVflag(nn)*PVt*cinit_tmp &
686 : - PVflag(nn)*min(c0,PVt)*ocean_bio(nn)) & ! LCOV_EXCL_LINE
687 0 : + f_meltn(nn)*cling(nn) - congel_alg(nn)
688 :
689 0 : if (cinit_tmp*sk_l < flux_bio_temp(nn)*dt) then
690 0 : flux_bio_temp(nn) = cinit_tmp*sk_l/dt*(c1-puny)
691 : endif
692 :
693 0 : cinit(nn) = cinit_tmp*sk_l - flux_bio_temp(nn)*dt
694 0 : flux_bio(nn) = flux_bio(nn) + flux_bio_temp(nn)*phi_sk
695 :
696 : ! Uncomment to update ocean concentration
697 : ! Currently not coupled with ocean biogeochemistry
698 : ! ocean_bio(nn) = ocean_bio(nn) + flux_bio(nn)/hmix*aicen
699 :
700 0 : if (.not. conserve_N) then
701 0 : write(warnstr,*) subname, 'N not conserved in skl_bgc, Nerror:',Nerror
702 0 : call icepack_warnings_add(warnstr)
703 0 : write(warnstr,*) subname, 'sk_bgc < 0 after algal fluxes, nn,cinit,flux_bio',&
704 0 : nn,cinit(nn),flux_bio(nn)
705 0 : call icepack_warnings_add(warnstr)
706 0 : write(warnstr,*) subname, 'cinit_tmp,flux_bio_temp,f_meltn,congel_alg,PVt,PVflag: '
707 0 : call icepack_warnings_add(warnstr)
708 0 : write(warnstr,*) subname, cinit_tmp,flux_bio_temp(nn),f_meltn(nn), &
709 0 : congel_alg(nn),PVt,PVflag(nn)
710 0 : call icepack_warnings_add(warnstr)
711 0 : write(warnstr,*) subname, 'congel, meltb: ',congel,meltb
712 0 : call icepack_warnings_add(warnstr)
713 0 : call icepack_warnings_add(subname//' N not conserved in skl_bgc')
714 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
715 0 : elseif (cinit(nn) < c0) then
716 0 : write(warnstr,*) subname, 'sk_bgc < 0 after algal fluxes, nn,cinit,flux_bio',&
717 0 : nn,cinit(nn),flux_bio(nn)
718 0 : call icepack_warnings_add(warnstr)
719 0 : write(warnstr,*) subname, 'cinit_tmp,flux_bio_temp,f_meltn,congel_alg,PVt,PVflag: '
720 0 : call icepack_warnings_add(warnstr)
721 0 : write(warnstr,*) subname, cinit_tmp,flux_bio_temp(nn),f_meltn(nn), &
722 0 : congel_alg(nn),PVt,PVflag(nn)
723 0 : call icepack_warnings_add(warnstr)
724 0 : write(warnstr,*) subname, 'congel, meltb: ',congel,meltb
725 0 : call icepack_warnings_add(warnstr)
726 0 : call icepack_warnings_add(subname//'sk_bgc < 0 after algal fluxes')
727 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
728 : endif
729 :
730 0 : if (icepack_warnings_aborted(subname)) return
731 :
732 : !-----------------------------------------------------------------
733 : ! reload tracer array: Bulk tracer concentration (mmol or mg per m^3)
734 : !-----------------------------------------------------------------
735 :
736 0 : trcrn(bio_index(nn)) = cinit(nn) * phi_sk/sk_l
737 :
738 : enddo !nbtrcr
739 :
740 : end subroutine skl_biogeochemistry
741 :
742 : !=======================================================================
743 : !
744 : ! Solve the scalar vertical diffusion equation implicitly using
745 : ! tridiag_solver. Calculate the diffusivity from temperature and salinity.
746 : !
747 : ! NOTE: In this subroutine, trcrn(nt_fbri) is the volume fraction of ice with
748 : ! dynamic salinity or the height ratio == hinS/vicen*aicen, where hinS is the
749 : ! height of the brine surface relative to the bottom of the ice. This volume fraction
750 : ! may be > 1 in which case there is brine above the ice surface (meltponds).
751 : !
752 :
753 0 : subroutine z_biogeochemistry (n_cat, dt, &
754 : nilyr, & ! LCOV_EXCL_LINE
755 : nblyr, nbtrcr, & ! LCOV_EXCL_LINE
756 : n_algae, n_doc, & ! LCOV_EXCL_LINE
757 : n_dic, n_don, & ! LCOV_EXCL_LINE
758 : n_fed, n_fep, & ! LCOV_EXCL_LINE
759 : n_zaero, first_ice, & ! LCOV_EXCL_LINE
760 : aicen, vicen, & ! LCOV_EXCL_LINE
761 : hice_old, ocean_bio, & ! LCOV_EXCL_LINE
762 : flux_bio, bphin, & ! LCOV_EXCL_LINE
763 : iphin, trcrn, & ! LCOV_EXCL_LINE
764 : iDin, & ! LCOV_EXCL_LINE
765 : fswthrul, grow_alg, & ! LCOV_EXCL_LINE
766 : upNOn, upNHn, & ! LCOV_EXCL_LINE
767 : dh_top, dh_bot, & ! LCOV_EXCL_LINE
768 : zfswin, hbri, & ! LCOV_EXCL_LINE
769 : hbri_old, darcy_V, & ! LCOV_EXCL_LINE
770 : ! darcy_V_chl, bgrid, & ! LCOV_EXCL_LINE
771 : bgrid, & ! LCOV_EXCL_LINE
772 : i_grid, ic_grid, & ! LCOV_EXCL_LINE
773 : bphi_min, zbgc_snow, & ! LCOV_EXCL_LINE
774 : zbgc_atm, & ! LCOV_EXCL_LINE
775 : iTin, dh_direct, & ! LCOV_EXCL_LINE
776 : Zoo, meltb, & ! LCOV_EXCL_LINE
777 : congel )
778 :
779 : integer (kind=int_kind), intent(in) :: &
780 : n_cat, & ! category number ! LCOV_EXCL_LINE
781 : nilyr, & ! number of ice layers ! LCOV_EXCL_LINE
782 : nblyr, & ! number of bio layers ! LCOV_EXCL_LINE
783 : nbtrcr, n_algae, & ! number of bgc tracers, number of autotrophs ! LCOV_EXCL_LINE
784 : n_zaero, & ! number of aerosols ! LCOV_EXCL_LINE
785 : n_doc, n_dic, n_don, n_fed, n_fep
786 :
787 : logical (kind=log_kind), intent(in) :: &
788 : first_ice ! initialized values should be used
789 :
790 : real (kind=dbl_kind), intent(in) :: &
791 : dt , & ! time step ! LCOV_EXCL_LINE
792 : hbri , & ! brine height (m) ! LCOV_EXCL_LINE
793 : bphi_min , & ! surface porosity ! LCOV_EXCL_LINE
794 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
795 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
796 : hice_old , & ! ice height (m) ! LCOV_EXCL_LINE
797 : meltb , & ! bottom melt in dt (m) ! LCOV_EXCL_LINE
798 : congel , & ! bottom growth in dt (m) ! LCOV_EXCL_LINE
799 : darcy_V , & ! darcy velocity ! LCOV_EXCL_LINE
800 : ! darcy_V_chl, & ! darcy velocity for algae ! LCOV_EXCL_LINE
801 : dh_bot , & ! change in brine bottom (m) ! LCOV_EXCL_LINE
802 : dh_top , & ! change in brine top (m) ! LCOV_EXCL_LINE
803 : dh_direct ! surface flooding or runoff (m)
804 :
805 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
806 : bgrid , & ! biology nondimensional vertical grid points ! LCOV_EXCL_LINE
807 : flux_bio , & ! total ocean tracer flux (mmol/m^2/s) ! LCOV_EXCL_LINE
808 : zfswin , & ! visible Short wave flux on igrid (W/m^2) ! LCOV_EXCL_LINE
809 : Zoo , & ! N losses to the system from reaction terms ! LCOV_EXCL_LINE
810 : ! (ie. zooplankton/bacteria) (mmol/m^3)
811 : trcrn ! bulk tracer concentration (mmol/m^3)
812 :
813 : real (kind=dbl_kind), dimension (:), intent(in) :: &
814 : i_grid , & ! biology vertical interface points ! LCOV_EXCL_LINE
815 : iTin , & ! salinity vertical interface points ! LCOV_EXCL_LINE
816 : iphin , & ! Porosity on the igrid ! LCOV_EXCL_LINE
817 : iDin , & ! Diffusivity/h on the igrid (1/s) ! LCOV_EXCL_LINE
818 : ic_grid , & ! CICE interface coordinate ! LCOV_EXCL_LINE
819 : fswthrul , & ! visible short wave radiation on icgrid (W/m^2) ! LCOV_EXCL_LINE
820 : zbgc_snow , & ! tracer input from snow (mmol/m^3*m) ! LCOV_EXCL_LINE
821 : zbgc_atm , & ! tracer input from atm (mmol/m^3 *m) ! LCOV_EXCL_LINE
822 : ocean_bio , & ! ocean concentrations (mmol/m^3) ! LCOV_EXCL_LINE
823 : bphin ! Porosity on the bgrid
824 :
825 : real (kind=dbl_kind), intent(inout) :: &
826 : hbri_old ! brine height (m)
827 :
828 : real (kind=dbl_kind), dimension (:,:), intent(out) :: &
829 : upNOn , & ! algal nitrate uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
830 : upNHn , & ! algal ammonium uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
831 : grow_alg ! algal growth rate (mmol/m^3/s)
832 :
833 : !-----------------------------------------------------------------------------
834 : ! algae absorption coefficient for 0.5 m thick layer
835 : ! Grenfell (1991): SA = specific absorption coefficient= 0.004 m^2/mg Chla
836 : ! generalizing kalg_bio(k) = SA*\sum R_chl2N(m)*trcrn(i,j,nt_bgc_N(m)+k-1)
837 : ! output kalg on the icgrid
838 : !-----------------------------------------------------------------------------
839 :
840 : ! local variables
841 :
842 : integer (kind=int_kind) :: &
843 : k, m, mm ! vertical biology layer index
844 :
845 : real (kind=dbl_kind) :: &
846 : hin , & ! ice thickness (m) ! LCOV_EXCL_LINE
847 : hin_old , & ! ice thickness before current melt/growth (m) ! LCOV_EXCL_LINE
848 : ice_conc , & ! algal concentration in ice above hin > hinS ! LCOV_EXCL_LINE
849 : sum_old , & ! ! LCOV_EXCL_LINE
850 : sum_new , & ! ! LCOV_EXCL_LINE
851 : sum_tot , & ! ! LCOV_EXCL_LINE
852 : zspace , & ! 1/nblyr ! LCOV_EXCL_LINE
853 : darcyV , & ! ! LCOV_EXCL_LINE
854 : dhtop , & ! ! LCOV_EXCL_LINE
855 : dhbot , & ! ! LCOV_EXCL_LINE
856 0 : dhflood ! >=0 (m) surface flooding from the ocean
857 :
858 : real (kind=dbl_kind), dimension (nblyr+2) :: &
859 0 : bphin_N ! porosity for tracer model has minimum
860 : ! bphin_N >= bphimin
861 :
862 : real (kind=dbl_kind), dimension (nblyr+1) :: &
863 : iphin_N , & ! tracer porosity on the igrid ! LCOV_EXCL_LINE
864 : sbdiagz , & ! sub-diagonal matrix elements ! LCOV_EXCL_LINE
865 : diagz , & ! diagonal matrix elements ! LCOV_EXCL_LINE
866 : spdiagz , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
867 : rhsz , & ! rhs of tri-diagonal matrix equation ! LCOV_EXCL_LINE
868 : ML_diag , & ! lumped mass matrix ! LCOV_EXCL_LINE
869 : D_spdiag , & ! artificial diffusion matrix ! LCOV_EXCL_LINE
870 : D_sbdiag , & ! artificial diffusion matrix ! LCOV_EXCL_LINE
871 : biomat_low , & ! Low order solution ! LCOV_EXCL_LINE
872 0 : Nerror ! Change in N after reactions
873 :
874 : real (kind=dbl_kind), dimension(nblyr+1,nbtrcr):: &
875 0 : react ! biological sources and sinks for equation matrix
876 :
877 : real (kind=dbl_kind), dimension(nblyr+1,nbtrcr):: &
878 : in_init_cons , & ! Initial bulk concentration*h (mmol/m^2) ! LCOV_EXCL_LINE
879 : biomat_cons , & ! Matrix output of (mmol/m^2) ! LCOV_EXCL_LINE
880 0 : biomat_brine ! brine concentration (mmol/m^3)
881 :
882 : real (kind=dbl_kind), dimension(nbtrcr):: &
883 : C_top, & ! bulk top tracer source: h phi C(meltwater) (mmol/m^2) ! LCOV_EXCL_LINE
884 : C_bot, & ! bulk bottom tracer source: h phi C(ocean) (mmol/m^2) ! LCOV_EXCL_LINE
885 : Source_top, & ! For cons: (+) top tracer source into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
886 : Source_bot, & ! For cons: (+) bottom tracer source into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
887 : Sink_bot, & ! For cons: (+ or -) remaining bottom flux into ice(mmol/m^2/s) ! LCOV_EXCL_LINE
888 : Sink_top, & ! For cons: (+ or -) remaining bottom flux into ice(mmol/m^2/s) ! LCOV_EXCL_LINE
889 : exp_ret, & ! exp dt/retention frequency ! LCOV_EXCL_LINE
890 : exp_rel, & ! exp dt/release frequency ! LCOV_EXCL_LINE
891 : atm_add_cons , & ! zbgc_snow+zbgc_atm (mmol/m^3*m) ! LCOV_EXCL_LINE
892 : dust_Fe , & ! contribution of dust surface flux to dFe (umol/m*3*m) ! LCOV_EXCL_LINE
893 0 : source ! mmol/m^2 surface input from snow/atmosphere
894 :
895 : real (kind=dbl_kind), dimension (ntrcr+2) :: &
896 : trtmp0 , & ! temporary, remapped tracers ! LCOV_EXCL_LINE
897 0 : trtmp ! temporary, remapped tracers
898 :
899 : logical (kind=log_kind), dimension(nblyr+1) :: &
900 0 : conserve_N
901 :
902 : real (kind=dbl_kind), dimension(nblyr+1):: & ! temporary variables for
903 : Diff , & ! diffusivity ! LCOV_EXCL_LINE
904 : initcons , & ! initial concentration ! LCOV_EXCL_LINE
905 : biocons , & ! new concentration ! LCOV_EXCL_LINE
906 : dmobile , & ! ! LCOV_EXCL_LINE
907 : initcons_mobile,&! ! LCOV_EXCL_LINE
908 0 : initcons_stationary
909 :
910 : real (kind=dbl_kind):: &
911 0 : top_conc ! 1% (min_bgc) of surface concentration
912 : ! when hin > hbri: just used in sw calculation
913 :
914 : real (kind=dbl_kind):: &
915 : bio_tmp, & ! temporary variable ! LCOV_EXCL_LINE
916 0 : exp_min ! temporary exp var
917 :
918 : real (kind=dbl_kind):: &
919 : Sat_conc , & ! adsorbing saturation concentration (mmols/m^3) ! LCOV_EXCL_LINE
920 : phi_max , & ! maximum porosity ! LCOV_EXCL_LINE
921 : S_col , & ! surface area of collector (um^2) ! LCOV_EXCL_LINE
922 : P_b , & ! projected area of diatoms (um^2) ! LCOV_EXCL_LINE
923 : V_c , & ! volume of collector (um^3) ! LCOV_EXCL_LINE
924 0 : V_alg ! volume of algae (um^3)
925 :
926 : real (kind=dbl_kind), dimension(nbtrcr) :: &
927 0 : mobile ! c1 if mobile, c0 otherwise
928 :
929 : ! local parameters
930 :
931 : real (kind=dbl_kind), parameter :: &
932 : accuracy = 1.0e-14_dbl_kind, & ! LCOV_EXCL_LINE
933 : r_c = 3.0e3_dbl_kind , & ! ice crystal radius (um) ! LCOV_EXCL_LINE
934 : r_bac= 15.0_dbl_kind , & ! diatom large radius (um) ! LCOV_EXCL_LINE
935 : r_alg= 10.0_dbl_kind , & ! diatom small radius (um) ! LCOV_EXCL_LINE
936 : N_vol = 0.04e-12_dbl_kind , & ! (g) Nitrogen per um^3 ! LCOV_EXCL_LINE
937 : Ng_to_mmol =0.0140067_dbl_kind , & ! (g/mmol) Nitrogen ! LCOV_EXCL_LINE
938 : f_s = c1 , & ! fracton of sites available for saturation ! LCOV_EXCL_LINE
939 : f_a = c1 , & ! fraction of collector available for attachment ! LCOV_EXCL_LINE
940 : f_v = 0.7854 ! fraction of algal coverage on area availabel for attachment
941 : ! 4(pi r^2)/(4r)^2 [Johnson et al, 1995, water res. research]
942 :
943 : integer, parameter :: &
944 : nt_zfswin = 1 ! for interpolation of short wave to bgrid
945 :
946 : character(len=*),parameter :: subname='(z_biogeochemistry)'
947 :
948 : !-------------------------------------
949 : ! Initialize
950 : !-----------------------------------
951 :
952 0 : zspace = c1/real(nblyr,kind=dbl_kind)
953 0 : in_init_cons(:,:) = c0
954 0 : atm_add_cons(:) = c0
955 0 : dhtop = c0
956 0 : dhbot = c0
957 0 : darcyV = c0
958 0 : C_top(:) = c0
959 0 : mobile(:) = c0
960 0 : conserve_N(:) = .true.
961 :
962 0 : do m = 1, nbtrcr
963 0 : do k = 1, nblyr+1
964 :
965 0 : bphin_N(nblyr+2) =c1
966 0 : bphin_N(k) = bphin(k)
967 0 : iphin_N(k) = iphin(k)
968 0 : bphin_N(1) = bphi_min
969 :
970 0 : if (first_ice) then
971 0 : trcrn(bio_index(m) + k-1) = ocean_bio(m)*zbgc_init_frac(m)
972 0 : in_init_cons(k,m) = trcrn(bio_index(m) + k-1)*hbri_old
973 0 : elseif (abs(trcrn(bio_index(m) + k-1)) < puny) then
974 0 : trcrn(bio_index(m) + k-1) = c0
975 0 : in_init_cons(k,m) = c0
976 : else
977 0 : in_init_cons(k,m) = trcrn(bio_index(m) + k-1)* hbri_old
978 : endif ! first_ice
979 :
980 0 : if (trcrn(bio_index(m) + k-1) < c0 ) then
981 0 : write(warnstr,*) subname,'zbgc initialization error, first ice = ', first_ice
982 0 : call icepack_warnings_add(warnstr)
983 0 : write(warnstr,*) subname,'Category,m:',n_cat,m
984 0 : call icepack_warnings_add(warnstr)
985 0 : write(warnstr,*) subname,'hbri,hbri_old'
986 0 : call icepack_warnings_add(warnstr)
987 0 : write(warnstr,*) subname, hbri,hbri_old
988 0 : call icepack_warnings_add(warnstr)
989 0 : write(warnstr,*) subname,'trcrn(bio_index(m) + k-1)'
990 0 : call icepack_warnings_add(warnstr)
991 0 : write(warnstr,*) subname, trcrn(bio_index(m) + k-1)
992 0 : call icepack_warnings_add(warnstr)
993 0 : call icepack_warnings_add(subname//' zbgc initialization error')
994 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
995 : endif
996 0 : if (icepack_warnings_aborted(subname)) return
997 : enddo !k
998 : enddo !m
999 :
1000 : !-----------------------------------------------------------------
1001 : ! boundary conditions
1002 : !-----------------------------------------------------------------
1003 :
1004 0 : ice_conc = c0
1005 0 : hin = vicen/aicen
1006 0 : hin_old = hice_old
1007 :
1008 : !-----------------------------------------------------------------
1009 : ! calculate the saturation concentration for attachment: Sat_conc
1010 : !-----------------------------------------------------------------
1011 :
1012 0 : phi_max = maxval(bphin_N(2:nblyr+1))
1013 0 : S_col = 4.0_dbl_kind*pi*r_c**2
1014 0 : P_b = pi*r_bac**2 !*10-6 for colloids
1015 0 : V_c = 4.0_dbl_kind*pi*r_c**3/3.0_dbl_kind*(1.0e-6_dbl_kind)**3 ! (m^3) sphere
1016 0 : V_alg = pi/6.0_dbl_kind*r_bac*r_alg**2 ! prolate spheroid (*10-9 for colloids)
1017 0 : Sat_conc= f_s*f_a*f_v*(c1-phi_max)/V_c*S_col/P_b*N_vol*V_alg/Ng_to_mmol
1018 : !mmol/m^3 (algae, don, hum...) and umols/m^3 for colloids
1019 :
1020 : !-----------------------------------------------------------------
1021 : ! convert surface dust flux (n_zaero > 2) to dFe(1) flux
1022 : !-----------------------------------------------------------------
1023 :
1024 0 : dust_Fe(:) = c0
1025 :
1026 0 : if (tr_zaero .and. n_zaero > 2 .and. tr_bgc_Fe) then
1027 0 : do m = 3,n_zaero
1028 0 : dust_Fe(nlt_bgc_Fed(1)) = dust_Fe(nlt_bgc_Fed(1)) + &
1029 : (zbgc_snow(nlt_zaero(m)) + zbgc_atm(nlt_zaero(m))) * & ! LCOV_EXCL_LINE
1030 0 : R_dFe2dust * dustFe_sol
1031 : ! dust_Fe(nlt_zaero(m)) = -(zbgc_snow(nlt_zaero(m)) + zbgc_atm(nlt_zaero(m))) * &
1032 : ! dustFe_sol
1033 : enddo
1034 : endif
1035 :
1036 0 : do m = 1,nbtrcr
1037 : !-----------------------------------------------------------------
1038 : ! time constants for mobile/stationary phase changes
1039 : !-----------------------------------------------------------------
1040 :
1041 0 : exp_rel(m) = c0
1042 0 : exp_ret(m) = c0
1043 0 : if (tau_ret(m) > c0) then
1044 0 : exp_min = min(dt/tau_ret(m),exp_argmax)
1045 0 : exp_ret(m) = exp(-exp_min)
1046 : endif
1047 0 : if (tau_rel(m) > c0) then
1048 0 : exp_min = min(dt/tau_rel(m),exp_argmax)
1049 0 : exp_rel(m) = exp(-exp_min)
1050 : endif
1051 0 : if (m .ne. nlt_bgc_N(1)) then
1052 0 : if (hin_old > hin) then !melting
1053 0 : exp_ret(m) = c1
1054 : else !not melting
1055 0 : exp_rel(m) = c1
1056 : endif
1057 0 : elseif (tr_bgc_N .and. hin_old > hin + algal_vel*dt) then
1058 0 : exp_ret(m) = c1
1059 0 : elseif (tr_bgc_N) then
1060 0 : exp_rel(m) = c1
1061 : endif
1062 :
1063 0 : dhtop = dh_top
1064 0 : dhbot = dh_bot
1065 0 : darcyV = darcy_V
1066 0 : C_top(m) = in_init_cons(1,m)*trcrn(nt_zbgc_frac+m-1)!mobile fraction
1067 0 : source(m) = abs(zbgc_snow(m) + zbgc_atm(m) + dust_Fe(m))
1068 0 : dhflood = max(c0,-dh_direct) ! ocean water flooding surface
1069 :
1070 0 : if (dhtop+darcyV/bphin_N(1)*dt < -puny) then !snow/top ice melt
1071 0 : C_top(m) = (zbgc_snow(m)+zbgc_atm(m) + dust_Fe(m))/abs(dhtop &
1072 0 : + darcyV/bphin_N(1)*dt + puny)*hbri_old
1073 0 : elseif (dhtop+darcyV/bphin_N(1)*dt >= -puny .and. &
1074 : abs((zbgc_snow(m)+zbgc_atm(m) + dust_Fe(m)) + & ! LCOV_EXCL_LINE
1075 0 : ocean_bio(m)*bphin_N(1)*dhflood) > puny) then
1076 0 : atm_add_cons(m) = abs(zbgc_snow(m) + zbgc_atm(m)+ dust_Fe(m)) + &
1077 0 : ocean_bio(m)*bphin_N(1)*dhflood
1078 : else ! only positive fluxes
1079 0 : atm_add_cons(m) = abs(zbgc_snow(m) + zbgc_atm(m)+ dust_Fe(m))
1080 : endif
1081 :
1082 0 : C_bot(m) = ocean_bio(m)*hbri_old*iphin_N(nblyr+1)
1083 :
1084 : enddo ! m
1085 :
1086 : !-----------------------------------------------------------------
1087 : ! Interpolate shortwave flux, fswthrul (defined at top to bottom with nilyr+1
1088 : ! evenly spaced with spacing = (1/nilyr) to grid variable zfswin:
1089 : !-----------------------------------------------------------------
1090 :
1091 0 : trtmp(:) = c0
1092 0 : trtmp0(:)= c0
1093 0 : zfswin(:) = c0
1094 :
1095 0 : do k = 1, nilyr+1
1096 : ! contains cice values (fswthrul(1) is surface value)
1097 : ! and fwsthrul(nilyr+1) is output
1098 0 : trtmp0(nt_zfswin+k-1) = fswthrul(k)
1099 : enddo !k
1100 :
1101 : call remap_zbgc(nilyr+1, &
1102 : nt_zfswin, & ! LCOV_EXCL_LINE
1103 : trtmp0(1:ntrcr), trtmp(1:ntrcr+2), & ! LCOV_EXCL_LINE
1104 : 0, nblyr+1, & ! LCOV_EXCL_LINE
1105 : hin, hbri, & ! LCOV_EXCL_LINE
1106 : ic_grid(1:nilyr+1), & ! LCOV_EXCL_LINE
1107 0 : i_grid(1:nblyr+1),ice_conc )
1108 :
1109 0 : if (icepack_warnings_aborted(subname)) return
1110 :
1111 0 : do k = 1,nblyr+1
1112 0 : zfswin(k) = trtmp(nt_zfswin+k-1)
1113 : enddo
1114 :
1115 : !-----------------------------------------------------------------
1116 : ! Initialze Biology
1117 : !-----------------------------------------------------------------
1118 :
1119 0 : do mm = 1, nbtrcr
1120 0 : mobile(mm) = c0
1121 0 : if (bgc_tracer_type(mm) .GE. c0) mobile(mm) = c1
1122 :
1123 0 : do k = 1, nblyr+1
1124 0 : biomat_cons(k,mm) = in_init_cons(k,mm)
1125 : enddo !k
1126 : enddo !mm
1127 :
1128 : !-----------------------------------------------------------------
1129 : ! Compute FCT
1130 : !-----------------------------------------------------------------
1131 :
1132 0 : do mm = 1, nbtrcr
1133 :
1134 0 : if (hbri_old > thinS .and. hbri > thinS) then
1135 0 : do k = 1,nblyr+1
1136 0 : initcons_mobile(k) = in_init_cons(k,mm)*trcrn(nt_zbgc_frac+mm-1)
1137 0 : initcons_stationary(k) = mobile(mm)*(in_init_cons(k,mm)-initcons_mobile(k))
1138 0 : dmobile(k) = mobile(mm)*(initcons_mobile(k)*(exp_ret(mm)-c1) + &
1139 0 : initcons_stationary(k)*(c1-exp_rel(mm)))
1140 0 : initcons_mobile(k) = max(c0,initcons_mobile(k) + dmobile(k))
1141 0 : initcons_stationary(k) = max(c0,initcons_stationary(k) - dmobile(k))
1142 0 : if (initcons_stationary(k)/hbri_old > Sat_conc) then
1143 0 : initcons_mobile(k) = initcons_mobile(k) + initcons_stationary(k) - Sat_conc*hbri_old
1144 0 : initcons_stationary(k) = Sat_conc*hbri_old
1145 : endif
1146 :
1147 0 : Diff(k) = iDin(k)
1148 0 : initcons(k) = initcons_mobile(k)
1149 0 : biocons(k) = initcons_mobile(k)
1150 : enddo
1151 :
1152 : call compute_FCT_matrix &
1153 : (initcons,sbdiagz, dt, nblyr, & ! LCOV_EXCL_LINE
1154 : diagz, spdiagz, rhsz, bgrid, & ! LCOV_EXCL_LINE
1155 : darcyV, dhtop, & ! LCOV_EXCL_LINE
1156 : dhbot, iphin_N, & ! LCOV_EXCL_LINE
1157 : Diff, hbri_old, & ! LCOV_EXCL_LINE
1158 : atm_add_cons(mm), bphin_N, & ! LCOV_EXCL_LINE
1159 : C_top(mm), C_bot(mm), & ! LCOV_EXCL_LINE
1160 : Source_bot(mm), Source_top(mm),& ! LCOV_EXCL_LINE
1161 : Sink_bot(mm),Sink_top(mm), & ! LCOV_EXCL_LINE
1162 0 : D_sbdiag, D_spdiag, ML_diag)
1163 0 : if (icepack_warnings_aborted(subname)) return
1164 :
1165 : call tridiag_solverz &
1166 : (nblyr+1, sbdiagz, & ! LCOV_EXCL_LINE
1167 : diagz, spdiagz, & ! LCOV_EXCL_LINE
1168 0 : rhsz, biocons)
1169 0 : if (icepack_warnings_aborted(subname)) return
1170 :
1171 : call check_conservation_FCT &
1172 : (initcons, & ! LCOV_EXCL_LINE
1173 : biocons, & ! LCOV_EXCL_LINE
1174 : biomat_low, & ! LCOV_EXCL_LINE
1175 : Source_top(mm), & ! LCOV_EXCL_LINE
1176 : Source_bot(mm), & ! LCOV_EXCL_LINE
1177 : Sink_bot(mm), & ! LCOV_EXCL_LINE
1178 : Sink_top(mm), & ! LCOV_EXCL_LINE
1179 : dt, flux_bio(mm), & ! LCOV_EXCL_LINE
1180 : nblyr, & ! LCOV_EXCL_LINE
1181 0 : source(mm))
1182 0 : if (icepack_warnings_aborted(subname)) return
1183 :
1184 : call compute_FCT_corr &
1185 : (initcons, & ! LCOV_EXCL_LINE
1186 : biocons, dt, nblyr, & ! LCOV_EXCL_LINE
1187 0 : D_sbdiag, D_spdiag, ML_diag)
1188 0 : if (icepack_warnings_aborted(subname)) return
1189 :
1190 0 : top_conc = c0 ! or frazil ice concentration
1191 :
1192 : ! assume diatoms actively maintain there relative position in the ice
1193 :
1194 0 : if (mm .ne. nlt_bgc_N(1)) then
1195 :
1196 : call regrid_stationary &
1197 : (initcons_stationary, hbri_old, & ! LCOV_EXCL_LINE
1198 : hbri, dt, & ! LCOV_EXCL_LINE
1199 : ntrcr, & ! LCOV_EXCL_LINE
1200 : nblyr, top_conc, & ! LCOV_EXCL_LINE
1201 : i_grid, flux_bio(mm),& ! LCOV_EXCL_LINE
1202 0 : meltb, congel)
1203 0 : if (icepack_warnings_aborted(subname)) return
1204 :
1205 0 : elseif (tr_bgc_N .and. mm .eq. nlt_bgc_N(1)) then
1206 0 : if (meltb > algal_vel*dt .or. aicen < 0.001_dbl_kind) then
1207 :
1208 : call regrid_stationary &
1209 : (initcons_stationary, hbri_old, & ! LCOV_EXCL_LINE
1210 : hbri, dt, & ! LCOV_EXCL_LINE
1211 : ntrcr, & ! LCOV_EXCL_LINE
1212 : nblyr, top_conc, & ! LCOV_EXCL_LINE
1213 : i_grid, flux_bio(mm),& ! LCOV_EXCL_LINE
1214 0 : meltb, congel)
1215 0 : if (icepack_warnings_aborted(subname)) return
1216 :
1217 : endif
1218 : endif
1219 :
1220 0 : biomat_cons(:,mm) = biocons(:) + initcons_stationary(:)
1221 :
1222 0 : sum_old = (biomat_low(1) + biomat_low(nblyr+1))*zspace/c2
1223 0 : sum_new = (biocons(1)+ biocons(nblyr+1))*zspace/c2
1224 0 : sum_tot = (biomat_cons(1,mm) + biomat_cons(nblyr+1,mm))*zspace/c2
1225 0 : do k = 2,nblyr
1226 0 : sum_old = sum_old + biomat_low(k)*zspace
1227 0 : sum_new = sum_new + biocons(k)*zspace
1228 0 : sum_tot = sum_tot + biomat_cons(k,mm)*zspace
1229 : enddo
1230 0 : trcrn(nt_zbgc_frac+mm-1) = zbgc_frac_init(mm)
1231 0 : if (sum_tot > c0 .and. mobile(mm) > c0) trcrn(nt_zbgc_frac+mm-1) = sum_new/sum_tot
1232 :
1233 : if (abs(sum_new-sum_old) > accuracy*sum_old .or. &
1234 : minval(biocons(:)) < c0 .or. minval(initcons_stationary(:)) < c0 & ! LCOV_EXCL_LINE
1235 0 : .or. icepack_warnings_aborted()) then
1236 0 : write(warnstr,*) subname,'zbgc FCT tracer solution failed'
1237 0 : call icepack_warnings_add(warnstr)
1238 0 : write(warnstr,*) subname,'sum_new,sum_old:',sum_new,sum_old
1239 0 : call icepack_warnings_add(warnstr)
1240 0 : write(warnstr,*) subname,'mm,biocons(:):',mm,biocons(:)
1241 0 : call icepack_warnings_add(warnstr)
1242 0 : write(warnstr,*) subname,'biomat_low:',biomat_low
1243 0 : call icepack_warnings_add(warnstr)
1244 0 : write(warnstr,*) subname,'Diff(:):',Diff(:)
1245 0 : call icepack_warnings_add(warnstr)
1246 0 : write(warnstr,*) subname,'dmobile(:):',dmobile(:)
1247 0 : call icepack_warnings_add(warnstr)
1248 0 : write(warnstr,*) subname,'mobile(mm):',mobile(mm)
1249 0 : call icepack_warnings_add(warnstr)
1250 0 : write(warnstr,*) subname,'initcons_stationary(:):',initcons_stationary(:)
1251 0 : call icepack_warnings_add(warnstr)
1252 0 : write(warnstr,*) subname, 'trcrn(nt_zbgc_frac+mm-1):',trcrn(nt_zbgc_frac+mm-1)
1253 0 : call icepack_warnings_add(warnstr)
1254 0 : write(warnstr,*) subname, 'in_init_cons(:,mm):',in_init_cons(:,mm)
1255 0 : call icepack_warnings_add(warnstr)
1256 0 : write(warnstr,*) subname, 'exp_ret( mm),exp_rel( mm)',exp_ret( mm),exp_rel( mm)
1257 0 : call icepack_warnings_add(warnstr)
1258 0 : write(warnstr,*) subname,'darcyV,dhtop,dhbot'
1259 0 : call icepack_warnings_add(warnstr)
1260 0 : write(warnstr,*) subname,darcyV,dhtop,dhbot
1261 0 : call icepack_warnings_add(warnstr)
1262 0 : write(warnstr,*) subname,'Category,mm:',n_cat,mm
1263 0 : call icepack_warnings_add(warnstr)
1264 0 : call icepack_warnings_add(subname//'zbgc FCT tracer solution failed')
1265 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1266 : endif
1267 0 : if (icepack_warnings_aborted(subname)) return
1268 :
1269 : else
1270 :
1271 0 : call thin_ice_flux(hbri,hbri_old,biomat_cons(:,mm), &
1272 : flux_bio(mm),source(mm), & ! LCOV_EXCL_LINE
1273 0 : dt, nblyr,ocean_bio(mm))
1274 0 : if (icepack_warnings_aborted(subname)) return
1275 :
1276 : endif ! thin or not
1277 :
1278 0 : do k = 1,nblyr+1
1279 0 : biomat_brine(k,mm) = biomat_cons(k,mm)/hbri/iphin_N(k)
1280 : enddo ! k
1281 : enddo ! mm
1282 :
1283 0 : react(:,:) = c0
1284 0 : grow_alg(:,:) = c0
1285 :
1286 0 : if (solve_zbgc) then
1287 0 : do k = 1, nblyr+1
1288 : call algal_dyn (dt, &
1289 : n_doc, n_dic, n_don, n_fed, n_fep, & ! LCOV_EXCL_LINE
1290 : dEdd_algae, & ! LCOV_EXCL_LINE
1291 : zfswin(k), react(k,:), & ! LCOV_EXCL_LINE
1292 : biomat_brine(k,:), & ! LCOV_EXCL_LINE
1293 : grow_alg(k,:), n_algae, & ! LCOV_EXCL_LINE
1294 : iTin(k), & ! LCOV_EXCL_LINE
1295 : upNOn(k,:), upNHn(k,:), & ! LCOV_EXCL_LINE
1296 : Zoo(k), & ! LCOV_EXCL_LINE
1297 0 : Nerror(k), conserve_N(k))
1298 0 : if (icepack_warnings_aborted(subname)) return
1299 :
1300 : enddo ! k
1301 : endif ! solve_zbgc
1302 :
1303 : !-----------------------------------------------------------------
1304 : ! Update the tracer variable
1305 : !-----------------------------------------------------------------
1306 :
1307 0 : do m = 1,nbtrcr
1308 0 : do k = 1,nblyr+1 ! back to bulk quantity
1309 0 : bio_tmp = (biomat_brine(k,m) + react(k,m))*iphin_N(k)
1310 :
1311 0 : if (.not. conserve_N(k)) then
1312 0 : write(warnstr,*) subname, 'N in algal_dyn not conserved'
1313 0 : call icepack_warnings_add(warnstr)
1314 0 : write(warnstr,*) subname, 'Nerror(k):', Nerror(k)
1315 0 : call icepack_warnings_add(warnstr)
1316 0 : write(warnstr,*) subname, 'k,m,hbri,hbri_old,bio_tmp,biomat_cons(k,m),ocean_bio(m)'
1317 0 : call icepack_warnings_add(warnstr)
1318 0 : write(warnstr,*) subname, k,m,hbri,hbri_old,bio_tmp,biomat_cons(k,m),ocean_bio(m)
1319 0 : call icepack_warnings_add(warnstr)
1320 0 : write(warnstr,*) subname, 'react(k,m),iphin_N(k),biomat_brine(k,m)'
1321 0 : call icepack_warnings_add(warnstr)
1322 0 : write(warnstr,*) subname, react(k,m),iphin_N(k),biomat_brine(k,m)
1323 0 : call icepack_warnings_add(warnstr)
1324 0 : call icepack_warnings_add(subname//' N in algal_dyn not conserved')
1325 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1326 0 : elseif (abs(bio_tmp) < puny) then
1327 0 : bio_tmp = c0
1328 0 : elseif (bio_tmp > 1.0e6_dbl_kind) then
1329 0 : write(warnstr,*) subname, 'very large bgc value'
1330 0 : call icepack_warnings_add(warnstr)
1331 0 : write(warnstr,*) subname, 'k,m,hbri,hbri_old,bio_tmp,biomat_cons(k,m),ocean_bio(m)'
1332 0 : call icepack_warnings_add(warnstr)
1333 0 : write(warnstr,*) subname, k,m,hbri,hbri_old,bio_tmp,biomat_cons(k,m),ocean_bio(m)
1334 0 : call icepack_warnings_add(warnstr)
1335 0 : write(warnstr,*) subname, 'react(k,m),iphin_N(k),biomat_brine(k,m)'
1336 0 : call icepack_warnings_add(warnstr)
1337 0 : write(warnstr,*) subname, react(k,m),iphin_N(k),biomat_brine(k,m)
1338 0 : call icepack_warnings_add(warnstr)
1339 0 : call icepack_warnings_add(subname//' very large bgc value')
1340 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1341 0 : elseif (bio_tmp < c0) then
1342 0 : write(warnstr,*) subname, 'negative bgc'
1343 0 : call icepack_warnings_add(warnstr)
1344 0 : write(warnstr,*) subname, 'k,m,nlt_bgc_Nit,hbri,hbri_old'
1345 0 : call icepack_warnings_add(warnstr)
1346 0 : write(warnstr,*) subname, k,m,nlt_bgc_Nit,hbri,hbri_old
1347 0 : call icepack_warnings_add(warnstr)
1348 0 : write(warnstr,*) subname, 'bio_tmp,biomat_cons(k,m),ocean_bio(m)'
1349 0 : call icepack_warnings_add(warnstr)
1350 0 : write(warnstr,*) subname, bio_tmp,biomat_cons(k,m),ocean_bio(m)
1351 0 : call icepack_warnings_add(warnstr)
1352 0 : write(warnstr,*) subname, 'react(k,m),iphin_N(k),biomat_brine(k,m)'
1353 0 : call icepack_warnings_add(warnstr)
1354 0 : write(warnstr,*) subname, react(k,m),iphin_N(k),biomat_brine(k,m)
1355 0 : call icepack_warnings_add(warnstr)
1356 0 : write(warnstr,*) subname, 'exp_ret( m),exp_ret( m)',exp_ret( m),exp_ret( m)
1357 0 : call icepack_warnings_add(warnstr)
1358 0 : call icepack_warnings_add(subname//'negative bgc')
1359 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1360 : endif
1361 0 : if (icepack_warnings_aborted()) then
1362 0 : write(warnstr,*) subname, 'trcrn(nt_zbgc_frac+m-1):',trcrn(nt_zbgc_frac+m-1)
1363 0 : call icepack_warnings_add(warnstr)
1364 0 : write(warnstr,*) subname, 'in_init_cons(k,m):',in_init_cons(k,m)
1365 0 : call icepack_warnings_add(warnstr)
1366 0 : write(warnstr,*) subname, 'trcrn(bio_index(m) + k-1)'
1367 0 : call icepack_warnings_add(warnstr)
1368 0 : write(warnstr,*) subname, trcrn(bio_index(m) + k-1)
1369 0 : call icepack_warnings_add(warnstr)
1370 0 : write(warnstr,*) subname, 'Category,m:',n_cat,m
1371 0 : call icepack_warnings_add(warnstr)
1372 0 : return
1373 : endif
1374 0 : trcrn(bio_index(m)+k-1) = max(c0, bio_tmp)
1375 0 : if (ocean_bio(m) .le. c0 .and. flux_bio(m) < c0) then
1376 : ! if (flux_bio(m) < -1.0e-12_dbl_kind) then
1377 : ! write(warnstr,*) subname, 'no ocean_bio but flux_bio < c0'
1378 : ! call icepack_warnings_add(warnstr)
1379 : ! write(warnstr,*) subname, 'm,ocean_bio(m),flux_bio(m)'
1380 : ! call icepack_warnings_add(warnstr)
1381 : ! write(warnstr,*) subname, m,ocean_bio(m),flux_bio(m)
1382 : ! call icepack_warnings_add(warnstr)
1383 : ! write(warnstr,*) subname, 'setting flux_bio(m) = c0'
1384 : ! call icepack_warnings_add(warnstr)
1385 : ! call icepack_warnings_add(subname//' flux_bio < 0 when ocean_bio = 0')
1386 : ! call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1387 : ! endif
1388 0 : flux_bio(m) = max(c0,flux_bio(m))
1389 : endif
1390 : enddo ! k
1391 : enddo ! m
1392 :
1393 : end subroutine z_biogeochemistry
1394 :
1395 : !=======================================================================
1396 : !
1397 : ! Do biogeochemistry from subroutine algal_dynamics
1398 : ! authors: Scott Elliott, LANL
1399 : ! Nicole Jeffery, LANL
1400 :
1401 0 : subroutine algal_dyn (dt, &
1402 : n_doc, n_dic, n_don, n_fed, n_fep, & ! LCOV_EXCL_LINE
1403 : dEdd_algae, & ! LCOV_EXCL_LINE
1404 : fswthru, reactb, & ! LCOV_EXCL_LINE
1405 : ltrcrn, & ! LCOV_EXCL_LINE
1406 : grow_alg, n_algae, & ! LCOV_EXCL_LINE
1407 : T_bot, & ! LCOV_EXCL_LINE
1408 : upNOn, upNHn, & ! LCOV_EXCL_LINE
1409 : Zoo, & ! LCOV_EXCL_LINE
1410 : Nerror, conserve_N)
1411 :
1412 : integer (kind=int_kind), intent(in) :: &
1413 : n_doc, n_dic, n_don, n_fed, n_fep, & ! LCOV_EXCL_LINE
1414 : n_algae ! number of autotrophic types
1415 :
1416 : real (kind=dbl_kind), intent(in) :: &
1417 : dt , & ! time step ! LCOV_EXCL_LINE
1418 : T_bot , & ! ice temperature (oC) ! LCOV_EXCL_LINE
1419 : fswthru ! average shortwave passing through current ice layer (W/m^2)
1420 :
1421 : real (kind=dbl_kind), intent(inout) :: &
1422 : Zoo, & ! N losses from zooplankton/bacteria... (mmol/m^3) ! LCOV_EXCL_LINE
1423 : Nerror ! Change in N after reactions (mmol/m^3)
1424 :
1425 : real (kind=dbl_kind), dimension (:), intent(out) :: &
1426 : grow_alg,& ! algal growth rate (mmol/m^3/s) ! LCOV_EXCL_LINE
1427 : upNOn, & ! algal NO uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
1428 : upNHn ! algal NH uptake rate (mmol/m^3/s)
1429 :
1430 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1431 : reactb ! biological reaction terms (mmol/m3)
1432 :
1433 : real (kind=dbl_kind), dimension(:), intent(in) :: &
1434 : ltrcrn ! brine concentrations in layer (mmol/m^3)
1435 :
1436 : logical (kind=log_kind), intent(inout) :: &
1437 : conserve_N
1438 :
1439 : logical (kind=log_kind), intent(in) :: &
1440 : dEdd_algae ! .true. chla impact on shortwave computed in dEdd
1441 :
1442 : ! local variables
1443 : !------------------------------------------------------------------------------------
1444 : ! 3 possible autotrophs nt_bgc_N(1:3): diatoms, flagellates, phaeocystis
1445 : ! 2 types of dissolved organic carbon nt_bgc_DOC(1:2):
1446 : ! polysaccharids, lipids
1447 : ! 1 DON (proteins)
1448 : ! 1 particulate iron (nt_bgc_Fe) n_fep
1449 : ! 1 dossp;ved orpm m+fed
1450 : ! Limiting macro/micro nutrients: nt_bgc_Nit -> nitrate, nt_bgc_NH -> ammonium,
1451 : ! nt_bgc_Sil -> silicate, nt_bgc_Fe -> dissolved iron
1452 : ! --------------------------------------------------------------------------------------
1453 :
1454 : ! real (kind=dbl_kind), parameter, dimension(max_algae) :: &
1455 : ! alpha2max_high = (/ 0.25_dbl_kind, 0.25_dbl_kind, 0.25_dbl_kind/) ! light limitation (1/(W/m^2))
1456 :
1457 : integer (kind=int_kind) :: k, n
1458 :
1459 : real (kind=dbl_kind), dimension(n_algae) :: &
1460 : Nin , & ! algal nitrogen concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1461 : ! Cin , & ! algal carbon concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1462 0 : chlin ! algal chlorophyll concentration on volume (mg/m^3)
1463 :
1464 : real (kind=dbl_kind), dimension(n_doc) :: &
1465 0 : DOCin ! dissolved organic carbon concentration on volume (mmolC/m^3)
1466 :
1467 : ! real (kind=dbl_kind), dimension(n_dic) :: &
1468 : ! DICin ! dissolved inorganic carbon concentration on volume (mmolC/m^3)
1469 :
1470 : real (kind=dbl_kind), dimension(n_don) :: & !proteins
1471 0 : DONin ! dissolved organic nitrogen concentration on volume (mmolN/m^3)
1472 :
1473 : real (kind=dbl_kind), dimension(n_fed) :: & !iron
1474 0 : Fedin ! dissolved iron concentration on volume (umol/m^3)
1475 :
1476 : real (kind=dbl_kind), dimension(n_fep) :: & !iron
1477 0 : Fepin ! algal nitrogen concentration on volume (umol/m^3)
1478 :
1479 : real (kind=dbl_kind) :: &
1480 : Nitin , & ! nitrate concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1481 : Amin , & ! ammonia/um concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1482 : Silin , & ! silicon concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1483 : ! DMSPpin , & ! DMSPp concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1484 : DMSPdin , & ! DMSPd concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1485 : DMSin , & ! DMS concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1486 : ! PONin , & ! PON concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1487 : op_dep , & ! bottom layer attenuation exponent (optical depth) ! LCOV_EXCL_LINE
1488 0 : Iavg_loc ! bottom layer attenuated Fswthru (W/m^2)
1489 :
1490 : real (kind=dbl_kind), dimension(n_algae) :: &
1491 : L_lim , & ! overall light limitation ! LCOV_EXCL_LINE
1492 : Nit_lim , & ! overall nitrate limitation ! LCOV_EXCL_LINE
1493 : Am_lim , & ! overall ammonium limitation ! LCOV_EXCL_LINE
1494 : N_lim , & ! overall nitrogen species limitation ! LCOV_EXCL_LINE
1495 : Sil_lim , & ! overall silicon limitation ! LCOV_EXCL_LINE
1496 : Fe_lim , & ! overall iron limitation ! LCOV_EXCL_LINE
1497 : fr_Nit , & ! fraction of local ecological growth as nitrate ! LCOV_EXCL_LINE
1498 : fr_Am , & ! fraction of local ecological growth as ammonia ! LCOV_EXCL_LINE
1499 : growmax_N, & ! maximum growth rate in N currency (mmol/m^3/s) ! LCOV_EXCL_LINE
1500 : grow_N , & ! true growth rate in N currency (mmol/m^3/s) ! LCOV_EXCL_LINE
1501 : ! potU_Nit , & ! potential nitrate uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1502 : potU_Am , & ! potential ammonium uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1503 : U_Nit , & ! actual nitrate uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1504 : U_Am , & ! actual ammonium uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1505 : U_Sil , & ! actual silicon uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1506 : U_Fe , & ! actual iron uptake (umol/m^3/s) ! LCOV_EXCL_LINE
1507 : U_Nit_f , & ! fraction of Nit uptake due to each algal species ! LCOV_EXCL_LINE
1508 : U_Am_f , & ! fraction of Am uptake due to each algal species ! LCOV_EXCL_LINE
1509 : U_Sil_f , & ! fraction of Sil uptake due to each algal species ! LCOV_EXCL_LINE
1510 0 : U_Fe_f ! fraction of Fe uptake due to each algal species
1511 :
1512 : real (kind=dbl_kind) :: &
1513 : dTemp , & ! sea ice temperature minus sst (oC) < 0 ! LCOV_EXCL_LINE
1514 : U_Nit_tot , & ! actual nitrate uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1515 : U_Am_tot , & ! actual ammonium uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1516 : U_Sil_tot , & ! actual silicon uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1517 : U_Fe_tot , & ! actual iron uptake (umol/m^3/s) ! LCOV_EXCL_LINE
1518 : nitrif , & ! nitrification (mmol/m^3/s) ! LCOV_EXCL_LINE
1519 : mort_N , & ! total algal mortality (mmol N/m^3) ! LCOV_EXCL_LINE
1520 : mort_C , & ! total algal mortality (mmol C/m^3) ! LCOV_EXCL_LINE
1521 : graze_N , & ! total algae grazed (mmol N/m^3) ! LCOV_EXCL_LINE
1522 : graze_C , & ! total algae grazed (mmol C/m^3) ! LCOV_EXCL_LINE
1523 : exude_C , & ! total carbon exuded by algae (mmol C/m^3) ! LCOV_EXCL_LINE
1524 : resp_N , & ! total N in respiration (mmol N/m^3) ! LCOV_EXCL_LINE
1525 0 : growth_N ! total algal growth (mmol N/m^3)
1526 : ! fr_graze_p , & ! fraction of N grazed that becomes protein
1527 : ! ! (rest is assimilated) < (1-fr_graze_a)
1528 : ! ! and fr_graze_a*fr_graze_e becomes ammonia
1529 : ! fr_mort_p ! fraction of N mortality that becomes protein
1530 : ! ! < (1-fr_mort2min)
1531 :
1532 : real (kind=dbl_kind), dimension(n_algae) :: &
1533 : resp , & ! respiration (mmol/m^3/s) ! LCOV_EXCL_LINE
1534 : graze , & ! grazing (mmol/m^3/s) ! LCOV_EXCL_LINE
1535 0 : mort ! sum of mortality and excretion (mmol/m^3/s)
1536 :
1537 : ! source terms underscore s, removal underscore r
1538 :
1539 : real (kind=dbl_kind), dimension(n_algae) :: &
1540 : N_s , & ! net algal nitrogen sources (mmol/m^3) ! LCOV_EXCL_LINE
1541 0 : N_r ! net algal nitrogen removal (mmol/m^3)
1542 :
1543 : real (kind=dbl_kind), dimension(n_doc) :: &
1544 : DOC_r , & ! net DOC removal (mmol/m^3) ! LCOV_EXCL_LINE
1545 0 : DOC_s ! net DOC sources (mmol/m^3)
1546 :
1547 : real (kind=dbl_kind), dimension(n_don) :: &
1548 : DON_r , & ! net DON removal (mmol/m^3) ! LCOV_EXCL_LINE
1549 0 : DON_s ! net DON sources (mmol/m^3)
1550 :
1551 : real (kind=dbl_kind), dimension(n_fed) :: &
1552 : Fed_r_l , & ! removal due to loss of binding saccharids (nM) ! LCOV_EXCL_LINE
1553 : Fed_r , & ! net Fed removal (nM) ! LCOV_EXCL_LINE
1554 : Fed_s , & ! net Fed sources (nM) ! LCOV_EXCL_LINE
1555 0 : rFed ! ratio of dissolved Fe to tot Fed
1556 :
1557 : real (kind=dbl_kind), dimension(n_fep) :: &
1558 : Fep_r , & ! net Fep removal (nM) ! LCOV_EXCL_LINE
1559 : Fep_s , & ! net Fep sources (nM) ! LCOV_EXCL_LINE
1560 0 : rFep ! ratio of particulate Fe to tot Fep
1561 :
1562 : real (kind=dbl_kind) :: &
1563 : dN , & ! change in N (mmol/m^3) ! LCOV_EXCL_LINE
1564 : ! N_s_p , & ! algal nitrogen photosynthesis (mmol/m^3) ! LCOV_EXCL_LINE
1565 : ! N_r_g , & ! algal nitrogen losses to grazing (mmol/m^3) ! LCOV_EXCL_LINE
1566 : ! N_r_r , & ! algal nitrogen losses to respiration (mmol/m^3) ! LCOV_EXCL_LINE
1567 : ! N_r_mo , & ! algal nitrogen losses to mortality (mmol/m^3) ! LCOV_EXCL_LINE
1568 : Nit_s_n , & ! nitrate from nitrification (mmol/m^3) ! LCOV_EXCL_LINE
1569 : ! Nit_s_r , & ! nitrate from respiration (mmol/m^3) ! LCOV_EXCL_LINE
1570 : Nit_r_p , & ! nitrate uptake by algae (mmol/m^3) ! LCOV_EXCL_LINE
1571 : Nit_s , & ! net nitrate sources (mmol/m^3) ! LCOV_EXCL_LINE
1572 : Nit_r , & ! net nitrate removal (mmol/m^3) ! LCOV_EXCL_LINE
1573 : Am_s_e , & ! ammonium source from excretion (mmol/m^3) ! LCOV_EXCL_LINE
1574 : Am_s_r , & ! ammonium source from respiration (mmol/m^3) ! LCOV_EXCL_LINE
1575 : Am_s_mo , & ! ammonium source from mort/remin (mmol/m^3) ! LCOV_EXCL_LINE
1576 : Am_r_p , & ! ammonium uptake by algae (mmol/m^3) ! LCOV_EXCL_LINE
1577 : Am_s , & ! net ammonium sources (mmol/m^3) ! LCOV_EXCL_LINE
1578 : Am_r , & ! net ammonium removal (mmol/m^3) ! LCOV_EXCL_LINE
1579 : Sil_r_p , & ! silicon uptake by algae (mmol/m^3) ! LCOV_EXCL_LINE
1580 : Sil_r , & ! net silicon removal (mmol/m^3) ! LCOV_EXCL_LINE
1581 0 : Fe_r_p ! iron uptake by algae (nM)
1582 : ! DOC_r_c , & ! net doc removal from bacterial consumption (mmol/m^3)
1583 : ! doc_s_m , & ! protein source due to algal mortality (mmol/m^3) ! LCOV_EXCL_LINE
1584 : ! doc_s_g ! protein source due to grazing (mmol/m^3)
1585 :
1586 : real (kind=dbl_kind) :: &
1587 : DMSPd_s_r , & ! skl dissolved DMSP from respiration (mmol/m^3) ! LCOV_EXCL_LINE
1588 : DMSPd_s_mo, & ! skl dissolved DMSP from MBJ algal mortexc (mmol/m^3) ! LCOV_EXCL_LINE
1589 : DMSPd_r , & ! skl dissolved DMSP conversion (mmol/m^3) DMSPD_sk_r ! LCOV_EXCL_LINE
1590 : DMSPd_s , & ! net skl dissolved DMSP sources (mmol/m^3) ! LCOV_EXCL_LINE
1591 : DMS_s_c , & ! skl DMS source via conversion (mmol/m^3) ! LCOV_EXCL_LINE
1592 : DMS_r_o , & ! skl DMS losses due to oxidation (mmol/m^3) ! LCOV_EXCL_LINE
1593 : DMS_s , & ! net skl DMS sources (mmol/m^3) ! LCOV_EXCL_LINE
1594 : DMS_r , & ! net skl DMS removal (mmol/m^3) ! LCOV_EXCL_LINE
1595 : Fed_tot , & ! total dissolved iron from all sources (nM) ! LCOV_EXCL_LINE
1596 : Fed_tot_r , & ! total dissolved iron losses (nM) ! LCOV_EXCL_LINE
1597 : Fed_tot_s , & ! total dissolved iron sources (nM) ! LCOV_EXCL_LINE
1598 : Fep_tot , & ! total particulate iron from all sources (nM) ! LCOV_EXCL_LINE
1599 : ! Fep_tot_r , & ! total particulate iron losses (nM) ! LCOV_EXCL_LINE
1600 : Fep_tot_s , & ! total particulate iron sources (nM) ! LCOV_EXCL_LINE
1601 : Zoo_s_a , & ! N Losses due to zooplankton assimilation (mmol/m^3) ! LCOV_EXCL_LINE
1602 : Zoo_s_s , & ! N Losses due to grazing spillage (mmol/m^3) ! LCOV_EXCL_LINE
1603 : Zoo_s_m , & ! N Losses due to algal mortality (mmol/m^3) ! LCOV_EXCL_LINE
1604 0 : Zoo_s_b ! N losses due to bacterial recycling of DON (mmol/m^3)
1605 :
1606 : character(len=*),parameter :: subname='(algal_dyn)'
1607 :
1608 : !-----------------------------------------------------------------------
1609 : ! Initialize
1610 : !-----------------------------------------------------------------------
1611 :
1612 0 : conserve_N = .true.
1613 0 : Nin(:) = c0
1614 : ! Cin(:) = c0
1615 0 : chlin(:) = c0
1616 0 : DOCin(:) = c0
1617 : ! DICin(:) = c0
1618 0 : DONin(:) = c0
1619 0 : Fedin(:) = c0
1620 0 : Fepin(:) = c0
1621 0 : Nitin = c0
1622 0 : Amin = c0
1623 0 : Silin = c0
1624 : ! DMSPpin = c0
1625 0 : DMSPdin = c0
1626 0 : DMSin = c0
1627 0 : U_Am_tot = c0
1628 0 : U_Nit_tot = c0
1629 0 : U_Sil_tot = c0
1630 0 : U_Fe_tot = c0
1631 0 : U_Am_f(:) = c0
1632 0 : U_Nit_f(:) = c0
1633 0 : U_Sil_f(:) = c0
1634 0 : U_Fe_f(:) = c0
1635 0 : DOC_s(:) = c0
1636 0 : DOC_r(:) = c0
1637 : ! DOC_r_c = c0
1638 0 : nitrif = c0
1639 0 : mort_N = c0
1640 0 : mort_C = c0
1641 0 : graze_N = c0
1642 0 : graze_C = c0
1643 0 : exude_C = c0
1644 0 : resp_N = c0
1645 0 : growth_N = c0
1646 0 : Nit_r = c0
1647 0 : Am_s = c0
1648 0 : Am_r = c0
1649 0 : Sil_r = c0
1650 0 : Fed_r(:) = c0
1651 0 : Fed_s(:) = c0
1652 0 : Fep_r(:) = c0
1653 0 : Fep_s(:) = c0
1654 0 : DMSPd_s = c0
1655 0 : dTemp = min(T_bot-T_max,c0)
1656 0 : Fed_tot = c0
1657 0 : Fed_tot_r = c0
1658 0 : Fed_tot_s = c0
1659 0 : rFed(:) = c0
1660 0 : Fep_tot = c0
1661 : ! Fep_tot_r = c0
1662 0 : Fep_tot_s = c0
1663 0 : rFep(:) = c0
1664 :
1665 0 : Nitin = ltrcrn(nlt_bgc_Nit)
1666 0 : op_dep = c0
1667 0 : do k = 1, n_algae
1668 0 : Nin(k) = ltrcrn(nlt_bgc_N(k))
1669 0 : chlin(k) = R_chl2N(k)* Nin(k)
1670 0 : op_dep = op_dep + chlabs(k)*chlin(k)
1671 : enddo
1672 0 : if (tr_bgc_C) then
1673 : ! do k = 1, n_algae
1674 : ! Cin(k)= ltrcrn(nlt_bgc_C(k))
1675 : ! enddo
1676 0 : do k = 1, n_doc
1677 0 : DOCin(k)= ltrcrn(nlt_bgc_DOC(k))
1678 : enddo
1679 : ! do k = 1, n_dic
1680 : ! DICin(k)= ltrcrn(nlt_bgc_DIC(k))
1681 : ! enddo
1682 : endif
1683 0 : if (tr_bgc_Am) Amin = ltrcrn(nlt_bgc_Am)
1684 0 : if (tr_bgc_Sil) Silin = ltrcrn(nlt_bgc_Sil)
1685 0 : if (tr_bgc_DMS) then
1686 : ! DMSPpin = ltrcrn(nlt_bgc_DMSPp)
1687 0 : DMSPdin = ltrcrn(nlt_bgc_DMSPd)
1688 0 : DMSin = ltrcrn(nlt_bgc_DMS)
1689 : endif
1690 : ! if (tr_bgc_PON) then
1691 : ! PONin = c0
1692 : ! PONin = ltrcrn(nlt_bgc_PON)
1693 : ! endif
1694 0 : if (tr_bgc_DON) then
1695 0 : do k = 1, n_don
1696 0 : DONin(k) = ltrcrn(nlt_bgc_DON(k))
1697 : enddo
1698 : endif
1699 0 : if (tr_bgc_Fe ) then
1700 0 : do k = 1, n_fed
1701 0 : Fedin(k) = ltrcrn(nlt_bgc_Fed(k))
1702 : enddo
1703 0 : do k = 1, n_fep
1704 0 : Fepin(k) = ltrcrn(nlt_bgc_Fep(k))
1705 : enddo
1706 : endif
1707 :
1708 : !-----------------------------------------------------------------------
1709 : ! Total iron from all pools
1710 : !-----------------------------------------------------------------------
1711 :
1712 0 : do k = 1,n_fed
1713 0 : Fed_tot = Fed_tot + Fedin(k)
1714 : enddo
1715 0 : do k = 1,n_fep
1716 0 : Fep_tot = Fep_tot + Fepin(k)
1717 : enddo
1718 0 : if (Fed_tot > puny) then
1719 0 : do k = 1,n_fed
1720 0 : rFed(k) = Fedin(k)/Fed_tot
1721 : enddo
1722 : endif
1723 0 : if (Fep_tot > puny) then
1724 0 : do k = 1,n_fep
1725 0 : rFep(k) = Fepin(k)/Fep_tot
1726 : enddo
1727 : endif
1728 :
1729 : !-----------------------------------------------------------------------
1730 : ! Light limitation (op_dep) defined above
1731 : !-----------------------------------------------------------------------
1732 :
1733 0 : if (op_dep > op_dep_min .and. .not. dEdd_algae) then
1734 0 : Iavg_loc = fswthru * (c1 - exp(-op_dep)) / op_dep
1735 : else
1736 0 : Iavg_loc = fswthru
1737 : endif
1738 :
1739 0 : do k = 1, n_algae
1740 : ! With light inhibition ! Maybe include light inhibition for diatoms but phaeocystis
1741 0 : L_lim = (c1 - exp(-alpha2max_low(k)*Iavg_loc)) * exp(-beta2max(k)*Iavg_loc)
1742 :
1743 : ! Without light inhibition
1744 : ! L_lim(k) = (c1 - exp(-alpha2max_low(k)*Iavg_loc))
1745 :
1746 : !-----------------------------------------------------------------------
1747 : ! Nutrient limitation
1748 : !-----------------------------------------------------------------------
1749 :
1750 0 : Nit_lim(k) = Nitin/(Nitin + K_Nit(k))
1751 0 : Am_lim(k) = c0
1752 0 : N_lim(k) = Nit_lim(k)
1753 0 : if (tr_bgc_Am) then
1754 0 : Am_lim(k) = Amin/(Amin + K_Am(k))
1755 0 : N_lim(k) = min(c1, Nit_lim(k) + Am_lim(k))
1756 : endif
1757 0 : Sil_lim(k) = c1
1758 0 : if (tr_bgc_Sil .and. K_Sil(k) > c0) Sil_lim(k) = Silin/(Silin + K_Sil(k))
1759 :
1760 : !-----------------------------------------------------------------------
1761 : ! Iron limitation
1762 : !-----------------------------------------------------------------------
1763 :
1764 0 : Fe_lim(k) = c1
1765 0 : if (tr_bgc_Fe .and. K_Fe (k) > c0) Fe_lim (k) = Fed_tot/(Fed_tot + K_Fe(k))
1766 :
1767 : !----------------------------------------------------------------------------
1768 : ! Growth and uptake computed within the bottom layer
1769 : ! Note here per A93 discussions and MBJ model, salinity is a universal
1770 : ! restriction. Comparison with available column nutrients inserted
1771 : ! but in tests had no effect.
1772 : ! Primary production reverts to SE form, see MBJ below and be careful
1773 : !----------------------------------------------------------------------------
1774 :
1775 0 : growmax_N(k) = mu_max(k) / secday * exp(grow_Tdep(k) * dTemp)* Nin(k) *fsal
1776 0 : grow_N(k) = min(L_lim(k), N_lim(k), Sil_lim(k), Fe_lim(k)) * growmax_N(k)
1777 : ! potU_Nit(k) = Nit_lim(k)* growmax_N(k)
1778 0 : potU_Am(k) = Am_lim(k)* growmax_N(k)
1779 0 : U_Am(k) = min(grow_N(k), potU_Am(k))
1780 0 : U_Nit(k) = grow_N(k) - U_Am(k)
1781 0 : U_Sil(k) = R_Si2N(k) * grow_N(k)
1782 0 : U_Fe (k) = R_Fe2N(k) * grow_N(k)
1783 :
1784 0 : U_Am_tot = U_Am_tot + U_Am(k)
1785 0 : U_Nit_tot = U_Nit_tot + U_Nit(k)
1786 0 : U_Sil_tot = U_Sil_tot + U_Sil(k)
1787 0 : U_Fe_tot = U_Fe_tot + U_Fe(k)
1788 : enddo
1789 0 : do k = 1, n_algae
1790 0 : if (U_Am_tot > c0) U_Am_f(k) = U_Am(k)/U_Am_tot
1791 0 : if (U_Nit_tot > c0) U_Nit_f(k) = U_Nit(k)/U_Nit_tot
1792 0 : if (U_Sil_tot > c0) U_Sil_f(k) = U_Sil(k)/U_Sil_tot
1793 0 : if (U_Fe_tot > c0) U_Fe_f(k) = U_Fe(k)/U_Fe_tot
1794 : enddo
1795 :
1796 0 : if (tr_bgc_Sil) U_Sil_tot = min(U_Sil_tot, max_loss * Silin/dt)
1797 0 : if (tr_bgc_Fe) U_Fe_tot = min(U_Fe_tot, max_loss * Fed_tot/dt)
1798 0 : U_Nit_tot = min(U_Nit_tot, max_loss * Nitin/dt)
1799 0 : U_Am_tot = min(U_Am_tot, max_loss * Amin/dt)
1800 :
1801 0 : do k = 1, n_algae
1802 0 : U_Am(k) = U_Am_f(k)*U_Am_tot
1803 0 : U_Nit(k) = U_Nit_f(k)*U_Nit_tot
1804 0 : U_Sil(k) = U_Sil_f(k)*U_Sil_tot
1805 0 : U_Fe(k) = U_Fe_f(k)*U_Fe_tot
1806 :
1807 0 : if (R_Si2N(k) > c0) then
1808 0 : grow_N(k) = min(U_Sil(k)/R_Si2N(k),U_Nit(k) + U_Am(k), U_Fe(k)/R_Fe2N(k))
1809 : else
1810 0 : grow_N(k) = min(U_Nit(k) + U_Am(k),U_Fe(k)/R_Fe2N(k))
1811 : endif
1812 :
1813 0 : fr_Am(k) = c0
1814 0 : if (tr_bgc_Am) then
1815 0 : fr_Am(k) = p5
1816 0 : if (grow_N(k) > c0) fr_Am(k) = min(U_Am(k)/grow_N(k), c1)
1817 : endif
1818 0 : fr_Nit(k) = c1 - fr_Am(k)
1819 0 : U_Nit(k) = fr_Nit(k) * grow_N(k)
1820 0 : U_Am(k) = fr_Am(k) * grow_N(k)
1821 0 : U_Sil(k) = R_Si2N(k) * grow_N(k)
1822 0 : U_Fe (k) = R_Fe2N(k) * grow_N(k)
1823 :
1824 : !-----------------------------------------------------------------------
1825 : ! Define reaction terms
1826 : !-----------------------------------------------------------------------
1827 :
1828 : ! Since the framework remains incomplete at this point,
1829 : ! it is assumed as a starting expedient that
1830 : ! DMSP loss to melting results in 10% conversion to DMS
1831 : ! which is then given a ten day removal constant.
1832 : ! Grazing losses are channeled into rough spillage and assimilation
1833 : ! then following ammonia there is some recycling.
1834 :
1835 : !--------------------------------------------------------------------
1836 : ! Algal reaction term
1837 : ! N_react = (grow_N*(c1 - fr_graze-fr_resp) - mort)*dt
1838 : !--------------------------------------------------------------------
1839 :
1840 0 : resp(k) = fr_resp * grow_N(k)
1841 0 : graze(k) = fr_graze(k) * grow_N(k)
1842 0 : mort(k) = min(max_loss * Nin(k)/dt, &
1843 0 : mort_pre(k)*exp(mort_Tdep(k)*dTemp) * Nin(k)/secday)
1844 :
1845 : ! history variables
1846 0 : grow_alg(k) = grow_N(k)
1847 0 : upNOn(k) = U_Nit(k)
1848 0 : upNHn(k) = U_Am(k)
1849 :
1850 : ! N_s_p = grow_N(k) * dt
1851 : ! N_r_g = graze(k) * dt
1852 : ! N_r_r = resp(k) * dt
1853 : ! N_r_mo = mort(k) * dt
1854 0 : N_s(k) = (c1- fr_resp - fr_graze(k)) * grow_N(k) *dt !N_s_p
1855 0 : N_r(k) = mort(k) * dt !N_r_g + N_r_mo + N_r_r
1856 :
1857 0 : graze_N = graze_N + graze(k)
1858 0 : graze_C = graze_C + R_C2N(k)*graze(k)
1859 0 : mort_N = mort_N + mort(k)
1860 0 : mort_C = mort_C + R_C2N(k)*mort(k)
1861 0 : resp_N = resp_N + resp(k)
1862 0 : growth_N = growth_N + grow_N(k)
1863 :
1864 : enddo ! n_algae
1865 : !--------------------------------------------------------------------
1866 : ! Ammonium source: algal grazing, respiration, and mortality
1867 : !--------------------------------------------------------------------
1868 :
1869 0 : Am_s_e = graze_N*(c1-fr_graze_s)*fr_graze_e*dt
1870 0 : Am_s_mo = mort_N*fr_mort2min*dt
1871 0 : Am_s_r = resp_N*dt
1872 0 : Am_s = Am_s_r + Am_s_e + Am_s_mo
1873 :
1874 : !--------------------------------------------------------------------
1875 : ! Nutrient net loss terms: algal uptake
1876 : !--------------------------------------------------------------------
1877 :
1878 0 : do k = 1, n_algae
1879 0 : Am_r_p = U_Am(k) * dt
1880 0 : Am_r = Am_r + Am_r_p
1881 0 : Nit_r_p = U_Nit(k) * dt
1882 0 : Nit_r = Nit_r + Nit_r_p
1883 0 : Sil_r_p = U_Sil(k) * dt
1884 0 : Sil_r = Sil_r + Sil_r_p
1885 0 : Fe_r_p = U_Fe (k) * dt
1886 0 : Fed_tot_r = Fed_tot_r + Fe_r_p
1887 0 : exude_C = exude_C + k_exude(k)* R_C2N(k)*Nin(k) / secday
1888 : enddo
1889 :
1890 : !--------------------------------------------------------------------
1891 : ! nitrification
1892 : !--------------------------------------------------------------------
1893 :
1894 0 : nitrif = k_nitrif /secday * Amin
1895 0 : Am_r = Am_r + nitrif*dt
1896 0 : Nit_s_n = nitrif * dt ! source from NH4
1897 0 : Nit_s = Nit_s_n
1898 :
1899 : !--------------------------------------------------------------------
1900 : ! PON: currently using PON to shadow nitrate
1901 : !
1902 : ! N Losses are counted in Zoo. These arise from mortality not
1903 : ! remineralized (Zoo_s_m), assimilated grazing not excreted (Zoo_s_a),
1904 : !spilled N not going to DON (Zoo_s_s) and bacterial recycling
1905 : ! of DON (Zoo_s_b).
1906 : !--------------------------------------------------------------------
1907 :
1908 0 : if (tr_bgc_Am) then
1909 0 : Zoo_s_a = graze_N*(c1-fr_graze_e)*(c1-fr_graze_s) *dt
1910 0 : Zoo_s_s = graze_N*fr_graze_s*dt
1911 0 : Zoo_s_m = mort_N*dt - Am_s_mo
1912 : else
1913 0 : Zoo_s_a = graze_N*dt*(c1-fr_graze_s)
1914 0 : Zoo_s_s = graze_N*fr_graze_s*dt
1915 0 : Zoo_s_m = mort_N*dt
1916 : endif
1917 :
1918 0 : Zoo_s_b = c0
1919 :
1920 : !--------------------------------------------------------------------
1921 : ! DON (n_don = 1)
1922 : ! Proteins
1923 : !--------------------------------------------------------------------
1924 :
1925 0 : DON_r(:) = c0
1926 0 : DON_s(:) = c0
1927 :
1928 0 : if (tr_bgc_DON) then
1929 0 : do n = 1, n_don
1930 0 : DON_r(n) = kn_bac(n)/secday * DONin(n) * dt
1931 0 : DON_s(n) = graze_N*f_don(n)*fr_graze_s * dt
1932 0 : Zoo_s_s = Zoo_s_s - DON_s(n)
1933 0 : Zoo_s_b = Zoo_s_b + DON_r(n)*(c1-f_don_Am(n))
1934 : !Am_s = Am_s + DON_r(n)*f_don_Am(n)
1935 : enddo
1936 : endif
1937 :
1938 0 : Zoo = Zoo_s_a + Zoo_s_s + Zoo_s_m + Zoo_s_b
1939 :
1940 : !--------------------------------------------------------------------
1941 : ! DOC
1942 : ! polysaccharids, lipids
1943 : !--------------------------------------------------------------------
1944 :
1945 0 : do n = 1, n_doc
1946 :
1947 0 : DOC_r(n) = k_bac(n)/secday * DOCin(n) * dt
1948 0 : DOC_s(n) = f_doc(n)*(fr_graze_s *graze_C + mort_C)*dt &
1949 0 : + f_exude(n)*exude_C
1950 : enddo
1951 :
1952 : !--------------------------------------------------------------------
1953 : ! Iron sources from remineralization (follows ammonium but reduced)
1954 : ! only Fed_s(1) has remineralized sources
1955 : !--------------------------------------------------------------------
1956 :
1957 0 : Fed_s(1) = Fed_s(1) + Am_s * R_Fe2N(1) * fr_dFe ! remineralization source
1958 :
1959 : !--------------------------------------------------------------------
1960 : ! Conversion to dissolved Fe from Particulate requires DOC(1)
1961 : ! Otherwise the only source of dFe is from remineralization
1962 : !--------------------------------------------------------------------
1963 :
1964 0 : if (tr_bgc_C .and. tr_bgc_Fe) then
1965 0 : if (DOCin(1) > c0) then
1966 0 : if (Fed_tot/DOCin(1) > max_dfe_doc1) then
1967 0 : do n = 1,n_fed ! low saccharid:dFe ratio leads to
1968 0 : Fed_r_l(n) = Fedin(n)/t_iron_conv*dt/secday ! loss of bioavailable Fe to particulate fraction
1969 0 : Fep_tot_s = Fep_tot_s + Fed_r_l(n)
1970 0 : Fed_r(n) = Fed_r_l(n) ! removal due to particulate scavenging
1971 : enddo
1972 0 : do n = 1,n_fep
1973 0 : Fep_s(n) = rFep(n)* Fep_tot_s ! source from dissolved Fe
1974 : enddo
1975 0 : elseif (Fed_tot/DOCin(1) < max_dfe_doc1) then
1976 0 : do n = 1,n_fep ! high saccharid:dFe ratio leads to
1977 0 : Fep_r(n) = Fepin(n)/t_iron_conv*dt/secday ! gain of bioavailable Fe from particulate fraction
1978 0 : Fed_tot_s = Fed_tot_s + Fep_r(n)
1979 : enddo
1980 0 : do n = 1,n_fed
1981 0 : Fed_s(n) = Fed_s(n) + rFed(n)* Fed_tot_s ! source from particulate Fe
1982 : enddo
1983 : endif
1984 : endif !Docin(1) > c0
1985 0 : elseif (tr_bgc_Fe) then
1986 0 : do n = 1,n_fed
1987 0 : Fed_r(n) = Fed_r(n) + rFed(n)*Fed_tot_r ! scavenging + uptake
1988 : enddo
1989 :
1990 : ! source from algal mortality/grazing and fraction of remineralized nitrogen that does
1991 : ! not become immediately bioavailable
1992 :
1993 0 : do n = 1,n_fep
1994 0 : Fep_s(n) = Fep_s(n) + rFep(n)* (Am_s * R_Fe2N(1) * (c1-fr_dFe))
1995 : enddo ! losses not direct to Fed
1996 : endif
1997 :
1998 : !--------------------------------------------------------------------
1999 : ! Sulfur cycle begins here
2000 : !--------------------------------------------------------------------
2001 : ! Grazing losses are channeled into rough spillage and assimilation
2002 : ! then onward and the MBJ mortality channel is included
2003 : ! It is assumed as a starting expedient that
2004 : ! DMSP loss to melting gives partial conversion to DMS in product layer
2005 : ! which then undergoes Stefels removal.
2006 :
2007 : !--------------------------------------------------------------------
2008 : ! DMSPd reaction term (DMSPd conversion is outside the algal loop)
2009 : ! DMSPd_react = R_S2N*((fr_graze_s+fr_excrt_2S*fr_graze_e*fr_graze_a)
2010 : ! *fr_graze*grow_N + fr_mort2min*mort)*dt
2011 : ! - [\DMSPd]/t_sk_conv*dt
2012 : !--------------------------------------------------------------------
2013 0 : do k = 1,n_algae
2014 0 : DMSPd_s_r = fr_resp_s * R_S2N(k) * resp(k) * dt !respiration fraction to DMSPd
2015 0 : DMSPd_s_mo= fr_mort2min * R_S2N(k)* mort(k) * dt !mortality and extracellular excretion
2016 0 : DMSPd_s = DMSPd_s + DMSPd_s_r + DMSPd_s_mo
2017 : enddo
2018 0 : DMSPd_r = (c1/t_sk_conv) * (c1/secday) * (DMSPdin) * dt
2019 :
2020 : !--------------------------------------------------------------------
2021 : ! DMS reaction term + DMSPd loss term
2022 : ! DMS_react = ([\DMSPd]*y_sk_DMS/t_sk_conv - c1/t_sk_ox *[\DMS])*dt
2023 : !--------------------------------------------------------------------
2024 :
2025 0 : DMS_s_c = y_sk_DMS * DMSPd_r
2026 0 : DMS_r_o = DMSin * dt / (t_sk_ox * secday)
2027 0 : DMS_s = DMS_s_c
2028 0 : DMS_r = DMS_r_o
2029 :
2030 : !-----------------------------------------------------------------------
2031 : ! Load reaction array
2032 : !-----------------------------------------------------------------------
2033 :
2034 0 : dN = c0
2035 0 : do k = 1,n_algae
2036 0 : reactb(nlt_bgc_N(k)) = N_s(k) - N_r(k)
2037 0 : dN = dN + reactb(nlt_bgc_N(k))
2038 : enddo
2039 0 : if (tr_bgc_C) then
2040 : ! do k = 1,n_algae
2041 : ! reactb(nlt_bgc_C(k)) = R_C2N(k)*reactb(nlt_bgc_N(k))
2042 : ! enddo
2043 0 : do k = 1,n_doc
2044 0 : reactb(nlt_bgc_DOC(k))= DOC_s(k) - DOC_r(k)
2045 : enddo
2046 : endif
2047 0 : reactb(nlt_bgc_Nit) = Nit_s - Nit_r
2048 0 : dN = dN + reactb(nlt_bgc_Nit)
2049 0 : if (tr_bgc_Am) then
2050 0 : reactb(nlt_bgc_Am) = Am_s - Am_r
2051 0 : dN = dN + reactb(nlt_bgc_Am)
2052 : endif
2053 0 : if (tr_bgc_Sil) then
2054 0 : reactb(nlt_bgc_Sil) = - Sil_r
2055 : endif
2056 0 : if (tr_bgc_DON) then
2057 0 : do k = 1,n_don
2058 0 : reactb(nlt_bgc_DON(k))= DON_s(k) - DON_r(k)
2059 0 : dN = dN + reactb(nlt_bgc_DON(k))
2060 : enddo
2061 : endif
2062 0 : if (tr_bgc_Fe ) then
2063 0 : do k = 1,n_fed
2064 0 : reactb(nlt_bgc_Fed(k))= Fed_s (k) - Fed_r (k)
2065 : enddo
2066 0 : do k = 1,n_fep
2067 0 : reactb(nlt_bgc_Fep(k))= Fep_s (k) - Fep_r (k)
2068 : enddo
2069 : endif
2070 0 : if (tr_bgc_DMS) then
2071 0 : reactb(nlt_bgc_DMSPd) = DMSPd_s - DMSPd_r
2072 0 : reactb(nlt_bgc_DMS) = DMS_s - DMS_r
2073 : endif
2074 0 : Nerror = dN + Zoo
2075 : ! if (abs(Nerror) > max(reactb(:))*1.0e-5) then
2076 : ! conserve_N = .false.
2077 : ! write(warnstr,*) subname, 'Conservation error!'
2078 : ! call icepack_warnings_add(warnstr)
2079 : ! write(warnstr,*) subname, 'Nerror,dN, DONin(1),kn_bac(1),secday,dt,n_doc'
2080 : ! call icepack_warnings_add(warnstr)
2081 : ! write(warnstr,*) subname, Nerror,dN, DONin(1),kn_bac(1),secday,dt,n_doc
2082 : ! call icepack_warnings_add(warnstr)
2083 : ! write(warnstr,*) subname, 'reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(1)),reactb(nlt_bgc_N(2)'
2084 : ! call icepack_warnings_add(warnstr)
2085 : ! write(warnstr,*) subname, reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(1)),reactb(nlt_bgc_N(2))
2086 : ! call icepack_warnings_add(warnstr)
2087 : ! write(warnstr,*) subname, 'reactb(nlt_bgc_Am),reactb(nlt_bgc_DON(1)), DON_r(1),DON_s(1)'
2088 : ! call icepack_warnings_add(warnstr)
2089 : ! write(warnstr,*) subname, reactb(nlt_bgc_Am),reactb(nlt_bgc_DON(1)),DON_r(1),DON_s(1)
2090 : ! call icepack_warnings_add(warnstr)
2091 : ! write(warnstr,*) subname, 'Zoo:',Zoo
2092 : ! endif
2093 :
2094 0 : end subroutine algal_dyn
2095 :
2096 : !=======================================================================
2097 : !
2098 : ! Find ice-ocean flux when ice is thin and internal dynamics/reactions are
2099 : ! assumed to be zero
2100 : !
2101 : ! authors Nicole Jeffery, LANL
2102 :
2103 0 : subroutine thin_ice_flux (hin, hin_old, Cin, flux_o_tot, &
2104 : source, dt, nblyr, & ! LCOV_EXCL_LINE
2105 : ocean_bio)
2106 :
2107 : integer (kind=int_kind), intent(in) :: &
2108 : nblyr ! number of bio layers
2109 :
2110 : real (kind=dbl_kind), dimension(nblyr+1), intent(inout) :: &
2111 : Cin ! initial concentration*hin_old*phin
2112 :
2113 : real (kind=dbl_kind), intent(in) :: &
2114 : hin_old , & ! brine thickness (m) ! LCOV_EXCL_LINE
2115 : hin , & ! new brine thickness (m) ! LCOV_EXCL_LINE
2116 : dt , & ! time step ! LCOV_EXCL_LINE
2117 : source , & ! atm, ocean, dust flux (mmol/m^2) ! LCOV_EXCL_LINE
2118 : ocean_bio ! ocean tracer concentration (mmol/m^3)
2119 :
2120 : real (kind=dbl_kind), intent(inout) :: &
2121 : flux_o_tot ! tracer flux, gravity+molecular drainage flux ,
2122 : ! and boundary flux to ocean (mmol/m^2/s)
2123 : ! positive into the ocean
2124 :
2125 : ! local variables
2126 :
2127 : integer (kind=int_kind) :: &
2128 : k ! vertical biology layer index
2129 :
2130 : real (kind=dbl_kind) :: &
2131 : sum_bio, & ! initial bio mass (mmol/m^2) ! LCOV_EXCL_LINE
2132 : zspace, & ! 1/nblyr ! LCOV_EXCL_LINE
2133 : dC, & ! added ocean bio mass (mmol/m^2) ! LCOV_EXCL_LINE
2134 0 : dh ! change in thickness (m)
2135 :
2136 : character(len=*),parameter :: subname='(thin_ice_flux)'
2137 :
2138 0 : zspace = c1/real(nblyr,kind=dbl_kind)
2139 :
2140 0 : dC = c0
2141 0 : sum_bio = c0
2142 0 : dh = hin-hin_old
2143 :
2144 0 : if (dh .le. c0) then ! keep the brine concentration fixed
2145 0 : sum_bio = (Cin(1)+Cin(nblyr+1))/hin_old*zspace*p5
2146 0 : Cin(1) = Cin(1)/hin_old*hin
2147 0 : Cin(nblyr+1) = Cin(nblyr+1)/hin_old*hin
2148 0 : do k = 2, nblyr
2149 0 : sum_bio = sum_bio + Cin(k)/hin_old*zspace
2150 0 : Cin(k) = Cin(k)/hin_old*hin + dC
2151 : enddo
2152 : else
2153 0 : dC = dh*ocean_bio
2154 0 : do k = 1, nblyr+1
2155 0 : Cin(k) = Cin(k) + dC
2156 : enddo
2157 : endif
2158 :
2159 0 : flux_o_tot = - dh*sum_bio/dt - dC/dt + source/dt
2160 :
2161 0 : end subroutine thin_ice_flux
2162 :
2163 : !=======================================================================
2164 : !
2165 : ! Compute matrix elements for the low order solution of FEM-FCT scheme
2166 : ! Predictor step
2167 : !
2168 : ! July, 2014 by N. Jeffery
2169 : !
2170 0 : subroutine compute_FCT_matrix (C_in, sbdiag, dt, nblyr, &
2171 : diag, spdiag, rhs, bgrid, & ! LCOV_EXCL_LINE
2172 : darcyV, dhtop, dhbot, & ! LCOV_EXCL_LINE
2173 : iphin_N, iDin, hbri_old, & ! LCOV_EXCL_LINE
2174 : atm_add, bphin_N, & ! LCOV_EXCL_LINE
2175 : C_top, C_bot, Qbot, Qtop, & ! LCOV_EXCL_LINE
2176 : Sink_bot, Sink_top, & ! LCOV_EXCL_LINE
2177 0 : D_sbdiag, D_spdiag, ML)
2178 :
2179 : integer (kind=int_kind), intent(in) :: &
2180 : nblyr ! number of bio layers
2181 :
2182 : real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
2183 : C_in ! Initial (bulk) concentration*hbri_old (mmol/m^2)
2184 : ! conserved quantity on i_grid
2185 :
2186 : real (kind=dbl_kind), intent(in) :: &
2187 : dt ! time step
2188 :
2189 : real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
2190 : iDin ! Diffusivity on the igrid (1/s)
2191 :
2192 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
2193 : iphin_N ! Porosity with min condition on igrid
2194 :
2195 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
2196 : bphin_N, & ! Porosity with min condition on igrid ! LCOV_EXCL_LINE
2197 : bgrid
2198 :
2199 : real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
2200 : sbdiag , & ! sub-diagonal matrix elements ! LCOV_EXCL_LINE
2201 : diag , & ! diagonal matrix elements ! LCOV_EXCL_LINE
2202 : spdiag , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
2203 : rhs , & ! rhs of tri-diagonal matrix eqn. ! LCOV_EXCL_LINE
2204 : ML, & ! lumped mass matrix ! LCOV_EXCL_LINE
2205 : D_spdiag, D_sbdiag ! artificial diffusion matrix
2206 :
2207 : real (kind=dbl_kind), intent(in) :: &
2208 : dhtop , & ! Change in brine top (m) ! LCOV_EXCL_LINE
2209 : dhbot , & ! Change in brine bottom (m) ! LCOV_EXCL_LINE
2210 : hbri_old , & ! brine height (m) ! LCOV_EXCL_LINE
2211 : atm_add , & ! atm-ice flux ! LCOV_EXCL_LINE
2212 : C_top , & ! bulk surface source (mmol/m^2) ! LCOV_EXCL_LINE
2213 : C_bot , & ! bulk bottom source (mmol/m^2) ! LCOV_EXCL_LINE
2214 : darcyV ! Darcy velocity (m/s
2215 :
2216 : real (kind=dbl_kind), intent(inout) :: & ! positive into ice
2217 : Qtop , & ! top flux source (mmol/m^2/s) ! LCOV_EXCL_LINE
2218 : Qbot , & ! bottom flux source (mmol/m^2/s) ! LCOV_EXCL_LINE
2219 : Sink_bot , & ! rest of bottom flux (sink or source) (mmol/m^2/s) ! LCOV_EXCL_LINE
2220 : Sink_top ! rest oftop flux (sink or source) (mmol/m^2/s)
2221 :
2222 : ! local variables
2223 :
2224 : real (kind=dbl_kind) :: &
2225 0 : vel, vel2, dphi_dx, zspace
2226 :
2227 : integer (kind=int_kind) :: &
2228 : k ! vertical index
2229 :
2230 : real (kind=dbl_kind), dimension (nblyr+1) :: &
2231 : Q_top, Q_bot, & ! surface and bottom source ! LCOV_EXCL_LINE
2232 : K_diag, K_spdiag, K_sbdiag, & ! advection matrix ! LCOV_EXCL_LINE
2233 : S_diag, S_spdiag, S_sbdiag, & ! diffusion matrix ! LCOV_EXCL_LINE
2234 0 : D_diag, iDin_phi
2235 :
2236 : real (kind=dbl_kind), dimension (nblyr) :: &
2237 0 : kvector1, kvectorn1
2238 :
2239 : character(len=*),parameter :: subname='(compute_FCT_matrix)'
2240 :
2241 : !---------------------------------------------------------------------
2242 : ! Diag (jj) solve for j = 1:nblyr+1
2243 : ! spdiag(j) == (j,j+1) solve for j = 1:nblyr otherwise 0
2244 : ! sbdiag(j) == (j,j-1) solve for j = 2:nblyr+1 otherwise 0
2245 : !---------------------------------------------------------------------
2246 0 : kvector1(:) = c0
2247 0 : kvector1(1) = c1
2248 0 : kvectorn1(:) = c1
2249 0 : kvectorn1(1) = c0
2250 :
2251 0 : zspace = c1/real(nblyr,kind=dbl_kind)
2252 0 : Qbot = c0
2253 0 : Qtop = c0
2254 0 : Sink_bot = c0
2255 0 : Sink_top = c0
2256 :
2257 : ! compute the lumped mass matrix
2258 :
2259 0 : ML(:) = zspace
2260 0 : ML(1) = zspace/c2
2261 0 : ML(nblyr+1) = zspace/c2
2262 :
2263 : ! compute matrix K: K_diag , K_sbdiag, K_spdiag
2264 : ! compute matrix S: S_diag, S_sbdiag, S_spdiag
2265 :
2266 0 : K_diag(:) = c0
2267 0 : D_diag(:) = c0
2268 0 : D_spdiag(:) = c0
2269 0 : D_sbdiag(:) = c0
2270 0 : K_spdiag(:) = c0
2271 0 : K_sbdiag(:) = c0
2272 0 : S_diag(:) = c0
2273 0 : S_spdiag(:) = c0
2274 0 : S_sbdiag(:) = c0
2275 0 : iDin_phi(:) = c0
2276 :
2277 0 : iDin_phi(1) = c0 !element 1
2278 0 : iDin_phi(nblyr+1) = iDin(nblyr+1)/iphin_N(nblyr+1) !outside ice at bottom
2279 0 : iDin_phi(nblyr) = p5*(iDin(nblyr+1)/iphin_N(nblyr+1)+iDin(nblyr)/iphin_N(nblyr))
2280 :
2281 0 : vel = (bgrid(2)*dhbot - (bgrid(2)-c1)*dhtop)/dt+darcyV/bphin_N(2)
2282 0 : K_diag(1) = p5*vel/hbri_old
2283 0 : dphi_dx = (iphin_N(nblyr+1) - iphin_N(nblyr))/(zspace)
2284 0 : vel = (bgrid(nblyr+1)*dhbot - (bgrid(nblyr+1)-c1)*dhtop)/dt +darcyV/bphin_N(nblyr+1)
2285 0 : vel = vel/hbri_old
2286 0 : vel2 = (dhbot/hbri_old/dt +darcyV/hbri_old)
2287 0 : K_diag(nblyr+1) = min(c0, vel2) - iDin_phi(nblyr+1)/(zspace+ grid_o/hbri_old) &
2288 0 : + p5*(-vel + iDin_phi(nblyr)/bphin_N(nblyr+1)*dphi_dx)
2289 :
2290 0 : do k = 1, nblyr-1
2291 0 : vel = (bgrid(k+1)*dhbot - (bgrid(k+1)-c1)*dhtop)/dt+darcyV/bphin_N(k+1)
2292 0 : iDin_phi(k+1) = p5*(iDin(k+2)/iphin_N(k+2) + iDin(k+1)/iphin_N(k+1))
2293 0 : dphi_dx = (iphin_N(k+1) - iphin_N(k))/(zspace)
2294 0 : K_spdiag(k)= p5*(vel/hbri_old - &
2295 0 : iDin_phi(k)/(bphin_N(k+1))*dphi_dx)
2296 :
2297 0 : vel = (bgrid(k+1)*dhbot - (bgrid(k+1)-c1)*dhtop)/dt +darcyV/bphin_N(k+1)
2298 0 : dphi_dx = c0
2299 0 : dphi_dx = kvectorn1(k)*(iphin_N(k+1) - iphin_N(k))/(zspace)
2300 0 : K_sbdiag(k+1)= -p5*(vel/hbri_old- &
2301 0 : iDin_phi(k)/bphin_N(k+1)*dphi_dx)
2302 0 : K_diag(k) = K_diag(1)*kvector1(k) + (K_spdiag(k) + K_sbdiag(k))*kvectorn1(k)
2303 :
2304 0 : S_diag(k+1) = -(iDin_phi(k)+ iDin_phi(k+1))/zspace
2305 0 : S_spdiag(k) = iDin_phi(k)/zspace
2306 0 : S_sbdiag(k+1) = iDin_phi(k)/zspace
2307 : enddo
2308 :
2309 : ! k = nblyr
2310 :
2311 0 : vel = (bgrid(nblyr+1)*dhbot - (bgrid(nblyr+1)-c1)*dhtop)/dt+darcyV/bphin_N(nblyr+1)
2312 0 : dphi_dx = (iphin_N(nblyr+1) - iphin_N(nblyr))/(zspace)
2313 0 : K_spdiag(nblyr)= p5*(vel/hbri_old - &
2314 0 : iDin_phi(nblyr)/(bphin_N(nblyr+1))*dphi_dx)
2315 0 : vel = (bgrid(nblyr+1)*dhbot - (bgrid(nblyr+1)-c1)*dhtop)/dt +darcyV/bphin_N(nblyr+1)
2316 0 : dphi_dx = kvectorn1(nblyr)*(iphin_N(nblyr+1) - iphin_N(nblyr))/(zspace)
2317 0 : K_sbdiag(nblyr+1)= -p5*(vel/hbri_old- &
2318 0 : iDin_phi(nblyr)/bphin_N(nblyr+1)*dphi_dx)
2319 0 : K_diag(nblyr) = K_spdiag(nblyr) + K_sbdiag(nblyr)
2320 0 : S_diag(nblyr+1) = -iDin_phi(nblyr)/zspace
2321 0 : S_spdiag(nblyr) = iDin_phi(nblyr)/zspace
2322 0 : S_sbdiag(nblyr+1) = iDin_phi(nblyr)/zspace
2323 :
2324 : ! compute matrix artificial D: D_spdiag, D_diag (D_spdiag(k) = D_sbdiag(k+1))
2325 :
2326 0 : do k = 1,nblyr
2327 0 : D_spdiag(k) = max(-K_spdiag(k), c0, -K_sbdiag(k+1))
2328 0 : D_sbdiag(k+1) = D_spdiag(k)
2329 : enddo
2330 0 : do k = 1,nblyr+1
2331 0 : D_diag(k) = D_diag(k) - D_spdiag(k) - D_sbdiag(k)
2332 : enddo
2333 :
2334 : ! compute Q_top and Q_bot: top and bottom sources
2335 :
2336 0 : vel2 = -(dhtop/hbri_old/dt +darcyV/bphin_N(1)/hbri_old)
2337 :
2338 0 : Q_top(:) = c0
2339 0 : Q_top(1) = max(c0,vel2*C_top + atm_add/dt)
2340 0 : Qtop = Q_top(1)
2341 :
2342 0 : vel = (dhbot/hbri_old/dt +darcyV/hbri_old) ! going from iphin_N(nblyr+1) to c1 makes a difference
2343 :
2344 0 : Q_bot(:) = c0
2345 0 : Q_bot(nblyr+1) = max(c0,vel*C_bot) + iDin_phi(nblyr+1)*C_bot&
2346 0 : /(zspace + grid_o/hbri_old)
2347 :
2348 0 : Qbot = Q_bot(nblyr+1)
2349 :
2350 0 : Sink_bot = K_diag(nblyr+1) + K_spdiag(nblyr)
2351 0 : Sink_top = K_diag(1) + K_sbdiag(2)
2352 :
2353 : !compute matrix elements (1 to nblyr+1)
2354 :
2355 0 : spdiag = -dt * (D_spdiag + K_spdiag + S_spdiag)
2356 0 : sbdiag = -dt * (D_sbdiag + K_sbdiag + S_sbdiag)
2357 0 : diag = ML - dt * (D_diag + K_diag + S_diag)
2358 0 : rhs = ML * C_in + dt * Q_top + dt* Q_bot
2359 :
2360 0 : end subroutine compute_FCT_matrix
2361 :
2362 : !=======================================================================
2363 : !
2364 : ! Compute matrices for final solution FCT for passive tracers
2365 : ! Corrector step
2366 : !
2367 : ! July, 2014 by N. Jeffery
2368 : !
2369 0 : subroutine compute_FCT_corr (C_in, C_low, dt, nblyr, &
2370 0 : D_sbdiag, D_spdiag, ML)
2371 :
2372 : integer (kind=int_kind), intent(in) :: &
2373 : nblyr ! number of bio layers
2374 :
2375 : real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
2376 : C_in ! Initial (bulk) concentration*hbri_old (mmol/m^2)
2377 : ! conserved quantity on igrid
2378 :
2379 : real (kind=dbl_kind), dimension(nblyr+1), intent(inout) :: &
2380 : C_low ! Low order solution (mmol/m^2) corrected
2381 :
2382 : real (kind=dbl_kind), intent(in) :: &
2383 : dt ! time step
2384 :
2385 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
2386 : D_sbdiag , & ! sub-diagonal artificial diffusion matrix elements ! LCOV_EXCL_LINE
2387 : ML , & ! Lumped mass diagonal matrix elements ! LCOV_EXCL_LINE
2388 : D_spdiag ! super-diagonal artificial diffusion matrix elements
2389 :
2390 : ! local variables
2391 :
2392 : real (kind=dbl_kind) :: &
2393 0 : zspace
2394 :
2395 : integer (kind=int_kind) :: &
2396 : k ! vertical index
2397 :
2398 : real (kind=dbl_kind), dimension (nblyr+1) :: &
2399 : M_spdiag, M_sbdiag, & ! mass matrix ! LCOV_EXCL_LINE
2400 : F_diag, F_spdiag, F_sbdiag, & ! anti-diffusive matrix ! LCOV_EXCL_LINE
2401 : Pplus, Pminus, & ! ! LCOV_EXCL_LINE
2402 : Qplus, Qminus, & ! ! LCOV_EXCL_LINE
2403 : Rplus, Rminus, & ! ! LCOV_EXCL_LINE
2404 0 : a_spdiag, a_sbdiag ! weightings of F
2405 :
2406 : character(len=*),parameter :: subname='(compute_FCT_corr)'
2407 :
2408 : !---------------------------------------------------------------------
2409 : ! Diag (jj) solve for j = 1:nblyr+1
2410 : ! spdiag(j) == (j,j+1) solve for j = 1:nblyr otherwise 0
2411 : ! sbdiag(j) == (j,j-1) solve for j = 2:nblyr+1 otherwise 0
2412 : !---------------------------------------------------------------------
2413 :
2414 0 : zspace = c1/real(nblyr,kind=dbl_kind)
2415 :
2416 : ! compute the mass matrix
2417 :
2418 0 : M_spdiag(:) = zspace/c6
2419 0 : M_spdiag(nblyr+1) = c0
2420 0 : M_sbdiag(:) = zspace/c6
2421 0 : M_sbdiag(1) = c0
2422 :
2423 : ! compute off matrix F:
2424 :
2425 0 : F_diag(:) = c0
2426 0 : F_spdiag(:) = c0
2427 0 : F_sbdiag(:) = c0
2428 :
2429 0 : do k = 1, nblyr
2430 0 : F_spdiag(k) = M_spdiag(k)*(C_low(k)-C_in(k) - C_low(k+1)+ C_in(k+1))/dt &
2431 0 : + D_spdiag(k)*(C_low(k)-C_low(k+1))
2432 0 : F_sbdiag(k+1) = M_sbdiag(k+1)*(C_low(k+1)-C_in(k+1) - C_low(k)+ C_in(k))/dt &
2433 0 : + D_sbdiag(k+1)*(C_low(k+1)-C_low(k))
2434 :
2435 0 : if (F_spdiag(k)*(C_low(k) - C_low(k+1)) > c0) F_spdiag(k) = c0
2436 0 : if (F_sbdiag(k+1)*(C_low(k+1) - C_low(k)) > c0) F_sbdiag(k+1) = c0
2437 : enddo
2438 :
2439 0 : if (maxval(abs(F_spdiag)) > c0) then
2440 :
2441 : ! compute the weighting factors: a_spdiag, a_sbdiag
2442 :
2443 0 : a_spdiag(:) = c0
2444 0 : a_sbdiag(:) = c0
2445 :
2446 0 : Pplus(1) = max(c0, F_spdiag(1))
2447 0 : Pminus(1) = min(c0, F_spdiag(1))
2448 0 : Pplus(nblyr+1) = max(c0, F_sbdiag(nblyr+1))
2449 0 : Pminus(nblyr+1) = min(c0, F_sbdiag(nblyr+1))
2450 0 : Qplus(1) = max(c0,C_low(2)-C_low(1))
2451 0 : Qminus(1)= min(c0,C_low(2)-C_low(1))
2452 0 : Qplus(nblyr+1) = max(c0,C_low(nblyr)-C_low(nblyr+1))
2453 0 : Qminus(nblyr+1)= min(c0,C_low(nblyr)-C_low(nblyr+1))
2454 0 : Rplus(1) = min(c1, ML(1)*Qplus(1)/dt/(Pplus(1)+puny))
2455 0 : Rminus(1) = min(c1, ML(1)*Qminus(1)/dt/(Pminus(1)-puny))
2456 0 : Rplus(nblyr+1) = min(c1, ML(nblyr+1)*Qplus(nblyr+1)/dt/(Pplus(nblyr+1)+puny))
2457 0 : Rminus(nblyr+1) = min(c1, ML(nblyr+1)*Qminus(nblyr+1)/dt/(Pminus(nblyr+1)-puny))
2458 0 : do k = 2,nblyr
2459 0 : Pplus(k) = max(c0,F_spdiag(k)) + max(c0,F_sbdiag(k))
2460 0 : Pminus(k) = min(c0,F_spdiag(k)) + min(c0,F_sbdiag(k))
2461 0 : Qplus(k) = max(c0, max(C_low(k+1)-C_low(k),C_low(k-1)-C_low(k)))
2462 0 : Qminus(k) = min(c0, min(C_low(k+1)-C_low(k),C_low(k-1)-C_low(k)))
2463 0 : Rplus(k) = min(c1, ML(k)*Qplus(k)/dt/(Pplus(k)+puny))
2464 0 : Rminus(k) = min(c1, ML(k)*Qminus(k)/dt/(Pminus(k)-puny))
2465 : enddo
2466 :
2467 0 : do k = 1, nblyr
2468 0 : a_spdiag(k) = min(Rminus(k),Rplus(k+1))
2469 0 : if (F_spdiag(k) > c0) a_spdiag(k) = min(Rplus(k),Rminus(k+1))
2470 0 : a_sbdiag(k+1) = min(Rminus(k+1),Rplus(k))
2471 0 : if (F_sbdiag(k+1) > c0) a_sbdiag(k+1) = min(Rplus(k+1),Rminus(k))
2472 : enddo
2473 :
2474 : !compute F_diag:
2475 :
2476 0 : F_diag(1) = a_spdiag(1)*F_spdiag(1)
2477 0 : F_diag(nblyr+1) = a_sbdiag(nblyr+1)* F_sbdiag(nblyr+1)
2478 0 : C_low(1) = C_low(1) + dt*F_diag(1)/ML(1)
2479 0 : C_low(nblyr+1) = C_low(nblyr+1) + dt*F_diag(nblyr+1)/ML(nblyr+1)
2480 :
2481 0 : do k = 2,nblyr
2482 0 : F_diag(k) = a_spdiag(k)*F_spdiag(k) + a_sbdiag(k)*F_sbdiag(k)
2483 0 : C_low(k) = C_low(k) + dt*F_diag(k)/ML(k)
2484 : enddo
2485 :
2486 : endif !F_spdiag is nonzero
2487 :
2488 0 : end subroutine compute_FCT_corr
2489 :
2490 : !=======================================================================
2491 : !
2492 : ! Tridiagonal matrix solver-- for salinity
2493 : !
2494 : ! authors William H. Lipscomb, LANL
2495 : ! C. M. Bitz, UW
2496 : !
2497 0 : subroutine tridiag_solverz (nmat, sbdiag, &
2498 : diag, spdiag, & ! LCOV_EXCL_LINE
2499 0 : rhs, xout)
2500 :
2501 : integer (kind=int_kind), intent(in) :: &
2502 : nmat ! matrix dimension
2503 :
2504 : real (kind=dbl_kind), dimension (nmat), intent(in) :: &
2505 : sbdiag , & ! sub-diagonal matrix elements ! LCOV_EXCL_LINE
2506 : diag , & ! diagonal matrix elements ! LCOV_EXCL_LINE
2507 : spdiag , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
2508 : rhs ! rhs of tri-diagonal matrix eqn.
2509 :
2510 : real (kind=dbl_kind), dimension (nmat), intent(inout) :: &
2511 : xout ! solution vector
2512 :
2513 : ! local variables
2514 :
2515 : integer (kind=int_kind) :: &
2516 : k ! row counter
2517 :
2518 : real (kind=dbl_kind) :: &
2519 0 : wbeta ! temporary matrix variable
2520 :
2521 : real (kind=dbl_kind), dimension(nmat):: &
2522 0 : wgamma ! temporary matrix variable
2523 :
2524 : character(len=*),parameter :: subname='(tridiag_solverz)'
2525 :
2526 0 : wbeta = diag(1)
2527 0 : xout(1) = rhs(1) / wbeta
2528 :
2529 0 : do k = 2, nmat
2530 0 : wgamma(k) = spdiag(k-1) / wbeta
2531 0 : wbeta = diag(k) - sbdiag(k)*wgamma(k)
2532 0 : xout(k) = (rhs(k) - sbdiag(k)*xout(k-1)) / wbeta
2533 : enddo ! k
2534 :
2535 0 : do k = nmat-1, 1, -1
2536 0 : xout(k) = xout(k) - wgamma(k+1)*xout(k+1)
2537 : enddo ! k
2538 :
2539 0 : end subroutine tridiag_solverz
2540 :
2541 : !=======================================================================
2542 : !
2543 : ! authors Nicole Jeffery, LANL
2544 :
2545 0 : subroutine check_conservation_FCT (C_init, C_new, C_low, S_top, &
2546 : S_bot, L_bot, L_top, dt, & ! LCOV_EXCL_LINE
2547 : fluxbio, nblyr, & ! LCOV_EXCL_LINE
2548 : source)
2549 :
2550 : integer (kind=int_kind), intent(in) :: &
2551 : nblyr ! number of bio layers
2552 :
2553 : real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
2554 : C_init , & ! initial bulk concentration * h_old (mmol/m^2) ! LCOV_EXCL_LINE
2555 : C_new ! new bulk concentration * h_new (mmol/m^2)
2556 :
2557 : real (kind=dbl_kind), dimension(nblyr+1), intent(out) :: &
2558 : C_low ! define low order solution = C_new
2559 :
2560 : real (kind=dbl_kind), intent(in) :: &
2561 : S_top , & ! surface flux into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
2562 : S_bot , & ! bottom flux into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
2563 : L_bot , & ! remaining bottom flux into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
2564 : L_top , & ! remaining top flux into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
2565 : dt , & ! LCOV_EXCL_LINE
2566 : source ! nutrient source from snow and atmosphere (mmol/m^2)
2567 :
2568 : real (kind=dbl_kind), intent(inout) :: &
2569 : fluxbio ! (mmol/m^2/s) positive down (into the ocean)
2570 :
2571 : ! local variables
2572 :
2573 : integer (kind=int_kind) :: &
2574 : k
2575 :
2576 : real (kind=dbl_kind) :: &
2577 : diff_dt , & ! LCOV_EXCL_LINE
2578 : C_init_tot , & ! LCOV_EXCL_LINE
2579 : C_new_tot , & ! LCOV_EXCL_LINE
2580 : zspace , & !1/nblyr ! LCOV_EXCL_LINE
2581 0 : accuracy ! centered difference is Order(zspace^2)
2582 :
2583 : character(len=*),parameter :: subname='(check_conservation_FCT)'
2584 :
2585 0 : zspace = c1/real(nblyr,kind=dbl_kind)
2586 :
2587 : !-------------------------------------
2588 : ! Ocean flux: positive into the ocean
2589 : !-------------------------------------
2590 0 : C_init_tot = (C_init(1) + C_init(nblyr+1))*zspace*p5
2591 0 : C_new_tot = (C_new(1) + C_new(nblyr+1))*zspace*p5
2592 0 : C_low(1) = C_new(1)
2593 0 : C_low(nblyr+1) = C_new(nblyr+1)
2594 :
2595 0 : do k = 2, nblyr
2596 0 : C_init_tot = C_init_tot + C_init(k)*zspace
2597 0 : C_new_tot = C_new_tot + C_new(k)*zspace
2598 0 : C_low(k) = C_new(k)
2599 : enddo
2600 :
2601 0 : accuracy = 1.0e-14_dbl_kind*max(c1, C_init_tot, C_new_tot)
2602 0 : fluxbio = (C_init_tot - C_new_tot + source)/dt
2603 0 : diff_dt =C_new_tot - C_init_tot - (S_top+S_bot+L_bot*C_new(nblyr+1)+L_top*C_new(1))*dt
2604 :
2605 0 : if (minval(C_low) < c0) then
2606 0 : write(warnstr,*) subname, 'Positivity of zbgc low order solution failed: C_low:',C_low
2607 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
2608 : endif
2609 :
2610 0 : if (abs(diff_dt) > accuracy ) then
2611 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
2612 0 : write(warnstr,*) subname, 'Conservation of zbgc low order solution failed: diff_dt:',&
2613 0 : diff_dt
2614 0 : write(warnstr,*) subname, 'Total initial tracer', C_init_tot
2615 0 : write(warnstr,*) subname, 'Total final1 tracer', C_new_tot
2616 0 : write(warnstr,*) subname, 'bottom final tracer', C_new(nblyr+1)
2617 0 : write(warnstr,*) subname, 'top final tracer', C_new(1)
2618 0 : write(warnstr,*) subname, 'Near bottom final tracer', C_new(nblyr)
2619 0 : write(warnstr,*) subname, 'Near top final tracer', C_new(2)
2620 0 : write(warnstr,*) subname, 'Top flux*dt into ice:', S_top*dt
2621 0 : write(warnstr,*) subname, 'Bottom flux*dt into ice:', S_bot*dt
2622 0 : write(warnstr,*) subname, 'Remaining bot flux*dt into ice:', L_bot*C_new(nblyr+1)*dt
2623 0 : write(warnstr,*) subname, 'S_bot*dt + L_bot*C_new(nblyr+1)*dt'
2624 0 : write(warnstr,*) subname, S_bot*dt + L_bot*C_new(nblyr+1)*dt
2625 0 : write(warnstr,*) subname, 'fluxbio*dt:', fluxbio*dt
2626 0 : write(warnstr,*) subname, 'fluxbio:', fluxbio
2627 0 : write(warnstr,*) subname, 'Remaining top flux*dt into ice:', L_top*C_new(1)*dt
2628 : endif
2629 :
2630 0 : end subroutine check_conservation_FCT
2631 :
2632 : !=======================================================================
2633 :
2634 : ! For each grid cell, sum field over all ice and snow layers
2635 : !
2636 : ! author: Nicole Jeffery, LANL
2637 :
2638 0 : subroutine bgc_column_sum (nblyr, nslyr, hsnow, hbrine, xin, xout)
2639 :
2640 : integer (kind=int_kind), intent(in) :: &
2641 : nblyr, & ! number of ice layers ! LCOV_EXCL_LINE
2642 : nslyr ! number of snow layers
2643 :
2644 : real (kind=dbl_kind), dimension(nblyr+3), intent(in) :: &
2645 : xin ! input field
2646 :
2647 : real (kind=dbl_kind), intent(in) :: &
2648 : hsnow, & ! snow thickness ! LCOV_EXCL_LINE
2649 : hbrine ! brine height
2650 :
2651 : real (kind=dbl_kind), intent(out) :: &
2652 : xout ! output field
2653 :
2654 : ! local variables
2655 :
2656 : real (kind=dbl_kind) :: &
2657 : dzssl, & ! snow surface layer (m) ! LCOV_EXCL_LINE
2658 : dzint, & ! snow interior depth (m) ! LCOV_EXCL_LINE
2659 : hslyr, & ! snow layer thickness (m) ! LCOV_EXCL_LINE
2660 0 : zspace ! brine layer thickness/hbrine
2661 :
2662 : integer (kind=int_kind) :: &
2663 : n ! category/layer index
2664 :
2665 : character(len=*),parameter :: subname='(bgc_column_sum)'
2666 :
2667 0 : hslyr = hsnow/real(nslyr,kind=dbl_kind)
2668 0 : dzssl = min(hslyr*p5, hs_ssl)
2669 0 : dzint = max(c0,hsnow - dzssl)
2670 0 : zspace = c1/real(nblyr,kind=dbl_kind)
2671 :
2672 0 : xout = c0
2673 0 : xout = (xin(1) + xin(nblyr+1))*hbrine*p5*zspace
2674 0 : do n = 2, nblyr
2675 0 : xout = xout + xin(n)*zspace*hbrine
2676 : enddo ! n
2677 0 : xout = xout + dzssl*xin(nblyr+2) + dzint*xin(nblyr+3)
2678 :
2679 0 : end subroutine bgc_column_sum
2680 :
2681 : !=======================================================================
2682 :
2683 : end module icepack_algae
2684 :
2685 : !=======================================================================
|
