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

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

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


Revision 1.11 - (hide annotations) (download)
Wed May 24 17:21:15 2017 UTC (7 years 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 mmazloff 1.11 C $Header: /u/gcmpack/MITgcm/pkg/bling/bling_remineralization.F,v 1.10 2017/03/29 15:51:19 mmazloff Exp $
2 mmazloff 1.1 C $Name: $
3    
4     #include "BLING_OPTIONS.h"
5    
6     CBOP
7 jmc 1.2 subroutine BLING_REMIN(
8 mmazloff 1.1 I PTR_NO3, PTR_FE, PTR_O2, irr_inst,
9 jmc 1.2 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 mmazloff 1.1 I bi, bj, imin, imax, jmin, jmax,
13     I myIter, myTime, myThid )
14    
15     C =================================================================
16     C | subroutine bling_remin
17 mmazloff 1.4 C | o Organic matter export and remineralization.
18 jmc 1.2 C | - Sinking particulate flux and diel migration contribute to
19     C | export.
20 mmazloff 1.4 C | - Benthic denitrification
21     C | - Iron source from sediments
22     C | - Iron scavenging
23 mmazloff 1.1 C =================================================================
24    
25     implicit none
26 jmc 1.2
27 mmazloff 1.1 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 mmazloff 1.4
65 mmazloff 1.1 C === Output ===
66 mmazloff 1.4 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 mmazloff 1.1 _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 jmc 1.2
77 mmazloff 1.1 #ifdef ALLOW_BLING
78     C === Local variables ===
79     C i,j,k :: loop indices
80 mmazloff 1.4
81 jmc 1.2 INTEGER i,j,k
82 mmazloff 1.1 INTEGER bttmlyr
83     _RL PONflux_u
84 mmazloff 1.4 _RL PONflux_l
85 mmazloff 1.1 _RL POPflux_u
86 mmazloff 1.4 _RL POPflux_l
87 mmazloff 1.1 _RL PFEflux_u
88 mmazloff 1.4 _RL PFEflux_l
89 mmazloff 1.1 _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 mmazloff 1.4 _RL Fe_ads_org(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104 mmazloff 1.1 _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 jmc 1.2
109 mmazloff 1.9 C ---------------------------------------------------------------------
110     C Initialize output and diagnostics
111 mmazloff 1.1
112     DO k=1,Nr
113     DO j=jmin,jmax
114     DO i=imin,imax
115 mmazloff 1.4 Fe_ads_org(i,j,k) = 0. _d 0
116 mmazloff 1.1 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 mmazloff 1.9 C ---------------------------------------------------------------------
136     C Remineralization
137 jmc 1.2
138 mmazloff 1.1 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 mmazloff 1.6
149     DO k=1,Nr
150    
151     C Initialization here helps taf
152     Fe_ads_org(i,j,k) = 0. _d 0
153 mmazloff 1.1
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 mmazloff 1.8 wsink = wsink0_2d(i,j,bi,bj)
167 mmazloff 1.1 ELSE
168 mmazloff 1.8 wsink = wsinkacc * (depth_l - wsink0z) + wsink0_2d(i,j,bi,bj)
169 mmazloff 1.1 ENDIF
170    
171     C Nutrient remineralization lengthscale
172 jmc 1.2 C Not an e-folding scale: this term increases with remineralization.
173 mmazloff 1.8 zremin = gamma_POM2d(i,j,bi,bj) * ( PTR_O2(i,j,k)**2 /
174 jmc 1.2 & (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min)
175 mmazloff 1.1 & + remin_min )/(wsink + epsln)
176    
177 jmc 1.2 C Calcium remineralization relaxed toward the inverse of the
178     C ca_remin_depth constant value as the calcite saturation approaches 0.
179 mmazloff 1.1 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 jmc 1.2 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 mmazloff 1.1 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 mmazloff 1.9 C NOW DO BOTTOM LAYER
213 mmazloff 1.1 ELSE
214 jmc 1.2 C If this layer is adjacent to bottom topography or it is the deepest
215 mmazloff 1.1 C cell of the domain, then remineralize/dissolve in this grid cell
216 jmc 1.2 C i.e. do not subtract off lower boundary fluxes when calculating remin
217 mmazloff 1.1
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 jmc 1.2
227 mmazloff 1.9 C Efflux Fed out of sediments
228 jmc 1.2 C The phosphate flux hitting the bottom boundary
229 mmazloff 1.1 C is used to scale the return of iron to the water column.
230 jmc 1.2 C Maximum value added for numerical stability.
231 mmazloff 1.1
232     POC_sed = PONflux_l * CtoN
233    
234 mmazloff 1.9 Fe_sed(i,j,k) = max(epsln, FetoC_sed * POC_sed * recip_drF(k)
235     & *recip_hFacC(i,j,k,bi,bj))
236 jmc 1.2
237 mmazloff 1.9 cav log_btm_flx = 0. _d 0
238     log_btm_flx = 1. _d -20
239 mmazloff 1.1
240 mmazloff 1.10 CMM: this is causing instability in the adjoint. Needs debugging
241     #ifndef BLING_ADJOINT_SAFE
242 mmazloff 1.1 IF (POC_sed .gt. 0. _d 0) THEN
243    
244 mmazloff 1.9 C Convert from mol N m-2 s-1 to umol C cm-2 d-1 and take the log
245 mmazloff 1.1
246 jmc 1.2 log_btm_flx = log10(min(43.0 _d 0, POC_sed *
247     & 86400. _d 0 * 100.0 _d 0))
248 mmazloff 1.1
249 mmazloff 1.9 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 mmazloff 1.1
252     N_den_benthic(i,j,k) = min (POC_sed * NO3toN / CtoN,
253 jmc 1.2 & (10 _d 0)**(-0.9543 _d 0 + 0.7662 _d 0 *
254 mmazloff 1.1 & log_btm_flx - 0.235 _d 0 * log_btm_flx * log_btm_flx)
255     & / (CtoN * 86400. _d 0 * 100.0 _d 0) * NO3toN *
256 jmc 1.2 & PTR_NO3(i,j,k) / (k_no3 + PTR_NO3(i,j,k)) ) *
257 mmazloff 1.1 & recip_drF(k)
258 jmc 1.2
259 mmazloff 1.9 ENDIF
260 mmazloff 1.10 #endif
261 mmazloff 1.1
262 mmazloff 1.9 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 mmazloff 1.1
280     NO3_sed(i,j) = PONflux_l*drF(k)*hFacC(i,j,k,bi,bj)
281     & - N_den_benthic(i,j,k) / NO3toN
282 jmc 1.2
283 mmazloff 1.1 PO4_sed(i,j) = POPflux_l*drF(k)*hFacC(i,j,k,bi,bj)
284    
285 mmazloff 1.9 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 mmazloff 1.1
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 jmc 1.2 ENDIF
294 mmazloff 1.1
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 jmc 1.2 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 mmazloff 1.1 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 mmazloff 1.9 C Two mechanisms for iron uptake, in addition to biological production:
322     C colloidal scavenging and scavenging by organic matter.
323 mmazloff 1.1
324 jmc 1.2 Fe_ads_inorg(i,j,k) =
325 mmazloff 1.1 & kFe_inorg*(max(1. _d -8,FreeFe))**(1.5)
326    
327     C Scavenging of iron by organic matter:
328 mmazloff 1.9 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 mmazloff 1.1 IF ( PONflux_l .GT. 0. _d 0 ) THEN
331 mmazloff 1.4 Fe_ads_org(i,j,k) =
332 jmc 1.2 & kFE_org*(PONflux_l/(epsln + wsink)
333 mmazloff 1.1 & * MasstoN)**(0.58)*FreeFe
334     ENDIF
335 jmc 1.2
336 mmazloff 1.1 C If water is oxic then the iron is remineralized normally. Otherwise
337 jmc 1.2 C it is completely remineralized (fe 2+ is soluble, but unstable
338 mmazloff 1.1 C in oxidizing environments).
339    
340     PFEflux_l = (PFEflux_u+(Fe_spm(i,j,k)+Fe_ads_inorg(i,j,k)
341 mmazloff 1.4 & +Fe_ads_org(i,j,k))*drF(k)
342 mmazloff 1.1 & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
343     & *hFacC(i,j,k,bi,bj))
344    
345 mmazloff 1.9 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 mmazloff 1.1
349 mmazloff 1.9 Fe_burial(i,j) = PFEflux_l
350 mmazloff 1.1
351     IF ( PTR_O2(i,j,k) .LT. oxic_min ) THEN
352 mmazloff 1.9 PFEflux_l = 0. _d 0
353 mmazloff 1.1 ENDIF
354    
355 mmazloff 1.9 Fe_reminp(i,j,k) = (PFEflux_u+(Fe_spm(i,j,k)
356 mmazloff 1.1 & +Fe_ads_inorg(i,j,k)
357 mmazloff 1.4 & +Fe_ads_org(i,j,k))*drF(k)
358 mmazloff 1.9 & *hFacC(i,j,k,bi,bj)-PFEflux_l)*recip_drF(k)
359 mmazloff 1.1 & *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 mmazloff 1.4 & - Fe_ads_org(i,j,k) - Fe_ads_inorg(i,j,k)
369 mmazloff 1.1
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 mmazloff 1.11 CALL DIAGNOSTICS_FILL(Fe_ads_org, 'BLGFEADO',0,Nr,2,bi,bj,
383 mmazloff 1.4 & myThid)
384 mmazloff 1.11 CALL DIAGNOSTICS_FILL(Fe_ads_inorg, 'BLGFEADI',0,Nr,2,bi,bj,
385 mmazloff 1.1 & 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