/[MITgcm]/MITgcm_contrib/bling/pkg/bling_remineralization.F
ViewVC logotype

Contents of /MITgcm_contrib/bling/pkg/bling_remineralization.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (show annotations) (download)
Sun May 15 00:30:35 2016 UTC (9 years, 2 months ago) by mmazloff
Branch: MAIN
Changes since 1.2: +16 -28 lines
Cleaning up TAF stores

1 C $Header: /u/gcmpack/MITgcm_contrib/bling/pkg/bling_remineralization.F,v 1.2 2016/02/28 21:49:24 mmazloff Exp $
2 C $Name: $
3
4 #include "BLING_OPTIONS.h"
5
6 CBOP
7 subroutine BLING_REMIN(
8 I PTR_NO3, PTR_FE, PTR_O2, irr_inst,
9 I N_spm, P_spm, Fe_spm, CaCO3_uptake,
10 O N_reminp, P_reminp, Fe_reminsum,
11 O N_den_benthic, CaCO3_diss,
12 I bi, bj, imin, imax, jmin, jmax,
13 I myIter, myTime, myThid )
14
15 C =================================================================
16 C | subroutine bling_remin
17 C | o Organic matter export and remineralization
18 C | - Sinking particulate flux and diel migration contribute to
19 C | export.
20 C | - Denitrification xxx
21 C | o Sediments
22 C =================================================================
23
24 implicit none
25
26 C === Global variables ===
27
28 #include "SIZE.h"
29 #include "DYNVARS.h"
30 #include "EEPARAMS.h"
31 #include "PARAMS.h"
32 #include "GRID.h"
33 #include "BLING_VARS.h"
34 #include "PTRACERS_SIZE.h"
35 #include "PTRACERS_PARAMS.h"
36 #ifdef ALLOW_AUTODIFF
37 # include "tamc.h"
38 #endif
39
40 C === Routine arguments ===
41 C bi,bj :: tile indices
42 C iMin,iMax :: computation domain: 1rst index range
43 C jMin,jMax :: computation domain: 2nd index range
44 C myTime :: current time
45 C myIter :: current timestep
46 C myThid :: thread Id. number
47 INTEGER bi, bj, imin, imax, jmin, jmax
48 _RL myTime
49 INTEGER myIter
50 INTEGER myThid
51 C === Input ===
52 C PTR_NO3 :: nitrate concentration
53 C PTR_FE :: iron concentration
54 C PTR_O2 :: oxygen concentration
55 _RL PTR_NO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
56 _RL PTR_FE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
57 _RL PTR_O2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
58 _RL irr_inst(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
59 _RL N_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
60 _RL P_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
61 _RL Fe_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
62 _RL CaCO3_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
63 C === Output ===
64 C
65 _RL N_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
66 _RL N_den_benthic(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
67 _RL P_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
68 _RL Fe_reminsum(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
69 _RL CaCO3_diss(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
70
71 #ifdef ALLOW_BLING
72 C === Local variables ===
73 C i,j,k :: loop indices
74 C irr_eff :: effective irradiance
75 C NO3_lim :: nitrate limitation
76 C PO4_lim :: phosphate limitation
77 C Fe_lim :: iron limitation for phytoplankton
78 C Fe_lim_diaz :: iron limitation for diazotrophs
79 C alpha_Fe :: initial slope of the P-I curve
80 C theta_Fe :: Chl:C ratio
81 C theta_Fe_max :: Fe-replete maximum Chl:C ratio
82 C irrk :: nut-limited efficiency of algal photosystems
83 C Pc_m :: light-saturated max photosynthesis rate for phyt
84 C Pc_m_diaz :: light-saturated max photosynthesis rate for diaz
85 C Pc_tot :: carbon-specific photosynthesis rate
86 C expkT :: temperature function
87 C mu :: net carbon-specific growth rate for phyt
88 C mu_diaz :: net carbon-specific growth rate for diaz
89 C biomass_sm :: nitrogen concentration in small phyto biomass
90 C biomass_lg :: nitrogen concentration in large phyto biomass
91 C N_uptake :: nitrogen uptake
92 C N_fix :: nitrogen fixation
93 C P_uptake :: phosphorus uptake
94 C POC_flux :: carbon export flux 3d field
95 C chl :: chlorophyll diagnostic
96 C PtoN :: variable ratio of phosphorus to nitrogen
97 C FetoN :: variable ratio of iron to nitrogen
98 C N_spm :: particulate sinking of nitrogen
99 C P_spm :: particulate sinking of phosphorus
100 C Fe_spm :: particulate sinking of iron
101 C N_dvm :: vertical transport of nitrogen by DVM
102 C P_dvm :: vertical transport of phosphorus by DVM
103 C Fe_dvm :: vertical transport of iron by DVM
104 C N_recycle :: recycling of newly-produced organic nitrogen
105 C P_recycle :: recycling of newly-produced organic phosphorus
106 C Fe_recycle :: recycling of newly-produced organic iron
107 c xxx to be completed
108 INTEGER i,j,k
109 _RL PONflux_u
110 _RL POPflux_u
111 _RL PFEflux_u
112 _RL CaCO3flux_u
113 _RL PONflux_l
114 _RL POPflux_l
115 _RL PFEflux_l
116 _RL CaCO3flux_l
117 _RL depth_l
118 _RL zremin
119 _RL zremin_caco3
120 _RL wsink
121 _RL POC_sed
122 _RL Fe_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
123 _RL NO3_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124 _RL PO4_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125 _RL O2_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126 _RL lig_stability
127 _RL FreeFe
128 _RL Fe_ads_inorg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
129 _RL Fe_ads_org
130 _RL log_btm_flx
131 _RL Fe_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
132 _RL o2_upper
133 _RL o2_lower
134 _RL dz_upper
135 _RL dz_lower
136 _RL temp_upper
137 _RL temp_lower
138 _RL z_dvm_regr
139 _RL frac_migr
140 _RL fdvm_migr
141 _RL fdvm_stat
142 _RL fdvmn_vint
143 _RL fdvmp_vint
144 _RL fdvmfe_vint
145 _RL z_dvm
146 _RL N_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
147 _RL P_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
148 _RL Fe_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
149 _RL dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
150 _RL mld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151 _RL Fe_burial(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152 _RL x_erfcc,z_erfcc,t_erfcc,erfcc
153 cxx order
154 CEOP
155
156 c ---------------------------------------------------------------------
157 c Initialize output and diagnostics
158
159 DO k=1,Nr
160 DO j=jmin,jmax
161 DO i=imin,imax
162 Fe_ads_inorg(i,j,k) = 0. _d 0
163 N_reminp(i,j,k) = 0. _d 0
164 P_reminp(i,j,k) = 0. _d 0
165 Fe_reminp(i,j,k) = 0. _d 0
166 Fe_reminsum(i,j,k) = 0. _d 0
167 N_remindvm(i,j,k) = 0. _d 0
168 P_remindvm(i,j,k) = 0. _d 0
169 Fe_remindvm(i,j,k) = 0. _d 0
170 N_den_benthic(i,j,k)= 0. _d 0
171 CaCO3_diss(i,j,k) = 0. _d 0
172 ENDDO
173 ENDDO
174 ENDDO
175 DO j=jmin,jmax
176 DO i=imin,imax
177 Fe_burial(i,j) = 0. _d 0
178 NO3_sed(i,j) = 0. _d 0
179 PO4_sed(i,j) = 0. _d 0
180 O2_sed(i,j) = 0. _d 0
181 ENDDO
182 ENDDO
183
184 c ---------------------------------------------------------------------
185 c Remineralization
186
187 C$TAF LOOP = parallel
188 DO j=jmin,jmax
189 C$TAF LOOP = parallel
190 DO i=imin,imax
191
192 C Initialize upper flux
193 PONflux_u = 0. _d 0
194 POPflux_u = 0. _d 0
195 PFEflux_u = 0. _d 0
196 CaCO3flux_u = 0. _d 0
197
198 c C$TAF init remin_stuff = static, Nr
199
200 DO k=1,Nr
201
202 Fe_ads_org = 0. _d 0
203 IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN
204
205 C Sinking speed is evaluated at the bottom of the cell
206 depth_l=-rF(k+1)
207 IF (depth_l .LE. wsink0z) THEN
208 wsink = wsink0
209 ELSE
210 wsink = wsinkacc * (depth_l - wsink0z) + wsink0
211 ENDIF
212
213 C Nutrient remineralization lengthscale
214 C Not an e-folding scale: this term increases with remineralization.
215 zremin = gamma_POM * ( PTR_O2(i,j,k)**2 /
216 & (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min)
217 & + remin_min )/(wsink + epsln)
218
219 C Calcium remineralization relaxed toward the inverse of the
220 C ca_remin_depth constant value as the calcite saturation approaches 0.
221 zremin_caco3 = 1. _d 0/ca_remin_depth*(1. _d 0 - min(1. _d 0,
222 & omegaC(i,j,k,bi,bj) + epsln ))
223
224
225 C POM flux leaving the cell
226 PONflux_l = (PONflux_u+N_spm(i,j,k)*drF(k)
227 & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
228 & *hFacC(i,j,k,bi,bj))
229 C!! multiply by intercept_frac ???
230
231 POPflux_l = (POPflux_u+P_spm(i,j,k)*drF(k)
232 & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
233 & *hFacC(i,j,k,bi,bj))
234 C!! multiply by intercept_frac ???
235
236 C CaCO3 flux leaving the cell
237 CaCO3flux_l = (caco3flux_u+CaCO3_uptake(i,j,k)*drF(k)
238 & *hFacC(i,j,k,bi,bj))/(1+zremin_caco3*drF(k)
239 & *hFacC(i,j,k,bi,bj))
240 C!! multiply by intercept_frac ???
241
242
243 C Start with cells that are not the deepest cells
244 IF ((k.LT.Nr) .AND. (hFacC(i,j,k+1,bi,bj).GT.0)) THEN
245
246 C Nutrient accumulation in a cell is given by the biological production
247 C (and instant remineralization) of particulate organic matter
248 C plus flux thought upper interface minus flux through lower interface.
249 C (Since not deepest cell: hFacC=1)
250 N_reminp(i,j,k) = (PONflux_u + N_spm(i,j,k)*drF(k)
251 & - PONflux_l)*recip_drF(k)
252
253 P_reminp(i,j,k) = (POPflux_u + P_spm(i,j,k)*drF(k)
254 & - POPflux_l)*recip_drF(k)
255
256 CaCO3_diss(i,j,k) = (CaCO3flux_u + CaCO3_uptake(i,j,k)
257 & *drF(k) - CaCO3flux_l)*recip_drF(k)
258
259 Fe_sed(i,j,k) = 0. _d 0
260
261
262 ELSE
263 C If this layer is adjacent to bottom topography or it is the deepest
264 C cell of the domain, then remineralize/dissolve in this grid cell
265 C i.e. don't subtract off lower boundary fluxes when calculating remin
266
267 N_reminp(i,j,k) = PONflux_u*recip_drF(k)
268 & *recip_hFacC(i,j,k,bi,bj) + N_spm(i,j,k)
269
270 P_reminp(i,j,k) = POPflux_u*recip_drF(k)
271 & *recip_hFacC(i,j,k,bi,bj) + P_spm(i,j,k)
272
273 CaCO3_diss(i,j,k) = CaCO3flux_u*recip_drF(k)
274 & *recip_hFacC(i,j,k,bi,bj) + CaCO3_uptake(i,j,k)
275
276
277 c Efflux Fed out of sediments
278 C The phosphate flux hitting the bottom boundary
279 C is used to scale the return of iron to the water column.
280 C Maximum value added for numerical stability.
281
282 POC_sed = PONflux_l * CtoN
283
284 Fe_sed(i,j,k) = min(1. _d -11,
285 & max(epsln, FetoC_sed * POC_sed * recip_drF(k)
286 & *recip_hFacC(i,j,k,bi,bj)))
287
288 #ifdef BLING_ADJOINT_SAFE
289 Fe_sed(i,j,k) = 0. _d 0
290 #endif
291
292
293 cav temporary until I figure out why this is problematic for adjoint
294 #ifndef BLING_ADJOINT_SAFE
295
296 #ifndef USE_SGS_SED
297 c Calculate benthic denitrification and Fe efflux here, if the subgridscale sediment module is not being used.
298
299 IF (POC_sed .gt. 0. _d 0) THEN
300
301 log_btm_flx = 0. _d 0
302
303 c Convert from mol N m-2 s-1 to umol C cm-2 d-1 and take the log
304
305 log_btm_flx = log10(min(43.0 _d 0, POC_sed *
306 & 86400. _d 0 * 100.0 _d 0))
307
308 c Metamodel gives units of umol C cm-2 d-1, convert to mol N m-2 s-1 and multiply by
309 c no3_2_n to give NO3 consumption rate
310
311 N_den_benthic(i,j,k) = min (POC_sed * NO3toN / CtoN,
312 & (10 _d 0)**(-0.9543 _d 0 + 0.7662 _d 0 *
313 & log_btm_flx - 0.235 _d 0 * log_btm_flx * log_btm_flx)
314 & / (CtoN * 86400. _d 0 * 100.0 _d 0) * NO3toN *
315 & PTR_NO3(i,j,k) / (k_no3 + PTR_NO3(i,j,k)) ) *
316 & recip_drF(k)
317
318 endif
319
320 #endif
321
322 #endif
323
324 c ---------------------------------------------------------------------
325 c Calculate external bottom fluxes for tracer_vertdiff. Positive fluxes are into the water
326 c column from the seafloor. For P, the bottom flux puts the sinking flux reaching the bottom
327 c cell into the water column through diffusion. For iron, the sinking flux disappears into the
328 c sediments if bottom waters are oxic (assumed adsorbed as oxides). If bottom waters are anoxic,
329 c the sinking flux of Fe is returned to the water column.
330 c
331 c For oxygen, the consumption of oxidant required to respire
332 c the settling flux of organic matter (in support of the
333 c no3 bottom flux) diffuses from the bottom water into the sediment.
334
335 c Assume all NO3 for benthic denitrification is supplied from the bottom water, and that
336 c all organic N is also consumed under denitrification (Complete Denitrification, sensu
337 c Paulmier, Biogeosciences 2009). Therefore, no NO3 is regenerated from organic matter
338 c respired by benthic denitrification (necessitating the second term in b_no3).
339
340 NO3_sed(i,j) = PONflux_l*drF(k)*hFacC(i,j,k,bi,bj)
341 & - N_den_benthic(i,j,k) / NO3toN
342
343 PO4_sed(i,j) = POPflux_l*drF(k)*hFacC(i,j,k,bi,bj)
344
345 c Oxygen flux into sediments is that required to support non-denitrification respiration,
346 c assuming a 4/5 oxidant ratio of O2 to NO3. Oxygen consumption is allowed to continue
347 c at negative oxygen concentrations, representing sulphate reduction.
348
349 O2_sed(i,j) = -(O2toN * PONflux_l*drF(k)*hFacC(i,j,k,bi,bj)
350 & - N_den_benthic(i,j,k)* 1.25)
351
352 ENDIF
353
354
355 C Begin iron uptake calculations by determining ligand bound and free iron.
356 C Both forms are available for biology, but only free iron is scavenged
357 C onto particles and forms colloids.
358
359 lig_stability = kFe_eq_lig_max-(KFe_eq_lig_max-kFe_eq_lig_min)
360 & *(irr_inst(i,j,k)**2
361 & /(kFe_eq_lig_irr**2+irr_inst(i,j,k)**2))
362 & *max(epsln,min(1. _d 0,(PTR_FE(i,j,k)
363 & -kFe_eq_lig_Femin)/
364 & (PTR_FE(i,j,k)+epsln)*1.2 _d 0))
365
366 C Use the quadratic equation to solve for binding between iron and ligands
367
368 FreeFe = (-(1+lig_stability*(ligand-PTR_FE(i,j,k)))
369 & +((1+lig_stability*(ligand-PTR_FE(i,j,k)))**2+4*
370 & lig_stability*PTR_FE(i,j,k))**(0.5))/(2*
371 & lig_stability)
372
373 C Iron scavenging doesn't occur in anoxic water (Fe2+ is soluble), so set
374 C FreeFe = 0 when anoxic. FreeFe should be interpreted the free iron that
375 C participates in scavenging.
376
377 #ifndef BLING_ADJOINT_SAFE
378 IF (PTR_O2(i,j,k) .LT. oxic_min) THEN
379 FreeFe = 0. _d 0
380 ENDIF
381 #endif
382
383 c Two mechanisms for iron uptake, in addition to biological production:
384 c colloidal scavenging and scavenging by organic matter.
385
386 c Colloidal scavenging:
387 c Minimum function for numerical stability
388 c Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+
389 c & min(0.5/PTRACERS_dTLev(1), kFe_inorg*FreeFe**(0.5))*FreeFe
390
391 Fe_ads_inorg(i,j,k) =
392 & kFe_inorg*(max(1. _d -8,FreeFe))**(1.5)
393
394 C Scavenging of iron by organic matter:
395 c The POM value used is the bottom boundary flux. This doesn't occur in
396 c oxic waters, but FreeFe is set to 0 in such waters earlier.
397 IF ( PONflux_l .GT. 0. _d 0 ) THEN
398
399 c Minimum function for numerical stability
400 c Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+
401 c & min(0.5/PTRACERS_dTLev(1), kFE_org*(POMflux_l
402 c & *CtoP/NUTfac*12.01/wsink)**(0.58)*FreeFe
403
404 #ifndef BLING_ADJOINT_SAFE
405 Fe_ads_org =
406 & kFE_org*(PONflux_l/(epsln + wsink)
407 & * MasstoN)**(0.58)*FreeFe
408 #else
409 Fe_ads_org =
410 & kFE_org*(PONflux_l/(epsln + wsink0)
411 & * MasstoN)**(0.58)*FreeFe
412 #endif
413 ENDIF
414
415
416
417 C If water is oxic then the iron is remineralized normally. Otherwise
418 C it is completely remineralized (fe 2+ is soluble, but unstable
419 C in oxidizing environments).
420
421 PFEflux_l = (PFEflux_u+(Fe_spm(i,j,k)+Fe_ads_inorg(i,j,k)
422 & +Fe_ads_org)*drF(k)
423 & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
424 & *hFacC(i,j,k,bi,bj))
425
426
427 c Added the burial flux of sinking particulate iron here as a
428 c diagnostic, needed to calculate mass balance of iron.
429 c this is calculated last for the deepest cell
430
431 Fe_burial(i,j) = PFeflux_l
432
433
434 #ifndef BLING_ADJOINT_SAFE
435 IF ( PTR_O2(i,j,k) .LT. oxic_min ) THEN
436 pfeflux_l = 0
437 ENDIF
438 #endif
439
440 Fe_reminp(i,j,k) = (pfeflux_u+(Fe_spm(i,j,k)
441 & +Fe_ads_inorg(i,j,k)
442 & +Fe_ads_org)*drF(k)
443 & *hFacC(i,j,k,bi,bj)-pfeflux_l)*recip_drF(k)
444 & *recip_hFacC(i,j,k,bi,bj)
445 C!! there's an intercept_frac here... need to add
446
447
448 C Prepare the tracers for the next layer down
449 PONflux_u = PONflux_l
450 POPflux_u = POPflux_l
451 PFEflux_u = PFEflux_l
452 CaCO3flux_u = CaCO3flux_l
453
454 c
455 Fe_reminsum(i,j,k) = Fe_reminp(i,j,k) + Fe_sed(i,j,k)
456 & - Fe_ads_org - Fe_ads_inorg(i,j,k)
457 cc Fe_reminsum(i,j,k) = 0. _d 0
458
459 ENDIF
460
461 Fe_ads_org = 0. _d 0
462
463 ENDDO
464 ENDDO
465 ENDDO
466
467 c ---------------------------------------------------------------------
468
469 #ifdef ALLOW_DIAGNOSTICS
470 IF ( useDiagnostics ) THEN
471
472 c 3d local variables
473 c CALL DIAGNOSTICS_FILL(POC_flux,'BLGPOCF ',0,Nr,2,bi,bj,myThid)
474 CALL DIAGNOSTICS_FILL(Fe_ads_inorg,'BLGFEAI',0,Nr,2,bi,bj,
475 & myThid)
476 CALL DIAGNOSTICS_FILL(Fe_sed,'BLGFESED',0,Nr,2,bi,bj,myThid)
477 CALL DIAGNOSTICS_FILL(Fe_reminp,'BLGFEREM',0,Nr,2,bi,bj,myThid)
478 CALL DIAGNOSTICS_FILL(N_den_benthic,'BLGNDENB',0,Nr,2,bi,bj,
479 & myThid)
480 c CALL DIAGNOSTICS_FILL(N_den_pelag,'BLGNDENP',0,Nr,2,bi,bj,myThid)
481 CALL DIAGNOSTICS_FILL(N_reminp,'BLGNREM ',0,Nr,2,bi,bj,myThid)
482 CALL DIAGNOSTICS_FILL(P_reminp,'BLGPREM ',0,Nr,2,bi,bj,myThid)
483 c 2d local variables
484 CALL DIAGNOSTICS_FILL(Fe_burial,'BLGFEBUR',0,1,2,bi,bj,myThid)
485 CALL DIAGNOSTICS_FILL(NO3_sed,'BLGNSED ',0,1,2,bi,bj,myThid)
486 CALL DIAGNOSTICS_FILL(PO4_sed,'BLGPSED ',0,1,2,bi,bj,myThid)
487 CALL DIAGNOSTICS_FILL(O2_sed,'BLGO2SED',0,1,2,bi,bj,myThid)
488 c these variables are currently 1d, could be 3d for diagnostics
489 c (or diag_fill could be called inside loop - which is faster?)
490 c CALL DIAGNOSTICS_FILL(zremin,'BLGZREM ',0,Nr,2,bi,bj,myThid)
491
492 ENDIF
493 #endif /* ALLOW_DIAGNOSTICS */
494
495
496
497 #endif /* ALLOW_BLING */
498
499 RETURN
500 END

  ViewVC Help
Powered by ViewVC 1.1.22