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
! 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