Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch
--- MITgcm_contrib/bling/pkg/bling_remineralization.F 2014/05/23 17:33:43 1.1 +++ MITgcm_contrib/bling/pkg/bling_remineralization.F 2016/02/28 21:49:24 1.2 @@ -1,314 +1,512 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/bling/pkg/bling_remineralization.F,v 1.1 2014/05/23 17:33:43 mmazloff Exp $ -C $Name: $ - -#include "BLING_OPTIONS.h" - -CBOP - subroutine BLING_REMIN( - I PTR_O2, PTR_FE, - O POM_prod, Fe_uptake, CaCO3_prod, - O POM_remin, POM_diss, Fe_remin, CaCO3_diss, - I bi, bj, imin, imax, jmin, jmax, - I myIter, myTime, myThid) - -C ================================================================= -C | subroutine bling_remin -C | o Calculate the nutrient flux to depth from bio activity. -C | Includes iron export and calcium carbonate (dissolution of -C | CaCO3 returns carbonate ions and changes alkalinity). -C | - Instant remineralization is assumed. -C | - A fraction of POM becomes DOM -C ================================================================= - - implicit none - -C === Global variables === -C irr_inst :: instantaneous irradiance - -#include "SIZE.h" -#include "DYNVARS.h" -#include "EEPARAMS.h" -#include "PARAMS.h" -#include "GRID.h" -#include "BLING_VARS.h" -#include "PTRACERS_SIZE.h" -#include "PTRACERS_PARAMS.h" -#ifdef ALLOW_AUTODIFF_TAMC -# include "tamc.h" -#endif - -C === Routine arguments === -C myTime :: current time -C myIter :: current timestep -C myThid :: thread number - _RL dt - _RL myTime - INTEGER myIter - INTEGER myThid -C === Input === -C POM_prod :: biological production of sinking particles -C Fe_uptake :: biological production of particulate iron -C CaCO3_prod :: biological production of CaCO3 shells - _RL POM_prod (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL Fe_uptake (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL CaCO3_prod (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL PTR_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL PTR_FE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - INTEGER imin, imax, jmin, jmax, bi, bj -C === Output === -C POM_remin :: remineralization of sinking particles -C Fe_remin :: remineralization of particulate iron -C CaCO3_diss :: dissolution of CaCO3 shells - _RL POM_remin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL POM_diss (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL Fe_remin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL CaCO3_diss (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - -C === Local variables === -C i,j,k :: loop indices -C depth_l :: depth of lower interface -C deltaPOM :: change in POM due to remin & dissolution -C *flux_u, *flux_l :: "*" flux through upper and lower interfaces -C *_export :: vertically-integrated export of "*" -C zremin :: remineralization lengthscale for nutrients -C zremin_caco3 :: remineralization lengthscale for CaCO3 -C wsink :: speed of sinking particles -C fe_sed_source :: iron source from sediments -C FreeFe :: ligand-free iron - INTEGER i,j,k - _RL depth_l - _RL deltaPOM - _RL POMflux_u - _RL POMflux_l - _RL PFEflux_u - _RL PFEflux_l - _RL CaCO3flux_u - _RL CaCO3flux_l - _RL POM_export (1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL PFE_export (1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL CaCO3_export(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL zremin - _RL zremin_caco3 - _RL wsink - _RL fe_sed_source - _RL lig_stability - _RL FreeFe -CEOP - -c --------------------------------------------------------------------- -c Initialize output and diagnostics - DO k=1,Nr - DO j=jmin,jmax - DO i=imin,imax - POM_remin(i,j,k) = 0. _d 0 - Fe_remin(i,j,k) = 0. _d 0 - CaCO3_diss(i,j,k) = 0. _d 0 - ENDDO - ENDDO - ENDDO - DO j=jmin,jmax - DO i=imin,imax - POM_export(i,j) = 0. _d 0 - PFE_export(i,j) = 0. _d 0 - CaCO3_export(i,j) = 0. _d 0 - ENDDO - ENDDO - POMflux_u = 0. _d 0 - PFEflux_u = 0. _d 0 - CaCO3flux_u = 0. _d 0 - -c --------------------------------------------------------------------- -c Nutrients export/remineralization, CaCO3 export/dissolution -c -c The flux at the bottom of a grid cell equals -C Fb = (Ft + prod*dz) / (1 + zremin*dz) -C where Ft is the flux at the top, and prod*dz is the integrated -C production of new sinking particles within the layer. -C Ft = 0 in the first layer. - -CADJ STORE Fe_uptake = comlev1, key = ikey_dynamics - -C$TAF LOOP = parallel - DO j=jmin,jmax -C$TAF LOOP = parallel - DO i=imin,imax -C$TAF init upper_flux = static, Nr - DO k=1,Nr -C$TAF STORE POMflux_u = upper_flux -C$TAF STORE PFEflux_u = upper_flux -C$TAF STORE CaCO3flux_u = upper_flux - - IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN - -C Sinking speed is evaluated at the bottom of the cell - depth_l=-rF(k+1) - IF (depth_l .LE. wsink0z) THEN - wsink = wsink0 - ELSE - wsink = wsinkacc * (depth_l - wsink0z) + wsink0 - ENDIF - -C Nutrient remineralization lengthscale -C Not an e-folding scale: this term increases with remineralization. - zremin = gamma_POM * ( PTR_O2(i,j,k)**2 / - & (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min) - & + remin_min )/(wsink + epsln) - -C Calcium remineralization relaxed toward the inverse of the -C ca_remin_depth constant value as the calcite saturation approaches 0. - zremin_caco3 = 1. _d 0/ca_remin_depth*(1. _d 0-min(1. _d 0, - & omegaC(i,j,k,bi,bj))) - -C POM flux leaving the cell - POMflux_l = (POMflux_u+POM_prod(i,j,k)*drF(k) - & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k) - & *hFacC(i,j,k,bi,bj)) - -C CaCO3 flux leaving the cell - CaCO3flux_l = (caco3flux_u+CaCO3_prod(i,j,k)*drF(k) - & *hFacC(i,j,k,bi,bj))/(1+zremin_caco3*drF(k) - & *hFacC(i,j,k,bi,bj)) - -C Start with cells that are not the deepest cells - IF ((k.LT.Nr) .AND. (hFacC(i,j,k+1,bi,bj).GT.0)) THEN - -C Nutrient accumulation in a cell is given by the biological production -C (and instant remineralization) of particulate organic matter -C plus flux thought upper interface minus flux through lower interface. -C (Since not deepest cell: hFacC=1) - deltaPOM = (POMflux_u + POM_prod(i,j,k)*drF(k) - & - POMflux_l)*recip_drF(k) - - CaCO3_diss(i,j,k) = (CaCO3flux_u + CaCO3_prod(i,j,k)*drF(k) - & - CaCO3flux_l)*recip_drF(k) - - fe_sed_source = 0. _d 0 - - ELSE -C If this layer is adjacent to bottom topography or it is the deepest -C cell of the domain, then remineralize/dissolve in this grid cell -C i.e. don't subtract off lower boundary fluxes when calculating remin - deltaPOM = POMflux_u*recip_drF(k) - & *recip_hFacC(i,j,k,bi,bj)+POM_prod(i,j,k) - - CaCO3_diss(i,j,k) = caco3flux_u*recip_drF(k) - & *recip_hFacC(i,j,k,bi,bj)+CaCO3_prod(i,j,k) - -C Iron from sediments: the phosphate flux hitting the bottom boundary -C is used to scale the return of iron to the water column. -C Maximum value added for numerical stability. - fe_sed_source = min(1. _d -11, - & max(0. _d 0,FetoPsed/NUTfac*POMflux_l*recip_drF(k) - & *recip_hFacC(i,j,k,bi,bj))) -#ifdef BLING_ADJOINT_SAFE - fe_sed_source = 0. _d 0 -#endif - ENDIF - -C A fraction of POM becomes DOM - POM_diss(i,j,k) = deltaPOM*phi_DOM - POM_remin(i,j,k) = deltaPOM*(1-phi_DOM) - -C Begin iron uptake calculations by determining ligand bound and free iron. -C Both forms are available for biology, but only free iron is scavenged -C onto particles and forms colloids. - lig_stability = KFeLeq_max-(KFeLeq_max-KFeLeq_min) - & *(irr_inst(i,j,k,bi,bj)**2 - & /(IFeL**2+irr_inst(i,j,k,bi,bj)**2)) - & *max(0. _d 0,min(1. _d 0,(PTR_FE(i,j,k)-Fe_min)/ - & (PTR_FE(i,j,k)+epsln)*b_const)) - -C Use the quadratic equation to solve for binding between iron and ligands - - FreeFe = (-(1+lig_stability*(ligand-PTR_FE(i,j,k))) - & +((1+lig_stability*(ligand-PTR_FE(i,j,k)))**2+4* - & lig_stability*PTR_FE(i,j,k))**(0.5))/(2* - & lig_stability) - -C Iron scavenging doesn't occur in anoxic water (Fe2+ is soluble), so set -C FreeFe = 0 when anoxic. FreeFe should be interpreted the free iron that -C participates in scavenging. -#ifndef BLING_ADJOINT_SAFE - IF (PTR_O2(i,j,k) .LT. O2_min) THEN - FreeFe = 0 - ENDIF -#endif - -c Two mechanisms for iron uptake, in addition to biological production: -c colloidal scavenging and scavenging by organic matter. - -c Colloidal scavenging: -c Minimum function for numerical stability -c Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+ -c & min(0.5/PTRACERS_dTLev(1), kFe_inorg*FreeFe**(0.5))*FreeFe - - Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+ - & kFe_inorg*FreeFe**(0.5)*FreeFe - -C Scavenging of iron by organic matter: -c The POM value used is the bottom boundary flux. This doesn't occur in -c oxic waters, but FreeFe is set to 0 in such waters earlier. - IF ( POMflux_l .GT. 0. _d 0 ) THEN - -c Minimum function for numerical stability -c Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+ -c & min(0.5/PTRACERS_dTLev(1), kFE_org*(POMflux_l -c & *CtoP/NUTfac*12.01/wsink)**(0.58)*FreeFe - -#ifndef BLING_ADJOINT_SAFE - Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+ - & kFE_org*(POMflux_l*CtoP/NUTfac - & *12.01/wsink)**(0.58)*FreeFe -#else - Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+ - & kFE_org*(POMflux_l*CtoP/NUTfac - & *12.01/wsink0)**(0.58)*FreeFe -#endif - ENDIF - -C If water is oxic then the iron is remineralized normally. Otherwise -C it is completely remineralized (fe 2+ is soluble, but unstable -C in oxidizing environments). - - pfeflux_l = (pfeflux_u+Fe_uptake(i,j,k)*drF(k) - & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k) - & *hFacC(i,j,k,bi,bj)) - -#ifndef BLING_ADJOINT_SAFE - IF ( PTR_O2(i,j,k) .LT. O2_min ) THEN - pfeflux_l = 0 - ENDIF -#endif - - Fe_remin(i,j,k) = (pfeflux_u+Fe_uptake(i,j,k)*drF(k) - & *hFacC(i,j,k,bi,bj)-pfeflux_l)*recip_drF(k) - & *recip_hFacC(i,j,k,bi,bj) - -C Add sediment source - Fe_remin(i,j,k) = Fe_remin(i,j,k) + fe_sed_source - -C Prepare the tracers for the next layer down - POMflux_u = POMflux_l - PFEflux_u = PFEflux_l - CaCO3flux_u = CaCO3flux_l - -C Depth-integrated export (through bottom of water column) -C This is calculated last for the deepest cell - POM_export(i,j) = POMflux_l - PFE_export(i,j) = PFEflux_l - CACO3_export(i,j) = CaCO3flux_l - - ENDIF - - ENDDO - -C Reset for next location (i,j) - POMflux_u = 0. _d 0 - PFEflux_u = 0. _d 0 - CaCO3flux_u = 0. _d 0 - - ENDDO - ENDDO - - RETURN - END +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/bling/pkg/bling_remineralization.F,v 1.2 2016/02/28 21:49:24 mmazloff Exp $ +C $Name: $ + +#include "BLING_OPTIONS.h" + +CBOP + subroutine BLING_REMIN( + I PTR_NO3, PTR_FE, PTR_O2, irr_inst, + I N_spm, P_spm, Fe_spm, CaCO3_uptake, + O N_reminp, P_reminp, Fe_reminsum, + O N_den_benthic, CaCO3_diss, + I bi, bj, imin, imax, jmin, jmax, + I myIter, myTime, myThid ) + +C ================================================================= +C | subroutine bling_remin +C | o Organic matter export and remineralization +C | - Sinking particulate flux and diel migration contribute to +C | export. +C | - Denitrification xxx +C | o Sediments +C ================================================================= + + implicit none + +C === Global variables === + +#include "SIZE.h" +#include "DYNVARS.h" +#include "EEPARAMS.h" +#include "PARAMS.h" +#include "GRID.h" +#include "BLING_VARS.h" +#include "PTRACERS_SIZE.h" +#include "PTRACERS_PARAMS.h" +#ifdef ALLOW_AUTODIFF +# include "tamc.h" +#endif + +C === Routine arguments === +C bi,bj :: tile indices +C iMin,iMax :: computation domain: 1rst index range +C jMin,jMax :: computation domain: 2nd index range +C myTime :: current time +C myIter :: current timestep +C myThid :: thread Id. number + INTEGER bi, bj, imin, imax, jmin, jmax + _RL myTime + INTEGER myIter + INTEGER myThid +C === Input === +C PTR_NO3 :: nitrate concentration +C PTR_FE :: iron concentration +C PTR_O2 :: oxygen concentration + _RL PTR_NO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL PTR_FE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL PTR_O2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL irr_inst(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL N_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL P_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL Fe_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL CaCO3_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) +C === Output === +C + _RL N_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL N_den_benthic(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL P_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL Fe_reminsum(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL CaCO3_diss(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + +#ifdef ALLOW_BLING +C === Local variables === +C i,j,k :: loop indices +C irr_eff :: effective irradiance +C NO3_lim :: nitrate limitation +C PO4_lim :: phosphate limitation +C Fe_lim :: iron limitation for phytoplankton +C Fe_lim_diaz :: iron limitation for diazotrophs +C alpha_Fe :: initial slope of the P-I curve +C theta_Fe :: Chl:C ratio +C theta_Fe_max :: Fe-replete maximum Chl:C ratio +C irrk :: nut-limited efficiency of algal photosystems +C Pc_m :: light-saturated max photosynthesis rate for phyt +C Pc_m_diaz :: light-saturated max photosynthesis rate for diaz +C Pc_tot :: carbon-specific photosynthesis rate +C expkT :: temperature function +C mu :: net carbon-specific growth rate for phyt +C mu_diaz :: net carbon-specific growth rate for diaz +C biomass_sm :: nitrogen concentration in small phyto biomass +C biomass_lg :: nitrogen concentration in large phyto biomass +C N_uptake :: nitrogen uptake +C N_fix :: nitrogen fixation +C P_uptake :: phosphorus uptake +C POC_flux :: carbon export flux 3d field +C chl :: chlorophyll diagnostic +C PtoN :: variable ratio of phosphorus to nitrogen +C FetoN :: variable ratio of iron to nitrogen +C N_spm :: particulate sinking of nitrogen +C P_spm :: particulate sinking of phosphorus +C Fe_spm :: particulate sinking of iron +C N_dvm :: vertical transport of nitrogen by DVM +C P_dvm :: vertical transport of phosphorus by DVM +C Fe_dvm :: vertical transport of iron by DVM +C N_recycle :: recycling of newly-produced organic nitrogen +C P_recycle :: recycling of newly-produced organic phosphorus +C Fe_recycle :: recycling of newly-produced organic iron +c xxx to be completed + INTEGER i,j,k + _RL PONflux_u + _RL POPflux_u + _RL PFEflux_u + _RL CaCO3flux_u + _RL PONflux_l + _RL POPflux_l + _RL PFEflux_l + _RL CaCO3flux_l + _RL depth_l + _RL zremin + _RL zremin_caco3 + _RL wsink + _RL POC_sed + _RL Fe_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL NO3_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL PO4_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL O2_sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL lig_stability + _RL FreeFe + _RL Fe_ads_inorg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL Fe_ads_org(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL log_btm_flx + _RL Fe_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL o2_upper + _RL o2_lower + _RL dz_upper + _RL dz_lower + _RL temp_upper + _RL temp_lower + _RL z_dvm_regr + _RL frac_migr + _RL fdvm_migr + _RL fdvm_stat + _RL fdvmn_vint + _RL fdvmp_vint + _RL fdvmfe_vint + _RL z_dvm + _RL N_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL P_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL Fe_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL mld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL Fe_burial(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL x_erfcc,z_erfcc,t_erfcc,erfcc +cxx order +CEOP + +c --------------------------------------------------------------------- +c Initialize output and diagnostics + DO k=1,Nr + DO j=jmin,jmax + DO i=imin,imax + Fe_ads_org(i,j,k) = 0. _d 0 + Fe_ads_inorg(i,j,k) = 0. _d 0 + N_reminp(i,j,k) = 0. _d 0 + P_reminp(i,j,k) = 0. _d 0 + Fe_reminp(i,j,k) = 0. _d 0 + Fe_reminsum(i,j,k) = 0. _d 0 + N_remindvm(i,j,k) = 0. _d 0 + P_remindvm(i,j,k) = 0. _d 0 + Fe_remindvm(i,j,k) = 0. _d 0 + N_den_benthic(i,j,k)= 0. _d 0 + CaCO3_diss(i,j,k) = 0. _d 0 + ENDDO + ENDDO + ENDDO + DO j=jmin,jmax + DO i=imin,imax + Fe_burial(i,j) = 0. _d 0 + NO3_sed(i,j) = 0. _d 0 + PO4_sed(i,j) = 0. _d 0 + O2_sed(i,j) = 0. _d 0 + ENDDO + ENDDO + +c --------------------------------------------------------------------- +c Remineralization + +CADJ STORE Fe_ads_org = comlev1, key = ikey_dynamics +cxx needed? + +C$TAF LOOP = parallel + DO j=jmin,jmax +C$TAF LOOP = parallel + DO i=imin,imax +cmm C$TAF init upper_flux = static, Nr + +C Initialize upper flux + PONflux_u = 0. _d 0 + POPflux_u = 0. _d 0 + PFEflux_u = 0. _d 0 + CaCO3flux_u = 0. _d 0 + + DO k=1,Nr +c C$TAF STORE PONflux_u = upper_flux +c C$TAF STORE POPflux_u = upper_flux +c C$TAF STORE PFEflux_u = upper_flux +c C$TAF STORE CaCO3flux_u = upper_flux +CADJ STORE PONflux_u, POPflux_u, PFEflux_u, CaCO3flux_u = +CADJ & comlev1, key = ikey_dynamics, kind = isbyte +CADJ STORE Fe_ads_org = +CADJ & comlev1, key = ikey_dynamics, kind = isbyte +CMM) + + IF ( hFacC(i,j,k,bi,bj).gt.0. _d 0 ) THEN + +C Sinking speed is evaluated at the bottom of the cell + depth_l=-rF(k+1) + IF (depth_l .LE. wsink0z) THEN + wsink = wsink0 + ELSE + wsink = wsinkacc * (depth_l - wsink0z) + wsink0 + ENDIF + +C Nutrient remineralization lengthscale +C Not an e-folding scale: this term increases with remineralization. + zremin = gamma_POM * ( PTR_O2(i,j,k)**2 / + & (k_O2**2 + PTR_O2(i,j,k)**2) * (1-remin_min) + & + remin_min )/(wsink + epsln) + +C Calcium remineralization relaxed toward the inverse of the +C ca_remin_depth constant value as the calcite saturation approaches 0. + zremin_caco3 = 1. _d 0/ca_remin_depth*(1. _d 0 - min(1. _d 0, + & omegaC(i,j,k,bi,bj) + epsln )) + + +C POM flux leaving the cell + PONflux_l = (PONflux_u+N_spm(i,j,k)*drF(k) + & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k) + & *hFacC(i,j,k,bi,bj)) +C!! multiply by intercept_frac ??? + + POPflux_l = (POPflux_u+P_spm(i,j,k)*drF(k) + & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k) + & *hFacC(i,j,k,bi,bj)) +C!! multiply by intercept_frac ??? + +C CaCO3 flux leaving the cell + CaCO3flux_l = (caco3flux_u+CaCO3_uptake(i,j,k)*drF(k) + & *hFacC(i,j,k,bi,bj))/(1+zremin_caco3*drF(k) + & *hFacC(i,j,k,bi,bj)) +C!! multiply by intercept_frac ??? + + +C Start with cells that are not the deepest cells + IF ((k.LT.Nr) .AND. (hFacC(i,j,k+1,bi,bj).GT.0)) THEN + +C Nutrient accumulation in a cell is given by the biological production +C (and instant remineralization) of particulate organic matter +C plus flux thought upper interface minus flux through lower interface. +C (Since not deepest cell: hFacC=1) + N_reminp(i,j,k) = (PONflux_u + N_spm(i,j,k)*drF(k) + & - PONflux_l)*recip_drF(k) + + P_reminp(i,j,k) = (POPflux_u + P_spm(i,j,k)*drF(k) + & - POPflux_l)*recip_drF(k) + + CaCO3_diss(i,j,k) = (CaCO3flux_u + CaCO3_uptake(i,j,k) + & *drF(k) - CaCO3flux_l)*recip_drF(k) + + Fe_sed(i,j,k) = 0. _d 0 + + + ELSE +C If this layer is adjacent to bottom topography or it is the deepest +C cell of the domain, then remineralize/dissolve in this grid cell +C i.e. don't subtract off lower boundary fluxes when calculating remin + + N_reminp(i,j,k) = PONflux_u*recip_drF(k) + & *recip_hFacC(i,j,k,bi,bj) + N_spm(i,j,k) + + P_reminp(i,j,k) = POPflux_u*recip_drF(k) + & *recip_hFacC(i,j,k,bi,bj) + P_spm(i,j,k) + + CaCO3_diss(i,j,k) = CaCO3flux_u*recip_drF(k) + & *recip_hFacC(i,j,k,bi,bj) + CaCO3_uptake(i,j,k) + + +c Efflux Fed out of sediments +C The phosphate flux hitting the bottom boundary +C is used to scale the return of iron to the water column. +C Maximum value added for numerical stability. + + POC_sed = PONflux_l * CtoN + + Fe_sed(i,j,k) = min(1. _d -11, + & max(epsln, FetoC_sed * POC_sed * recip_drF(k) + & *recip_hFacC(i,j,k,bi,bj))) + +#ifdef BLING_ADJOINT_SAFE + Fe_sed(i,j,k) = 0. _d 0 +#endif + + +cav temporary until I figure out why this is problematic for adjoint +#ifndef BLING_ADJOINT_SAFE + +#ifndef USE_SGS_SED +c Calculate benthic denitrification and Fe efflux here, if the subgridscale sediment module is not being used. + + IF (POC_sed .gt. 0. _d 0) THEN + + log_btm_flx = 0. _d 0 + +c Convert from mol N m-2 s-1 to umol C cm-2 d-1 and take the log + + log_btm_flx = log10(min(43.0 _d 0, POC_sed * + & 86400. _d 0 * 100.0 _d 0)) + +c Metamodel gives units of umol C cm-2 d-1, convert to mol N m-2 s-1 and multiply by +c no3_2_n to give NO3 consumption rate + + N_den_benthic(i,j,k) = min (POC_sed * NO3toN / CtoN, + & (10 _d 0)**(-0.9543 _d 0 + 0.7662 _d 0 * + & log_btm_flx - 0.235 _d 0 * log_btm_flx * log_btm_flx) + & / (CtoN * 86400. _d 0 * 100.0 _d 0) * NO3toN * + & PTR_NO3(i,j,k) / (k_no3 + PTR_NO3(i,j,k)) ) * + & recip_drF(k) + + endif + +#endif + +#endif + +c --------------------------------------------------------------------- +c Calculate external bottom fluxes for tracer_vertdiff. Positive fluxes are into the water +c column from the seafloor. For P, the bottom flux puts the sinking flux reaching the bottom +c cell into the water column through diffusion. For iron, the sinking flux disappears into the +c sediments if bottom waters are oxic (assumed adsorbed as oxides). If bottom waters are anoxic, +c the sinking flux of Fe is returned to the water column. +c +c For oxygen, the consumption of oxidant required to respire +c the settling flux of organic matter (in support of the +c no3 bottom flux) diffuses from the bottom water into the sediment. + +c Assume all NO3 for benthic denitrification is supplied from the bottom water, and that +c all organic N is also consumed under denitrification (Complete Denitrification, sensu +c Paulmier, Biogeosciences 2009). Therefore, no NO3 is regenerated from organic matter +c respired by benthic denitrification (necessitating the second term in b_no3). + + NO3_sed(i,j) = PONflux_l*drF(k)*hFacC(i,j,k,bi,bj) + & - N_den_benthic(i,j,k) / NO3toN + + PO4_sed(i,j) = POPflux_l*drF(k)*hFacC(i,j,k,bi,bj) + +c Oxygen flux into sediments is that required to support non-denitrification respiration, +c assuming a 4/5 oxidant ratio of O2 to NO3. Oxygen consumption is allowed to continue +c at negative oxygen concentrations, representing sulphate reduction. + + O2_sed(i,j) = -(O2toN * PONflux_l*drF(k)*hFacC(i,j,k,bi,bj) + & - N_den_benthic(i,j,k)* 1.25) + + ENDIF + + +C Begin iron uptake calculations by determining ligand bound and free iron. +C Both forms are available for biology, but only free iron is scavenged +C onto particles and forms colloids. + + lig_stability = kFe_eq_lig_max-(KFe_eq_lig_max-kFe_eq_lig_min) + & *(irr_inst(i,j,k)**2 + & /(kFe_eq_lig_irr**2+irr_inst(i,j,k)**2)) + & *max(epsln,min(1. _d 0,(PTR_FE(i,j,k) + & -kFe_eq_lig_Femin)/ + & (PTR_FE(i,j,k)+epsln)*1.2 _d 0)) + +C Use the quadratic equation to solve for binding between iron and ligands + + FreeFe = (-(1+lig_stability*(ligand-PTR_FE(i,j,k))) + & +((1+lig_stability*(ligand-PTR_FE(i,j,k)))**2+4* + & lig_stability*PTR_FE(i,j,k))**(0.5))/(2* + & lig_stability) + +C Iron scavenging doesn't occur in anoxic water (Fe2+ is soluble), so set +C FreeFe = 0 when anoxic. FreeFe should be interpreted the free iron that +C participates in scavenging. + +#ifndef BLING_ADJOINT_SAFE + IF (PTR_O2(i,j,k) .LT. oxic_min) THEN + FreeFe = 0. _d 0 + ENDIF +#endif + +c Two mechanisms for iron uptake, in addition to biological production: +c colloidal scavenging and scavenging by organic matter. + +c Colloidal scavenging: +c Minimum function for numerical stability +c Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+ +c & min(0.5/PTRACERS_dTLev(1), kFe_inorg*FreeFe**(0.5))*FreeFe + + Fe_ads_inorg(i,j,k) = + & kFe_inorg*(max(1. _d -8,FreeFe))**(1.5) + +C Scavenging of iron by organic matter: +c The POM value used is the bottom boundary flux. This doesn't occur in +c oxic waters, but FreeFe is set to 0 in such waters earlier. + IF ( PONflux_l .GT. 0. _d 0 ) THEN + +c Minimum function for numerical stability +c Fe_uptake(i,j,k) = Fe_uptake(i,j,k)+ +c & min(0.5/PTRACERS_dTLev(1), kFE_org*(POMflux_l +c & *CtoP/NUTfac*12.01/wsink)**(0.58)*FreeFe + +#ifndef BLING_ADJOINT_SAFE + Fe_ads_org(i,j,k) = + & kFE_org*(PONflux_l/(epsln + wsink) + & * MasstoN)**(0.58)*FreeFe +#else + Fe_ads_org(i,j,k) = + & kFE_org*(PONflux_l/(epsln + wsink0) + & * MasstoN)**(0.58)*FreeFe +#endif + ENDIF + + + +C If water is oxic then the iron is remineralized normally. Otherwise +C it is completely remineralized (fe 2+ is soluble, but unstable +C in oxidizing environments). + + PFEflux_l = (PFEflux_u+(Fe_spm(i,j,k)+Fe_ads_inorg(i,j,k) + & +Fe_ads_org(i,j,k))*drF(k) + & *hFacC(i,j,k,bi,bj))/(1+zremin*drF(k) + & *hFacC(i,j,k,bi,bj)) + + +c Added the burial flux of sinking particulate iron here as a +c diagnostic, needed to calculate mass balance of iron. +c this is calculated last for the deepest cell + + Fe_burial(i,j) = PFeflux_l + + +#ifndef BLING_ADJOINT_SAFE + IF ( PTR_O2(i,j,k) .LT. oxic_min ) THEN + pfeflux_l = 0 + ENDIF +#endif + + Fe_reminp(i,j,k) = (pfeflux_u+(Fe_spm(i,j,k) + & +Fe_ads_inorg(i,j,k) + & +Fe_ads_org(i,j,k))*drF(k) + & *hFacC(i,j,k,bi,bj)-pfeflux_l)*recip_drF(k) + & *recip_hFacC(i,j,k,bi,bj) +C!! there's an intercept_frac here... need to add + + +C Prepare the tracers for the next layer down + PONflux_u = PONflux_l + POPflux_u = POPflux_l + PFEflux_u = PFEflux_l + CaCO3flux_u = CaCO3flux_l + +c + Fe_reminsum(i,j,k) = Fe_reminp(i,j,k) + Fe_sed(i,j,k) + & - Fe_ads_org(i,j,k) - Fe_ads_inorg(i,j,k) +cc Fe_reminsum(i,j,k) = 0. _d 0 + + ENDIF + + ENDDO + ENDDO + ENDDO + +CADJ STORE Fe_ads_org = comlev1, key = ikey_dynamics +cxx needed? + + +c --------------------------------------------------------------------- + +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + +c 3d local variables +c CALL DIAGNOSTICS_FILL(POC_flux,'BLGPOCF ',0,Nr,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(Fe_ads_inorg,'BLGFEAI',0,Nr,2,bi,bj, + & myThid) + CALL DIAGNOSTICS_FILL(Fe_sed,'BLGFESED',0,Nr,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(Fe_reminp,'BLGFEREM',0,Nr,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(N_den_benthic,'BLGNDENB',0,Nr,2,bi,bj, + & myThid) +c CALL DIAGNOSTICS_FILL(N_den_pelag,'BLGNDENP',0,Nr,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(N_reminp,'BLGNREM ',0,Nr,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(P_reminp,'BLGPREM ',0,Nr,2,bi,bj,myThid) +c 2d local variables + CALL DIAGNOSTICS_FILL(Fe_burial,'BLGFEBUR',0,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(NO3_sed,'BLGNSED ',0,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(PO4_sed,'BLGPSED ',0,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(O2_sed,'BLGO2SED',0,1,2,bi,bj,myThid) +c these variables are currently 1d, could be 3d for diagnostics +c (or diag_fill could be called inside loop - which is faster?) +c CALL DIAGNOSTICS_FILL(zremin,'BLGZREM ',0,Nr,2,bi,bj,myThid) + + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + + + +#endif /* ALLOW_BLING */ + + RETURN + END
| ViewVC Help | |
| Powered by ViewVC 1.1.22 |