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 ! LCOV_EXCL_LINE
75 : accuracy = 1.0e-13_dbl_kind ! accuracy parameter for bgc tracers
76 : !=======================================================================
77 :
78 : contains
79 :
80 : !=======================================================================
81 :
82 33835340 : subroutine zbio (dt, &
83 : meltt, melts, & ! LCOV_EXCL_LINE
84 : meltb, congel, & ! LCOV_EXCL_LINE
85 : snoice, fsnow, & ! LCOV_EXCL_LINE
86 33835340 : trcrn, bio_index, & ! LCOV_EXCL_LINE
87 33835340 : bio_index_o, aice_old, & ! LCOV_EXCL_LINE
88 : vice_old, vsno_old, & ! LCOV_EXCL_LINE
89 : vicen, vsnon, & ! LCOV_EXCL_LINE
90 33835340 : aicen, flux_bio_atm,& ! LCOV_EXCL_LINE
91 : n_cat, first_ice, & ! LCOV_EXCL_LINE
92 33835340 : hice_old, ocean_bio, & ! LCOV_EXCL_LINE
93 33835340 : ocean_bio_dh, & ! LCOV_EXCL_LINE
94 33835340 : bphin, iphin, & ! LCOV_EXCL_LINE
95 33835340 : iDin, & ! LCOV_EXCL_LINE
96 33835340 : fswthrul, & ! LCOV_EXCL_LINE
97 : dh_top, dh_bot, & ! LCOV_EXCL_LINE
98 33835340 : zfswin, & ! LCOV_EXCL_LINE
99 : hbri, hbri_old, & ! LCOV_EXCL_LINE
100 : darcy_V, & ! LCOV_EXCL_LINE
101 : bphi_min, & ! LCOV_EXCL_LINE
102 33835340 : iTin, & ! LCOV_EXCL_LINE
103 33835340 : Zoo, & ! LCOV_EXCL_LINE
104 33835340 : flux_bio, dh_direct, & ! LCOV_EXCL_LINE
105 : upNO, upNH, & ! LCOV_EXCL_LINE
106 33835340 : fbio_snoice, fbio_atmice, & ! LCOV_EXCL_LINE
107 33835340 : PP_net, ice_bio_net, & ! LCOV_EXCL_LINE
108 33835340 : snow_bio_net, grow_net, & ! LCOV_EXCL_LINE
109 : totalChla, & ! LCOV_EXCL_LINE
110 33835340 : flux_bion, iSin, & ! LCOV_EXCL_LINE
111 33835340 : bioPorosityIceCell, & ! LCOV_EXCL_LINE
112 33835340 : bioSalinityIceCell, & ! LCOV_EXCL_LINE
113 33835340 : 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) ! LCOV_EXCL_LINE
120 : bio_index_o ! references index of data arrays (eg. kscavz)
121 :
122 : real (kind=dbl_kind), intent(in) :: &
123 : dt, & ! time step ! LCOV_EXCL_LINE
124 : hbri, & ! brine height (m) ! LCOV_EXCL_LINE
125 : bphi_min, & ! surface porosity ! LCOV_EXCL_LINE
126 : meltt, & ! thermodynamic melt/growth rates in dt (m) ! LCOV_EXCL_LINE
127 : melts, & ! LCOV_EXCL_LINE
128 : meltb, & ! LCOV_EXCL_LINE
129 : congel, & ! LCOV_EXCL_LINE
130 : snoice, & ! LCOV_EXCL_LINE
131 : fsnow, & ! snowfall rate (kg/m^2 s) ! LCOV_EXCL_LINE
132 : hice_old, & ! ice height (m) ! LCOV_EXCL_LINE
133 : vicen, & ! ice volume (m) ! LCOV_EXCL_LINE
134 : vsnon, & ! snow volume (m) ! LCOV_EXCL_LINE
135 : aicen, & ! ice area fraction ! LCOV_EXCL_LINE
136 : aice_old, & ! values prior to thermodynamic changes ! LCOV_EXCL_LINE
137 : vice_old, & ! LCOV_EXCL_LINE
138 : vsno_old, & ! LCOV_EXCL_LINE
139 : darcy_V, & ! darcy velocity ! LCOV_EXCL_LINE
140 : ! darcy_V_chl,& ! darcy velocity for algae ! LCOV_EXCL_LINE
141 : dh_bot, & ! change in brine bottom (m) ! LCOV_EXCL_LINE
142 : dh_top, & ! change in brine top (m) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
147 : ice_bio_net, & ! net bio tracer in ice (mmol/m^2) ! LCOV_EXCL_LINE
148 : fbio_atmice, & ! bio flux from atm to ice (mmol/m^2/s) ! LCOV_EXCL_LINE
149 : fbio_snoice, & ! bio flux from snow to ice (mmol/m^2/s) ! LCOV_EXCL_LINE
150 : flux_bio, & ! total ocean tracer flux (mmol/m^2/s) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
158 : iphin , & ! Porosity on the igrid ! LCOV_EXCL_LINE
159 : iDin , & ! Diffusivity/h on the igrid (1/s) ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
191 : grow_net , & ! net specific growth (m/d) times vice ! LCOV_EXCL_LINE
192 : upNO , & ! tot nitrate uptake rate (mmol/m^2/d) times aice ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
200 : bioSalinityIceCell , & ! salinity on vertical interface points (ppt) ! LCOV_EXCL_LINE
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 67670680 : upNOn , & ! algal nitrate uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
213 67670680 : upNHn , & ! algal ammonium uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
214 67670680 : grow_alg ! algal growth rate (mmol/m^3/s)
215 :
216 : real (kind=dbl_kind),dimension(nbtrcr) :: &
217 67670680 : zbgc_snown, & ! aerosol contribution from snow to ice ! LCOV_EXCL_LINE
218 67670680 : zbgc_atmn ! and atm to ice concentration * volume (mmol/m^3*m)
219 :
220 : real (kind=dbl_kind), dimension(nbtrcr) :: &
221 67670680 : Tot_BGC_i, & ! initial column sum, ice + snow, of BGC tracer (mmol/m^2) ! LCOV_EXCL_LINE
222 67670680 : Tot_BGC_f, & ! final column sum ! LCOV_EXCL_LINE
223 33835340 : flux_bio_sno !
224 :
225 : real (kind=dbl_kind) :: &
226 : hsnow_i, & ! initial snow thickness (m) ! LCOV_EXCL_LINE
227 : hsnow_f, & ! final snow thickness (m) ! LCOV_EXCL_LINE
228 : carbonError ! carbon conservation error (mmol/m2)
229 :
230 : real (kind=dbl_kind) :: &
231 : carbonInitial, & ! initial carbon content (mmol/m2) ! LCOV_EXCL_LINE
232 : carbonFinal, & ! final carbon content (mmol/m2) ! LCOV_EXCL_LINE
233 : carbonFlux ! carbon flux (mmol/m2/s)
234 :
235 : logical (kind=log_kind) :: &
236 : write_flux_diag, & ! LCOV_EXCL_LINE
237 : write_carbon_errors
238 :
239 : character(len=*),parameter :: subname='(zbio)'
240 :
241 595877548 : zbgc_snown(:) = c0
242 595877548 : zbgc_atmn (:) = c0
243 595877548 : flux_bion (:) = c0
244 595877548 : flux_bio_sno(:) = c0
245 595877548 : Tot_BGC_i (:) = c0
246 595877548 : Tot_BGC_f (:) = c0
247 304518060 : Zoo (:) = c0
248 33835340 : hsnow_i = c0
249 33835340 : hsnow_f = c0
250 33835340 : write_flux_diag = .false.
251 33835340 : write_carbon_errors = .true.
252 33835340 : if (.not. tr_bgc_C) write_carbon_errors = .false.
253 :
254 33835340 : call bgc_carbon_sum(hbri_old, trcrn(:), carbonInitial)
255 33835340 : if (icepack_warnings_aborted(subname)) return
256 :
257 33835340 : if (aice_old > puny) then
258 33835340 : hsnow_i = vsno_old/aice_old
259 595877548 : do mm = 1,nbtrcr
260 : call bgc_column_sum (hsnow_i, hbri_old, &
261 : trcrn(bio_index(mm):bio_index(mm)+nblyr+2), & ! LCOV_EXCL_LINE
262 562042208 : Tot_BGC_i(mm))
263 595877548 : if (icepack_warnings_aborted(subname)) return
264 : enddo
265 : endif
266 :
267 : call update_snow_bgc (dt, &
268 : meltt, melts, & ! LCOV_EXCL_LINE
269 : meltb, congel, & ! LCOV_EXCL_LINE
270 : snoice, fsnow, & ! LCOV_EXCL_LINE
271 : trcrn, bio_index, & ! LCOV_EXCL_LINE
272 : aice_old, zbgc_snown, & ! LCOV_EXCL_LINE
273 : vice_old, vsno_old, & ! LCOV_EXCL_LINE
274 : vicen, vsnon, & ! LCOV_EXCL_LINE
275 : aicen, flux_bio_atm, & ! LCOV_EXCL_LINE
276 : zbgc_atmn, flux_bio_sno, & ! LCOV_EXCL_LINE
277 33835340 : bio_index_o)
278 :
279 33835340 : if (icepack_warnings_aborted(subname)) return
280 :
281 : call z_biogeochemistry (n_cat, dt, &
282 : first_ice, & ! LCOV_EXCL_LINE
283 : aicen, vicen, & ! LCOV_EXCL_LINE
284 : hice_old, ocean_bio, & ! LCOV_EXCL_LINE
285 : ocean_bio_dh, & ! LCOV_EXCL_LINE
286 : flux_bion, bphin, & ! LCOV_EXCL_LINE
287 : iphin, trcrn, & ! LCOV_EXCL_LINE
288 : iDin, & ! LCOV_EXCL_LINE
289 : fswthrul, grow_alg, & ! LCOV_EXCL_LINE
290 : upNOn, upNHn, & ! LCOV_EXCL_LINE
291 : dh_top, dh_bot, & ! LCOV_EXCL_LINE
292 : zfswin, hbri, & ! LCOV_EXCL_LINE
293 : hbri_old, darcy_V, & ! LCOV_EXCL_LINE
294 : bphi_min, zbgc_snown,& ! LCOV_EXCL_LINE
295 : zbgc_atmn, & ! LCOV_EXCL_LINE
296 : iTin, dh_direct, & ! LCOV_EXCL_LINE
297 : Zoo, meltb, & ! LCOV_EXCL_LINE
298 33835340 : congel )
299 33835340 : if (icepack_warnings_aborted(subname)) return
300 :
301 595877548 : do mm = 1,nbtrcr
302 595877548 : flux_bion(mm) = flux_bion(mm) + flux_bio_sno(mm)
303 : enddo
304 :
305 33835340 : call bgc_carbon_sum(hbri, trcrn(:), carbonFinal)
306 33835340 : if (icepack_warnings_aborted(subname)) return
307 33835340 : call bgc_carbon_flux(flux_bio_atm,flux_bion,carbonFlux)
308 33835340 : if (icepack_warnings_aborted(subname)) return
309 :
310 33835340 : carbonError = carbonInitial-carbonFlux*dt-carbonFinal
311 :
312 101506020 : 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), & ! LCOV_EXCL_LINE
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 33835340 : if (icepack_warnings_aborted(subname)) return
364 :
365 33835340 : 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), & ! LCOV_EXCL_LINE
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 33835340 : if (icepack_warnings_aborted(subname)) return
392 :
393 : call merge_bgc_fluxes (dt, &
394 : bio_index, & ! LCOV_EXCL_LINE
395 : aicen, & ! LCOV_EXCL_LINE
396 : vicen, vsnon, & ! LCOV_EXCL_LINE
397 : iphin, &! ntrcr ! LCOV_EXCL_LINE
398 : trcrn, aice_old, &!aice_old ! LCOV_EXCL_LINE
399 : flux_bion, flux_bio, & ! LCOV_EXCL_LINE
400 : upNOn, upNHn, & ! LCOV_EXCL_LINE
401 : upNO, upNH, & ! LCOV_EXCL_LINE
402 : zbgc_snown, zbgc_atmn, & ! LCOV_EXCL_LINE
403 : fbio_snoice, fbio_atmice,& ! LCOV_EXCL_LINE
404 : PP_net, ice_bio_net,& ! LCOV_EXCL_LINE
405 : snow_bio_net, grow_alg, & ! LCOV_EXCL_LINE
406 : grow_net, totalChla, & ! LCOV_EXCL_LINE
407 : iTin, iSin, & ! LCOV_EXCL_LINE
408 : bioPorosityIceCell, & ! LCOV_EXCL_LINE
409 : bioSalinityIceCell, & ! LCOV_EXCL_LINE
410 33835340 : bioTemperatureIceCell)
411 33835340 : if (icepack_warnings_aborted(subname)) return
412 :
413 33835340 : 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 33835340 : 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, & ! LCOV_EXCL_LINE
439 : aicen, & ! LCOV_EXCL_LINE
440 : meltb, congel, & ! LCOV_EXCL_LINE
441 : fswthru, first_ice, & ! LCOV_EXCL_LINE
442 0 : trcrn, & ! LCOV_EXCL_LINE
443 : PP_net, upNO, & ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
451 : Tf, & ! basal freezing temperature (C) ! LCOV_EXCL_LINE
452 : ! hmix, & ! mixed layer depth (m) ! LCOV_EXCL_LINE
453 : aicen, & ! ice area fraction ! LCOV_EXCL_LINE
454 : meltb, & ! bottom melt (m) ! LCOV_EXCL_LINE
455 : congel, & ! bottom growth (m) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
470 : grow_net, & ! net specific growth (/s) ! LCOV_EXCL_LINE
471 : upNO , & ! tot nitrate uptake rate (mmol/m^2/s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
478 0 : upNHn , & ! algal ammonium uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
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, & ! LCOV_EXCL_LINE
493 : ! hmix, aicen, & ! LCOV_EXCL_LINE
494 : meltb, congel, & ! LCOV_EXCL_LINE
495 : fswthru, first_ice, & ! LCOV_EXCL_LINE
496 : trcrn, upNOn, & ! LCOV_EXCL_LINE
497 : upNHn, grow_alg, & ! LCOV_EXCL_LINE
498 0 : Tf)
499 0 : if (icepack_warnings_aborted(subname)) return
500 :
501 : call merge_bgc_fluxes_skl ( &
502 : aicen, trcrn, & ! LCOV_EXCL_LINE
503 : flux_bion, flux_bio, & ! LCOV_EXCL_LINE
504 : PP_net, upNOn, & ! LCOV_EXCL_LINE
505 : upNHn, upNO, & ! LCOV_EXCL_LINE
506 : upNH, grow_net, & ! LCOV_EXCL_LINE
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, & ! LCOV_EXCL_LINE
518 : ! hmix, aicen, & ! LCOV_EXCL_LINE
519 : meltb, congel, & ! LCOV_EXCL_LINE
520 : fswthru, first_ice, & ! LCOV_EXCL_LINE
521 0 : trcrn, upNOn, & ! LCOV_EXCL_LINE
522 0 : upNHn, grow_alg_skl, & ! LCOV_EXCL_LINE
523 : Tf)
524 :
525 : real (kind=dbl_kind), intent(in) :: &
526 : dt , & ! time step ! LCOV_EXCL_LINE
527 : ! hmix , & ! mixed layer depth ! LCOV_EXCL_LINE
528 : ! aicen , & ! ice area ! LCOV_EXCL_LINE
529 : meltb , & ! bottom ice melt ! LCOV_EXCL_LINE
530 : congel , & ! bottom ice growth ! LCOV_EXCL_LINE
531 : Tf , & ! bottom freezing temperature ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
550 : upNOn , & ! algal NO uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
559 0 : cinit , & ! initial brine concentration*sk_l (mmol/m^2) ! LCOV_EXCL_LINE
560 0 : cinit_v , & ! initial brine concentration (mmol/m^3) ! LCOV_EXCL_LINE
561 0 : congel_alg , & ! congelation flux contribution to ice algae (mmol/m^2 s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
565 0 : PVflag , & ! 1 for tracers that flow with the brine, 0 otherwise ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
570 : iTin , & ! LCOV_EXCL_LINE
571 : PVt , & ! type 'Jin2006' piston velocity (m/s) ! LCOV_EXCL_LINE
572 : ice_growth , & ! Jin2006 definition: either congel rate or bottom melt rate (m/s) ! LCOV_EXCL_LINE
573 : grow_val , & ! (m/x) ! LCOV_EXCL_LINE
574 : rphi_sk , & ! 1 / skeletal layer porosity ! LCOV_EXCL_LINE
575 : cinit_tmp , & ! temporary variable for concentration (mmol/m^2) ! LCOV_EXCL_LINE
576 : Cerror , & ! change in total carbon from reactions (mmol/m^3) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
581 : PV_scale_growth = p5 , & ! scale factor in Jin code PV during ice growth ! LCOV_EXCL_LINE
582 : PV_scale_melt = p05 , & ! scale factor in Jin code PV during ice melt ! LCOV_EXCL_LINE
583 : growth_max = 1.85e-5_dbl_kind , & ! PVt function reaches maximum here. (m/s) ! LCOV_EXCL_LINE
584 : MJ1 = 9.667e-9_dbl_kind , & ! (m/s) coefficients in Jin2008 ! LCOV_EXCL_LINE
585 : MJ2 = 38.8_dbl_kind , & ! (1) from:4.49e-4_dbl_kind*secday ! LCOV_EXCL_LINE
586 : MJ3 = 1.04e7_dbl_kind , & ! 1/(m/s) from: 1.39e-3_dbl_kind*secday^2 ! LCOV_EXCL_LINE
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)), & ! LCOV_EXCL_LINE
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)), & ! LCOV_EXCL_LINE
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, & ! LCOV_EXCL_LINE
702 : fswthru, react, & ! LCOV_EXCL_LINE
703 : cinit_v, & ! LCOV_EXCL_LINE
704 : grow_alg_skl(:), & ! LCOV_EXCL_LINE
705 : iTin, & ! LCOV_EXCL_LINE
706 : upNOn(:), upNHn(:), & ! LCOV_EXCL_LINE
707 : Zoo_skl, & ! LCOV_EXCL_LINE
708 : Cerror, conserve_C,& ! LCOV_EXCL_LINE
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)) & ! LCOV_EXCL_LINE
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 33835340 : subroutine z_biogeochemistry (n_cat, dt, &
794 : first_ice, & ! LCOV_EXCL_LINE
795 : aicen, vicen, & ! LCOV_EXCL_LINE
796 0 : hice_old, ocean_bio, & ! LCOV_EXCL_LINE
797 33835340 : ocean_bio_dh, & ! LCOV_EXCL_LINE
798 67670680 : flux_bio, bphin, & ! LCOV_EXCL_LINE
799 33835340 : iphin, trcrn, & ! LCOV_EXCL_LINE
800 33835340 : iDin, & ! LCOV_EXCL_LINE
801 67670680 : fswthrul, grow_alg, & ! LCOV_EXCL_LINE
802 33835340 : upNOn, upNHn, & ! LCOV_EXCL_LINE
803 : dh_top, dh_bot, & ! LCOV_EXCL_LINE
804 33835340 : zfswin, hbri, & ! LCOV_EXCL_LINE
805 : hbri_old, darcy_V, & ! LCOV_EXCL_LINE
806 33835340 : bphi_min, zbgc_snow, & ! LCOV_EXCL_LINE
807 33835340 : zbgc_atm, & ! LCOV_EXCL_LINE
808 33835340 : iTin, dh_direct, & ! LCOV_EXCL_LINE
809 33835340 : Zoo, meltb, & ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
820 : hbri , & ! brine height (m) ! LCOV_EXCL_LINE
821 : bphi_min , & ! surface porosity ! LCOV_EXCL_LINE
822 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
823 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
824 : hice_old , & ! ice height (m) ! LCOV_EXCL_LINE
825 : meltb , & ! bottom melt in dt (m) ! LCOV_EXCL_LINE
826 : congel , & ! bottom growth in dt (m) ! LCOV_EXCL_LINE
827 : darcy_V , & ! darcy velocity ! LCOV_EXCL_LINE
828 : dh_bot , & ! change in brine bottom (m) ! LCOV_EXCL_LINE
829 : dh_top , & ! change in brine top (m) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
834 : zfswin , & ! visible Short wave flux on igrid (W/m^2) ! LCOV_EXCL_LINE
835 : Zoo , & ! N losses to the system from reaction terms ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
841 : iphin , & ! Porosity on the igrid ! LCOV_EXCL_LINE
842 : iDin , & ! Diffusivity/h on the igrid (1/s) ! LCOV_EXCL_LINE
843 : fswthrul , & ! visible short wave radiation on icgrid (W/m^2) ! LCOV_EXCL_LINE
844 : zbgc_snow , & ! tracer input from snow (mmol/m^3*m) ! LCOV_EXCL_LINE
845 : zbgc_atm , & ! tracer input from atm (mmol/m^3 *m) ! LCOV_EXCL_LINE
846 : ocean_bio , & ! ocean concentrations (mmol/m^3) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
857 : upNHn , & ! algal ammonium uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
874 : hin_old , & ! ice thickness before current melt/growth (m) ! LCOV_EXCL_LINE
875 : ice_conc , & ! algal concentration in ice above hin > hinS ! LCOV_EXCL_LINE
876 : sum_initial , & ! ! LCOV_EXCL_LINE
877 : sum_old , & ! ! LCOV_EXCL_LINE
878 : sum_new , & ! ! LCOV_EXCL_LINE
879 : sum_tot , & ! ! LCOV_EXCL_LINE
880 : zspace , & ! 1/nblyr ! LCOV_EXCL_LINE
881 : darcyV , & ! ! LCOV_EXCL_LINE
882 : dhtop , & ! ! LCOV_EXCL_LINE
883 : dhbot , & ! ! LCOV_EXCL_LINE
884 : dhflood ! >=0 (m) surface flooding from the ocean
885 :
886 : real (kind=dbl_kind), dimension (nblyr+2) :: &
887 67670680 : bphin_N ! porosity for tracer model has minimum
888 : ! bphin_N >= bphimin
889 :
890 : real (kind=dbl_kind), dimension (nblyr+1) :: &
891 67670680 : iphin_N , & ! tracer porosity on the igrid ! LCOV_EXCL_LINE
892 67670680 : sbdiagz , & ! sub-diagonal matrix elements ! LCOV_EXCL_LINE
893 67670680 : diagz , & ! diagonal matrix elements ! LCOV_EXCL_LINE
894 67670680 : spdiagz , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
895 67670680 : rhsz , & ! rhs of tri-diagonal matrix equation ! LCOV_EXCL_LINE
896 67670680 : ML_diag , & ! lumped mass matrix ! LCOV_EXCL_LINE
897 67670680 : D_spdiag , & ! artificial diffusion matrix ! LCOV_EXCL_LINE
898 67670680 : D_sbdiag , & ! artificial diffusion matrix ! LCOV_EXCL_LINE
899 67670680 : biomat_low , & ! Low order solution ! LCOV_EXCL_LINE
900 67670680 : Cerror ! Change in N after reactions
901 :
902 : real (kind=dbl_kind), dimension(nblyr+1,nbtrcr):: &
903 67670680 : react ! biological sources and sinks for equation matrix
904 :
905 : real (kind=dbl_kind), dimension(nblyr+1,nbtrcr):: &
906 67670680 : in_init_cons , & ! Initial bulk concentration*h (mmol/m^2) ! LCOV_EXCL_LINE
907 67670680 : biomat_cons , & ! Matrix output of (mmol/m^2) ! LCOV_EXCL_LINE
908 67670680 : biomat_brine ! brine concentration (mmol/m^3)
909 :
910 : real (kind=dbl_kind), dimension(nbtrcr):: &
911 67670680 : C_top, & ! bulk top tracer source: h phi C(meltwater) (mmol/m^2) ! LCOV_EXCL_LINE
912 67670680 : C_bot, & ! bulk bottom tracer source: h phi C(ocean) (mmol/m^2) ! LCOV_EXCL_LINE
913 67670680 : Source_top, & ! For cons: (+) top tracer source into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
914 67670680 : Source_bot, & ! For cons: (+) bottom tracer source into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
915 67670680 : Sink_bot, & ! For cons: (+ or -) remaining bottom flux into ice(mmol/m^2/s) ! LCOV_EXCL_LINE
916 67670680 : Sink_top, & ! For cons: (+ or -) remaining bottom flux into ice(mmol/m^2/s) ! LCOV_EXCL_LINE
917 67670680 : exp_ret, & ! exp dt/retention frequency ! LCOV_EXCL_LINE
918 67670680 : exp_rel, & ! exp dt/release frequency ! LCOV_EXCL_LINE
919 33835340 : atm_add_cons , & ! zbgc_snow+zbgc_atm (mmol/m^3*m) ! LCOV_EXCL_LINE
920 67670680 : dust_Fe , & ! contribution of dust surface flux to dFe (umol/m*3*m) ! LCOV_EXCL_LINE
921 67670680 : source ! mmol/m^2 surface input from snow/atmosphere
922 :
923 : real (kind=dbl_kind), dimension (ntrcr+2) :: &
924 67670680 : trtmp0 , & ! temporary, remapped tracers ! LCOV_EXCL_LINE
925 67670680 : trtmp ! temporary, remapped tracers
926 :
927 : logical (kind=log_kind), dimension(nblyr+1) :: &
928 67670680 : conserve_C
929 :
930 : real (kind=dbl_kind), dimension(nblyr+1):: & ! temporary variables for
931 67670680 : Diff , & ! diffusivity ! LCOV_EXCL_LINE
932 67670680 : initcons , & ! initial concentration ! LCOV_EXCL_LINE
933 67670680 : biocons , & ! new concentration ! LCOV_EXCL_LINE
934 67670680 : dmobile , & ! ! LCOV_EXCL_LINE
935 67670680 : initcons_mobile,&! ! LCOV_EXCL_LINE
936 67670680 : initcons_stationary, & ! LCOV_EXCL_LINE
937 67670680 : dz , & ! normalized vertical grid spacing ! LCOV_EXCL_LINE
938 67670680 : 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 ! LCOV_EXCL_LINE
946 : exp_min ! temporary exp var
947 :
948 : real (kind=dbl_kind):: &
949 : Sat_conc , & ! adsorbing saturation concentration (mmols/m^3) ! LCOV_EXCL_LINE
950 : phi_max , & ! maximum porosity ! LCOV_EXCL_LINE
951 : S_col , & ! surface area of collector (um^2) ! LCOV_EXCL_LINE
952 : P_b , & ! projected area of diatoms (um^2) ! LCOV_EXCL_LINE
953 : V_c , & ! volume of collector (um^3) ! LCOV_EXCL_LINE
954 : V_alg ! volume of algae (um^3)
955 :
956 : real (kind=dbl_kind), dimension(nbtrcr) :: &
957 67670680 : mobile , & ! c1 if mobile, c0 otherwise ! LCOV_EXCL_LINE
958 33835340 : 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) ! LCOV_EXCL_LINE
964 : r_bac= 4.7_dbl_kind , & ! diatom large radius (um) ! LCOV_EXCL_LINE
965 : r_alg= 10.0_dbl_kind , & ! diatom small radius (um) ! LCOV_EXCL_LINE
966 : Nquota_A = 0.88_dbl_kind, & ! slope in Nitrogen quota to cell volume fit ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
970 : f_a = 0.3_dbl_kind, & !c1 , & ! fraction of collector available for attachment ! LCOV_EXCL_LINE
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 33835340 : write_carbon_errors = .true.
987 33835340 : if (.not. tr_bgc_C) write_carbon_errors = .false.
988 :
989 33835340 : zspace = c1/real(nblyr,kind=dbl_kind)
990 304518060 : dz(:) = zspace
991 33835340 : dz(1) = zspace/c2
992 33835340 : dz(nblyr+1) = zspace/c2
993 5092215212 : in_init_cons(:,:) = c0
994 595877548 : atm_add_cons(:) = c0
995 33835340 : dhtop = c0
996 33835340 : dhbot = c0
997 33835340 : darcyV = c0
998 595877548 : C_top(:) = c0
999 595877548 : C_bot(:) = c0
1000 33835340 : if (present(ocean_bio_dh)) ocean_bio_dh(:) = c0
1001 595877548 : mobile(:) = c0
1002 304518060 : conserve_C(:) = .true.
1003 304518060 : nitrification(:) = c0
1004 595877548 : flux_bio_tmp(:) = c0
1005 :
1006 595877548 : do m = 1, nbtrcr
1007 5092215212 : do k = 1, nblyr+1
1008 :
1009 4496337664 : bphin_N(nblyr+2) =c1
1010 4496337664 : bphin_N(k) = bphin(k)
1011 4496337664 : iphin_N(k) = iphin(k)
1012 4496337664 : bphin_N(1) = bphi_min
1013 :
1014 4496337664 : if (abs(trcrn(bio_index(m) + k-1)) < accuracy) then
1015 283825628 : flux_bio_tmp(m) = trcrn(bio_index(m) + k-1)* hbri_old * dz(k)/dt
1016 283825628 : trcrn(bio_index(m) + k-1) = c0
1017 283825628 : in_init_cons(k,m) = c0
1018 : else
1019 4212512036 : in_init_cons(k,m) = trcrn(bio_index(m) + k-1)* hbri_old
1020 : endif
1021 :
1022 4496337664 : 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 5058379872 : if (icepack_warnings_aborted(subname)) return
1039 : enddo !k
1040 : enddo !m
1041 :
1042 : !-----------------------------------------------------------------
1043 : ! boundary conditions
1044 : !-----------------------------------------------------------------
1045 :
1046 33835340 : ice_conc = c0
1047 33835340 : hin = vicen/aicen
1048 33835340 : hin_old = hice_old
1049 :
1050 : !-----------------------------------------------------------------
1051 : ! calculate the saturation concentration for attachment: Sat_conc
1052 : !-----------------------------------------------------------------
1053 :
1054 304518060 : phi_max = maxval(bphin_N(2:nblyr+1))
1055 33835340 : S_col = 4.0_dbl_kind*pi*r_c**2
1056 33835340 : P_b = pi*r_bac**2 !*10-6 for colloids
1057 33835340 : V_c = 4.0_dbl_kind*pi*r_c**3/3.0_dbl_kind ! (m^3) sphere
1058 33835340 : V_alg = pi/6.0_dbl_kind*r_bac*r_alg**2 ! prolate spheroid (*10-9 for colloids)
1059 33835340 : 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 595877548 : dust_Fe(:) = c0
1066 :
1067 33835340 : 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))) * & ! LCOV_EXCL_LINE
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 595877548 : do m = 1,nbtrcr
1078 : !-----------------------------------------------------------------
1079 : ! time constants for mobile/stationary phase changes
1080 : !-----------------------------------------------------------------
1081 :
1082 562042208 : exp_rel(m) = c0
1083 562042208 : exp_ret(m) = c0
1084 562042208 : if (tau_ret(m) > c0) then
1085 562042208 : exp_min = min(dt/tau_ret(m),exp_argmax)
1086 562042208 : exp_ret(m) = exp(-exp_min)
1087 : endif
1088 562042208 : if (tau_rel(m) > c0) then
1089 562042208 : exp_min = min(dt/tau_rel(m),exp_argmax)
1090 562042208 : exp_rel(m) = exp(-exp_min)
1091 : endif
1092 562042208 : if (m .ne. nlt_bgc_N(1)) then
1093 534951844 : if (hin_old > hin) then !melting
1094 369263707 : exp_ret(m) = c1
1095 : else !not melting
1096 165688137 : exp_rel(m) = c1
1097 : endif
1098 27090364 : elseif (tr_bgc_N .and. hin_old > hin + algal_vel*dt) then
1099 3294536 : exp_ret(m) = c1
1100 23795828 : elseif (tr_bgc_N) then
1101 23795828 : exp_rel(m) = c1
1102 : endif
1103 :
1104 562042208 : dhtop = dh_top
1105 562042208 : dhbot = dh_bot
1106 562042208 : darcyV = darcy_V
1107 562042208 : C_top(m) = in_init_cons(1,m)*trcrn(nt_zbgc_frac+m-1)!mobile fraction
1108 562042208 : source(m) = abs(zbgc_snow(m) + zbgc_atm(m) + dust_Fe(m))
1109 562042208 : dhflood = max(c0,-dh_direct) ! ocean water flooding surface
1110 :
1111 562042208 : 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 130867301 : + darcyV/bphin_N(1)*dt + puny)*hbri_old
1114 431174907 : elseif (dhtop+darcyV/bphin_N(1)*dt >= -puny .and. &
1115 : abs((zbgc_snow(m)+zbgc_atm(m) + dust_Fe(m)) + & ! LCOV_EXCL_LINE
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 32400950 : ocean_bio(m)*bphin_N(1)*dhflood
1119 : else ! only positive fluxes
1120 398773957 : atm_add_cons(m) = abs(zbgc_snow(m) + zbgc_atm(m)+ dust_Fe(m))
1121 : endif
1122 :
1123 562042208 : C_bot(m) = ocean_bio(m)*hbri_old*iphin_N(nblyr+1)
1124 595877548 : 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 7859589096 : trtmp(:) = c0
1134 7859589096 : trtmp0(:)= c0
1135 304518060 : zfswin(:) = c0
1136 :
1137 304518060 : do k = 1, nilyr+1
1138 : ! contains cice values (fswthrul(1) is surface value)
1139 : ! and fwsthrul(nilyr+1) is output
1140 304518060 : trtmp0(nt_zfswin+k-1) = fswthrul(k)
1141 : enddo !k
1142 :
1143 : call remap_zbgc(nilyr+1, &
1144 : nt_zfswin, & ! LCOV_EXCL_LINE
1145 : trtmp0(1:ntrcr), trtmp(1:ntrcr+2), & ! LCOV_EXCL_LINE
1146 : 0, nblyr+1, & ! LCOV_EXCL_LINE
1147 : hin, hbri, & ! LCOV_EXCL_LINE
1148 : icgrid(1:nilyr+1), & ! LCOV_EXCL_LINE
1149 33835340 : igrid(1:nblyr+1),ice_conc )
1150 :
1151 33835340 : if (icepack_warnings_aborted(subname)) return
1152 :
1153 304518060 : do k = 1,nblyr+1
1154 304518060 : zfswin(k) = trtmp(nt_zfswin+k-1)
1155 : enddo
1156 :
1157 : !-----------------------------------------------------------------
1158 : ! Initialze Biology
1159 : !-----------------------------------------------------------------
1160 :
1161 595877548 : do mm = 1, nbtrcr
1162 562042208 : mobile(mm) = c0
1163 562042208 : if (bgc_tracer_type(mm) .GE. c0) mobile(mm) = c1
1164 :
1165 5092215212 : do k = 1, nblyr+1
1166 5058379872 : biomat_cons(k,mm) = in_init_cons(k,mm)
1167 : enddo !k
1168 : enddo !mm
1169 :
1170 : !-----------------------------------------------------------------
1171 : ! Compute FCT
1172 : !-----------------------------------------------------------------
1173 :
1174 595877548 : do mm = 1, nbtrcr
1175 :
1176 562042208 : if (hbri_old > thinS .and. hbri > thinS) then
1177 4881852153 : do k = 1,nblyr+1
1178 4339424136 : initcons_mobile(k) = in_init_cons(k,mm)*trcrn(nt_zbgc_frac+mm-1)
1179 4339424136 : 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))) + & ! LCOV_EXCL_LINE
1183 4339424136 : (1-mobile(mm))*initcons_stationary(k)*(c1-exp_rel(mm))
1184 4339424136 : initcons_mobile(k) = max(c0,initcons_mobile(k) + dmobile(k))
1185 4339424136 : initcons_stationary(k) = max(c0,initcons_stationary(k) - dmobile(k))
1186 4339424136 : if (initcons_stationary(k)/hbri_old > Sat_conc) then
1187 426144 : initcons_mobile(k) = initcons_mobile(k) + initcons_stationary(k) - Sat_conc*hbri_old
1188 426144 : initcons_stationary(k) = Sat_conc*hbri_old
1189 : endif
1190 :
1191 4339424136 : Diff(k) = iDin(k)
1192 4339424136 : initcons(k) = initcons_mobile(k)
1193 4881852153 : biocons(k) = initcons_mobile(k)
1194 : enddo
1195 :
1196 : call compute_FCT_matrix &
1197 : (initcons,sbdiagz, dt, & ! LCOV_EXCL_LINE
1198 : diagz, spdiagz, rhsz, & ! LCOV_EXCL_LINE
1199 : darcyV, dhtop, & ! LCOV_EXCL_LINE
1200 : dhbot, iphin_N, & ! LCOV_EXCL_LINE
1201 : Diff, hbri_old, & ! LCOV_EXCL_LINE
1202 : atm_add_cons(mm), bphin_N, & ! LCOV_EXCL_LINE
1203 : C_top(mm), C_bot(mm), & ! LCOV_EXCL_LINE
1204 : Source_bot(mm), Source_top(mm),& ! LCOV_EXCL_LINE
1205 : Sink_bot(mm),Sink_top(mm), & ! LCOV_EXCL_LINE
1206 542428017 : D_sbdiag, D_spdiag, ML_diag)
1207 542428017 : if (icepack_warnings_aborted(subname)) return
1208 :
1209 : call tridiag_solverz &
1210 : (nblyr+1, sbdiagz, & ! LCOV_EXCL_LINE
1211 : diagz, spdiagz, & ! LCOV_EXCL_LINE
1212 542428017 : rhsz, biocons)
1213 542428017 : if (icepack_warnings_aborted(subname)) return
1214 :
1215 : call check_conservation_FCT &
1216 : (initcons, & ! LCOV_EXCL_LINE
1217 : biocons, & ! LCOV_EXCL_LINE
1218 : biomat_low, & ! LCOV_EXCL_LINE
1219 : Source_top(mm), & ! LCOV_EXCL_LINE
1220 : Source_bot(mm), & ! LCOV_EXCL_LINE
1221 : Sink_bot(mm), & ! LCOV_EXCL_LINE
1222 : Sink_top(mm), & ! LCOV_EXCL_LINE
1223 : dt, flux_bio(mm), & ! LCOV_EXCL_LINE
1224 542428017 : source(mm))
1225 542428017 : if (icepack_warnings_aborted(subname)) return
1226 :
1227 : call compute_FCT_corr &
1228 : (initcons, biocons, dt, & ! LCOV_EXCL_LINE
1229 542428017 : D_sbdiag, D_spdiag, ML_diag)
1230 542428017 : if (icepack_warnings_aborted(subname)) return
1231 :
1232 542428017 : top_conc = c0 ! or frazil ice concentration
1233 :
1234 : ! assume diatoms actively maintain there relative position in the ice
1235 :
1236 542428017 : if (mm .ne. nlt_bgc_N(1)) then
1237 :
1238 : call regrid_stationary &
1239 : (initcons_stationary, hbri_old, & ! LCOV_EXCL_LINE
1240 : hbri, dt, & ! LCOV_EXCL_LINE
1241 : top_conc, & ! LCOV_EXCL_LINE
1242 : igrid, flux_bio(mm),& ! LCOV_EXCL_LINE
1243 516306855 : meltb, congel)
1244 516306855 : if (icepack_warnings_aborted(subname)) return
1245 :
1246 26121162 : elseif (tr_bgc_N .and. mm .eq. nlt_bgc_N(1)) then
1247 26121162 : if (meltb > algal_vel*dt .or. aicen < 0.001_dbl_kind) then
1248 :
1249 : call regrid_stationary &
1250 : (initcons_stationary, hbri_old, & ! LCOV_EXCL_LINE
1251 : hbri, dt, & ! LCOV_EXCL_LINE
1252 : top_conc, & ! LCOV_EXCL_LINE
1253 : igrid, flux_bio(mm),& ! LCOV_EXCL_LINE
1254 13226444 : meltb, congel)
1255 13226444 : if (icepack_warnings_aborted(subname)) return
1256 :
1257 : endif
1258 : endif
1259 :
1260 4881852153 : biomat_cons(:,mm) = biocons(:) + initcons_stationary(:)
1261 :
1262 542428017 : sum_initial = (in_init_cons(1,mm) + in_init_cons(nblyr+1,mm))*zspace/c2
1263 542428017 : sum_old = (biomat_low(1) + biomat_low(nblyr+1))*zspace/c2
1264 542428017 : sum_new = (biocons(1)+ biocons(nblyr+1))*zspace/c2
1265 542428017 : sum_tot = (biomat_cons(1,mm) + biomat_cons(nblyr+1,mm))*zspace/c2
1266 3796996119 : do k = 2,nblyr
1267 3254568102 : sum_initial = sum_initial + in_init_cons(k,mm)*zspace
1268 3254568102 : sum_old = sum_old + biomat_low(k)*zspace
1269 3254568102 : sum_new = sum_new + biocons(k)*zspace
1270 3796996119 : sum_tot = sum_tot + biomat_cons(k,mm)*zspace
1271 : enddo
1272 542428017 : trcrn(nt_zbgc_frac+mm-1) = zbgc_frac_init(mm)
1273 542428017 : 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) & ! LCOV_EXCL_LINE
1277 10306132323 : .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 542428017 : 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), & ! LCOV_EXCL_LINE
1311 19614191 : dt, ocean_bio(mm))
1312 19614191 : if (icepack_warnings_aborted(subname)) return
1313 :
1314 : endif ! thin or not
1315 :
1316 5092215212 : do k = 1,nblyr+1
1317 5058379872 : biomat_brine(k,mm) = biomat_cons(k,mm)/hbri/iphin_N(k)
1318 : enddo ! k
1319 : enddo ! mm
1320 :
1321 5092215212 : react(:,:) = c0
1322 765275168 : grow_alg(:,:) = c0
1323 :
1324 33835340 : if (solve_zbgc) then
1325 243813276 : do k = 1, nblyr+1
1326 : call algal_dyn (dt, &
1327 : dEdd_algae, & ! LCOV_EXCL_LINE
1328 : zfswin(k), react(k,:), & ! LCOV_EXCL_LINE
1329 : biomat_brine(k,:), & ! LCOV_EXCL_LINE
1330 : grow_alg(k,:), & ! LCOV_EXCL_LINE
1331 : iTin(k), & ! LCOV_EXCL_LINE
1332 : upNOn(k,:), upNHn(k,:), & ! LCOV_EXCL_LINE
1333 : Zoo(k), & ! LCOV_EXCL_LINE
1334 : Cerror(k), conserve_C(k), & ! LCOV_EXCL_LINE
1335 216722912 : nitrification(k))
1336 243813276 : if (icepack_warnings_aborted(subname)) return
1337 :
1338 : enddo ! k
1339 : endif ! solve_zbgc
1340 :
1341 : !-----------------------------------------------------------------
1342 : ! Update the tracer variable
1343 : !-----------------------------------------------------------------
1344 :
1345 33835340 : sum_new = c0
1346 33835340 : sum_tot = c0
1347 :
1348 595877548 : do m = 1,nbtrcr
1349 5058379872 : do k = 1,nblyr+1 ! back to bulk quantity
1350 4496337664 : bio_tmp = (biomat_brine(k,m) + react(k,m))*iphin_N(k)
1351 :
1352 4496337664 : 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 4496337664 : if (m .eq. nlt_bgc_Nit) then
1372 : initcons_mobile(k) = max(c0,(biomat_brine(k,m)-nitrification(k) + &
1373 216722912 : 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 216722912 : nitrification(k) + react(k,m)) + nitrification(k))*iphin_N(k))
1376 :
1377 216722912 : sum_new = sum_new + initcons_mobile(k)*dz(k)
1378 216722912 : sum_tot = sum_tot + (initcons_mobile(k) + initcons_stationary(k))*dz(k)
1379 :
1380 : end if ! m .eq. nlt_bgc_Nit
1381 4496337664 : 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 4496337664 : elseif (abs(bio_tmp) < accuracy) then
1397 281965903 : flux_bio(m) = flux_bio(m) + bio_tmp*dz(k)*hbri/dt
1398 281965903 : bio_tmp = c0
1399 4214371761 : 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 4214371761 : elseif (bio_tmp < c0) then
1413 115555 : write(warnstr,*) subname, 'negative bgc'
1414 115555 : call icepack_warnings_add(warnstr)
1415 115555 : write(warnstr,*) subname, 'k,m,nlt_bgc_Nit,hbri,hbri_old'
1416 115555 : call icepack_warnings_add(warnstr)
1417 115555 : write(warnstr,*) subname, k,m,nlt_bgc_Nit,hbri,hbri_old
1418 115555 : call icepack_warnings_add(warnstr)
1419 115555 : write(warnstr,*) subname, 'bio_tmp,biomat_cons(k,m),ocean_bio(m)'
1420 115555 : call icepack_warnings_add(warnstr)
1421 115555 : write(warnstr,*) subname, bio_tmp,biomat_cons(k,m),ocean_bio(m)
1422 115555 : call icepack_warnings_add(warnstr)
1423 115555 : write(warnstr,*) subname, 'react(k,m),iphin_N(k),biomat_brine(k,m)'
1424 115555 : call icepack_warnings_add(warnstr)
1425 115555 : write(warnstr,*) subname, react(k,m),iphin_N(k),biomat_brine(k,m)
1426 115555 : call icepack_warnings_add(warnstr)
1427 115555 : write(warnstr,*) subname, 'exp_ret( m),exp_ret( m)',exp_ret( m),exp_ret( m)
1428 115555 : call icepack_warnings_add(warnstr)
1429 115555 : call icepack_warnings_add(subname//'negative bgc')
1430 : !call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1431 : endif
1432 4496337664 : trcrn(bio_index(m)+k-1) = max(c0, bio_tmp)
1433 5058379872 : 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 5620422080 : if (m .eq. nlt_bgc_Nit .and. MAXVAL(nitrification) > c0) then
1448 27090364 : trcrn(nt_zbgc_frac+m-1) = zbgc_frac_init(m)
1449 27090364 : if (sum_tot > c0) trcrn(nt_zbgc_frac+m-1) = sum_new/sum_tot
1450 : end if
1451 595877548 : 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 216722912 : subroutine algal_dyn (dt, &
1463 : dEdd_algae, & ! LCOV_EXCL_LINE
1464 0 : fswthru, reactb, & ! LCOV_EXCL_LINE
1465 216722912 : ltrcrn, & ! LCOV_EXCL_LINE
1466 216722912 : grow_alg, & ! LCOV_EXCL_LINE
1467 : T_bot, & ! LCOV_EXCL_LINE
1468 216722912 : upNOn, upNHn, & ! LCOV_EXCL_LINE
1469 : Zoo, & ! LCOV_EXCL_LINE
1470 : Cerror, conserve_C, & ! LCOV_EXCL_LINE
1471 : nitrification)
1472 :
1473 : real (kind=dbl_kind), intent(in) :: &
1474 : dt , & ! time step ! LCOV_EXCL_LINE
1475 : T_bot , & ! ice temperature (oC) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
1480 : Cerror, & ! Change in C after reactions (mmol/m^3) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
1485 : upNOn, & ! algal NO uptake rate (mmol/m^3/s) ! LCOV_EXCL_LINE
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 433445824 : Nin , & ! algal nitrogen concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1516 433445824 : Cin , & ! algal carbon concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1517 433445824 : chlin ! algal chlorophyll concentration on volume (mg/m^3)
1518 :
1519 : real (kind=dbl_kind), dimension(n_doc) :: &
1520 433445824 : DOCin ! dissolved organic carbon concentration on volume (mmolC/m^3)
1521 :
1522 : real (kind=dbl_kind), dimension(n_dic) :: &
1523 433445824 : DICin ! dissolved inorganic carbon concentration on volume (mmolC/m^3)
1524 :
1525 : real (kind=dbl_kind), dimension(n_don) :: & !proteins
1526 433445824 : DONin ! dissolved organic nitrogen concentration on volume (mmolN/m^3)
1527 :
1528 : real (kind=dbl_kind), dimension(n_fed) :: & !iron
1529 433445824 : Fedin ! dissolved iron concentration on volume (umol/m^3)
1530 :
1531 : real (kind=dbl_kind), dimension(n_fep) :: & !iron
1532 433445824 : Fepin ! algal nitrogen concentration on volume (umol/m^3)
1533 :
1534 : real (kind=dbl_kind) :: &
1535 : Nitin , & ! nitrate concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1536 : Amin , & ! ammonia/um concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1537 : Silin , & ! silicon concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1538 : ! DMSPpin , & ! DMSPp concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1539 : DMSPdin , & ! DMSPd concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1540 : DMSin , & ! DMS concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1541 : ! PONin , & ! PON concentration on volume (mmol/m^3) ! LCOV_EXCL_LINE
1542 : op_dep , & ! bottom layer attenuation exponent (optical depth) ! LCOV_EXCL_LINE
1543 : Iavg_loc ! bottom layer attenuated Fswthru (W/m^2)
1544 :
1545 : real (kind=dbl_kind), dimension(n_algae) :: &
1546 433445824 : L_lim , & ! overall light limitation ! LCOV_EXCL_LINE
1547 433445824 : Nit_lim , & ! overall nitrate limitation ! LCOV_EXCL_LINE
1548 433445824 : Am_lim , & ! overall ammonium limitation ! LCOV_EXCL_LINE
1549 433445824 : N_lim , & ! overall nitrogen species limitation ! LCOV_EXCL_LINE
1550 433445824 : Sil_lim , & ! overall silicon limitation ! LCOV_EXCL_LINE
1551 433445824 : Fe_lim , & ! overall iron limitation ! LCOV_EXCL_LINE
1552 433445824 : fr_Nit , & ! fraction of local ecological growth as nitrate ! LCOV_EXCL_LINE
1553 433445824 : fr_Am , & ! fraction of local ecological growth as ammonia ! LCOV_EXCL_LINE
1554 433445824 : growmax_N, & ! maximum growth rate in N currency (mmol/m^3/s) ! LCOV_EXCL_LINE
1555 433445824 : grow_N , & ! true growth rate in N currency (mmol/m^3/s) ! LCOV_EXCL_LINE
1556 433445824 : potU_Am , & ! potential ammonium uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1557 433445824 : U_Nit , & ! actual nitrate uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1558 433445824 : U_Am , & ! actual ammonium uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1559 433445824 : U_Sil , & ! actual silicon uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1560 433445824 : U_Fe , & ! actual iron uptake (umol/m^3/s) ! LCOV_EXCL_LINE
1561 433445824 : U_Nit_f , & ! fraction of Nit uptake due to each algal species ! LCOV_EXCL_LINE
1562 433445824 : U_Am_f , & ! fraction of Am uptake due to each algal species ! LCOV_EXCL_LINE
1563 433445824 : U_Sil_f , & ! fraction of Sil uptake due to each algal species ! LCOV_EXCL_LINE
1564 433445824 : 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 ! LCOV_EXCL_LINE
1568 : U_Nit_tot , & ! actual nitrate uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1569 : U_Am_tot , & ! actual ammonium uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1570 : U_Sil_tot , & ! actual silicon uptake (mmol/m^3/s) ! LCOV_EXCL_LINE
1571 : U_Fe_tot , & ! actual iron uptake (umol/m^3/s) ! LCOV_EXCL_LINE
1572 : nitrif , & ! nitrification (mmol/m^3/s) ! LCOV_EXCL_LINE
1573 : mort_N , & ! total algal mortality (mmol N/m^3) ! LCOV_EXCL_LINE
1574 : mort_C , & ! total algal mortality (mmol C/m^3) ! LCOV_EXCL_LINE
1575 : graze_N , & ! total algae grazed (mmol N/m^3) ! LCOV_EXCL_LINE
1576 : graze_C , & ! total algae grazed (mmol C/m^3) ! LCOV_EXCL_LINE
1577 : exude_C , & ! total carbon exuded by algae (mmol C/m^3) ! LCOV_EXCL_LINE
1578 : resp_N , & ! total N in respiration (mmol N/m^3) ! LCOV_EXCL_LINE
1579 : growth_N ! total algal growth (mmol N/m^3)
1580 :
1581 : real (kind=dbl_kind), dimension(n_algae) :: &
1582 433445824 : resp , & ! respiration (mmol/m^3/s) ! LCOV_EXCL_LINE
1583 433445824 : graze , & ! grazing (mmol/m^3/s) ! LCOV_EXCL_LINE
1584 433445824 : 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 433445824 : N_s , & ! net algal nitrogen sources (mmol/m^3) ! LCOV_EXCL_LINE
1590 433445824 : N_r ! net algal nitrogen removal (mmol/m^3)
1591 :
1592 : real (kind=dbl_kind), dimension(n_doc) :: &
1593 433445824 : DOC_r , & ! net DOC removal (mmol/m^3) ! LCOV_EXCL_LINE
1594 433445824 : DOC_s ! net DOC sources (mmol/m^3)
1595 :
1596 : real (kind=dbl_kind), dimension(n_dic) :: &
1597 433445824 : DIC_r , & ! net DIC removal (mmol/m^3) ! LCOV_EXCL_LINE
1598 433445824 : DIC_s ! net DIC sources (mmol/m^3)
1599 :
1600 : real (kind=dbl_kind), dimension(n_don) :: &
1601 433445824 : DON_r , & ! net DON removal (mmol/m^3) ! LCOV_EXCL_LINE
1602 433445824 : 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) ! LCOV_EXCL_LINE
1606 433445824 : Fed_r , & ! net Fed removal (nM) ! LCOV_EXCL_LINE
1607 433445824 : Fed_s , & ! net Fed sources (nM) ! LCOV_EXCL_LINE
1608 433445824 : rFed ! ratio of dissolved Fe to tot Fed
1609 :
1610 : real (kind=dbl_kind), dimension(n_fep) :: &
1611 433445824 : Fep_r , & ! net Fep removal (nM) ! LCOV_EXCL_LINE
1612 433445824 : Fep_s , & ! net Fep sources (nM) ! LCOV_EXCL_LINE
1613 433445824 : rFep ! ratio of particulate Fe to tot Fep
1614 :
1615 : real (kind=dbl_kind) :: &
1616 : dN , & ! change in N (mmol/m^3) ! LCOV_EXCL_LINE
1617 : dC , & ! change in Carbon (mmol C/m^3) ! LCOV_EXCL_LINE
1618 : N_s_p , & ! algal nitrogen photosynthesis (mmol/m^3) ! LCOV_EXCL_LINE
1619 : N_r_g , & ! algal nitrogen losses to grazing (mmol/m^3) ! LCOV_EXCL_LINE
1620 : N_r_r , & ! algal nitrogen losses to respiration (mmol/m^3) ! LCOV_EXCL_LINE
1621 : N_r_mo , & ! algal nitrogen losses to mortality (mmol/m^3) ! LCOV_EXCL_LINE
1622 : Nit_s_n , & ! nitrate from nitrification (mmol/m^3) ! LCOV_EXCL_LINE
1623 : Nit_r_p , & ! nitrate uptake by algae (mmol/m^3) ! LCOV_EXCL_LINE
1624 : Nit_s , & ! net nitrate sources (mmol/m^3) ! LCOV_EXCL_LINE
1625 : Nit_r , & ! net nitrate removal (mmol/m^3) ! LCOV_EXCL_LINE
1626 : Am_s_e , & ! ammonium source from excretion (mmol/m^3) ! LCOV_EXCL_LINE
1627 : Am_s_r , & ! ammonium source from respiration (mmol/m^3) ! LCOV_EXCL_LINE
1628 : Am_s_mo , & ! ammonium source from mort/remin (mmol/m^3) ! LCOV_EXCL_LINE
1629 : Am_r_p , & ! ammonium uptake by algae (mmol/m^3) ! LCOV_EXCL_LINE
1630 : Am_s , & ! net ammonium sources (mmol/m^3) ! LCOV_EXCL_LINE
1631 : Am_r , & ! net ammonium removal (mmol/m^3) ! LCOV_EXCL_LINE
1632 : Sil_r_p , & ! silicon uptake by algae (mmol/m^3) ! LCOV_EXCL_LINE
1633 : Sil_r , & ! net silicon removal (mmol/m^3) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
1641 : DMSPd_s_mo, & ! skl dissolved DMSP from MBJ algal mortexc (mmol/m^3) ! LCOV_EXCL_LINE
1642 : DMSPd_r , & ! skl dissolved DMSP conversion (mmol/m^3) DMSPD_sk_r ! LCOV_EXCL_LINE
1643 : DMSPd_s , & ! net skl dissolved DMSP sources (mmol/m^3) ! LCOV_EXCL_LINE
1644 : DMS_s_c , & ! skl DMS source via conversion (mmol/m^3) ! LCOV_EXCL_LINE
1645 : DMS_r_o , & ! skl DMS losses due to oxidation (mmol/m^3) ! LCOV_EXCL_LINE
1646 : DMS_s , & ! net skl DMS sources (mmol/m^3) ! LCOV_EXCL_LINE
1647 : DMS_r , & ! net skl DMS removal (mmol/m^3) ! LCOV_EXCL_LINE
1648 : Fed_tot , & ! total dissolved iron from all sources (nM) ! LCOV_EXCL_LINE
1649 : Fed_tot_r , & ! total dissolved iron losses (nM) ! LCOV_EXCL_LINE
1650 : Fed_tot_s , & ! total dissolved iron sources (nM) ! LCOV_EXCL_LINE
1651 : Fep_tot , & ! total particulate iron from all sources (nM) ! LCOV_EXCL_LINE
1652 : ! Fep_tot_r , & ! total particulate iron losses (nM) ! LCOV_EXCL_LINE
1653 : Fep_tot_s , & ! total particulate iron sources (nM) ! LCOV_EXCL_LINE
1654 : Zoo_s_a , & ! N Losses due to zooplankton assimilation (mmol/m^3) ! LCOV_EXCL_LINE
1655 : Zoo_s_s , & ! N Losses due to grazing spillage (mmol/m^3) ! LCOV_EXCL_LINE
1656 : Zoo_s_m , & ! N Losses due to algal mortality (mmol/m^3) ! LCOV_EXCL_LINE
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 216722912 : conserve_C = .true.
1666 866891648 : Nin(:) = c0
1667 866891648 : Cin(:) = c0
1668 866891648 : chlin(:) = c0
1669 650168736 : DOCin(:) = c0
1670 433445824 : DICin(:) = c0
1671 433445824 : DONin(:) = c0
1672 433445824 : Fedin(:) = c0
1673 433445824 : Fepin(:) = c0
1674 216722912 : Nitin = c0
1675 216722912 : Amin = c0
1676 216722912 : Silin = c0
1677 216722912 : DMSPdin = c0
1678 216722912 : DMSin = c0
1679 216722912 : U_Am_tot = c0
1680 216722912 : U_Nit_tot = c0
1681 216722912 : U_Sil_tot = c0
1682 216722912 : U_Fe_tot = c0
1683 866891648 : U_Am_f(:) = c0
1684 866891648 : U_Nit_f(:) = c0
1685 866891648 : U_Sil_f(:) = c0
1686 866891648 : U_Fe_f(:) = c0
1687 650168736 : DOC_s(:) = c0
1688 650168736 : DOC_r(:) = c0
1689 216722912 : nitrif = c0
1690 216722912 : mort_N = c0
1691 216722912 : mort_C = c0
1692 216722912 : graze_N = c0
1693 216722912 : graze_C = c0
1694 216722912 : exude_C = c0
1695 216722912 : resp_N = c0
1696 216722912 : growth_N = c0
1697 216722912 : Nit_r = c0
1698 216722912 : Am_s = c0
1699 216722912 : Am_r = c0
1700 216722912 : Sil_r = c0
1701 433445824 : Fed_r(:) = c0
1702 433445824 : Fed_s(:) = c0
1703 433445824 : Fep_r(:) = c0
1704 433445824 : Fep_s(:) = c0
1705 216722912 : DMSPd_s = c0
1706 216722912 : dTemp = min(T_bot-T_max,c0)
1707 216722912 : Fed_tot = c0
1708 216722912 : Fed_tot_r = c0
1709 216722912 : Fed_tot_s = c0
1710 433445824 : rFed(:) = c0
1711 216722912 : Fep_tot = c0
1712 216722912 : Fep_tot_s = c0
1713 433445824 : rFep(:) = c0
1714 433445824 : DIC_r(:) = c0
1715 433445824 : DIC_s(:) = c0
1716 216722912 : Zoo = c0
1717 :
1718 216722912 : Nitin = ltrcrn(nlt_bgc_Nit)
1719 216722912 : op_dep = c0
1720 866891648 : do k = 1, n_algae
1721 650168736 : Nin(k) = ltrcrn(nlt_bgc_N(k))
1722 650168736 : chlin(k) = R_chl2N(k)* Nin(k)
1723 866891648 : op_dep = op_dep + chlabs(k)*chlin(k)
1724 : enddo
1725 216722912 : if (tr_bgc_C) then
1726 : ! do k = 1, n_algae
1727 : ! Cin(k)= ltrcrn(nlt_bgc_C(k))
1728 : ! enddo
1729 650168736 : do k = 1, n_doc
1730 650168736 : DOCin(k)= ltrcrn(nlt_bgc_DOC(k))
1731 : enddo
1732 433445824 : do k = 1, n_dic
1733 433445824 : DICin(k)= ltrcrn(nlt_bgc_DIC(k))
1734 : enddo
1735 : endif
1736 216722912 : if (tr_bgc_Am) Amin = ltrcrn(nlt_bgc_Am)
1737 216722912 : if (tr_bgc_Sil) Silin = ltrcrn(nlt_bgc_Sil)
1738 216722912 : if (tr_bgc_DMS) then
1739 : ! DMSPpin = ltrcrn(nlt_bgc_DMSPp)
1740 216722912 : DMSPdin = ltrcrn(nlt_bgc_DMSPd)
1741 216722912 : 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 216722912 : if (tr_bgc_DON) then
1748 433445824 : do k = 1, n_don
1749 433445824 : DONin(k) = ltrcrn(nlt_bgc_DON(k))
1750 : enddo
1751 : endif
1752 216722912 : if (tr_bgc_Fe ) then
1753 433445824 : do k = 1, n_fed
1754 433445824 : Fedin(k) = ltrcrn(nlt_bgc_Fed(k))
1755 : enddo
1756 433445824 : do k = 1, n_fep
1757 433445824 : Fepin(k) = ltrcrn(nlt_bgc_Fep(k))
1758 : enddo
1759 : endif
1760 :
1761 : !-----------------------------------------------------------------------
1762 : ! Total iron from all pools
1763 : !-----------------------------------------------------------------------
1764 :
1765 433445824 : do k = 1,n_fed
1766 433445824 : Fed_tot = Fed_tot + Fedin(k)
1767 : enddo
1768 433445824 : do k = 1,n_fep
1769 433445824 : Fep_tot = Fep_tot + Fepin(k)
1770 : enddo
1771 216722912 : if (Fed_tot > puny) then
1772 433445824 : do k = 1,n_fed
1773 433445824 : rFed(k) = Fedin(k)/Fed_tot
1774 : enddo
1775 : endif
1776 216722912 : if (Fep_tot > puny) then
1777 433445824 : do k = 1,n_fep
1778 433445824 : rFep(k) = Fepin(k)/Fep_tot
1779 : enddo
1780 : endif
1781 :
1782 : !-----------------------------------------------------------------------
1783 : ! Light limitation (op_dep) defined above
1784 : !-----------------------------------------------------------------------
1785 :
1786 216722912 : if (op_dep > op_dep_min .and. .not. dEdd_algae) then
1787 212239968 : Iavg_loc = fswthru * (c1 - exp(-op_dep)) / op_dep
1788 : else
1789 4482944 : Iavg_loc = fswthru
1790 : endif
1791 :
1792 866891648 : do k = 1, n_algae
1793 : ! With light inhibition ! Maybe include light inhibition for diatoms but phaeocystis
1794 2600674944 : 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 650168736 : Nit_lim(k) = Nitin/(Nitin + K_Nit(k))
1804 650168736 : Am_lim(k) = c0
1805 650168736 : N_lim(k) = Nit_lim(k)
1806 650168736 : if (tr_bgc_Am) then
1807 650168736 : Am_lim(k) = Amin/(Amin + K_Am(k))
1808 650168736 : N_lim(k) = min(c1, Nit_lim(k) + Am_lim(k))
1809 : endif
1810 650168736 : Sil_lim(k) = c1
1811 650168736 : 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 650168736 : Fe_lim(k) = c1
1818 650168736 : 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 650168736 : growmax_N(k) = mu_max(k) / secday * exp(grow_Tdep(k) * dTemp)* Nin(k) *fsal
1829 650168736 : 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 650168736 : 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 650168736 : potU_Am(k) = Am_lim(k)* growmax_N(k)
1836 650168736 : U_Am(k) = min(grow_N(k), potU_Am(k))
1837 650168736 : U_Nit(k) = grow_N(k) - U_Am(k)
1838 650168736 : U_Sil(k) = R_Si2N(k) * grow_N(k)
1839 650168736 : U_Fe (k) = R_Fe2N(k) * grow_N(k)
1840 :
1841 650168736 : U_Am_tot = U_Am_tot + U_Am(k)
1842 650168736 : U_Nit_tot = U_Nit_tot + U_Nit(k)
1843 650168736 : U_Sil_tot = U_Sil_tot + U_Sil(k)
1844 866891648 : U_Fe_tot = U_Fe_tot + U_Fe(k)
1845 : enddo
1846 866891648 : do k = 1, n_algae
1847 650168736 : if (U_Am_tot > c0) U_Am_f(k) = U_Am(k)/U_Am_tot
1848 650168736 : if (U_Nit_tot > c0) U_Nit_f(k) = U_Nit(k)/U_Nit_tot
1849 650168736 : if (U_Sil_tot > c0) U_Sil_f(k) = U_Sil(k)/U_Sil_tot
1850 866891648 : if (U_Fe_tot > c0) U_Fe_f(k) = U_Fe(k)/U_Fe_tot
1851 : enddo
1852 :
1853 216722912 : if (tr_bgc_Sil) U_Sil_tot = min(U_Sil_tot, max_loss * Silin/dt)
1854 216722912 : if (tr_bgc_Fe) U_Fe_tot = min(U_Fe_tot, max_loss * Fed_tot/dt)
1855 216722912 : U_Nit_tot = min(U_Nit_tot, max_loss * Nitin/dt)
1856 216722912 : U_Am_tot = min(U_Am_tot, max_loss * Amin/dt)
1857 :
1858 866891648 : do k = 1, n_algae
1859 650168736 : U_Am(k) = U_Am_f(k)*U_Am_tot
1860 650168736 : U_Nit(k) = U_Nit_f(k)*U_Nit_tot
1861 650168736 : U_Sil(k) = U_Sil_f(k)*U_Sil_tot
1862 650168736 : U_Fe(k) = U_Fe_f(k)*U_Fe_tot
1863 :
1864 650168736 : 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 650168736 : fr_Am(k) = c0
1873 650168736 : if (tr_bgc_Am) then
1874 650168736 : fr_Am(k) = p5
1875 650168736 : if (grow_N(k) > c0) fr_Am(k) = min(U_Am(k)/grow_N(k), c1)
1876 : endif
1877 650168736 : fr_Nit(k) = c1 - fr_Am(k)
1878 650168736 : U_Nit(k) = fr_Nit(k) * grow_N(k)
1879 650168736 : U_Am(k) = fr_Am(k) * grow_N(k)
1880 650168736 : U_Sil(k) = R_Si2N(k) * grow_N(k)
1881 650168736 : 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 650168736 : resp(k) = fr_resp * grow_N(k)
1902 650168736 : 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 650168736 : mort_pre(k)*exp(mort_Tdep(k)*dTemp) * Nin(k)/secday)
1905 :
1906 : ! history variables
1907 650168736 : grow_alg(k) = grow_N(k)
1908 650168736 : upNOn(k) = U_Nit(k)
1909 650168736 : upNHn(k) = U_Am(k)
1910 :
1911 650168736 : N_s_p = grow_N(k) * dt
1912 650168736 : N_r_g = graze(k) * dt
1913 650168736 : N_r_r = resp(k) * dt
1914 650168736 : N_r_mo = mort(k) * dt
1915 650168736 : N_s(k) = N_s_p
1916 650168736 : N_r(k) = N_r_g + N_r_mo + N_r_r
1917 650168736 : graze_N = graze_N + graze(k)
1918 650168736 : graze_C = graze_C + R_C2N(k)*graze(k)
1919 650168736 : mort_N = mort_N + mort(k)
1920 650168736 : mort_C = mort_C + R_C2N(k)*mort(k)
1921 650168736 : resp_N = resp_N + resp(k)
1922 866891648 : growth_N = growth_N + grow_N(k)
1923 :
1924 : enddo ! n_algae
1925 : !--------------------------------------------------------------------
1926 : ! Ammonium source: algal grazing, respiration, and mortality
1927 : !--------------------------------------------------------------------
1928 :
1929 216722912 : Am_s_e = graze_N*(c1-fr_graze_s)*fr_graze_e*dt
1930 216722912 : Am_s_mo = mort_N*fr_mort2min*dt
1931 216722912 : Am_s_r = resp_N*dt
1932 216722912 : Am_s = Am_s_r + Am_s_e + Am_s_mo
1933 :
1934 : !--------------------------------------------------------------------
1935 : ! Nutrient net loss terms: algal uptake
1936 : !--------------------------------------------------------------------
1937 :
1938 866891648 : do k = 1, n_algae
1939 650168736 : Am_r_p = U_Am(k) * dt
1940 650168736 : Am_r = Am_r + Am_r_p
1941 650168736 : Nit_r_p = U_Nit(k) * dt
1942 650168736 : Nit_r = Nit_r + Nit_r_p
1943 650168736 : Sil_r_p = U_Sil(k) * dt
1944 650168736 : Sil_r = Sil_r + Sil_r_p
1945 650168736 : Fe_r_p = U_Fe (k) * dt
1946 650168736 : Fed_tot_r = Fed_tot_r + Fe_r_p
1947 650168736 : exude_C = exude_C + k_exude(k)* R_C2N(k)*Nin(k) / secday
1948 866891648 : 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 216722912 : nitrification = c0
1956 216722912 : nitrif = k_nitrif /secday * Amin
1957 216722912 : Am_r = Am_r + nitrif*dt
1958 216722912 : Nit_s_n = nitrif * dt ! source from NH4
1959 216722912 : 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 216722912 : if (tr_bgc_Am) then
1971 216722912 : Zoo_s_a = graze_N*(c1-fr_graze_e)*(c1-fr_graze_s) *dt
1972 216722912 : Zoo_s_s = graze_N*fr_graze_s*dt
1973 216722912 : 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 216722912 : Zoo_s_b = c0
1981 :
1982 : !--------------------------------------------------------------------
1983 : ! DON (n_don = 1)
1984 : ! Proteins
1985 : !--------------------------------------------------------------------
1986 :
1987 433445824 : DON_r(:) = c0
1988 433445824 : DON_s(:) = c0
1989 :
1990 216722912 : if (tr_bgc_DON) then
1991 433445824 : do n = 1, n_don
1992 216722912 : DON_r(n) = kn_bac(n)/secday * DONin(n) * dt
1993 216722912 : DON_s(n) = graze_N*dt - Am_s_e + mort_N*dt - Am_s_mo
1994 216722912 : Zoo_s_s = Zoo_s_s - DON_s(n)
1995 216722912 : Zoo_s_b = Zoo_s_b + DON_r(n)*(c1-f_don_Am(n))
1996 216722912 : Am_s = Am_s + DON_r(n)*f_don_Am(n)
1997 433445824 : DIC_s(1) = DIC_s(1) + DON_r(n) * R_C2N_DON(n)
1998 : enddo
1999 : endif
2000 :
2001 216722912 : Zoo = Zoo_s_a + Zoo_s_s + Zoo_s_m + Zoo_s_b
2002 :
2003 : !--------------------------------------------------------------------
2004 : ! DOC
2005 : ! polysaccharids, lipids
2006 : !--------------------------------------------------------------------
2007 :
2008 650168736 : do n = 1, n_doc
2009 :
2010 433445824 : 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 433445824 : DOC_s(n) = f_doc(n) * (graze_C*dt + mort_C*dt - DON_s(1) * R_C2N_DON(1))
2014 650168736 : 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 216722912 : 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 216722912 : if (tr_bgc_C .and. tr_bgc_Fe) then
2030 216722912 : 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 216722912 : if (Fed_tot/DOCin(1) < max_dfe_doc1) then
2041 433425986 : do n = 1,n_fep ! high saccharid:dFe ratio leads to
2042 216712993 : Fep_r(n) = Fepin(n)/t_iron_conv*dt/secday ! gain of bioavailable Fe from particulate fraction
2043 433425986 : Fed_tot_s = Fed_tot_s + Fep_r(n)
2044 : enddo
2045 433425986 : do n = 1,n_fed
2046 433425986 : 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 216722912 : if (tr_bgc_Fe) then
2052 433445824 : do n = 1,n_fed
2053 433445824 : 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 433445824 : do n = 1,n_fep
2060 433445824 : 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 866891648 : do k = 1,n_algae
2080 650168736 : DMSPd_s_r = fr_resp_s * R_S2N(k) * resp(k) * dt !respiration fraction to DMSPd
2081 650168736 : DMSPd_s_mo= fr_mort2min * R_S2N(k)* mort(k) * dt !mortality and extracellular excretion
2082 866891648 : DMSPd_s = DMSPd_s + DMSPd_s_r + DMSPd_s_mo
2083 : enddo
2084 216722912 : 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 216722912 : DMS_s_c = y_sk_DMS * DMSPd_r
2092 216722912 : DMS_r_o = DMSin * dt / (t_sk_ox * secday)
2093 216722912 : DMS_s = DMS_s_c
2094 216722912 : DMS_r = DMS_r_o
2095 :
2096 : !-----------------------------------------------------------------------
2097 : ! Load reaction array
2098 : !-----------------------------------------------------------------------
2099 :
2100 216722912 : dN = c0
2101 216722912 : dC = c0
2102 866891648 : do k = 1,n_algae
2103 650168736 : reactb(nlt_bgc_N(k)) = N_s(k) - N_r(k)
2104 650168736 : dN = dN + reactb(nlt_bgc_N(k))
2105 866891648 : dC = dC + reactb(nlt_bgc_N(k)) * R_C2N(k)
2106 : enddo
2107 216722912 : 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 650168736 : do k = 1,n_doc
2112 433445824 : reactb(nlt_bgc_DOC(k))= DOC_s(k) - DOC_r(k)
2113 650168736 : dC = dC + reactb(nlt_bgc_DOC(k))
2114 : enddo
2115 433445824 : do k = 1,n_dic
2116 216722912 : reactb(nlt_bgc_DIC(k))= DIC_s(k) - DIC_r(k)
2117 433445824 : dC = dC + reactb(nlt_bgc_DIC(k))
2118 : enddo
2119 : endif
2120 216722912 : reactb(nlt_bgc_Nit) = Nit_s - Nit_r
2121 216722912 : nitrification = Nit_s_n
2122 216722912 : dN = dN + reactb(nlt_bgc_Nit)
2123 216722912 : if (tr_bgc_Am) then
2124 216722912 : reactb(nlt_bgc_Am) = Am_s - Am_r
2125 216722912 : dN = dN + reactb(nlt_bgc_Am)
2126 : endif
2127 216722912 : if (tr_bgc_Sil) then
2128 216722912 : reactb(nlt_bgc_Sil) = - Sil_r
2129 : endif
2130 216722912 : if (tr_bgc_DON) then
2131 433445824 : do k = 1,n_don
2132 216722912 : reactb(nlt_bgc_DON(k))= DON_s(k) - DON_r(k)
2133 216722912 : dN = dN + reactb(nlt_bgc_DON(k))
2134 433445824 : dC = dC + reactb(nlt_bgc_DON(k)) * R_C2N_DON(k)
2135 : enddo
2136 : endif
2137 216722912 : Cerror = dC
2138 216722912 : if (tr_bgc_Fe ) then
2139 433445824 : do k = 1,n_fed
2140 433445824 : reactb(nlt_bgc_Fed(k))= Fed_s (k) - Fed_r (k)
2141 : enddo
2142 433445824 : do k = 1,n_fep
2143 433445824 : reactb(nlt_bgc_Fep(k))= Fep_s (k) - Fep_r (k)
2144 : enddo
2145 : endif
2146 216722912 : if (tr_bgc_DMS) then
2147 216722912 : reactb(nlt_bgc_DMSPd) = DMSPd_s - DMSPd_r
2148 216722912 : reactb(nlt_bgc_DMS) = DMS_s - DMS_r
2149 : endif
2150 :
2151 216722912 : if (tr_bgc_C) then
2152 9319085216 : 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 216722912 : 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 19614191 : 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) ! LCOV_EXCL_LINE
2237 : hin , & ! new brine thickness (m) ! LCOV_EXCL_LINE
2238 : dt , & ! time step ! LCOV_EXCL_LINE
2239 : source , & ! atm, ocean, dust flux (mmol/m^2) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2254 : zspace, & ! 1/nblyr ! LCOV_EXCL_LINE
2255 : dC, & ! added ocean bio mass (mmol/m^2) ! LCOV_EXCL_LINE
2256 : dh ! change in thickness (m)
2257 :
2258 : character(len=*),parameter :: subname='(thin_ice_flux)'
2259 :
2260 19614191 : zspace = c1/real(nblyr,kind=dbl_kind)
2261 :
2262 19614191 : dC = c0
2263 19614191 : sum_bio = c0
2264 19614191 : dh = hin-hin_old
2265 :
2266 19614191 : if (dh .le. c0 .and. hin_old > puny) then ! keep the brine concentration fixed
2267 16557703 : sum_bio = (Cin(1)+Cin(nblyr+1))/hin_old*zspace*p5
2268 16557703 : Cin(1) = Cin(1)/hin_old*hin
2269 16557703 : Cin(nblyr+1) = Cin(nblyr+1)/hin_old*hin
2270 115903921 : do k = 2, nblyr
2271 99346218 : sum_bio = sum_bio + Cin(k)/hin_old*zspace
2272 115903921 : Cin(k) = Cin(k)/hin_old*hin + dC
2273 : enddo
2274 : else ! spread evenly in ice layers
2275 3056488 : dC = dh*ocean_bio
2276 27508392 : do k = 1, nblyr+1
2277 27508392 : Cin(k) = Cin(k) + dC
2278 : enddo
2279 : endif
2280 :
2281 19614191 : flux_o_tot = - dh*sum_bio/dt - dC/dt + source/dt
2282 :
2283 19614191 : 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 542428017 : subroutine compute_FCT_matrix (C_in, sbdiag, dt, &
2293 542428017 : diag, spdiag, rhs, & ! LCOV_EXCL_LINE
2294 : darcyV, dhtop, dhbot, & ! LCOV_EXCL_LINE
2295 542428017 : iphin_N, iDin, hbri_old, & ! LCOV_EXCL_LINE
2296 542428017 : atm_add, bphin_N, & ! LCOV_EXCL_LINE
2297 : C_top, C_bot, Qbot, Qtop, & ! LCOV_EXCL_LINE
2298 : Sink_bot, Sink_top, & ! LCOV_EXCL_LINE
2299 542428017 : 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 ! LCOV_EXCL_LINE
2319 : diag , & ! diagonal matrix elements ! LCOV_EXCL_LINE
2320 : spdiag , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
2321 : rhs , & ! rhs of tri-diagonal matrix eqn. ! LCOV_EXCL_LINE
2322 : ML, & ! lumped mass matrix ! LCOV_EXCL_LINE
2323 : D_spdiag, D_sbdiag ! artificial diffusion matrix
2324 :
2325 : real (kind=dbl_kind), intent(in) :: &
2326 : dhtop , & ! Change in brine top (m) ! LCOV_EXCL_LINE
2327 : dhbot , & ! Change in brine bottom (m) ! LCOV_EXCL_LINE
2328 : hbri_old , & ! brine height (m) ! LCOV_EXCL_LINE
2329 : atm_add , & ! atm-ice flux ! LCOV_EXCL_LINE
2330 : C_top , & ! bulk surface source (mmol/m^2) ! LCOV_EXCL_LINE
2331 : C_bot , & ! bulk bottom source (mmol/m^2) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2336 : Qbot , & ! bottom flux source (mmol/m^2/s) ! LCOV_EXCL_LINE
2337 : Sink_bot , & ! rest of bottom flux (sink or source) (mmol/m^2/s) ! LCOV_EXCL_LINE
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 1084856034 : Q_top, Q_bot, & ! surface and bottom source ! LCOV_EXCL_LINE
2350 1084856034 : K_diag, K_spdiag, K_sbdiag, & ! advection matrix ! LCOV_EXCL_LINE
2351 1084856034 : S_diag, S_spdiag, S_sbdiag, & ! diffusion matrix ! LCOV_EXCL_LINE
2352 1084856034 : D_diag, iDin_phi
2353 :
2354 : real (kind=dbl_kind), dimension (nblyr) :: &
2355 1084856034 : 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 4339424136 : kvector1(:) = c0
2365 542428017 : kvector1(1) = c1
2366 4339424136 : kvectorn1(:) = c1
2367 542428017 : kvectorn1(1) = c0
2368 :
2369 542428017 : zspace = c1/real(nblyr,kind=dbl_kind)
2370 542428017 : Qbot = c0
2371 542428017 : Qtop = c0
2372 542428017 : Sink_bot = c0
2373 542428017 : Sink_top = c0
2374 :
2375 : ! compute the lumped mass matrix
2376 :
2377 4881852153 : ML(:) = zspace
2378 542428017 : ML(1) = zspace/c2
2379 542428017 : 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 4881852153 : K_diag(:) = c0
2385 4881852153 : D_diag(:) = c0
2386 4881852153 : D_spdiag(:) = c0
2387 4881852153 : D_sbdiag(:) = c0
2388 4881852153 : K_spdiag(:) = c0
2389 4881852153 : K_sbdiag(:) = c0
2390 4881852153 : S_diag(:) = c0
2391 4881852153 : S_spdiag(:) = c0
2392 4881852153 : S_sbdiag(:) = c0
2393 4881852153 : iDin_phi(:) = c0
2394 :
2395 542428017 : iDin_phi(1) = c0 !element 1
2396 542428017 : iDin_phi(nblyr+1) = iDin(nblyr+1)/iphin_N(nblyr+1) !outside ice at bottom
2397 542428017 : iDin_phi(nblyr) = p5*(iDin(nblyr+1)/iphin_N(nblyr+1)+iDin(nblyr)/iphin_N(nblyr))
2398 :
2399 542428017 : vel = (bgrid(2)*dhbot - (bgrid(2)-c1)*dhtop)/dt+darcyV/bphin_N(2)
2400 542428017 : K_diag(1) = p5*vel/hbri_old
2401 542428017 : dphi_dx = (iphin_N(nblyr+1) - iphin_N(nblyr))/(zspace)
2402 542428017 : vel = (bgrid(nblyr+1)*dhbot - (bgrid(nblyr+1)-c1)*dhtop)/dt +darcyV/bphin_N(nblyr+1)
2403 542428017 : vel = vel/hbri_old
2404 542428017 : 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 542428017 : + p5*(-vel + iDin_phi(nblyr)/bphin_N(nblyr+1)*dphi_dx)
2407 :
2408 3796996119 : do k = 1, nblyr-1
2409 3254568102 : vel = (bgrid(k+1)*dhbot - (bgrid(k+1)-c1)*dhtop)/dt+darcyV/bphin_N(k+1)
2410 3254568102 : iDin_phi(k+1) = p5*(iDin(k+2)/iphin_N(k+2) + iDin(k+1)/iphin_N(k+1))
2411 3254568102 : dphi_dx = (iphin_N(k+1) - iphin_N(k))/(zspace)
2412 : K_spdiag(k)= p5*(vel/hbri_old - &
2413 3254568102 : iDin_phi(k)/(bphin_N(k+1))*dphi_dx)
2414 :
2415 3254568102 : vel = (bgrid(k+1)*dhbot - (bgrid(k+1)-c1)*dhtop)/dt +darcyV/bphin_N(k+1)
2416 3254568102 : dphi_dx = c0
2417 3254568102 : dphi_dx = kvectorn1(k)*(iphin_N(k+1) - iphin_N(k))/(zspace)
2418 : K_sbdiag(k+1)= -p5*(vel/hbri_old- &
2419 3254568102 : iDin_phi(k)/bphin_N(k+1)*dphi_dx)
2420 3254568102 : K_diag(k) = K_diag(1)*kvector1(k) + (K_spdiag(k) + K_sbdiag(k))*kvectorn1(k)
2421 :
2422 3254568102 : S_diag(k+1) = -(iDin_phi(k)+ iDin_phi(k+1))/zspace
2423 3254568102 : S_spdiag(k) = iDin_phi(k)/zspace
2424 3796996119 : S_sbdiag(k+1) = iDin_phi(k)/zspace
2425 : enddo
2426 :
2427 : ! k = nblyr
2428 :
2429 542428017 : vel = (bgrid(nblyr+1)*dhbot - (bgrid(nblyr+1)-c1)*dhtop)/dt+darcyV/bphin_N(nblyr+1)
2430 542428017 : dphi_dx = (iphin_N(nblyr+1) - iphin_N(nblyr))/(zspace)
2431 : K_spdiag(nblyr)= p5*(vel/hbri_old - &
2432 542428017 : iDin_phi(nblyr)/(bphin_N(nblyr+1))*dphi_dx)
2433 542428017 : vel = (bgrid(nblyr+1)*dhbot - (bgrid(nblyr+1)-c1)*dhtop)/dt +darcyV/bphin_N(nblyr+1)
2434 542428017 : dphi_dx = kvectorn1(nblyr)*(iphin_N(nblyr+1) - iphin_N(nblyr))/(zspace)
2435 : K_sbdiag(nblyr+1)= -p5*(vel/hbri_old- &
2436 542428017 : iDin_phi(nblyr)/bphin_N(nblyr+1)*dphi_dx)
2437 542428017 : K_diag(nblyr) = K_spdiag(nblyr) + K_sbdiag(nblyr)
2438 542428017 : S_diag(nblyr+1) = -iDin_phi(nblyr)/zspace
2439 542428017 : S_spdiag(nblyr) = iDin_phi(nblyr)/zspace
2440 542428017 : 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 4339424136 : do k = 1,nblyr
2445 3796996119 : D_spdiag(k) = max(-K_spdiag(k), c0, -K_sbdiag(k+1))
2446 4339424136 : D_sbdiag(k+1) = D_spdiag(k)
2447 : enddo
2448 4881852153 : do k = 1,nblyr+1
2449 4881852153 : 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 542428017 : vel2 = -(dhtop/hbri_old/dt +darcyV/bphin_N(1)/hbri_old)
2455 :
2456 4881852153 : Q_top(:) = c0
2457 542428017 : Q_top(1) = max(c0,vel2*C_top + atm_add/dt)
2458 542428017 : Qtop = Q_top(1)
2459 :
2460 542428017 : vel = (dhbot/hbri_old/dt +darcyV/hbri_old) ! going from iphin_N(nblyr+1) to c1 makes a difference
2461 :
2462 4881852153 : Q_bot(:) = c0
2463 : Q_bot(nblyr+1) = max(c0,vel*C_bot) + iDin_phi(nblyr+1)*C_bot&
2464 542428017 : /(zspace + grid_o/hbri_old)
2465 :
2466 542428017 : Qbot = Q_bot(nblyr+1)
2467 :
2468 542428017 : Sink_bot = K_diag(nblyr+1) + K_spdiag(nblyr)
2469 542428017 : Sink_top = K_diag(1) + K_sbdiag(2)
2470 :
2471 : !compute matrix elements (1 to nblyr+1)
2472 :
2473 4881852153 : spdiag = -dt * (D_spdiag + K_spdiag + S_spdiag)
2474 4881852153 : sbdiag = -dt * (D_sbdiag + K_sbdiag + S_sbdiag)
2475 4881852153 : diag = ML - dt * (D_diag + K_diag + S_diag)
2476 4881852153 : rhs = ML * C_in + dt * Q_top + dt* Q_bot
2477 :
2478 542428017 : 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 542428017 : subroutine compute_FCT_corr(C_in, C_low, dt, &
2488 542428017 : 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 ! LCOV_EXCL_LINE
2502 : ML , & ! Lumped mass diagonal matrix elements ! LCOV_EXCL_LINE
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 1084856034 : M_spdiag, M_sbdiag, & ! mass matrix ! LCOV_EXCL_LINE
2515 1084856034 : F_diag, F_spdiag, F_sbdiag, & ! anti-diffusive matrix ! LCOV_EXCL_LINE
2516 1084856034 : Pplus, Pminus, & ! ! LCOV_EXCL_LINE
2517 1084856034 : Qplus, Qminus, & ! ! LCOV_EXCL_LINE
2518 1084856034 : Rplus, Rminus, & ! ! LCOV_EXCL_LINE
2519 1084856034 : 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 542428017 : zspace = c1/real(nblyr,kind=dbl_kind)
2530 :
2531 : ! compute the mass matrix
2532 :
2533 4881852153 : M_spdiag(:) = zspace/c6
2534 542428017 : M_spdiag(nblyr+1) = c0
2535 4881852153 : M_sbdiag(:) = zspace/c6
2536 542428017 : M_sbdiag(1) = c0
2537 :
2538 : ! compute off matrix F:
2539 :
2540 4881852153 : F_diag(:) = c0
2541 4881852153 : F_spdiag(:) = c0
2542 4881852153 : F_sbdiag(:) = c0
2543 :
2544 4339424136 : 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 3796996119 : + 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 3796996119 : + D_sbdiag(k+1)*(C_low(k+1)-C_low(k))
2549 :
2550 3796996119 : if (F_spdiag(k)*(C_low(k) - C_low(k+1)) > c0) F_spdiag(k) = c0
2551 4339424136 : if (F_sbdiag(k+1)*(C_low(k+1) - C_low(k)) > c0) F_sbdiag(k+1) = c0
2552 : enddo
2553 :
2554 4881852153 : if (maxval(abs(F_spdiag)) > c0) then
2555 :
2556 : ! compute the weighting factors: a_spdiag, a_sbdiag
2557 :
2558 3144374280 : a_spdiag(:) = c0
2559 3144374280 : a_sbdiag(:) = c0
2560 :
2561 349374920 : Pplus(1) = max(c0, F_spdiag(1))
2562 349374920 : Pminus(1) = min(c0, F_spdiag(1))
2563 349374920 : Pplus(nblyr+1) = max(c0, F_sbdiag(nblyr+1))
2564 349374920 : Pminus(nblyr+1) = min(c0, F_sbdiag(nblyr+1))
2565 349374920 : Qplus(1) = max(c0,C_low(2)-C_low(1))
2566 349374920 : Qminus(1)= min(c0,C_low(2)-C_low(1))
2567 349374920 : Qplus(nblyr+1) = max(c0,C_low(nblyr)-C_low(nblyr+1))
2568 349374920 : Qminus(nblyr+1)= min(c0,C_low(nblyr)-C_low(nblyr+1))
2569 349374920 : Rplus(1) = min(c1, ML(1)*Qplus(1)/dt/(Pplus(1)+puny))
2570 349374920 : Rminus(1) = min(c1, ML(1)*Qminus(1)/dt/(Pminus(1)-puny))
2571 349374920 : Rplus(nblyr+1) = min(c1, ML(nblyr+1)*Qplus(nblyr+1)/dt/(Pplus(nblyr+1)+puny))
2572 349374920 : Rminus(nblyr+1) = min(c1, ML(nblyr+1)*Qminus(nblyr+1)/dt/(Pminus(nblyr+1)-puny))
2573 2445624440 : do k = 2,nblyr
2574 2096249520 : Pplus(k) = max(c0,F_spdiag(k)) + max(c0,F_sbdiag(k))
2575 2096249520 : Pminus(k) = min(c0,F_spdiag(k)) + min(c0,F_sbdiag(k))
2576 2096249520 : Qplus(k) = max(c0, max(C_low(k+1)-C_low(k),C_low(k-1)-C_low(k)))
2577 2096249520 : Qminus(k) = min(c0, min(C_low(k+1)-C_low(k),C_low(k-1)-C_low(k)))
2578 2096249520 : Rplus(k) = min(c1, ML(k)*Qplus(k)/dt/(Pplus(k)+puny))
2579 2445624440 : Rminus(k) = min(c1, ML(k)*Qminus(k)/dt/(Pminus(k)-puny))
2580 : enddo
2581 :
2582 2794999360 : do k = 1, nblyr
2583 2445624440 : a_spdiag(k) = min(Rminus(k),Rplus(k+1))
2584 2445624440 : if (F_spdiag(k) > c0) a_spdiag(k) = min(Rplus(k),Rminus(k+1))
2585 2445624440 : a_sbdiag(k+1) = min(Rminus(k+1),Rplus(k))
2586 2794999360 : 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 349374920 : F_diag(1) = a_spdiag(1)*F_spdiag(1)
2592 349374920 : F_diag(nblyr+1) = a_sbdiag(nblyr+1)* F_sbdiag(nblyr+1)
2593 349374920 : C_low(1) = C_low(1) + dt*F_diag(1)/ML(1)
2594 349374920 : C_low(nblyr+1) = C_low(nblyr+1) + dt*F_diag(nblyr+1)/ML(nblyr+1)
2595 :
2596 2445624440 : do k = 2,nblyr
2597 2096249520 : F_diag(k) = a_spdiag(k)*F_spdiag(k) + a_sbdiag(k)*F_sbdiag(k)
2598 2445624440 : C_low(k) = C_low(k) + dt*F_diag(k)/ML(k)
2599 : enddo
2600 :
2601 : endif !F_spdiag is nonzero
2602 :
2603 542428017 : 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 542428017 : subroutine tridiag_solverz(nmat, sbdiag, &
2613 542428017 : diag, spdiag, & ! LCOV_EXCL_LINE
2614 542428017 : 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 ! LCOV_EXCL_LINE
2621 : diag , & ! diagonal matrix elements ! LCOV_EXCL_LINE
2622 : spdiag , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
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 1084856034 : wgamma ! temporary matrix variable
2638 :
2639 : character(len=*),parameter :: subname='(tridiag_solverz)'
2640 :
2641 542428017 : wbeta = diag(1)
2642 542428017 : xout(1) = rhs(1) / wbeta
2643 :
2644 4339424136 : do k = 2, nmat
2645 3796996119 : wgamma(k) = spdiag(k-1) / wbeta
2646 3796996119 : wbeta = diag(k) - sbdiag(k)*wgamma(k)
2647 4339424136 : xout(k) = (rhs(k) - sbdiag(k)*xout(k-1)) / wbeta
2648 : enddo ! k
2649 :
2650 4339424136 : do k = nmat-1, 1, -1
2651 4339424136 : xout(k) = xout(k) - wgamma(k+1)*xout(k+1)
2652 : enddo ! k
2653 :
2654 542428017 : end subroutine tridiag_solverz
2655 :
2656 : !=======================================================================
2657 : !
2658 : ! authors Nicole Jeffery, LANL
2659 :
2660 542428017 : subroutine check_conservation_FCT (C_init, C_new, C_low, S_top, &
2661 : S_bot, L_bot, L_top, dt, & ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2673 : S_bot , & ! bottom flux into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
2674 : L_bot , & ! remaining bottom flux into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
2675 : L_top , & ! remaining top flux into ice (mmol/m^2/s) ! LCOV_EXCL_LINE
2676 : dt , & ! LCOV_EXCL_LINE
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 , & ! LCOV_EXCL_LINE
2689 : C_init_tot , & ! LCOV_EXCL_LINE
2690 : C_new_tot , & ! LCOV_EXCL_LINE
2691 : zspace , & !1/nblyr ! LCOV_EXCL_LINE
2692 : accuracyC , & ! centered difference is Order(zspace^2) ! LCOV_EXCL_LINE
2693 : var_tmp ! temporary variable
2694 :
2695 : character(len=*),parameter :: subname='(check_conservation_FCT)'
2696 :
2697 542428017 : zspace = c1/real(nblyr,kind=dbl_kind)
2698 :
2699 : !-------------------------------------
2700 : ! Ocean flux: positive into the ocean
2701 : !-------------------------------------
2702 542428017 : C_init_tot = (C_init(1) + C_init(nblyr+1))*zspace*p5
2703 542428017 : C_new_tot = (C_new(1) + C_new(nblyr+1))*zspace*p5
2704 542428017 : C_low(1) = C_new(1)
2705 542428017 : C_low(nblyr+1) = C_new(nblyr+1)
2706 :
2707 3796996119 : do k = 2, nblyr
2708 3254568102 : C_init_tot = C_init_tot + C_init(k)*zspace
2709 3254568102 : C_new_tot = C_new_tot + C_new(k)*zspace
2710 3796996119 : C_low(k) = C_new(k)
2711 : enddo
2712 :
2713 542428017 : accuracyC = 1.0e-11_dbl_kind*max(c1, C_init_tot, C_new_tot)
2714 542428017 : fluxbio = fluxbio + (C_init_tot - C_new_tot + source)/dt
2715 542428017 : 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 4881852153 : 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 542428017 : 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 542428017 : 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 562042208 : 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 ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
2815 : dzint, & ! snow interior depth (m) ! LCOV_EXCL_LINE
2816 : hslyr, & ! snow layer thickness (m) ! LCOV_EXCL_LINE
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 562042208 : hslyr = hsnow/real(nslyr,kind=dbl_kind)
2825 562042208 : dzssl = hslyr*p5
2826 562042208 : dzint = max(c0,hsnow - dzssl)
2827 562042208 : zspace = c1/real(nblyr,kind=dbl_kind)
2828 :
2829 562042208 : xout = c0
2830 562042208 : xout = (xin(1) + xin(nblyr+1))*hbrine*p5*zspace
2831 3934295456 : do n = 2, nblyr
2832 3934295456 : xout = xout + xin(n)*zspace*hbrine
2833 : enddo ! n
2834 562042208 : xout = xout + dzssl*xin(nblyr+2) + dzint*xin(nblyr+3)
2835 :
2836 562042208 : 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 67670680 : 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 135341360 : 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 609036120 : zspace(:) = c1/real(nblyr,kind=dbl_kind)
2867 67670680 : zspace(1) = p5*zspace(1)
2868 67670680 : zspace(nblyr+1) = zspace(1)
2869 :
2870 67670680 : xout = c0
2871 :
2872 67670680 : if (tr_bgc_N) then
2873 54180728 : iBioCount = c0
2874 216722912 : do m = 1, n_algae
2875 162542184 : nBGC = nt_bgc_N(1)
2876 1462879656 : do n = 1, nblyr+1
2877 1300337472 : iLayer = iBioCount + n-1
2878 1462879656 : xout = xout + xin(nBGC+iLayer)*zspace(n)*hbrine*R_C2N(m)
2879 : enddo
2880 216722912 : iBioCount = iBioCount + nblyr+3
2881 : enddo
2882 : endif
2883 67670680 : if (tr_bgc_C) then
2884 54180728 : iBioCount = c0
2885 54180728 : nBGC = nt_bgc_DOC(1)
2886 162542184 : do m = 1, n_doc
2887 975253104 : do n = 1, nblyr+1
2888 866891648 : iLayer = iBioCount + n-1
2889 975253104 : xout = xout + xin(nBGC+iLayer)*zspace(n)*hbrine
2890 : enddo
2891 162542184 : iBioCount = iBioCount + nblyr+3
2892 : enddo
2893 108361456 : do m = 1, n_dic
2894 487626552 : do n = 1, nblyr+1
2895 433445824 : iLayer = iBioCount + n-1
2896 487626552 : xout = xout + xin(nBGC+iLayer)*zspace(n)*hbrine
2897 : enddo
2898 108361456 : iBioCount = iBioCount + nblyr+3
2899 : enddo
2900 : endif
2901 :
2902 67670680 : if (tr_bgc_DON) then
2903 54180728 : iBioCount = c0
2904 108361456 : do m = 1, n_don
2905 54180728 : nBGC = nt_bgc_DON(1)
2906 487626552 : do n = 1, nblyr+1
2907 433445824 : iLayer = iBioCount + n-1
2908 487626552 : xout = xout + xin(nBGC+iLayer)*zspace(n)*hbrine*R_C2N_DON(m)
2909 : enddo
2910 108361456 : iBioCount = iBioCount + nblyr+3
2911 : enddo
2912 : endif
2913 67670680 : if (tr_bgc_hum) then
2914 54180728 : nBGC = nt_bgc_hum
2915 487626552 : do n = 1, nblyr+1
2916 433445824 : iLayer = n-1
2917 487626552 : xout = xout + xin(nBGC+iLayer)*zspace(n)*hbrine
2918 : enddo
2919 : endif
2920 :
2921 67670680 : 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 33835340 : 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 ! LCOV_EXCL_LINE
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 33835340 : Tot_Carbon_flux = c0
2946 :
2947 33835340 : if (tr_bgc_N) then
2948 108361456 : do m = 1, n_algae
2949 108361456 : 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 33835340 : if (tr_bgc_C) then
2953 81271092 : do m = 1, n_doc
2954 81271092 : Tot_Carbon_flux = Tot_Carbon_flux - flux_bio_atm(nlt_bgc_DOC(m)) + flux_bion(nlt_bgc_DOC(m))
2955 : enddo
2956 54180728 : do m = 1, n_dic
2957 54180728 : Tot_Carbon_flux = Tot_Carbon_flux - flux_bio_atm(nlt_bgc_DIC(m)) + flux_bion(nlt_bgc_DIC(m))
2958 : enddo
2959 : endif
2960 33835340 : if (tr_bgc_DON) then
2961 54180728 : do m = 1, n_don
2962 54180728 : 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 33835340 : if (tr_bgc_hum) &
2966 27090364 : Tot_Carbon_flux = Tot_Carbon_flux - flux_bio_atm(nlt_bgc_hum) + flux_bion(nlt_bgc_hum)
2967 :
2968 33835340 : end subroutine bgc_carbon_flux
2969 :
2970 : !=======================================================================
2971 :
2972 : end module icepack_algae
2973 :
2974 : !=======================================================================
|
