diff --git a/cfgs/SHARED/field_def_nemo-pisces.xml b/cfgs/SHARED/field_def_nemo-pisces.xml
index 32ed08bc8f7070696025c03ac60bbac623606cad..c0181b557c16c428a70e7644a54b5905dc98c5fb 100644
--- a/cfgs/SHARED/field_def_nemo-pisces.xml
+++ b/cfgs/SHARED/field_def_nemo-pisces.xml
@@ -274,27 +274,35 @@
<field id="INTdtDIP" long_name="Vertically int. of change of phophate " unit="mol/m2/s" />
<field id="INTdtSil" long_name="Vertically int. of change of silicon " unit="mol/m2/s" />
-
- <!-- dbio_T on T grid : variables available with diaar5 -->
<field id="TPP" long_name="Total Primary production of phyto" unit="mol/m3/s" grid_ref="grid_T_3D_inner" />
<field id="TPNEW" long_name="New Primary production of phyto" unit="mol/m3/s" grid_ref="grid_T_3D_inner" />
<field id="TPBFE" long_name="Total biogenic iron production" unit="mol/m3/s" grid_ref="grid_T_3D_inner" />
<field id="INTDIC" long_name="DIC content" unit="kg/m2" />
<field id="O2MIN" long_name="Oxygen minimum concentration" unit="mol/m3" />
<field id="ZO2MIN" long_name="Depth of oxygen minimum concentration" unit="m" />
- <field id="INTNFIX" long_name="Nitrogen fixation rate : vert. integrated" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" > Nfix * e3t </field >
- <field id="INTPPPHYN" long_name="Vertically integrated primary production by nanophy" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" > PPPHYN * e3t </field >
- <field id="INTPPPHYD" long_name="Vertically integrated primary production by diatom" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" > PPPHYD * e3t </field >
- <field id="INTPPPHYP" long_name="Vertically integrated primary production by picophy" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" > PPPHYP * e3t </field >
- <field id="INTPP" long_name="Vertically integrated primary production by phyto" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" > TPP * e3t </field >
- <field id="INTPNEW" long_name="Vertically integrated new primary production" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" > TPNEW * e3t </field >
- <field id="INTPBFE" long_name="Vertically integrated of biogenic iron production" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" > TPBFE * e3t </field >
- <field id="INTPBSI" long_name="Vertically integrated of biogenic Si production" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" > PBSi * e3t </field >
+
+ <field id="Nfix_e3t" long_name="Nfix * e3t" unit="molN/m2/s" grid_ref="grid_T_3D_inner" > Nfix * e3t </field >
+ <field id="PPPHYN_e3t" long_name="PPPHYN * e3t" unit="mol/m2/s" grid_ref="grid_T_3D_inner" > PPPHYN * e3t </field >
+ <field id="PPPHYD_e3t" long_name="PPPHYD * e3t" unit="mol/m2/s" grid_ref="grid_T_3D_inner" > PPPHYD * e3t </field >
+ <field id="PPPHYP_e3t" long_name="PPPHYP * e3t" unit="mol/m2/s" grid_ref="grid_T_3D_inner" > PPPHYP * e3t </field >
+ <field id="PBSi_e3t" long_name="PBSi * e3t" unit="mol/m2/s" grid_ref="grid_T_3D_inner" > PBSi * e3t </field >
+ <field id="TPP_e3t" long_name="TPP * e3t" unit="mol/m2/s" grid_ref="grid_T_3D_inner" > TPP * e3t </field >
+ <field id="TPNEW_e3t" long_name="TPNEW * e3t" unit="mol/m2/s" grid_ref="grid_T_3D_inner" > TPNEW * e3t </field >
+ <field id="TPBFE_e3t" long_name="TPBFE * e3t" unit="mol/m2/s" grid_ref="grid_T_3D_inner" > TPBFE * e3t </field >
+
+ <field id="INTNFIX" long_name="Nitrogen fixation rate : vert. integrated" field_ref="Nfix_e3t" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" />
+ <field id="INTPPPHYN" long_name="Vertically integrated primary production by nanophy" field_ref="PPPHYN_e3t" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" />
+ <field id="INTPPPHYD" long_name="Vertically integrated primary production by diatom" field_ref="PPPHYD_e3t" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" />
+ <field id="INTPPPHYP" long_name="Vertically integrated primary production by picophy" field_ref="PPPHYP_e3t" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" />
+ <field id="INTPP" long_name="Vertically integrated primary production by phyto" field_ref="TPP_e3t" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" />
+ <field id="INTPNEW" long_name="Vertically integrated new primary production" field_ref="TPNEW_e3t" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" />
+ <field id="INTPBFE" long_name="Vertically integrated of biogenic iron production" field_ref="TPBFE_e3t" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" />
+ <field id="INTPBSI" long_name="Vertically integrated of biogenic Si production" field_ref="PBSi_e3t" unit="mol/m2/s" grid_ref="grid_T_vsum" detect_missing_value="true" />
<!-- PISCES light : variables available with key_pisces_reduced -->
- <field id="FNO3PHY" long_name="FNO3PHY" unit="" grid_ref="grid_T_3D_inner" />
- <field id="FNH4PHY" long_name="FNH4PHY" unit="" grid_ref="grid_T_3D_inner" />
- <field id="FNH4NO3" long_name="FNH4NO3" unit="" grid_ref="grid_T_3D_inner" />
+ <field id="FNO3PHY" long_name="FNO3PHY" unit="" grid_ref="grid_T_3D_inner" />
+ <field id="FNH4PHY" long_name="FNH4PHY" unit="" grid_ref="grid_T_3D_inner" />
+ <field id="FNH4NO3" long_name="FNH4NO3" unit="" grid_ref="grid_T_3D_inner" />
<field id="TNO3PHY" long_name="TNO3PHY" unit="" />
<field id="TNH4PHY" long_name="TNH4PHY" unit="" />
<field id="TPHYDOM" long_name="TPHYDOM" unit="" />
diff --git a/mk/bldxag.cfg b/mk/bldxag.cfg
index 56b707aaae4efaa19fc210bd83c4e77d914ee553..19ea847577665a653b1db9b57c60ec2ef623e33b 100644
--- a/mk/bldxag.cfg
+++ b/mk/bldxag.cfg
@@ -14,6 +14,7 @@ search_src 1
src::ioipsl $MAIN_DIR/ext/IOIPSL/src
src::nemo $CONFIG_DIR/$NEW_CONF/WORK
+src::ppr_1d $MAIN_DIR/ext/PPR/src
bld::target nemo.exe
bld::exe_dep
@@ -35,8 +36,10 @@ bld::tool::make %MK
# Pre-process code before analysing dependencies
bld::pp::ioipsl 1
bld::pp::nemo 1
+bld::pp::ppr_1d 1
bld::tool::fppflags::nemo %FPPFLAGS
bld::tool::fppflags::ioipsl %FPPFLAGS
+bld::tool::fppflags::ppr_1d %FPPFLAGS
# Ignore the following dependencies
bld::excl_dep inc::netcdf.inc
diff --git a/mk/bldxagxcdf.cfg b/mk/bldxagxcdf.cfg
index c3e85802ea35d1866e1dbf58228f1e4c1c985d7d..d562b61b6f1b686d562889216a606792f10a9c97 100644
--- a/mk/bldxagxcdf.cfg
+++ b/mk/bldxagxcdf.cfg
@@ -15,6 +15,7 @@ search_src 1
src::nocdf $MAIN_DIR/ext/DUMMY_NETCDF
src::ioipsl $MAIN_DIR/ext/IOIPSL/src
src::nemo $CONFIG_DIR/$NEW_CONF/WORK
+src::ppr_1d $MAIN_DIR/ext/PPR/src
bld::target nemo.exe
bld::exe_dep
@@ -37,9 +38,11 @@ bld::tool::make %MK
bld::pp::nocdf 1
bld::pp::ioipsl 1
bld::pp::nemo 1
+bld::pp::ppr_1d 1
bld::tool::fppflags::nemo %FPPFLAGS
bld::tool::fppflags::nocdf %FPPFLAGS
bld::tool::fppflags::ioipsl %FPPFLAGS
+bld::tool::fppflags::ppr_1d %FPPFLAGS
# Ignore the following dependencies
bld::excl_dep inc::VT.inc
diff --git a/src/TOP/PISCES/P4Z/p4zsed.F90 b/src/TOP/PISCES/P4Z/p4zsed.F90
index d978d72acf03f1181752c4959300329cef04f4dc..45d5dafba939d609e8ea8258c62e5aba83f4368a 100644
--- a/src/TOP/PISCES/P4Z/p4zsed.F90
+++ b/src/TOP/PISCES/P4Z/p4zsed.F90
@@ -15,6 +15,7 @@ MODULE p4zsed
USE sms_pisces ! PISCES Source Minus Sink variables
USE p4zlim ! Co-limitations of differents nutrients
USE p4zint ! interpolation and computation of various fields
+ USE p4zsink ! Sinking fluxes
USE sed ! Sediment module
USE iom ! I/O manager
USE prtctl ! print control for debugging
@@ -62,19 +63,17 @@ CONTAINS
INTEGER, INTENT(in) :: kt, knt ! ocean time step
INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices
INTEGER :: ji, jj, jk, ikt
- REAL(wp) :: zrivalk, zrivsil, zrivno3
+ REAL(wp) :: zbureff, zrivsil
REAL(wp) :: zlim, zfact, zfactcal
REAL(wp) :: zo2, zno3, zflx, zpdenit, z1pdenit, zolimit
- REAL(wp) :: zsiloss, zcaloss, zws3, zws4, zwsc, zdep
+ REAL(wp) :: zsiloss, zcaloss, zdep
REAL(wp) :: zwstpoc, zwstpon, zwstpop
REAL(wp) :: ztrfer, ztrpo4s, ztrdp, zwdust, zmudia, ztemp
REAL(wp) :: zsoufer, zlight, ztrpo4, ztrdop
REAL(wp) :: xdiano3, xdianh4
!
CHARACTER (len=25) :: charout
- REAL(wp), DIMENSION(A2D(0)) :: zdenit2d, zbureff
- REAL(wp), DIMENSION(A2D(0)) :: zwsbio3, zwsbio4
- REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsedcal, zsedsi, zsedc
+ REAL(wp), DIMENSION(A2D(0)) :: zdenit2d, zrivno3, zrivalk
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
!!---------------------------------------------------------------------
!
@@ -86,17 +85,10 @@ CONTAINS
l_dia_sed = .NOT.lk_sed .AND. ( iom_use( "SedC" ) .OR. iom_use( "SedCal" ) .OR. iom_use( "SedSi" ) )
ENDIF
- ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments
- ! --------------------------------------------------------------------
- DO_2D( 0, 0, 0, 0 )
- ikt = mbkt(ji,jj)
- zdep = e3t(ji,jj,ikt,Kmm) / xstep
- zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
- zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) )
- END_2D
!
zdenit2d(:,:) = 0.e0
- zbureff (:,:) = 0.e0
+ zrivno3 (:,:) = 0.e0
+ zrivalk (:,:) = 0.e0
!
IF( .NOT.lk_sed ) THEN
! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used
@@ -105,8 +97,7 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
IF( tmask(ji,jj,1) == 1 ) THEN
ikt = mbkt(ji,jj)
- zflx = ( tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj) &
- & + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) ) * 1E3 * 1E6 / 1E4
+ zflx = sinkpocb(ji,jj) / xstep * 1E3 * 1E6 / 1E4
zflx = LOG10( MAX( 1E-3, zflx ) )
zo2 = LOG10( MAX( 10. , tr(ji,jj,ikt,jpoxy,Kbb) * 1E6 ) )
zno3 = LOG10( MAX( 1. , tr(ji,jj,ikt,jpno3,Kbb) * 1E6 * rno3 ) )
@@ -115,9 +106,9 @@ CONTAINS
& + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2
zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) )
!
- zflx = ( tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj) &
- & + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) ) * 1E6
- zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2
+ zflx = sinkpocb(ji,jj) / xstep * 1E6
+ zbureff = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2
+ zrivno3(ji,jj) = 1. - zbureff
ENDIF
END_2D
!
@@ -126,73 +117,33 @@ CONTAINS
! This loss is scaled at each bottom grid cell for equilibrating the total budget of silica in the ocean.
! Thus, the amount of silica lost in the sediments equal the supply at the surface (dust+rivers)
! ------------------------------------------------------
- IF( .NOT.lk_sed ) zrivsil = 1._wp - sedsilfrac
-
- DO_2D( 0, 0, 0, 0 )
- ikt = mbkt(ji,jj)
- zdep = xstep / e3t(ji,jj,ikt,Kmm)
- zwsc = zwsbio4(ji,jj) * zdep
- zsiloss = tr(ji,jj,ikt,jpgsi,Kbb) * zwsc
- zcaloss = tr(ji,jj,ikt,jpcal,Kbb) * zwsc
- !
- tr(ji,jj,ikt,jpgsi,Krhs) = tr(ji,jj,ikt,jpgsi,Krhs) - zsiloss
- tr(ji,jj,ikt,jpcal,Krhs) = tr(ji,jj,ikt,jpcal,Krhs) - zcaloss
- END_2D
!
IF( .NOT.lk_sed ) THEN
+ zrivsil = 1._wp - sedsilfrac
DO_2D( 0, 0, 0, 0 )
ikt = mbkt(ji,jj)
- zdep = xstep / e3t(ji,jj,ikt,Kmm)
- zwsc = zwsbio4(ji,jj) * zdep
- zsiloss = tr(ji,jj,ikt,jpgsi,Kbb) * zwsc
- zcaloss = tr(ji,jj,ikt,jpcal,Kbb) * zwsc
+ zdep = 1._wp / e3t(ji,jj,ikt,Kmm)
+ zsiloss = sinksilb(ji,jj) * zdep
+ zcaloss = sinkcalb(ji,jj) * zdep
tr(ji,jj,ikt,jpsil,Krhs) = tr(ji,jj,ikt,jpsil,Krhs) + zsiloss * zrivsil
!
zfactcal = MAX(-0.1, MIN( excess(ji,jj,ikt), 0.2 ) )
zfactcal = 0.3 + 0.7 * MIN( 1., (0.1 + zfactcal) / ( 0.5 - zfactcal ) )
- zrivalk = sedcalfrac * zfactcal
- tr(ji,jj,ikt,jptal,Krhs) = tr(ji,jj,ikt,jptal,Krhs) + zcaloss * zrivalk * 2.0
- tr(ji,jj,ikt,jpdic,Krhs) = tr(ji,jj,ikt,jpdic,Krhs) + zcaloss * zrivalk
+ zrivalk(ji,jj) = sedcalfrac * zfactcal
+ tr(ji,jj,ikt,jptal,Krhs) = tr(ji,jj,ikt,jptal,Krhs) + zcaloss * zrivalk(ji,jj) * 2.0
+ tr(ji,jj,ikt,jpdic,Krhs) = tr(ji,jj,ikt,jpdic,Krhs) + zcaloss * zrivalk(ji,jj)
END_2D
ENDIF
!
- !
- DO_2D( 0, 0, 0, 0 )
- ikt = mbkt(ji,jj)
- zdep = xstep / e3t(ji,jj,ikt,Kmm)
- zws4 = zwsbio4(ji,jj) * zdep
- zws3 = zwsbio3(ji,jj) * zdep
- tr(ji,jj,ikt,jpgoc,Krhs) = tr(ji,jj,ikt,jpgoc,Krhs) - tr(ji,jj,ikt,jpgoc,Kbb) * zws4
- tr(ji,jj,ikt,jppoc,Krhs) = tr(ji,jj,ikt,jppoc,Krhs) - tr(ji,jj,ikt,jppoc,Kbb) * zws3
- tr(ji,jj,ikt,jpbfe,Krhs) = tr(ji,jj,ikt,jpbfe,Krhs) - tr(ji,jj,ikt,jpbfe,Kbb) * zws4
- tr(ji,jj,ikt,jpsfe,Krhs) = tr(ji,jj,ikt,jpsfe,Krhs) - tr(ji,jj,ikt,jpsfe,Kbb) * zws3
- END_2D
- !
- IF( ln_p5z ) THEN
- DO_2D( 0, 0, 0, 0 )
- ikt = mbkt(ji,jj)
- zdep = xstep / e3t(ji,jj,ikt,Kmm)
- zws4 = zwsbio4(ji,jj) * zdep
- zws3 = zwsbio3(ji,jj) * zdep
- tr(ji,jj,ikt,jpgon,Krhs) = tr(ji,jj,ikt,jpgon,Krhs) - tr(ji,jj,ikt,jpgon,Kbb) * zws4
- tr(ji,jj,ikt,jppon,Krhs) = tr(ji,jj,ikt,jppon,Krhs) - tr(ji,jj,ikt,jppon,Kbb) * zws3
- tr(ji,jj,ikt,jpgop,Krhs) = tr(ji,jj,ikt,jpgop,Krhs) - tr(ji,jj,ikt,jpgop,Kbb) * zws4
- tr(ji,jj,ikt,jppop,Krhs) = tr(ji,jj,ikt,jppop,Krhs) - tr(ji,jj,ikt,jppop,Kbb) * zws3
- END_2D
- ENDIF
! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
! denitrification in the sediments. Not very clever, but simpliest option.
IF( .NOT.lk_sed ) THEN
DO_2D( 0, 0, 0, 0 )
ikt = mbkt(ji,jj)
- zdep = xstep / e3t(ji,jj,ikt,Kmm)
- zws4 = zwsbio4(ji,jj) * zdep
- zws3 = zwsbio3(ji,jj) * zdep
- zrivno3 = 1. - zbureff(ji,jj)
- zwstpoc = tr(ji,jj,ikt,jpgoc,Kbb) * zws4 + tr(ji,jj,ikt,jppoc,Kbb) * zws3
- zpdenit = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
- z1pdenit = zwstpoc * zrivno3 - zpdenit
+ zwstpoc = sinkpocb(ji,jj) / e3t(ji,jj,ikt,Kmm)
+ zpdenit = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3(ji,jj) )
+ z1pdenit = zwstpoc * zrivno3(ji,jj) - zpdenit
zolimit = MIN( ( tr(ji,jj,ikt,jpoxy,Kbb) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
tr(ji,jj,ikt,jpdoc,Krhs) = tr(ji,jj,ikt,jpdoc,Krhs) + z1pdenit - zolimit
tr(ji,jj,ikt,jppo4,Krhs) = tr(ji,jj,ikt,jppo4,Krhs) + zpdenit + zolimit
@@ -206,16 +157,13 @@ CONTAINS
IF( ln_p5z ) THEN
DO_2D( 0, 0, 0, 0 )
ikt = mbkt(ji,jj)
- zdep = xstep / e3t(ji,jj,ikt,Kmm)
- zws4 = zwsbio4(ji,jj) * zdep
- zws3 = zwsbio3(ji,jj) * zdep
- zrivno3 = 1. - zbureff(ji,jj)
- zwstpoc = tr(ji,jj,ikt,jpgoc,Kbb) * zws4 + tr(ji,jj,ikt,jppoc,Kbb) * zws3
- zpdenit = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
- z1pdenit = zwstpoc * zrivno3 - zpdenit
+ zdep = 1._wp / e3t(ji,jj,ikt,Kmm)
+ zwstpoc = sinkpocb(ji,jj) * zdep
+ zwstpop = sinkpopb(ji,jj) * zdep
+ zwstpon = sinkponb(ji,jj) * zdep
+ zpdenit = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3(ji,jj) )
+ z1pdenit = zwstpoc * zrivno3(ji,jj) - zpdenit
zolimit = MIN( ( tr(ji,jj,ikt,jpoxy,Kbb) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
- zwstpop = tr(ji,jj,ikt,jpgop,Kbb) * zws4 + tr(ji,jj,ikt,jppop,Kbb) * zws3
- zwstpon = tr(ji,jj,ikt,jpgon,Kbb) * zws4 + tr(ji,jj,ikt,jppon,Kbb) * zws3
tr(ji,jj,ikt,jpdon,Krhs) = tr(ji,jj,ikt,jpdon,Krhs) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn)
tr(ji,jj,ikt,jpdop,Krhs) = tr(ji,jj,ikt,jpdop,Krhs) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn)
END_2D
@@ -319,31 +267,17 @@ CONTAINS
ENDIF
!
IF( l_dia_sed ) THEN
- ALLOCATE( zsedcal(A2D(0)), zsedsi(A2D(0) ), zsedc(A2D(0) ) )
- zfact = 1.e+3 * rfact2r ! conversion from molC/l/kt to molC/m3/s
- DO_2D( 0, 0, 0, 0 )
- ikt = mbkt(ji,jj)
- zsedsi(ji,jj) = (1.0 - zrivsil) * tr(ji,jj,ikt,jpgsi,Kbb) * zwsbio4(ji,jj) * xstep
- !
- zfactcal = MAX(-0.1, MIN( excess(ji,jj,ikt), 0.2 ) )
- zfactcal = 0.3 + 0.7 * MIN( 1., (0.1 + zfactcal) / ( 0.5 - zfactcal ) )
- zrivalk = sedcalfrac * zfactcal
- zsedcal(ji,jj) = (1.0 - zrivalk) * tr(ji,jj,ikt,jpcal,Kbb) * zwsbio4(ji,jj) * xstep
- !
- zrivno3 = 1. - zbureff(ji,jj)
- zsedc(ji,jj) = (1. - zrivno3) * &
- & ( tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj) + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) ) * xstep
- END_2D
- CALL iom_put( "SedCal", zsedcal(:,:) * zfact )
- CALL iom_put( "SedSi" , zsedsi (:,:) * zfact )
- CALL iom_put( "SedC" , zsedc (:,:) * zfact )
- DEALLOCATE( zsedcal, zsedsi, zsedc )
- ENDIF
- IF( l_dia_sdenit ) THEN
+ zfact = 1.e+3 * rfact2r ! conversion from molC/l/kt to molC/m3/s
+ CALL iom_put( "SedCal", ( 1.0 - zrivalk(:,:) ) * sinkcalb(:,:) * zfact )
+ CALL iom_put( "SedSi" , ( 1.0 - zrivsil ) * sinksilb(:,:) * zfact )
+ CALL iom_put( "SedC" , ( 1.0 - zrivno3(:,:) ) * sinkpocb(:,:) * zfact )
+ ENDIF
+ !
+ IF( l_dia_sdenit ) THEN
zfact = rno3 * 1.e+3 * rfact2r ! conversion from molC/l/kt to molN/m3/s
CALL iom_put( "Sdenit", sdenit(:,:) * zfact )
- ENDIF
- !
+ ENDIF
+ !
ENDIF
!
IF(sn_cfctl%l_prttrc) THEN ! print mean trneds (USEd for debugging)
diff --git a/src/TOP/PISCES/P4Z/p4zsink.F90 b/src/TOP/PISCES/P4Z/p4zsink.F90
index 7b6346706bdb329ac6f8a1dd2d09de063f9a729f..e794e8a7d9dea887fd70324a228d821a42e44085 100644
--- a/src/TOP/PISCES/P4Z/p4zsink.F90
+++ b/src/TOP/PISCES/P4Z/p4zsink.F90
@@ -31,16 +31,15 @@ MODULE p4zsink
PUBLIC p4z_sink_init ! called in trcini_pisces.F90
PUBLIC p4z_sink_alloc ! called in trcini_pisces.F90
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinking, sinking2 !: POC sinking fluxes
- ! ! (different meanings depending on the parameterization)
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkingn, sinking2n !: PON sinking fluxes
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkingp, sinking2p !: POP sinking fluxes
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkcal, sinksil !: CaCO3 and BSi sinking fluxes
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer !: Small BFe sinking fluxes
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer2 !: Big iron sinking fluxes
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sinkpocb !: POC sinking fluxes at bottom
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sinkcalb !: CaCO3 sinking fluxes at bottom
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sinksilb !: BSi sinking fluxes at bottom
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sinkponb !: POC sinking fluxes at bottom
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sinkpopb !: POC sinking fluxes at bottom
INTEGER :: ik100
- LOGICAL :: l_dia_sink2d, l_dia_sink3d, l_dia_tcexp
+ REAL(wp) :: xfact
+ LOGICAL :: l_dia_sink, l_diag
!! * Substitutions
# include "do_loop_substitute.h90"
@@ -69,19 +68,21 @@ CONTAINS
INTEGER :: ji, jj, jk
CHARACTER (len=25) :: charout
REAL(wp) :: zmax, zfact
+ REAL(wp), DIMENSION(A2D(0),jpk) :: zsinking
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
- REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zw2d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_sink')
IF( kt == nittrc000 ) THEN
- l_dia_sink2d = iom_use( "EPC100" ) .OR. iom_use( "EPFE100" ) &
- & .OR. iom_use( "EPCAL100" ) .OR. iom_use( "EPSI100" )
- l_dia_sink3d = iom_use( "EXPC" ) .OR. iom_use( "EXPFE" ) &
- & .OR. iom_use( "EXPCAL" ) .OR. iom_use( "EXPSI" )
- l_dia_tcexp = iom_use( "tcexp" )
+ l_dia_sink = iom_use( "EPC100" ) .OR. iom_use( "EPFE100" ) .OR. iom_use( "EPCAL100" ) .OR. iom_use( "EPSI100" ) &
+ & .OR. iom_use( "EXPC" ) .OR. iom_use( "EXPFE" ) .OR. iom_use( "EXPCAL" ) .OR. iom_use( "EXPSI" ) &
+ & .OR. iom_use( "tcexp" )
+ !
+ xfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s
ENDIF
+ l_diag = l_dia_sink .OR. ( ln_check_mass .AND. kt == nitend )
+ l_diag = l_diag .AND. knt == nrdttrc
! Initialization of some global variables
! ---------------------------------------
@@ -106,91 +107,88 @@ CONTAINS
! Sinking speed of the small particles is always constant
wsbio3(:,:,:) = wsbio
- ! Initialize to zero all the sinking arrays
- ! -----------------------------------------
- sinking (:,:,:) = 0.e0
- sinking2(:,:,:) = 0.e0
- sinkcal (:,:,:) = 0.e0
- sinkfer (:,:,:) = 0.e0
- sinksil (:,:,:) = 0.e0
- sinkfer2(:,:,:) = 0.e0
-
! Compute the sedimentation term using trc_sink for all the sinking particles
! ---------------------------------------------------------------------------
- CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinking , jppoc, rfact2 )
- CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkfer , jpsfe, rfact2 )
- CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2, jpgoc, rfact2 )
- CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinkfer2, jpbfe, rfact2 )
- CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinksil , jpgsi, rfact2 )
- CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinkcal , jpcal, rfact2 )
-
+ zsinking(:,:,:) = 0.e0
+ !
+ CALL trc_sink( kt, Kbb, Kmm, wsbio3, zsinking , jppoc, rfact2 )
+ DO_2D( 0, 0, 0, 0 )
+ sinkpocb(ji,jj) = zsinking(ji,jj,mbkt(ji,jj)+1)
+ END_2D
+ IF( l_diag ) THEN
+ ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),:) = 0._wp
+ zw3d(A2D(0),1:jpkm1) = zsinking(A2D(0),1:jpkm1) * xfact * tmask(A2D(0),1:jpkm1)
+ ENDIF
+ !
+ CALL trc_sink( kt, Kbb, Kmm, wsbio4, zsinking, jpgoc, rfact2 )
+ DO_2D( 0, 0, 0, 0 )
+ sinkpocb(ji,jj) = sinkpocb(ji,jj) + zsinking(ji,jj,mbkt(ji,jj)+1)
+ END_2D
+ IF( l_diag ) THEN
+ zw3d(A2D(0),1:jpkm1) = zw3d(A2D(0),1:jpkm1) + zsinking(A2D(0),1:jpkm1) * xfact * tmask(A2D(0),1:jpkm1)
+ t_oce_co2_exp = glob_sum( 'p4zsink', zw3d(A2D(0),ik100) * e1e2t(A2D(0)) * tmask(A2D(0),1) )
+ CALL iom_put( "EPC100", zw3d(:,:,ik100) ) ! Export of carbon at 100m
+ CALL iom_put( "EXPC" , zw3d ) ! Export of carbon in the water column
+ CALL iom_put( "tcexp", t_oce_co2_exp ) ! Total cabon exort
+ ENDIF
+ !
+ CALL trc_sink( kt, Kbb, Kmm, wsbio4, zsinking, jpgsi, rfact2 )
+ DO_2D( 0, 0, 0, 0 )
+ sinksilb(ji,jj) = zsinking(ji,jj,mbkt(ji,jj)+1)
+ END_2D
+ IF( l_diag ) THEN
+ zw3d(A2D(0),1:jpkm1) = zsinking(A2D(0),1:jpkm1) * xfact * tmask(A2D(0),1:jpkm1)
+ CALL iom_put( "EPSI100", zw3d(:,:,ik100) ) ! Export of Silicate at 100m
+ CALL iom_put( "EXPSI" , zw3d ) ! Export of Silicate in the water column
+ ENDIF
+ !
+ CALL trc_sink( kt, Kbb, Kmm, wsbio4, zsinking, jpcal, rfact2 )
+ DO_2D( 0, 0, 0, 0 )
+ sinkcalb(ji,jj) = zsinking(ji,jj,mbkt(ji,jj)+1)
+ END_2D
+ IF( l_diag ) THEN
+ zw3d(A2D(0),1:jpkm1) = zsinking(A2D(0),1:jpkm1) * xfact * tmask(A2D(0),1:jpkm1)
+ CALL iom_put( "EPCAL100", zw3d(:,:,ik100) ) ! Export of calcite at 100m
+ CALL iom_put( "EXPCAL" , zw3d ) ! Export of calcite in the water column
+ ENDIF
+ !
+ CALL trc_sink( kt, Kbb, Kmm, wsbio3, zsinking, jpsfe, rfact2 )
+ IF( l_diag ) THEN
+ zw3d(A2D(0),1:jpkm1) = zsinking(A2D(0),1:jpkm1) * xfact * tmask(A2D(0),1:jpkm1)
+ ENDIF
+ CALL trc_sink( kt, Kbb, Kmm, wsbio4, zsinking, jpbfe, rfact2 )
+ IF( l_diag ) THEN
+ zw3d(A2D(0),1:jpkm1) = zw3d(A2D(0),1:jpkm1) + zsinking(A2D(0),1:jpkm1) * xfact * tmask(A2D(0),1:jpkm1)
+ CALL iom_put( "EPFE100", zw3d(:,:,ik100) ) ! Export of iron at 100m
+ CALL iom_put( "EXPFE" , zw3d ) ! Export of iron in the water column
+ ENDIF
+ !
! PISCES-QUOTA part
IF( ln_p5z ) THEN
- sinkingn (:,:,:) = 0.e0
- sinking2n(:,:,:) = 0.e0
- sinkingp (:,:,:) = 0.e0
- sinking2p(:,:,:) = 0.e0
-
- ! Compute the sedimentation term using trc_sink for all the sinking particles
- ! ---------------------------------------------------------------------------
- CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkingn , jppon, rfact2 )
- CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkingp , jppop, rfact2 )
- CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2n, jpgon, rfact2 )
- CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2p, jpgop, rfact2 )
- ENDIF
-
- ! Total carbon export per year
- IF( l_dia_tcexp .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc ) ) THEN
- ALLOCATE( zw2d(A2D(0)) )
- zw2d(A2D(0)) = ( sinking(A2D(0),ik100) + sinking2(A2D(0),ik100) ) * e1e2t(A2D(0)) * tmask(A2D(0),1)
- t_oce_co2_exp = glob_sum( 'p4zsink', zw2d(:,:) )
- DEALLOCATE( zw2d )
- ENDIF
- !
- IF( lk_iomput .AND. knt == nrdttrc ) THEN
- !
- IF( l_dia_sink2d ) THEN
- zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s
- ALLOCATE( zw2d(A2D(0)) )
- ! Export of carbon at 100m
- zw2d(A2D(0)) = ( sinking(A2D(0),ik100) + sinking2(A2D(0),ik100) ) * zfact * tmask(A2D(0),1)
- CALL iom_put( "EPC100", zw2d )
- ! Export of iron at 100m
- zw2d(A2D(0)) = ( sinkfer(A2D(0),ik100) + sinkfer2(A2D(0),ik100) ) * zfact * tmask(A2D(0),1)
- CALL iom_put( "EPFE100", zw2d )
- ! Export of calcite at 100m
- zw2d(A2D(0)) = sinkcal(A2D(0),ik100) * zfact * tmask(A2D(0),1)
- CALL iom_put( "EPCAL100", zw2d )
- ! Export of bigenic silica at 100m
- zw2d(A2D(0)) = sinksil(A2D(0),ik100) * zfact * tmask(A2D(0),1)
- CALL iom_put( "EPSI100", zw2d )
- !
- DEALLOCATE( zw2d )
- ENDIF
- !
- IF( l_dia_sink3d ) THEN
- zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s
- ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
- ! Export of carbon in the water column
- zw3d(A2D(0),1:jpkm1) = ( sinking(A2D(0),1:jpkm1) + sinking2(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1)
- CALL iom_put( "EXPC", zw3d )
- ! Export of iron
- zw3d(A2D(0),1:jpkm1) = ( sinkfer(A2D(0),1:jpkm1) + sinkfer2(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1)
- CALL iom_put( "EXPFE", zw3d )
- ! Export of calcite
- zw3d(A2D(0),1:jpkm1) = sinkcal(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
- CALL iom_put( "EXPCAL", zw3d )
- ! Export of bigenic silica
- zw3d(A2D(0),1:jpkm1) = sinksil(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
- CALL iom_put( "EXPSI", zw3d )
- !
- DEALLOCATE( zw3d )
- ENDIF
- !
- IF( l_dia_tcexp ) CALL iom_put( "tcexp", t_oce_co2_exp * 1.e+3 * rfact2r )
- !
+ !
+ CALL trc_sink( kt, Kbb, Kmm, wsbio3, zsinking, jppon, rfact2 )
+ DO_2D( 0, 0, 0, 0 )
+ sinkponb(ji,jj) = zsinking(ji,jj,mbkt(ji,jj)+1)
+ END_2D
+ !
+ CALL trc_sink( kt, Kbb, Kmm, wsbio4, zsinking, jpgon, rfact2 )
+ DO_2D( 0, 0, 0, 0 )
+ sinkponb(ji,jj) = sinkponb(ji,jj) + zsinking(ji,jj,mbkt(ji,jj)+1)
+ END_2D
+ !
+ CALL trc_sink( kt, Kbb, Kmm, wsbio3, zsinking, jppop, rfact2 )
+ DO_2D( 0, 0, 0, 0 )
+ sinkpopb(ji,jj) = zsinking(ji,jj,mbkt(ji,jj)+1)
+ END_2D
+ !
+ CALL trc_sink( kt, Kbb, Kmm, wsbio4, zsinking, jpgop, rfact2 )
+ DO_2D( 0, 0, 0, 0 )
+ sinkpopb(ji,jj) = sinkpopb(ji,jj) + zsinking(ji,jj,mbkt(ji,jj)+1)
+ END_2D
ENDIF
!
+ IF( l_diag ) DEALLOCATE( zw3d )
+ !
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
WRITE(charout, FMT="('sink')")
CALL prt_ctl_info( charout, cdcomp = 'top' )
@@ -223,8 +221,6 @@ CONTAINS
IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ', ik100 + 1
IF (lwp) WRITE(numout,*)
!
- t_oce_co2_exp = 0._wp
- !
END SUBROUTINE p4z_sink_init
INTEGER FUNCTION p4z_sink_alloc()
@@ -236,13 +232,10 @@ CONTAINS
!
ierr(:) = 0
!
- ALLOCATE( sinking(A2D(0),jpk) , sinking2(A2D(0),jpk) , &
- & sinkcal(A2D(0),jpk) , sinksil (A2D(0),jpk) , &
- & sinkfer2(A2D(0),jpk) , &
- & sinkfer(A2D(0),jpk) , STAT=ierr(1) )
- !
- IF( ln_p5z ) ALLOCATE( sinkingn(A2D(0),jpk), sinking2n(A2D(0),jpk) , &
- & sinkingp(A2D(0),jpk), sinking2p(A2D(0),jpk) , STAT=ierr(2) )
+ ALLOCATE( sinkpocb(A2D(0)), sinkcalb(A2D(0)), sinksilb(A2D(0)), STAT=ierr(1) )
+ !
+ IF( ln_p5z ) ALLOCATE( sinkponb(A2D(0)), sinkpopb(A2D(0)), STAT=ierr(2) )
+
!
p4z_sink_alloc = MAXVAL( ierr )
IF( p4z_sink_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_sink_alloc : failed to allocate arrays.' )
diff --git a/src/TOP/TRP/trcsink.F90 b/src/TOP/TRP/trcsink.F90
index 572f501b7b620a40622ec95b8a190469fe0a1dc7..d7735a8e182e7a3839c3cc0b3c1ac898ad250b95 100644
--- a/src/TOP/TRP/trcsink.F90
+++ b/src/TOP/TRP/trcsink.F90
@@ -86,13 +86,8 @@ CONTAINS
ENDIF
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- IF( tmask(ji,jj,jk) == 1.0 ) THEN
- zwsmax = 0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact
- zwsink(ji,jj,jk) = MIN( pwsink(ji,jj,jk), zwsmax * REAL( iiter(ji,jj), wp ) )
- ELSE
- ! provide a default value so there is no use of undefinite value in trc_sink2 for zwsink2 initialization
- zwsink(ji,jj,jk) = 0.
- ENDIF
+ zwsmax = 0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact
+ zwsink(ji,jj,jk+1) = -MIN( pwsink(ji,jj,jk), zwsmax * REAL( iiter(ji,jj), wp ) ) / rday
END_3D
! Initializa to zero all the sinking arrays
@@ -126,19 +121,12 @@ CONTAINS
REAL(wp), INTENT(inout), DIMENSION(A2D(0),jpk) :: psinkflx ! sinking fluxe
!
INTEGER :: ji, jj, jk, jn, jt
- REAL(wp) :: zigma,zew,zign, zflx, zstep
- REAL(wp), DIMENSION(A2D(0),jpk) :: ztraz, zakz, zwsink2, ztrb, psinking
+ REAL(wp) :: zigma,z0w,zign, zflx, zstep, zzwx, zzwy, zalpha
+ REAL(wp), DIMENSION(A2D(0),jpk) :: ztraz, zakz, ztrb, zsinking
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('trc_sink2')
!
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- zwsink2(ji,jj,jk+1) = -pwsink(ji,jj,jk) / rday * tmask(ji,jj,jk+1)
- END_3D
- DO_2D( 0, 0, 0, 0 )
- zwsink2(ji,jj,1) = 0.e0
- END_2D
-
DO_2D( 0, 0, 0, 0 )
! Vertical advective flux
zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2.
@@ -168,27 +156,28 @@ CONTAINS
! vertical advective flux
DO jk = 1, jpkm1
- zigma = zwsink2(ji,jj,jk+1) * zstep / e3w(ji,jj,jk+1,Kmm)
- zew = zwsink2(ji,jj,jk+1)
- psinking(ji,jj,jk+1) = -zew * ( tr(ji,jj,jk,jp_tra,Kbb) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
+ z0w = SIGN( 0.5_wp, pwsink(ji,jj,jk+1) )
+ zalpha = 0.5 + z0w
+ zigma = z0w - 0.5 * pwsink(ji,jj,jk+1) * zstep / e3w(ji,jj,jk+1,Kmm)
+ zzwx = tr(ji,jj,jk+1,jp_tra,Kbb) + zigma * zakz(ji,jj,jk+1)
+ zzwy = tr(ji,jj,jk,jp_tra,Kbb) + zigma * zakz(ji,jj,jk)
+ zsinking(ji,jj,jk+1) = -pwsink(ji,jj,jk+1) * ( zalpha * zzwx + (1.0 - zalpha) * zzwy ) * zstep
END DO
!
! Boundary conditions
- psinking(ji,jj,1 ) = 0.e0
- psinking(ji,jj,jpk) = 0.e0
-
+ zsinking(ji,jj,1 ) = 0.e0
DO jk = 1, jpkm1
- zflx = ( psinking(ji,jj,jk) - psinking(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
- tr(ji,jj,jk,jp_tra,Kbb) = tr(ji,jj,jk,jp_tra,Kbb) + zflx
+ zflx = ( zsinking(ji,jj,jk) - zsinking(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
+ tr(ji,jj,jk,jp_tra,Kbb) = tr(ji,jj,jk,jp_tra,Kbb) + zflx * tmask(ji,jj,jk)
END DO
END DO
DO jk = 1, jpkm1
- zflx = ( psinking(ji,jj,jk) - psinking(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
- ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
+ zflx = ( zsinking(ji,jj,jk) - zsinking(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
+ ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx * tmask(ji,jj,jk)
END DO
tr(ji,jj,:,jp_tra,Kbb) = ztrb(ji,jj,:)
- psinkflx(ji,jj,:) = psinkflx(ji,jj,:) + 2. * psinking(ji,jj,:)
+ psinkflx(ji,jj,:) = psinkflx(ji,jj,:) + 2. * zsinking(ji,jj,:)
END DO
END_2D
!