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

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

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

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