/[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.1 - (show annotations) (download)
Fri May 23 17:33:43 2014 UTC (11 years, 2 months ago) by mmazloff
Branch: MAIN
Adding package BLING

1 C $Header: $
2 C $Name: $
3
4 #include "BLING_OPTIONS.h"
5
6 CBOP
7 subroutine BLING_REMIN(
8 I PTR_O2, PTR_FE,
9 O POM_prod, Fe_uptake, CaCO3_prod,
10 O POM_remin, POM_diss, Fe_remin, CaCO3_diss,
11 I bi, bj, imin, imax, jmin, jmax,
12 I myIter, myTime, myThid)
13
14 C =================================================================
15 C | subroutine bling_remin
16 C | o Calculate the nutrient flux to depth from bio activity.
17 C | Includes iron export and calcium carbonate (dissolution of
18 C | CaCO3 returns carbonate ions and changes alkalinity).
19 C | - Instant remineralization is assumed.
20 C | - A fraction of POM becomes DOM
21 C =================================================================
22
23 implicit none
24
25 C === Global variables ===
26 C irr_inst :: instantaneous irradiance
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_TAMC
37 # include "tamc.h"
38 #endif
39
40 C === Routine arguments ===
41 C myTime :: current time
42 C myIter :: current timestep
43 C myThid :: thread number
44 _RL dt
45 _RL myTime
46 INTEGER myIter
47 INTEGER myThid
48 C === Input ===
49 C POM_prod :: biological production of sinking particles
50 C Fe_uptake :: biological production of particulate iron
51 C CaCO3_prod :: biological production of CaCO3 shells
52 _RL POM_prod (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
53 _RL Fe_uptake (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
54 _RL CaCO3_prod (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
55 _RL PTR_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
56 _RL PTR_FE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
57 INTEGER imin, imax, jmin, jmax, bi, bj
58 C === Output ===
59 C POM_remin :: remineralization of sinking particles
60 C Fe_remin :: remineralization of particulate iron
61 C CaCO3_diss :: dissolution of CaCO3 shells
62 _RL POM_remin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
63 _RL POM_diss (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
64 _RL Fe_remin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
65 _RL CaCO3_diss (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
66
67 C === Local variables ===
68 C i,j,k :: loop indices
69 C depth_l :: depth of lower interface
70 C deltaPOM :: change in POM due to remin & dissolution
71 C *flux_u, *flux_l :: "*" flux through upper and lower interfaces
72 C *_export :: vertically-integrated export of "*"
73 C zremin :: remineralization lengthscale for nutrients
74 C zremin_caco3 :: remineralization lengthscale for CaCO3
75 C wsink :: speed of sinking particles
76 C fe_sed_source :: iron source from sediments
77 C FreeFe :: ligand-free iron
78 INTEGER i,j,k
79 _RL depth_l
80 _RL deltaPOM
81 _RL POMflux_u
82 _RL POMflux_l
83 _RL PFEflux_u
84 _RL PFEflux_l
85 _RL CaCO3flux_u
86 _RL CaCO3flux_l
87 _RL POM_export (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 _RL PFE_export (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL CaCO3_export(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 _RL zremin
91 _RL zremin_caco3
92 _RL wsink
93 _RL fe_sed_source
94 _RL lig_stability
95 _RL FreeFe
96 CEOP
97
98 c ---------------------------------------------------------------------
99 c Initialize output and diagnostics
100 DO k=1,Nr
101 DO j=jmin,jmax
102 DO i=imin,imax
103 POM_remin(i,j,k) = 0. _d 0
104 Fe_remin(i,j,k) = 0. _d 0
105 CaCO3_diss(i,j,k) = 0. _d 0
106 ENDDO
107 ENDDO
108 ENDDO
109 DO j=jmin,jmax
110 DO i=imin,imax
111 POM_export(i,j) = 0. _d 0
112 PFE_export(i,j) = 0. _d 0
113 CaCO3_export(i,j) = 0. _d 0
114 ENDDO
115 ENDDO
116 POMflux_u = 0. _d 0
117 PFEflux_u = 0. _d 0
118 CaCO3flux_u = 0. _d 0
119
120 c ---------------------------------------------------------------------
121 c Nutrients export/remineralization, CaCO3 export/dissolution
122 c
123 c The flux at the bottom of a grid cell equals
124 C Fb = (Ft + prod*dz) / (1 + zremin*dz)
125 C where Ft is the flux at the top, and prod*dz is the integrated
126 C production of new sinking particles within the layer.
127 C Ft = 0 in the first layer.
128
129 CADJ STORE Fe_uptake = comlev1, key = ikey_dynamics
130
131 C$TAF LOOP = parallel
132 DO j=jmin,jmax
133 C$TAF LOOP = parallel
134 DO i=imin,imax
135 C$TAF init upper_flux = static, Nr
136 DO k=1,Nr
137 C$TAF STORE POMflux_u = upper_flux
138 C$TAF STORE PFEflux_u = upper_flux
139 C$TAF STORE CaCO3flux_u = upper_flux
140
141 IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN
142
143 C Sinking speed is evaluated at the bottom of the cell
144 depth_l=-rF(k+1)
145 IF (depth_l .LE. wsink0z) THEN
146 wsink = wsink0
147 ELSE
148 wsink = wsinkacc * (depth_l - wsink0z) + wsink0
149 ENDIF
150
151 C Nutrient remineralization lengthscale
152 C Not an e-folding scale: this term increases with remineralization.
153 zremin = gamma_POM * ( PTR_O2(i,j,k)**2 /
154 & (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min)
155 & + remin_min )/(wsink + epsln)
156
157 C Calcium remineralization relaxed toward the inverse of the
158 C ca_remin_depth constant value as the calcite saturation approaches 0.
159 zremin_caco3 = 1. _d 0/ca_remin_depth*(1. _d 0-min(1. _d 0,
160 & omegaC(i,j,k,bi,bj)))
161
162 C POM flux leaving the cell
163 POMflux_l = (POMflux_u+POM_prod(i,j,k)*drF(k)
164 & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
165 & *hFacC(i,j,k,bi,bj))
166
167 C CaCO3 flux leaving the cell
168 CaCO3flux_l = (caco3flux_u+CaCO3_prod(i,j,k)*drF(k)
169 & *hFacC(i,j,k,bi,bj))/(1+zremin_caco3*drF(k)
170 & *hFacC(i,j,k,bi,bj))
171
172 C Start with cells that are not the deepest cells
173 IF ((k.LT.Nr) .AND. (hFacC(i,j,k+1,bi,bj).GT.0)) THEN
174
175 C Nutrient accumulation in a cell is given by the biological production
176 C (and instant remineralization) of particulate organic matter
177 C plus flux thought upper interface minus flux through lower interface.
178 C (Since not deepest cell: hFacC=1)
179 deltaPOM = (POMflux_u + POM_prod(i,j,k)*drF(k)
180 & - POMflux_l)*recip_drF(k)
181
182 CaCO3_diss(i,j,k) = (CaCO3flux_u + CaCO3_prod(i,j,k)*drF(k)
183 & - CaCO3flux_l)*recip_drF(k)
184
185 fe_sed_source = 0. _d 0
186
187 ELSE
188 C If this layer is adjacent to bottom topography or it is the deepest
189 C cell of the domain, then remineralize/dissolve in this grid cell
190 C i.e. don't subtract off lower boundary fluxes when calculating remin
191 deltaPOM = POMflux_u*recip_drF(k)
192 & *recip_hFacC(i,j,k,bi,bj)+POM_prod(i,j,k)
193
194 CaCO3_diss(i,j,k) = caco3flux_u*recip_drF(k)
195 & *recip_hFacC(i,j,k,bi,bj)+CaCO3_prod(i,j,k)
196
197 C Iron from sediments: the phosphate flux hitting the bottom boundary
198 C is used to scale the return of iron to the water column.
199 C Maximum value added for numerical stability.
200 fe_sed_source = min(1. _d -11,
201 & max(0. _d 0,FetoPsed/NUTfac*POMflux_l*recip_drF(k)
202 & *recip_hFacC(i,j,k,bi,bj)))
203 #ifdef BLING_ADJOINT_SAFE
204 fe_sed_source = 0. _d 0
205 #endif
206 ENDIF
207
208 C A fraction of POM becomes DOM
209 POM_diss(i,j,k) = deltaPOM*phi_DOM
210 POM_remin(i,j,k) = deltaPOM*(1-phi_DOM)
211
212 C Begin iron uptake calculations by determining ligand bound and free iron.
213 C Both forms are available for biology, but only free iron is scavenged
214 C onto particles and forms colloids.
215 lig_stability = KFeLeq_max-(KFeLeq_max-KFeLeq_min)
216 & *(irr_inst(i,j,k,bi,bj)**2
217 & /(IFeL**2+irr_inst(i,j,k,bi,bj)**2))
218 & *max(0. _d 0,min(1. _d 0,(PTR_FE(i,j,k)-Fe_min)/
219 & (PTR_FE(i,j,k)+epsln)*b_const))
220
221 C Use the quadratic equation to solve for binding between iron and ligands
222
223 FreeFe = (-(1+lig_stability*(ligand-PTR_FE(i,j,k)))
224 & +((1+lig_stability*(ligand-PTR_FE(i,j,k)))**2+4*
225 & lig_stability*PTR_FE(i,j,k))**(0.5))/(2*
226 & lig_stability)
227
228 C Iron scavenging doesn't occur in anoxic water (Fe2+ is soluble), so set
229 C FreeFe = 0 when anoxic. FreeFe should be interpreted the free iron that
230 C participates in scavenging.
231 #ifndef BLING_ADJOINT_SAFE
232 IF (PTR_O2(i,j,k) .LT. O2_min) THEN
233 FreeFe = 0
234 ENDIF
235 #endif
236
237 c Two mechanisms for iron uptake, in addition to biological production:
238 c colloidal scavenging and scavenging by organic matter.
239
240 c Colloidal scavenging:
241 c Minimum function for numerical stability
242 c Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+
243 c & min(0.5/PTRACERS_dTLev(1), kFe_inorg*FreeFe**(0.5))*FreeFe
244
245 Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+
246 & kFe_inorg*FreeFe**(0.5)*FreeFe
247
248 C Scavenging of iron by organic matter:
249 c The POM value used is the bottom boundary flux. This doesn't occur in
250 c oxic waters, but FreeFe is set to 0 in such waters earlier.
251 IF ( POMflux_l .GT. 0. _d 0 ) THEN
252
253 c Minimum function for numerical stability
254 c Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+
255 c & min(0.5/PTRACERS_dTLev(1), kFE_org*(POMflux_l
256 c & *CtoP/NUTfac*12.01/wsink)**(0.58)*FreeFe
257
258 #ifndef BLING_ADJOINT_SAFE
259 Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+
260 & kFE_org*(POMflux_l*CtoP/NUTfac
261 & *12.01/wsink)**(0.58)*FreeFe
262 #else
263 Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+
264 & kFE_org*(POMflux_l*CtoP/NUTfac
265 & *12.01/wsink0)**(0.58)*FreeFe
266 #endif
267 ENDIF
268
269 C If water is oxic then the iron is remineralized normally. Otherwise
270 C it is completely remineralized (fe 2+ is soluble, but unstable
271 C in oxidizing environments).
272
273 pfeflux_l = (pfeflux_u+Fe_uptake(i,j,k)*drF(k)
274 & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
275 & *hFacC(i,j,k,bi,bj))
276
277 #ifndef BLING_ADJOINT_SAFE
278 IF ( PTR_O2(i,j,k) .LT. O2_min ) THEN
279 pfeflux_l = 0
280 ENDIF
281 #endif
282
283 Fe_remin(i,j,k) = (pfeflux_u+Fe_uptake(i,j,k)*drF(k)
284 & *hFacC(i,j,k,bi,bj)-pfeflux_l)*recip_drF(k)
285 & *recip_hFacC(i,j,k,bi,bj)
286
287 C Add sediment source
288 Fe_remin(i,j,k) = Fe_remin(i,j,k) + fe_sed_source
289
290 C Prepare the tracers for the next layer down
291 POMflux_u = POMflux_l
292 PFEflux_u = PFEflux_l
293 CaCO3flux_u = CaCO3flux_l
294
295 C Depth-integrated export (through bottom of water column)
296 C This is calculated last for the deepest cell
297 POM_export(i,j) = POMflux_l
298 PFE_export(i,j) = PFEflux_l
299 CACO3_export(i,j) = CaCO3flux_l
300
301 ENDIF
302
303 ENDDO
304
305 C Reset for next location (i,j)
306 POMflux_u = 0. _d 0
307 PFEflux_u = 0. _d 0
308 CaCO3flux_u = 0. _d 0
309
310 ENDDO
311 ENDDO
312
313 RETURN
314 END

  ViewVC Help
Powered by ViewVC 1.1.22