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

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

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


Revision 1.1 - (hide annotations) (download)
Fri May 23 17:33:43 2014 UTC (11 years, 2 months ago) by mmazloff
Branch: MAIN
Adding package BLING

1 mmazloff 1.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