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