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