ecearth3/sources/nemo-3.6/CONFIG/ORCA1L75_LIM3_PISCES_standalone/MY_SRC/p4zfechem.F90

MODULE p4zfechem
   !!======================================================================
   !!                         ***  MODULE p4zfechem  ***
   !! TOP :   PISCES Compute iron chemistry and scavenging
   !!======================================================================
   !! History :   3.5  !  2012-07 (O. Aumont, A. Tagliabue, C. Ethe) Original code
   !!----------------------------------------------------------------------
#if defined key_pisces
   !!----------------------------------------------------------------------
   !!   'key_top'       and                                      TOP models
   !!   'key_pisces'                                       PISCES bio-model
   !!----------------------------------------------------------------------
   !!   p4z_fechem       :  Compute remineralization/scavenging of iron
   !!   p4z_fechem_init  :  Initialisation of parameters for remineralisation
   !!   p4z_fechem_alloc :  Allocate remineralisation variables
   !!----------------------------------------------------------------------
   USE oce_trc         !  shared variables between ocean and passive tracers
   USE trc             !  passive tracers common variables 
   USE sms_pisces      !  PISCES Source Minus Sink variables
   USE p4zopt          !  optical model
   USE p4zche          !  chemical model
   USE p4zsbc          !  Boundary conditions from sediments
   USE prtctl_trc      !  print control for debugging
   USE iom             !  I/O manager

   IMPLICIT NONE
   PRIVATE

   PUBLIC   p4z_fechem      ! called in p4zbio.F90
   PUBLIC   p4z_fechem_init ! called in trcsms_pisces.F90

   !! * Shared module variables
   LOGICAL          ::  ln_fechem    !: boolean for complex iron chemistry following Tagliabue and voelker
   LOGICAL          ::  ln_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker
   REAL(wp), PUBLIC ::  xlam1        !: scavenging rate of Iron 
   REAL(wp), PUBLIC ::  xlamdust     !: scavenging rate of Iron by dust 
   REAL(wp), PUBLIC ::  ligand       !: ligand concentration in the ocean 

   REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth

   !!* Substitution
#  include "top_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
   !! $Id: p4zrem.F90 3160 2011-11-20 14:27:18Z cetlod $ 
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE p4z_fechem( kt, knt )
      !!---------------------------------------------------------------------
      !!                     ***  ROUTINE p4z_fechem  ***
      !!
      !! ** Purpose :   Compute remineralization/scavenging of iron
      !!
      !! ** Method  :   2 different chemistry models are available for iron
      !!                (1) The simple chemistry model of Aumont and Bopp (2006)
      !!                    based on one ligand and one inorganic form
      !!                (2) The complex chemistry model of Tagliabue and 
      !!                    Voelker (2009) based on 2 ligands, 2 inorganic forms
      !!                    and one particulate form (ln_fechem)
      !!---------------------------------------------------------------------
      !
      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
      !
      INTEGER  ::   ji, jj, jk, jic
      REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll
      REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag
      REAL(wp) ::   ztrc, zdust
#if ! defined key_kriest
      REAL(wp) ::   zdenom, zdenom2
#endif
      REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig
      REAL(wp), POINTER, DIMENSION(:,:,:) :: zFeL1, zFeL2, zTL2, zFe2, zFeP
      REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand
      REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2
      REAL(wp) :: zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2
      REAL(wp) :: ztfe, zoxy
      REAL(wp) :: zstep
      CHARACTER (len=25) :: charout
      !!---------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('p4z_fechem')
      !
      ! Allocate temporary workspace
      CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )
      zFe3 (:,:,:) = 0.
      zFeL1(:,:,:) = 0.
      zTL1 (:,:,:) = 0.
      IF( ln_fechem ) THEN
         CALL wrk_alloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
         zFe2 (:,:,:) = 0.
         zFeL2(:,:,:) = 0.
         zTL2 (:,:,:) = 0.
         zFeP (:,:,:) = 0.
      ENDIF
      !Initialization of zxs to zero to avoid SEG FAULT if it's not assigned any value at first time step - Raffa July 2018
      zxs = 0.
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! Total ligand concentration : Ligands can be chosen to be constant or variable
      ! Parameterization from Tagliabue and Voelker (2011)
      ! -------------------------------------------------
      IF( ln_ligvar ) THEN
         ztotlig(:,:,:) =  0.09 * trb(:,:,:,jpdoc) * 1E6 + ligand * 1E9
         ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
      ELSE
         ztotlig(:,:,:) = ligand * 1E9
      ENDIF

      IF( ln_fechem ) THEN
         ! ------------------------------------------------------------
         ! NEW FE CHEMISTRY ROUTINE from Tagliabue and Volker (2009)
         ! This model is based on two ligands, Fe2+, Fe3+ and Fep
         ! Chemistry is supposed to be fast enough to be at equilibrium
         ! ------------------------------------------------------------
!CDIR NOVERRCHK
         DO jk = 1, jpkm1
!CDIR NOVERRCHK
            DO jj = 1, jpj
!CDIR NOVERRCHK
               DO ji = 1, jpi
                  ! Calculate ligand concentrations : assume 2/3rd of excess goes to
                  ! strong ligands (L1) and 1/3rd to weak ligands (L2)
                  ztligand       = ztotlig(ji,jj,jk) - ligand * 1E9
                  zTL1(ji,jj,jk) =                0.000001 + 0.67 * ztligand
                  zTL2(ji,jj,jk) = ligand * 1E9 - 0.000001 + 0.33 * ztligand
                  ! ionic strength from Millero et al. 1987
                  zionic = 19.9201 * tsn(ji,jj,jk,jp_sal) / ( 1000. - 1.00488 * tsn(ji,jj,jk,jp_sal) + rtrn )
                  zph    = -LOG10( MAX( hi(ji,jj,jk), rtrn) )
                  zoxy   = trb(ji,jj,jk,jpoxy) * ( rhop(ji,jj,jk) / 1.e3 )
                  ! Fe2+ oxydation rate from Santana-Casiano et al. (2005)
                  zkox   = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tsn(ji,jj,jk,jp_tem) + 273.15 )  &
                    &    - 0.04406 * SQRT( tsn(ji,jj,jk,jp_sal) ) - 0.002847 * tsn(ji,jj,jk,jp_sal)
                  zkox   = ( 10.** zkox ) * spd
                  zkox   = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn )
                  ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006)
                  zkph2 = MAX( 0., 15. * etot(ji,jj,jk) / ( etot(ji,jj,jk) + 2. ) )
                  zkph1 = zkph2 / 5.
                  ! pass the dfe concentration from PISCES
                  ztfe = trb(ji,jj,jk,jpfer) * 1e9
                  ! ----------------------------------------------------------
                  ! ANALYTICAL SOLUTION OF ROOTS OF THE FE3+ EQUATION
                  ! As shown in Tagliabue and Voelker (2009), Fe3+ is the root of a 3rd order polynom. 
                  ! ----------------------------------------------------------
                  ! calculate some parameters
                  za = 1 + ks / kpr
                  zb = 1 + ( zkph1 + kth ) / ( zkox + rtrn )
                  zc = 1 + zkph2 / ( zkox + rtrn )
                  zkappa1 = ( kb1 + zkph1 + kth ) / kl1
                  zkappa2 = ( kb2 + zkph2 ) / kl2
                  za2 = zTL1(ji,jj,jk) * zb / za + zTL2(ji,jj,jk) * zc / za + zkappa1 + zkappa2 - ztfe / za
                  za1 = zkappa2 * zTL1(ji,jj,jk) * zb / za + zkappa1 * zTL2(ji,jj,jk) * zc / za &
                      & + zkappa1 * zkappa2 - ( zkappa1 + zkappa2 ) * ztfe / za
                  za0 = -zkappa1 * zkappa2 * ztfe / za
                  zp  = za1 - za2 * za2 / 3.
                  zq  = za2 * za2 * za2 * 2. / 27. - za2 * za1 / 3. + za0
                  zp3 = zp / 3.
                  zq2 = zq / 2.
                  zd  = zp3 * zp3 * zp3 + zq2 * zq2
                  zr  = zq / ABS( zq ) * SQRT( ABS( zp ) / 3. )
                  ! compute the roots
                  IF( zp > 0.) THEN
                     ! zphi = ASINH( zq / ( 2. * zr * zr * zr ) )
                     zphi =  zq / ( 2. * zr * zr * zr ) 
                     zphi = LOG( zphi + SQRT( zphi * zphi + 1 ) )  ! asinh(x) = log(x + sqrt(x^2+1))
                     zxs  = -2. * zr * SINH( zphi / 3. ) - za1 / 3.
                  ELSE
                     IF( zd > 0. ) THEN
                        zfff = MAX( 1., zq / ( 2. * zr * zr * zr ) )
                        ! zphi = ACOSH( zfff )
                        zphi = LOG( zfff + SQRT( zfff * zfff - 1 ) )  ! acosh(x) = log(x + sqrt(x^2-1))
                        zxs = -2. * zr * COSH( zphi / 3. ) - za1 / 3.
                     ELSE
                        zfff = MIN( 1., zq / ( 2. * zr * zr * zr ) )
                        zphi = ACOS( zfff )
                        DO jic = 1, 3
                           zfunc = -2 * zr * COS( zphi / 3. + 2. * FLOAT( jic - 1 ) * rpi / 3. ) - za2 / 3.
                           IF( zfunc > 0. .AND. zfunc <= ztfe)  zxs = zfunc
                        END DO
                     ENDIF
                  ENDIF
                  ! solve for the other Fe species
                  zFe3(ji,jj,jk) = MAX( 0., zxs ) 
                  zFep(ji,jj,jk) = MAX( 0., ( ks * zFe3(ji,jj,jk) / kpr ) )
                  zkappa2 = ( kb2 + zkph2 ) / kl2
                  zFeL2(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * zTL2(ji,jj,jk) ) / ( zkappa2 + zFe3(ji,jj,jk) ) )
                  zFeL1(ji,jj,jk) = MAX( 0., ( ztfe / zb - za / zb * zFe3(ji,jj,jk) - zc / zb * zFeL2(ji,jj,jk) ) )
                  zFe2 (ji,jj,jk) = MAX( 0., ( ( zkph1 * zFeL1(ji,jj,jk) + zkph2 * zFeL2(ji,jj,jk) ) / zkox ) )
               END DO
            END DO
         END DO
      ELSE
         ! ------------------------------------------------------------
         ! OLD FE CHEMISTRY ROUTINE from Aumont and Bopp (2006)
         ! This model is based on one ligand and Fe' 
         ! Chemistry is supposed to be fast enough to be at equilibrium
         ! ------------------------------------------------------------
!CDIR NOVERRCHK
         DO jk = 1, jpkm1
!CDIR NOVERRCHK
            DO jj = 1, jpj
!CDIR NOVERRCHK
               DO ji = 1, jpi
                  zTL1(ji,jj,jk) = ztotlig(ji,jj,jk)
                  zkeq           = fekeq(ji,jj,jk)
                  zfesatur       = zTL1(ji,jj,jk) * 1E-9
                  ztfe           = trb(ji,jj,jk,jpfer) 
                  ! Fe' is the root of a 2nd order polynom
                  zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
                     &             + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       &
                     &               + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
                  zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
                  zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
              END DO
            END DO
         END DO
         !
      ENDIF

      zdust = 0.         ! if no dust available
!CDIR NOVERRCHK
      DO jk = 1, jpkm1
!CDIR NOVERRCHK
         DO jj = 1, jpj
!CDIR NOVERRCHK
            DO ji = 1, jpi
               zstep = xstep
# if defined key_degrad
               zstep = zstep * facvol(ji,jj,jk)
# endif
               ! Scavenging rate of iron. This scavenging rate depends on the load of particles of sea water. 
               ! This parameterization assumes a simple second order kinetics (k[Particles][Fe]).
               ! Scavenging onto dust is also included as evidenced from the DUNE experiments.
               ! --------------------------------------------------------------------------------------
               IF( ln_fechem ) THEN
                  zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) + zFeP(ji,jj,jk) ) * 1E-9
                  zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
               ELSE
                  zfeequi = zFe3(ji,jj,jk) * 1E-9 
                  zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
               ENDIF
#if defined key_kriest
               ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
#else
               ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
#endif
               IF( ln_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk) ! dust in kg/m2/s
               zlam1b = 3.e-5 + xlamdust * zdust + xlam1 * ztrc
               zscave = zfeequi * zlam1b * zstep

               ! Compute the different ratios for scavenging of iron
               ! to later allocate scavenged iron to the different organic pools
               ! ---------------------------------------------------------
               zdenom1 = xlam1 * trb(ji,jj,jk,jppoc) / zlam1b
#if ! defined key_kriest
               zdenom2 = xlam1 * trb(ji,jj,jk,jpgoc) / zlam1b
#endif

               !  Increased scavenging for very high iron concentrations found near the coasts 
               !  due to increased lithogenic particles and let say it is unknown processes (precipitation, ...)
               !  -----------------------------------------------------------
               zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
               zlamfac = MIN( 1.  , zlamfac )
               zdep    = MIN( 1., 1000. / fsdept(ji,jj,jk) )
               zlam1b  = xlam1 * MAX( 0.e0, ( trb(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) )
               zcoag   = zfeequi * zlam1b * zstep + 1E-4 * ( 1. - zlamfac ) * zdep * zstep * trb(ji,jj,jk,jpfer)

               !  Compute the coagulation of colloidal iron. This parameterization 
               !  could be thought as an equivalent of colloidal pumping.
               !  It requires certainly some more work as it is very poorly constrained.
               !  ----------------------------------------------------------------
               zlam1a  = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
                   &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) + 5.09E3 * trb(ji,jj,jk,jppoc) )
               zaggdfea = zlam1a * zstep * zfecoll
#if defined key_kriest
               zaggdfeb = 0.
               !
               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag
               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea + zaggdfeb
#else
               zlam1b = 3.53E3 *   trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
               zaggdfeb = zlam1b * zstep * zfecoll
               !
               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag
               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea
               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
#endif
            END DO
         END DO
      END DO
      !
      !  Define the bioavailable fraction of iron
      !  ----------------------------------------
      IF( ln_fechem ) THEN
          biron(:,:,:) = MAX( 0., trb(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 )
      ELSE
          biron(:,:,:) = trb(:,:,:,jpfer) 
      ENDIF

      !  Output of some diagnostics variables
      !     ---------------------------------
      IF( lk_iomput .AND. knt == nrdttrc ) THEN
         IF( iom_use("Fe3")    )  CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
         IF( iom_use("FeL1")   )  CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
         IF( iom_use("TL1")    )  CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
         IF( iom_use("Totlig") )  CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
         IF( iom_use("Biron")  )  CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! biron
         IF( ln_fechem ) THEN
            IF( iom_use("Fe2")  ) CALL iom_put("Fe2"    , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
            IF( iom_use("FeL2") ) CALL iom_put("FeL2"   , zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
            IF( iom_use("FeP")  ) CALL iom_put("FeP"    , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
            IF( iom_use("TL2")  ) CALL iom_put("TL2"    , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
         ENDIF
      ENDIF

      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
         WRITE(charout, FMT="('fechem')")
         CALL prt_ctl_trc_info(charout)
         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
      ENDIF
      !
                       CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )
      IF( ln_fechem )  CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
      !
      IF( nn_timing == 1 )  CALL timing_stop('p4z_fechem')
      !
   END SUBROUTINE p4z_fechem


   SUBROUTINE p4z_fechem_init
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE p4z_fechem_init  ***
      !!
      !! ** Purpose :   Initialization of iron chemistry parameters
      !!
      !! ** Method  :   Read the nampisfer namelist and check the parameters
      !!      called at the first timestep
      !!
      !! ** input   :   Namelist nampisfer
      !!
      !!----------------------------------------------------------------------
      NAMELIST/nampisfer/ ln_fechem, ln_ligvar, xlam1, xlamdust, ligand 
      INTEGER :: ios                 ! Local integer output status for namelist read

      REWIND( numnatp_ref )              ! Namelist nampisfer in reference namelist : Pisces iron chemistry
      READ  ( numnatp_ref, nampisfer, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisfer in reference namelist', lwp )

      REWIND( numnatp_cfg )              ! Namelist nampisfer in configuration namelist : Pisces iron chemistry
      READ  ( numnatp_cfg, nampisfer, IOSTAT = ios, ERR = 902 )
902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisfer in configuration namelist', lwp )
      IF(lwm) WRITE ( numonp, nampisfer )

      IF(lwp) THEN                         ! control print
         WRITE(numout,*) ' '
         WRITE(numout,*) ' Namelist parameters for Iron chemistry, nampisfer'
         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
         WRITE(numout,*) '    enable complex iron chemistry scheme      ln_fechem =', ln_fechem
         WRITE(numout,*) '    variable concentration of ligand          ln_ligvar =', ln_ligvar
         WRITE(numout,*) '    scavenging rate of Iron                   xlam1     =', xlam1
         WRITE(numout,*) '    scavenging rate of Iron by dust           xlamdust  =', xlamdust
         WRITE(numout,*) '    ligand concentration in the ocean         ligand    =', ligand
      ENDIF
      !
      IF( ln_fechem ) THEN
         ! initialization of some constants used by the complexe chemistry scheme
         ! ----------------------------------------------------------------------
         spd = 3600. * 24.
         con = 1.E9
         ! LIGAND KINETICS (values from Witter et al. 2000)
         ! Weak (L2) ligands
         ! Phaeophytin
         kl2 = 12.2E5  * spd / con
         kb2 = 12.3E-6 * spd
         ! Strong (L1) ligands
         ! Saccharides
         ! kl1 = 12.2E5  * spd / con
         ! kb1 = 12.3E-6 * spd
         ! DFOB-
         kl1 = 19.6e5  * spd / con
         kb1 = 1.5e-6  * spd
         ! pcp and remin of Fe3p
         ks  = 0.075
         kpr = 0.05
         ! thermal reduction of Fe3
         kth = 0.0048 * 24.
         !
      ENDIF
      !
   END SUBROUTINE p4z_fechem_init

#else
   !!======================================================================
   !!  Dummy module :                                   No PISCES bio-model
   !!======================================================================
CONTAINS
   SUBROUTINE p4z_fechem                    ! Empty routine
   END SUBROUTINE p4z_fechem
#endif 

   !!======================================================================
END MODULE p4zfechem