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