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

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

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


Revision 1.11 - (show annotations) (download)
Wed May 24 17:21:15 2017 UTC (6 years, 11 months ago) by mmazloff
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, HEAD
Changes since 1.10: +3 -3 lines
Add additional DVM-remin CPP flag. Some diagnostics cleaning.

1 C $Header: /u/gcmpack/MITgcm/pkg/bling/bling_remineralization.F,v 1.10 2017/03/29 15:51:19 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 | - Benthic denitrification
21 C | - Iron source from sediments
22 C | - Iron scavenging
23 C =================================================================
24
25 implicit none
26
27 C === Global variables ===
28
29 #include "SIZE.h"
30 #include "DYNVARS.h"
31 #include "EEPARAMS.h"
32 #include "PARAMS.h"
33 #include "GRID.h"
34 #include "BLING_VARS.h"
35 #include "PTRACERS_SIZE.h"
36 #include "PTRACERS_PARAMS.h"
37 #ifdef ALLOW_AUTODIFF
38 # include "tamc.h"
39 #endif
40
41 C === Routine arguments ===
42 C bi,bj :: tile indices
43 C iMin,iMax :: computation domain: 1rst index range
44 C jMin,jMax :: computation domain: 2nd index range
45 C myTime :: current time
46 C myIter :: current timestep
47 C myThid :: thread Id. number
48 INTEGER bi, bj, imin, imax, jmin, jmax
49 _RL myTime
50 INTEGER myIter
51 INTEGER myThid
52 C === Input ===
53 C PTR_NO3 :: nitrate concentration
54 C PTR_FE :: iron concentration
55 C PTR_O2 :: oxygen concentration
56 _RL PTR_NO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
57 _RL PTR_FE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
58 _RL PTR_O2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
59 _RL irr_inst(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
60 _RL N_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
61 _RL P_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
62 _RL Fe_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
63 _RL CaCO3_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
64
65 C === Output ===
66 C N_reminp :: remineralization of particulate organic nitrogen
67 C N_den_benthic :: Benthic denitrification
68 C P_reminp :: remineralization of particulate organic nitrogen
69 C Fe_reminsum :: iron remineralization and adsorption
70 C CaCO3_diss :: Calcium carbonate dissolution
71 _RL N_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
72 _RL N_den_benthic(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73 _RL P_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
74 _RL Fe_reminsum(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
75 _RL CaCO3_diss(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
76
77 #ifdef ALLOW_BLING
78 C === Local variables ===
79 C i,j,k :: loop indices
80
81 INTEGER i,j,k
82 INTEGER bttmlyr
83 _RL PONflux_u
84 _RL PONflux_l
85 _RL POPflux_u
86 _RL POPflux_l
87 _RL PFEflux_u
88 _RL PFEflux_l
89 _RL CaCO3flux_u
90 _RL CaCO3flux_l
91 _RL depth_l
92 _RL zremin
93 _RL zremin_caco3
94 _RL wsink
95 _RL POC_sed
96 _RL Fe_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97 _RL NO3_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 _RL PO4_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99 _RL O2_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 _RL lig_stability
101 _RL FreeFe
102 _RL Fe_ads_inorg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103 _RL Fe_ads_org(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104 _RL log_btm_flx
105 _RL Fe_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106 _RL Fe_burial(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107 CEOP
108
109 C ---------------------------------------------------------------------
110 C Initialize output and diagnostics
111
112 DO k=1,Nr
113 DO j=jmin,jmax
114 DO i=imin,imax
115 Fe_ads_org(i,j,k) = 0. _d 0
116 Fe_ads_inorg(i,j,k) = 0. _d 0
117 N_reminp(i,j,k) = 0. _d 0
118 P_reminp(i,j,k) = 0. _d 0
119 Fe_reminp(i,j,k) = 0. _d 0
120 Fe_reminsum(i,j,k) = 0. _d 0
121 N_den_benthic(i,j,k)= 0. _d 0
122 CaCO3_diss(i,j,k) = 0. _d 0
123 ENDDO
124 ENDDO
125 ENDDO
126 DO j=jmin,jmax
127 DO i=imin,imax
128 Fe_burial(i,j) = 0. _d 0
129 NO3_sed(i,j) = 0. _d 0
130 PO4_sed(i,j) = 0. _d 0
131 O2_sed(i,j) = 0. _d 0
132 ENDDO
133 ENDDO
134
135 C ---------------------------------------------------------------------
136 C Remineralization
137
138 C$TAF LOOP = parallel
139 DO j=jmin,jmax
140 C$TAF LOOP = parallel
141 DO i=imin,imax
142
143 C Initialize upper flux
144 PONflux_u = 0. _d 0
145 POPflux_u = 0. _d 0
146 PFEflux_u = 0. _d 0
147 CaCO3flux_u = 0. _d 0
148
149 DO k=1,Nr
150
151 C Initialization here helps taf
152 Fe_ads_org(i,j,k) = 0. _d 0
153
154 C ARE WE ON THE BOTTOM
155 bttmlyr = 1
156 IF (k.LT.Nr) THEN
157 IF (hFacC(i,j,k+1,bi,bj).GT.0) bttmlyr = 0
158 C we are not yet at the bottom
159 ENDIF
160
161 IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN
162
163 C Sinking speed is evaluated at the bottom of the cell
164 depth_l=-rF(k+1)
165 IF (depth_l .LE. wsink0z) THEN
166 wsink = wsink0_2d(i,j,bi,bj)
167 ELSE
168 wsink = wsinkacc * (depth_l - wsink0z) + wsink0_2d(i,j,bi,bj)
169 ENDIF
170
171 C Nutrient remineralization lengthscale
172 C Not an e-folding scale: this term increases with remineralization.
173 zremin = gamma_POM2d(i,j,bi,bj) * ( PTR_O2(i,j,k)**2 /
174 & (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min)
175 & + remin_min )/(wsink + epsln)
176
177 C Calcium remineralization relaxed toward the inverse of the
178 C ca_remin_depth constant value as the calcite saturation approaches 0.
179 zremin_caco3 = 1. _d 0/ca_remin_depth*(1. _d 0 - min(1. _d 0,
180 & omegaC(i,j,k,bi,bj) + epsln ))
181
182 C POM flux leaving the cell
183 PONflux_l = (PONflux_u+N_spm(i,j,k)*drF(k)
184 & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
185 & *hFacC(i,j,k,bi,bj))
186
187 POPflux_l = (POPflux_u+P_spm(i,j,k)*drF(k)
188 & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
189 & *hFacC(i,j,k,bi,bj))
190
191 C CaCO3 flux leaving the cell
192 CaCO3flux_l = (caco3flux_u+CaCO3_uptake(i,j,k)*drF(k)
193 & *hFacC(i,j,k,bi,bj))/(1+zremin_caco3*drF(k)
194 & *hFacC(i,j,k,bi,bj))
195
196 C Start with cells that are not the deepest cells
197 IF (bttmlyr.EQ.0) THEN
198 C Nutrient accumulation in a cell is given by the biological production
199 C (and instant remineralization) of particulate organic matter
200 C plus flux thought upper interface minus flux through lower interface.
201 C (Since not deepest cell: hFacC=1)
202 N_reminp(i,j,k) = (PONflux_u + N_spm(i,j,k)*drF(k)
203 & - PONflux_l)*recip_drF(k)
204
205 P_reminp(i,j,k) = (POPflux_u + P_spm(i,j,k)*drF(k)
206 & - POPflux_l)*recip_drF(k)
207
208 CaCO3_diss(i,j,k) = (CaCO3flux_u + CaCO3_uptake(i,j,k)
209 & *drF(k) - CaCO3flux_l)*recip_drF(k)
210
211 Fe_sed(i,j,k) = 0. _d 0
212 C NOW DO BOTTOM LAYER
213 ELSE
214 C If this layer is adjacent to bottom topography or it is the deepest
215 C cell of the domain, then remineralize/dissolve in this grid cell
216 C i.e. do not subtract off lower boundary fluxes when calculating remin
217
218 N_reminp(i,j,k) = PONflux_u*recip_drF(k)
219 & *recip_hFacC(i,j,k,bi,bj) + N_spm(i,j,k)
220
221 P_reminp(i,j,k) = POPflux_u*recip_drF(k)
222 & *recip_hFacC(i,j,k,bi,bj) + P_spm(i,j,k)
223
224 CaCO3_diss(i,j,k) = CaCO3flux_u*recip_drF(k)
225 & *recip_hFacC(i,j,k,bi,bj) + CaCO3_uptake(i,j,k)
226
227 C Efflux Fed out of sediments
228 C The phosphate flux hitting the bottom boundary
229 C is used to scale the return of iron to the water column.
230 C Maximum value added for numerical stability.
231
232 POC_sed = PONflux_l * CtoN
233
234 Fe_sed(i,j,k) = max(epsln, FetoC_sed * POC_sed * recip_drF(k)
235 & *recip_hFacC(i,j,k,bi,bj))
236
237 cav log_btm_flx = 0. _d 0
238 log_btm_flx = 1. _d -20
239
240 CMM: this is causing instability in the adjoint. Needs debugging
241 #ifndef BLING_ADJOINT_SAFE
242 IF (POC_sed .gt. 0. _d 0) THEN
243
244 C Convert from mol N m-2 s-1 to umol C cm-2 d-1 and take the log
245
246 log_btm_flx = log10(min(43.0 _d 0, POC_sed *
247 & 86400. _d 0 * 100.0 _d 0))
248
249 C Metamodel gives units of umol C cm-2 d-1, convert to mol N m-2 s-1 and
250 C multiply by no3_2_n to give NO3 consumption rate
251
252 N_den_benthic(i,j,k) = min (POC_sed * NO3toN / CtoN,
253 & (10 _d 0)**(-0.9543 _d 0 + 0.7662 _d 0 *
254 & log_btm_flx - 0.235 _d 0 * log_btm_flx * log_btm_flx)
255 & / (CtoN * 86400. _d 0 * 100.0 _d 0) * NO3toN *
256 & PTR_NO3(i,j,k) / (k_no3 + PTR_NO3(i,j,k)) ) *
257 & recip_drF(k)
258
259 ENDIF
260 #endif
261
262 C ---------------------------------------------------------------------
263 C Calculate external bottom fluxes for tracer_vertdiff. Positive fluxes
264 C are into the water column from the seafloor. For P, the bottom flux puts
265 C the sinking flux reaching the bottom cell into the water column through
266 C diffusion. For iron, the sinking flux disappears into the sediments if
267 C bottom waters are oxic (assumed adsorbed as oxides). If bottom waters are
268 C anoxic, the sinking flux of Fe is returned to the water column.
269 C
270 C For oxygen, the consumption of oxidant required to respire the settling flux
271 C of organic matter (in support of the no3 bottom flux) diffuses from the
272 C bottom water into the sediment.
273
274 C Assume all NO3 for benthic denitrification is supplied from the bottom water,
275 C and that all organic N is also consumed under denitrification (Complete
276 C Denitrification, sensu Paulmier, Biogeosciences 2009). Therefore, no NO3 is
277 C regenerated from organic matter respired by benthic denitrification
278 C (necessitating the second term in b_no3).
279
280 NO3_sed(i,j) = PONflux_l*drF(k)*hFacC(i,j,k,bi,bj)
281 & - N_den_benthic(i,j,k) / NO3toN
282
283 PO4_sed(i,j) = POPflux_l*drF(k)*hFacC(i,j,k,bi,bj)
284
285 C Oxygen flux into sediments is that required to support non-denitrification
286 C respiration, assuming a 4/5 oxidant ratio of O2 to NO3. Oxygen consumption
287 C is allowed to continue at negative oxygen concentrations, representing
288 C sulphate reduction.
289
290 O2_sed(i,j) = -(O2toN * PONflux_l*drF(k)*hFacC(i,j,k,bi,bj)
291 & - N_den_benthic(i,j,k)* 1.25)
292
293 ENDIF
294
295 C Begin iron uptake calculations by determining ligand bound and free iron.
296 C Both forms are available for biology, but only free iron is scavenged
297 C onto particles and forms colloids.
298
299 lig_stability = kFe_eq_lig_max-(KFe_eq_lig_max-kFe_eq_lig_min)
300 & *(irr_inst(i,j,k)**2
301 & /(kFe_eq_lig_irr**2+irr_inst(i,j,k)**2))
302 & *max(epsln,min(1. _d 0,(PTR_FE(i,j,k)
303 & -kFe_eq_lig_Femin)/
304 & (PTR_FE(i,j,k)+epsln)*1.2 _d 0))
305
306 C Use the quadratic equation to solve for binding between iron and ligands
307
308 FreeFe = (-(1+lig_stability*(ligand-PTR_FE(i,j,k)))
309 & +((1+lig_stability*(ligand-PTR_FE(i,j,k)))**2+4*
310 & lig_stability*PTR_FE(i,j,k))**(0.5))/(2*
311 & lig_stability)
312
313 C Iron scavenging does not occur in anoxic water (Fe2+ is soluble), so set
314 C FreeFe = 0 when anoxic. FreeFe should be interpreted the free iron that
315 C participates in scavenging.
316
317 IF (PTR_O2(i,j,k) .LT. oxic_min) THEN
318 FreeFe = 0. _d 0
319 ENDIF
320
321 C Two mechanisms for iron uptake, in addition to biological production:
322 C colloidal scavenging and scavenging by organic matter.
323
324 Fe_ads_inorg(i,j,k) =
325 & kFe_inorg*(max(1. _d -8,FreeFe))**(1.5)
326
327 C Scavenging of iron by organic matter:
328 C The POM value used is the bottom boundary flux. This does not occur in
329 C oxic waters, but FreeFe is set to 0 in such waters earlier.
330 IF ( PONflux_l .GT. 0. _d 0 ) THEN
331 Fe_ads_org(i,j,k) =
332 & kFE_org*(PONflux_l/(epsln + wsink)
333 & * MasstoN)**(0.58)*FreeFe
334 ENDIF
335
336 C If water is oxic then the iron is remineralized normally. Otherwise
337 C it is completely remineralized (fe 2+ is soluble, but unstable
338 C in oxidizing environments).
339
340 PFEflux_l = (PFEflux_u+(Fe_spm(i,j,k)+Fe_ads_inorg(i,j,k)
341 & +Fe_ads_org(i,j,k))*drF(k)
342 & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
343 & *hFacC(i,j,k,bi,bj))
344
345 C Added the burial flux of sinking particulate iron here as a
346 C diagnostic, needed to calculate mass balance of iron.
347 C this is calculated last for the deepest cell
348
349 Fe_burial(i,j) = PFEflux_l
350
351 IF ( PTR_O2(i,j,k) .LT. oxic_min ) THEN
352 PFEflux_l = 0. _d 0
353 ENDIF
354
355 Fe_reminp(i,j,k) = (PFEflux_u+(Fe_spm(i,j,k)
356 & +Fe_ads_inorg(i,j,k)
357 & +Fe_ads_org(i,j,k))*drF(k)
358 & *hFacC(i,j,k,bi,bj)-PFEflux_l)*recip_drF(k)
359 & *recip_hFacC(i,j,k,bi,bj)
360
361 C Prepare the tracers for the next layer down
362 PONflux_u = PONflux_l
363 POPflux_u = POPflux_l
364 PFEflux_u = PFEflux_l
365 CaCO3flux_u = CaCO3flux_l
366
367 Fe_reminsum(i,j,k) = Fe_reminp(i,j,k) + Fe_sed(i,j,k)
368 & - Fe_ads_org(i,j,k) - Fe_ads_inorg(i,j,k)
369
370 ENDIF
371
372 ENDDO
373 ENDDO
374 ENDDO
375
376 c ---------------------------------------------------------------------
377
378 #ifdef ALLOW_DIAGNOSTICS
379 IF ( useDiagnostics ) THEN
380
381 c 3d local variables
382 CALL DIAGNOSTICS_FILL(Fe_ads_org, 'BLGFEADO',0,Nr,2,bi,bj,
383 & myThid)
384 CALL DIAGNOSTICS_FILL(Fe_ads_inorg, 'BLGFEADI',0,Nr,2,bi,bj,
385 & myThid)
386 CALL DIAGNOSTICS_FILL(Fe_sed, 'BLGFESED',0,Nr,2,bi,bj,myThid)
387 CALL DIAGNOSTICS_FILL(Fe_reminp,'BLGFEREM',0,Nr,2,bi,bj,myThid)
388 CALL DIAGNOSTICS_FILL(N_reminp, 'BLGNREM ',0,Nr,2,bi,bj,myThid)
389 CALL DIAGNOSTICS_FILL(P_reminp, 'BLGPREM ',0,Nr,2,bi,bj,myThid)
390 c 2d local variables
391 CALL DIAGNOSTICS_FILL(Fe_burial,'BLGFEBUR',0,1,2,bi,bj,myThid)
392 CALL DIAGNOSTICS_FILL(NO3_sed, 'BLGNSED ',0,1,2,bi,bj,myThid)
393 CALL DIAGNOSTICS_FILL(PO4_sed, 'BLGPSED ',0,1,2,bi,bj,myThid)
394 CALL DIAGNOSTICS_FILL(O2_sed, 'BLGO2SED',0,1,2,bi,bj,myThid)
395
396 ENDIF
397 #endif /* ALLOW_DIAGNOSTICS */
398
399 #endif /* ALLOW_BLING */
400
401 RETURN
402 END

  ViewVC Help
Powered by ViewVC 1.1.22