/[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.2 by mmazloff, Sun Feb 28 21:49:24 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        _RL PONflux_u
110         DO i=imin,imax        _RL POPflux_u
111          POM_export(i,j)    = 0. _d 0        _RL PFEflux_u
112          PFE_export(i,j)    = 0. _d 0        _RL CaCO3flux_u
113          CaCO3_export(i,j)  = 0. _d 0        _RL PONflux_l
114         ENDDO        _RL POPflux_l
115        ENDDO        _RL PFEflux_l
116        POMflux_u            = 0. _d 0        _RL CaCO3flux_l
117        PFEflux_u            = 0. _d 0        _RL depth_l
118        CaCO3flux_u          = 0. _d 0        _RL zremin
119          _RL zremin_caco3
120  c ---------------------------------------------------------------------        _RL wsink
121  c  Nutrients export/remineralization, CaCO3 export/dissolution        _RL POC_sed
122  c        _RL Fe_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
123  c  The flux at the bottom of a grid cell equals        _RL NO3_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124  C  Fb = (Ft + prod*dz) / (1 + zremin*dz)        _RL PO4_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 O2_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126  C  production of new sinking particles within the layer.        _RL lig_stability
127  C  Ft = 0 in the first layer.        _RL FreeFe
128          _RL Fe_ads_inorg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
129  CADJ STORE Fe_uptake = comlev1, key = ikey_dynamics        _RL Fe_ads_org(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
130          _RL log_btm_flx
131  C$TAF LOOP = parallel        _RL Fe_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
132        DO j=jmin,jmax        _RL o2_upper
133  C$TAF LOOP = parallel        _RL o2_lower
134         DO i=imin,imax        _RL dz_upper
135  C$TAF init upper_flux = static, Nr        _RL dz_lower
136          DO k=1,Nr        _RL temp_upper
137  C$TAF STORE POMflux_u = upper_flux        _RL temp_lower
138  C$TAF STORE PFEflux_u = upper_flux        _RL z_dvm_regr
139  C$TAF STORE CaCO3flux_u = upper_flux        _RL frac_migr
140          _RL fdvm_migr
141           IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN        _RL fdvm_stat
142          _RL fdvmn_vint
143  C  Sinking speed is evaluated at the bottom of the cell        _RL fdvmp_vint
144            depth_l=-rF(k+1)        _RL fdvmfe_vint
145            IF (depth_l .LE. wsink0z)  THEN        _RL z_dvm
146             wsink = wsink0        _RL N_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
147            ELSE        _RL P_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
148             wsink = wsinkacc * (depth_l - wsink0z) + wsink0        _RL Fe_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
149            ENDIF        _RL dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
150          _RL mld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151  C  Nutrient remineralization lengthscale        _RL Fe_burial(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152  C  Not an e-folding scale: this term increases with remineralization.        _RL x_erfcc,z_erfcc,t_erfcc,erfcc
153            zremin = gamma_POM * ( PTR_O2(i,j,k)**2 /  cxx order
154       &               (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min)  CEOP
155       &               + remin_min )/(wsink + epsln)    
156    c ---------------------------------------------------------------------
157  C  Calcium remineralization relaxed toward the inverse of the  c  Initialize output and diagnostics
158  C  ca_remin_depth constant value as the calcite saturation approaches 0.           DO k=1,Nr
159            zremin_caco3 = 1. _d 0/ca_remin_depth*(1. _d 0-min(1. _d 0,          DO j=jmin,jmax
160       &               omegaC(i,j,k,bi,bj)))            DO i=imin,imax
161                  Fe_ads_org(i,j,k)   = 0. _d 0
162  C  POM flux leaving the cell                Fe_ads_inorg(i,j,k) = 0. _d 0
163            POMflux_l = (POMflux_u+POM_prod(i,j,k)*drF(k)                N_reminp(i,j,k)     = 0. _d 0
164       &           *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)                P_reminp(i,j,k)     = 0. _d 0
165       &           *hFacC(i,j,k,bi,bj))                Fe_reminp(i,j,k)    = 0. _d 0
166                  Fe_reminsum(i,j,k)  = 0. _d 0
167  C  CaCO3 flux leaving the cell                N_remindvm(i,j,k)   = 0. _d 0
168            CaCO3flux_l = (caco3flux_u+CaCO3_prod(i,j,k)*drF(k)                P_remindvm(i,j,k)   = 0. _d 0
169       &           *hFacC(i,j,k,bi,bj))/(1+zremin_caco3*drF(k)                Fe_remindvm(i,j,k)  = 0. _d 0
170       &           *hFacC(i,j,k,bi,bj))                N_den_benthic(i,j,k)= 0. _d 0
171                  CaCO3_diss(i,j,k)   = 0. _d 0
172  C  Start with cells that are not the deepest cells            ENDDO
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          DO j=jmin,jmax
176  C  (and instant remineralization) of particulate organic matter            DO i=imin,imax
177  C  plus flux thought upper interface minus flux through lower interface.                Fe_burial(i,j)       = 0. _d 0
178  C  (Since not deepest cell: hFacC=1)                NO3_sed(i,j)         = 0. _d 0
179             deltaPOM = (POMflux_u + POM_prod(i,j,k)*drF(k)                PO4_sed(i,j)         = 0. _d 0
180       &                    - POMflux_l)*recip_drF(k)                O2_sed(i,j)          = 0. _d 0
181              ENDDO
182             CaCO3_diss(i,j,k) = (CaCO3flux_u + CaCO3_prod(i,j,k)*drF(k)          ENDDO
183       &                    - CaCO3flux_l)*recip_drF(k)  
184    c ---------------------------------------------------------------------
185             fe_sed_source = 0. _d 0  c  Remineralization
186            
187            ELSE  CADJ STORE Fe_ads_org = comlev1, key = ikey_dynamics
188  C  If this layer is adjacent to bottom topography or it is the deepest  cxx needed?
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  C$TAF LOOP = parallel
191             deltaPOM = POMflux_u*recip_drF(k)        DO j=jmin,jmax
192       &            *recip_hFacC(i,j,k,bi,bj)+POM_prod(i,j,k)  C$TAF LOOP = parallel
193           DO i=imin,imax
194             CaCO3_diss(i,j,k) = caco3flux_u*recip_drF(k)  cmm C$TAF init upper_flux = static, Nr
195       &            *recip_hFacC(i,j,k,bi,bj)+CaCO3_prod(i,j,k)  
196    C  Initialize upper flux
197  C  Iron from sediments: the phosphate flux hitting the bottom boundary          PONflux_u            = 0. _d 0
198  C  is used to scale the return of iron to the water column.          POPflux_u            = 0. _d 0
199  C  Maximum value added for numerical stability.          PFEflux_u            = 0. _d 0
200             fe_sed_source = min(1. _d -11,          CaCO3flux_u          = 0. _d 0
201       &            max(0. _d 0,FetoPsed/NUTfac*POMflux_l*recip_drF(k)  
202       &            *recip_hFacC(i,j,k,bi,bj)))          DO k=1,Nr
203  #ifdef BLING_ADJOINT_SAFE  c C$TAF STORE PONflux_u = upper_flux
204             fe_sed_source = 0. _d 0  c C$TAF STORE POPflux_u = upper_flux
205  #endif            c C$TAF STORE PFEflux_u = upper_flux
206            ENDIF  c C$TAF STORE CaCO3flux_u = upper_flux
207    CADJ STORE PONflux_u, POPflux_u, PFEflux_u, CaCO3flux_u =
208  C  A fraction of POM becomes DOM  CADJ &     comlev1, key = ikey_dynamics, kind = isbyte
209            POM_diss(i,j,k) = deltaPOM*phi_DOM  CADJ STORE Fe_ads_org =
210            POM_remin(i,j,k) = deltaPOM*(1-phi_DOM)  CADJ &     comlev1, key = ikey_dynamics, kind = isbyte
211    CMM)
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           IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN
214  C  onto particles and forms colloids.  
215            lig_stability = KFeLeq_max-(KFeLeq_max-KFeLeq_min)  C  Sinking speed is evaluated at the bottom of the cell
216       &             *(irr_inst(i,j,k,bi,bj)**2            depth_l=-rF(k+1)
217       &             /(IFeL**2+irr_inst(i,j,k,bi,bj)**2))            IF (depth_l .LE. wsink0z)  THEN
218       &             *max(0. _d 0,min(1. _d 0,(PTR_FE(i,j,k)-Fe_min)/             wsink = wsink0
219       &             (PTR_FE(i,j,k)+epsln)*b_const))            ELSE
220               wsink = wsinkacc * (depth_l - wsink0z) + wsink0
221  C  Use the quadratic equation to solve for binding between iron and ligands            ENDIF
222    
223            FreeFe = (-(1+lig_stability*(ligand-PTR_FE(i,j,k)))  C  Nutrient remineralization lengthscale
224       &            +((1+lig_stability*(ligand-PTR_FE(i,j,k)))**2+4*  C  Not an e-folding scale: this term increases with remineralization.
225       &            lig_stability*PTR_FE(i,j,k))**(0.5))/(2*            zremin = gamma_POM * ( PTR_O2(i,j,k)**2 /
226       &            lig_stability)       &               (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min)
227         &               + remin_min )/(wsink + epsln)
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  C  Calcium remineralization relaxed toward the inverse of the
230  C  participates in scavenging.  C  ca_remin_depth constant value as the calcite saturation approaches 0.  
231  #ifndef BLING_ADJOINT_SAFE            zremin_caco3 = 1. _d 0/ca_remin_depth*(1. _d 0 - min(1. _d 0,
232            IF (PTR_O2(i,j,k) .LT. O2_min)  THEN       &               omegaC(i,j,k,bi,bj) + epsln ))
233             FreeFe = 0  
234            ENDIF  
235  #endif  C  POM flux leaving the cell
236              PONflux_l = (PONflux_u+N_spm(i,j,k)*drF(k)
237  c  Two mechanisms for iron uptake, in addition to biological production:       &           *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
238  c  colloidal scavenging and scavenging by organic matter.       &           *hFacC(i,j,k,bi,bj))
239    C!! multiply by intercept_frac ???
240  c  Colloidal scavenging:  
241  c  Minimum function for numerical stability            POPflux_l = (POPflux_u+P_spm(i,j,k)*drF(k)
242  c           Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+       &           *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
243  c     &       min(0.5/PTRACERS_dTLev(1), kFe_inorg*FreeFe**(0.5))*FreeFe       &           *hFacC(i,j,k,bi,bj))
244    C!! multiply by intercept_frac ???
245             Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+  
246       &       kFe_inorg*FreeFe**(0.5)*FreeFe  C  CaCO3 flux leaving the cell
247              CaCO3flux_l = (caco3flux_u+CaCO3_uptake(i,j,k)*drF(k)
248  C  Scavenging of iron by organic matter:       &           *hFacC(i,j,k,bi,bj))/(1+zremin_caco3*drF(k)
249  c  The POM value used is the bottom boundary flux. This doesn't occur in       &           *hFacC(i,j,k,bi,bj))
250  c  oxic waters, but FreeFe is set to 0 in such waters earlier.  C!! multiply by intercept_frac ???
251             IF ( POMflux_l .GT. 0. _d 0 ) THEN  
252              
253  c  Minimum function for numerical stability  C  Start with cells that are not the deepest cells
254  c            Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+            IF ((k.LT.Nr) .AND. (hFacC(i,j,k+1,bi,bj).GT.0)) THEN
255  c     &           min(0.5/PTRACERS_dTLev(1), kFE_org*(POMflux_l  
256  c     &           *CtoP/NUTfac*12.01/wsink)**(0.58)*FreeFe  C  Nutrient accumulation in a cell is given by the biological production
257    C  (and instant remineralization) of particulate organic matter
258  #ifndef BLING_ADJOINT_SAFE  C  plus flux thought upper interface minus flux through lower interface.
259              Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+  C  (Since not deepest cell: hFacC=1)
260       &           kFE_org*(POMflux_l*CtoP/NUTfac             N_reminp(i,j,k) = (PONflux_u + N_spm(i,j,k)*drF(k)
261       &            *12.01/wsink)**(0.58)*FreeFe       &                    - PONflux_l)*recip_drF(k)
262  #else  
263              Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+             P_reminp(i,j,k) = (POPflux_u + P_spm(i,j,k)*drF(k)
264       &           kFE_org*(POMflux_l*CtoP/NUTfac       &                    - POPflux_l)*recip_drF(k)
265       &            *12.01/wsink0)**(0.58)*FreeFe  
266  #endif             CaCO3_diss(i,j,k) = (CaCO3flux_u + CaCO3_uptake(i,j,k)
267             ENDIF       &                    *drF(k) - CaCO3flux_l)*recip_drF(k)
268              
269  C  If water is oxic then the iron is remineralized normally. Otherwise             Fe_sed(i,j,k) = 0. _d 0
270  C  it is completely remineralized (fe 2+ is soluble, but unstable  
271  C  in oxidizing environments).  
272              ELSE
273             pfeflux_l = (pfeflux_u+Fe_uptake(i,j,k)*drF(k)  C  If this layer is adjacent to bottom topography or it is the deepest
274       &            *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)  C  cell of the domain, then remineralize/dissolve in this grid cell
275       &            *hFacC(i,j,k,bi,bj))  C  i.e. don't subtract off lower boundary fluxes when calculating remin
276    
277  #ifndef BLING_ADJOINT_SAFE             N_reminp(i,j,k) = PONflux_u*recip_drF(k)
278            IF ( PTR_O2(i,j,k) .LT. O2_min ) THEN       &                    *recip_hFacC(i,j,k,bi,bj) + N_spm(i,j,k)
279              pfeflux_l = 0  
280            ENDIF             P_reminp(i,j,k) = POPflux_u*recip_drF(k)
281  #endif       &                    *recip_hFacC(i,j,k,bi,bj) + P_spm(i,j,k)
282    
283            Fe_remin(i,j,k) = (pfeflux_u+Fe_uptake(i,j,k)*drF(k)             CaCO3_diss(i,j,k) = CaCO3flux_u*recip_drF(k)
284       &            *hFacC(i,j,k,bi,bj)-pfeflux_l)*recip_drF(k)       &                  *recip_hFacC(i,j,k,bi,bj) + CaCO3_uptake(i,j,k)
285       &            *recip_hFacC(i,j,k,bi,bj)      
286        
287  C  Add sediment source  c  Efflux Fed out of sediments
288            Fe_remin(i,j,k) = Fe_remin(i,j,k) + fe_sed_source  C  The phosphate flux hitting the bottom boundary
289    C  is used to scale the return of iron to the water column.
290  C  Prepare the tracers for the next layer down  C  Maximum value added for numerical stability.
291            POMflux_u   = POMflux_l  
292            PFEflux_u   = PFEflux_l             POC_sed = PONflux_l * CtoN
293            CaCO3flux_u = CaCO3flux_l  
294               Fe_sed(i,j,k) = min(1. _d -11,
295  C  Depth-integrated export (through bottom of water column)       &            max(epsln, FetoC_sed * POC_sed * recip_drF(k)
296  C  This is calculated last for the deepest cell       &            *recip_hFacC(i,j,k,bi,bj)))
297            POM_export(i,j)   = POMflux_l      
298            PFE_export(i,j)   = PFEflux_l  #ifdef BLING_ADJOINT_SAFE
299            CACO3_export(i,j) = CaCO3flux_l             Fe_sed(i,j,k) = 0. _d 0
300    #endif
301           ENDIF  
302    
303          ENDDO  cav temporary until I figure out why this is problematic for adjoint          
304    #ifndef BLING_ADJOINT_SAFE
305  C  Reset for next location (i,j)            
306          POMflux_u   = 0. _d 0  #ifndef USE_SGS_SED
307          PFEflux_u   = 0. _d 0  c Calculate benthic denitrification and Fe efflux here, if the subgridscale sediment module is not being used.
308          CaCO3flux_u = 0. _d 0  
309               IF (POC_sed .gt. 0. _d 0) THEN
310         ENDDO  
311        ENDDO              log_btm_flx = 0. _d 0
312    
313        RETURN  c Convert from mol N m-2 s-1 to umol C cm-2 d-1 and take the log
314        END  
315                log_btm_flx = log10(min(43.0 _d 0, POC_sed *
316         &           86400. _d 0 * 100.0 _d 0))
317    
318    c Metamodel gives units of umol C cm-2 d-1, convert to mol N m-2 s-1 and multiply by
319    c no3_2_n to give NO3 consumption rate
320    
321                 N_den_benthic(i,j,k) = min (POC_sed * NO3toN / CtoN,
322         &         (10 _d 0)**(-0.9543 _d 0 + 0.7662 _d 0 *
323         &         log_btm_flx - 0.235 _d 0 * log_btm_flx * log_btm_flx)
324         &         / (CtoN * 86400. _d 0 * 100.0 _d 0) * NO3toN *
325         &         PTR_NO3(i,j,k) / (k_no3 + PTR_NO3(i,j,k)) ) *
326         &         recip_drF(k)
327    
328              endif
329    
330    #endif
331    
332    #endif
333    
334    c ---------------------------------------------------------------------
335    c Calculate external bottom fluxes for tracer_vertdiff. Positive fluxes are into the water
336    c column from the seafloor. For P, the bottom flux puts the sinking flux reaching the bottom
337    c cell into the water column through diffusion. For iron, the sinking flux disappears into the
338    c sediments if bottom waters are oxic (assumed adsorbed as oxides). If bottom waters are anoxic,
339    c the sinking flux of Fe is returned to the water column.
340    c
341    c For oxygen, the consumption of oxidant required to respire  
342    c the settling flux of organic matter (in support of the
343    c no3 bottom flux) diffuses from the bottom water into the sediment.
344    
345    c Assume all NO3 for benthic denitrification is supplied from the bottom water, and that
346    c all organic N is also consumed under denitrification (Complete Denitrification, sensu
347    c Paulmier, Biogeosciences 2009). Therefore, no NO3 is regenerated from organic matter
348    c respired by benthic denitrification (necessitating the second term in b_no3).
349    
350              NO3_sed(i,j) = PONflux_l*drF(k)*hFacC(i,j,k,bi,bj)
351         &                   - N_den_benthic(i,j,k) / NO3toN
352        
353              PO4_sed(i,j) = POPflux_l*drF(k)*hFacC(i,j,k,bi,bj)
354    
355    c Oxygen flux into sediments is that required to support non-denitrification respiration,
356    c assuming a 4/5 oxidant ratio of O2 to NO3. Oxygen consumption is allowed to continue
357    c at negative oxygen concentrations, representing sulphate reduction.
358    
359              O2_sed(i,j) = -(O2toN * PONflux_l*drF(k)*hFacC(i,j,k,bi,bj)
360         &                  - N_den_benthic(i,j,k)* 1.25)
361    
362              ENDIF
363    
364    
365    C  Begin iron uptake calculations by determining ligand bound and free iron.
366    C  Both forms are available for biology, but only free iron is scavenged
367    C  onto particles and forms colloids.
368    
369              lig_stability = kFe_eq_lig_max-(KFe_eq_lig_max-kFe_eq_lig_min)
370         &             *(irr_inst(i,j,k)**2
371         &             /(kFe_eq_lig_irr**2+irr_inst(i,j,k)**2))
372         &             *max(epsln,min(1. _d 0,(PTR_FE(i,j,k)
373         &             -kFe_eq_lig_Femin)/
374         &             (PTR_FE(i,j,k)+epsln)*1.2  _d 0))
375    
376    C  Use the quadratic equation to solve for binding between iron and ligands
377    
378              FreeFe = (-(1+lig_stability*(ligand-PTR_FE(i,j,k)))
379         &            +((1+lig_stability*(ligand-PTR_FE(i,j,k)))**2+4*
380         &            lig_stability*PTR_FE(i,j,k))**(0.5))/(2*
381         &            lig_stability)
382    
383    C  Iron scavenging doesn't occur in anoxic water (Fe2+ is soluble), so set  
384    C  FreeFe = 0 when anoxic.  FreeFe should be interpreted the free iron that
385    C  participates in scavenging.
386    
387    #ifndef BLING_ADJOINT_SAFE
388              IF (PTR_O2(i,j,k) .LT. oxic_min)  THEN
389               FreeFe = 0. _d 0
390              ENDIF
391    #endif
392    
393    c  Two mechanisms for iron uptake, in addition to biological production:
394    c  colloidal scavenging and scavenging by organic matter.
395    
396    c  Colloidal scavenging:
397    c  Minimum function for numerical stability
398    c           Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+
399    c     &       min(0.5/PTRACERS_dTLev(1), kFe_inorg*FreeFe**(0.5))*FreeFe
400    
401               Fe_ads_inorg(i,j,k) =
402         &       kFe_inorg*(max(1. _d -8,FreeFe))**(1.5)
403    
404    C  Scavenging of iron by organic matter:
405    c  The POM value used is the bottom boundary flux. This doesn't occur in
406    c  oxic waters, but FreeFe is set to 0 in such waters earlier.
407               IF ( PONflux_l .GT. 0. _d 0 ) THEN
408              
409    c  Minimum function for numerical stability
410    c            Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+
411    c     &           min(0.5/PTRACERS_dTLev(1), kFE_org*(POMflux_l
412    c     &           *CtoP/NUTfac*12.01/wsink)**(0.58)*FreeFe
413    
414    #ifndef BLING_ADJOINT_SAFE
415                Fe_ads_org(i,j,k) =
416         &           kFE_org*(PONflux_l/(epsln + wsink)
417         &             * MasstoN)**(0.58)*FreeFe
418    #else
419                Fe_ads_org(i,j,k) =
420         &           kFE_org*(PONflux_l/(epsln + wsink0)
421         &             * MasstoN)**(0.58)*FreeFe
422    #endif
423               ENDIF
424    
425    
426    
427    C  If water is oxic then the iron is remineralized normally. Otherwise
428    C  it is completely remineralized (fe 2+ is soluble, but unstable
429    C  in oxidizing environments).
430    
431               PFEflux_l = (PFEflux_u+(Fe_spm(i,j,k)+Fe_ads_inorg(i,j,k)
432         &            +Fe_ads_org(i,j,k))*drF(k)
433         &            *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k)
434         &            *hFacC(i,j,k,bi,bj))
435    
436    
437    c Added the burial flux of sinking particulate iron here as a
438    c diagnostic, needed to calculate mass balance of iron.
439    c this is calculated last for the deepest cell
440    
441               Fe_burial(i,j) = PFeflux_l
442    
443    
444    #ifndef BLING_ADJOINT_SAFE
445               IF ( PTR_O2(i,j,k) .LT. oxic_min ) THEN
446                pfeflux_l = 0
447               ENDIF
448    #endif
449    
450               Fe_reminp(i,j,k) = (pfeflux_u+(Fe_spm(i,j,k)
451         &            +Fe_ads_inorg(i,j,k)
452         &            +Fe_ads_org(i,j,k))*drF(k)
453         &            *hFacC(i,j,k,bi,bj)-pfeflux_l)*recip_drF(k)
454         &            *recip_hFacC(i,j,k,bi,bj)
455    C!! there's an intercept_frac here... need to add
456    
457    
458    C  Prepare the tracers for the next layer down
459              PONflux_u   = PONflux_l
460              POPflux_u   = POPflux_l
461              PFEflux_u   = PFEflux_l
462              CaCO3flux_u = CaCO3flux_l
463    
464    c
465              Fe_reminsum(i,j,k) = Fe_reminp(i,j,k) + Fe_sed(i,j,k)
466         &                       - Fe_ads_org(i,j,k) - Fe_ads_inorg(i,j,k)
467    cc             Fe_reminsum(i,j,k) = 0. _d 0
468    
469             ENDIF
470    
471            ENDDO
472           ENDDO
473          ENDDO
474    
475    CADJ STORE Fe_ads_org = comlev1, key = ikey_dynamics
476    cxx needed?
477    
478    
479    c ---------------------------------------------------------------------
480    
481    #ifdef ALLOW_DIAGNOSTICS
482          IF ( useDiagnostics ) THEN
483    
484    c 3d local variables
485    c        CALL DIAGNOSTICS_FILL(POC_flux,'BLGPOCF ',0,Nr,2,bi,bj,myThid)
486            CALL DIAGNOSTICS_FILL(Fe_ads_inorg,'BLGFEAI',0,Nr,2,bi,bj,
487         &       myThid)
488            CALL DIAGNOSTICS_FILL(Fe_sed,'BLGFESED',0,Nr,2,bi,bj,myThid)
489            CALL DIAGNOSTICS_FILL(Fe_reminp,'BLGFEREM',0,Nr,2,bi,bj,myThid)
490            CALL DIAGNOSTICS_FILL(N_den_benthic,'BLGNDENB',0,Nr,2,bi,bj,
491         &       myThid)
492    c        CALL DIAGNOSTICS_FILL(N_den_pelag,'BLGNDENP',0,Nr,2,bi,bj,myThid)
493            CALL DIAGNOSTICS_FILL(N_reminp,'BLGNREM ',0,Nr,2,bi,bj,myThid)
494            CALL DIAGNOSTICS_FILL(P_reminp,'BLGPREM ',0,Nr,2,bi,bj,myThid)
495    c 2d local variables
496            CALL DIAGNOSTICS_FILL(Fe_burial,'BLGFEBUR',0,1,2,bi,bj,myThid)
497            CALL DIAGNOSTICS_FILL(NO3_sed,'BLGNSED ',0,1,2,bi,bj,myThid)
498            CALL DIAGNOSTICS_FILL(PO4_sed,'BLGPSED ',0,1,2,bi,bj,myThid)
499            CALL DIAGNOSTICS_FILL(O2_sed,'BLGO2SED',0,1,2,bi,bj,myThid)
500    c these variables are currently 1d, could be 3d for diagnostics
501    c (or diag_fill could be called inside loop - which is faster?)
502    c        CALL DIAGNOSTICS_FILL(zremin,'BLGZREM ',0,Nr,2,bi,bj,myThid)
503    
504          ENDIF
505    #endif /* ALLOW_DIAGNOSTICS */
506    
507    
508    
509    #endif /* ALLOW_BLING */
510    
511          RETURN
512          END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22