Lon-lat variable sediment porosity (#189)

Introducing a static 3D sediment porosity field that can be optionally read in with effects on molecular pore water diffusion and shifting.

cime_config/buildnml

@@ -237,6 +237,9 @@ set OALKSCEN = "''"
 set OALKFILE = "''"
 set WITH_DMSPH = .false.
 set PI_PH_FILE = "''"
+set L_3DVARSEDPOR = .false.
+set SEDPORFILE = "''"
+
 
 # set DIAPHY defaults 
 set GLB_FNAMETAG = "'hd','hm','hy'"  
@@ -706,6 +709,7 @@ if      ($OCN_GRID == tnx2v1) then
   set SCFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/sss_clim_core_tnx2v1_20130927.nc'"
   set FEDEPFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/dustdep_mhw2006_tnx2v1_20130506.nc'"
   set SWACLIMFILE = "''"
+  set SEDPORFILE = "''"
   if ($BLOM_RIVER_NUTRIENTS == TRUE) then
     set RIVINFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/river_nutrients_GNEWS2000_tnx2v1_20170915.nc'"
   else
@@ -726,6 +730,7 @@ else if ($OCN_GRID == tnx1v4) then
   set CCFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/chlorophyll_concentration_tnx1v4_20170608.nc'"
   set SCFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/sss_clim_core_tnx1v4_20170604.nc'"
   set FEDEPFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/dustdep_mhw2006_tnx1v4_20171107.nc'"
+  set SEDPORFILE = "''"
   if ($HAMOCC_VSLS == TRUE) then
     set SWACLIMFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/Annual_clim_swa_tnx1v4_20210415.nc'"
   else
@@ -752,6 +757,7 @@ else if ($OCN_GRID == tnx0.25v4) then
   set SCFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/sss_clim_core_tnx0.25v4_20170623.nc'"
   set FEDEPFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/dustdep_mhw2006_tnx0.25v4_20181004.nc'"
   set SWACLIMFILE = "''"
+  set SEDPORFILE = "''"
   if ($BLOM_RIVER_NUTRIENTS == TRUE) then
     set RIVINFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/river_nutrients_GNEWS2000_tnx0.25v4_20170821.nc'"
   else
@@ -773,6 +779,7 @@ else if ($OCN_GRID == tnx0.125v4) then
   set SCFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/sss_clim_core_tnx0.125v4_20200722.nc'"
   set FEDEPFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/dustdep_mhw2006_tnx0.125v4_20200722.nc'"
   set SWACLIMFILE = "''"
+  set SEDPORFILE = "''"
   if ($BLOM_RIVER_NUTRIENTS == TRUE) then
     set RIVINFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/river_nutrients_GNEWS2000_tnx0.125v4_20170821.nc'"
   else
@@ -1445,6 +1452,8 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
 !                OXY, NO3, SIL, D13C, and D14C
 ! WITH_DMSPH  : Logical switch to activate DMS calculation as function of pH
 ! PI_PH_FILE  : File name (incl. full path) for surface PI pH input data.
+! L_3DVARSEDPOR : Logical switch to enable lon-lat-depth variable sediment porosity (as opposed to default: only depth) 
+! SEDPORFILE  : File name (incl. full path) for sediment porosity
 &BGCNML
   ATM_CO2      = $CCSM_CO2_PPMV
   FEDEPFILE    = $FEDEPFILE
@@ -1470,6 +1479,8 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
   INID14C      = $INID14C
   WITH_DMSPH   = $WITH_DMSPH
   PI_PH_FILE   = $PI_PH_FILE
+  L_3DVARSEDPOR = $L_3DVARSEDPOR
+  SEDPORFILE   = $SEDPORFILE
 /
 
 ! IO-NAMELIST FOR DIAGNOSTIC iHAMOCC OUTPUT
@@ -1871,6 +1882,11 @@ EOF
   if ($HAMOCC_VSLS == TRUE) then
 cat >> $CASEBUILD/blom.input_data_list << EOF
 swa_clim_file = `echo $SWACLIMFILE | tr -d '"' | tr -d "'"`
+EOF
+  endif
+  if ($L_3DVARSEDPOR == TRUE) then
+cat >> $CASEBUILD/blom.input_data_list << EOF
+sed_porosity_file = `echo $SEDPORFILE | tr -d '"' | tr -d "'"`
 EOF
   endif
 endif

hamocc/bodensed.F90

@@ -17,7 +17,7 @@
 ! along with BLOM. If not, see https://www.gnu.org/licenses/.
 
 
-subroutine bodensed(kpie,kpje,kpke,pddpo)
+subroutine bodensed(kpie,kpje,kpke,pddpo,omask,sed_por)
 !**********************************************************************
 !
 !**** *BODENSED* - .
@@ -44,15 +44,17 @@ subroutine bodensed(kpie,kpje,kpke,pddpo)
 !**********************************************************************
 
   use mo_sedmnt,      only: calcwei,calfa,clafa,claydens,calcdens,opaldens,opalwei,oplfa,orgdens,orgfa,seddzi,porwat,porwah,       &
-                          & porsol,dzs,seddw,sedict,solfu,orgwei
-  use mo_control_bgc, only: dtbgc,io_stdo_bgc
+                          & porsol,dzs,seddw,sedict,solfu,orgwei,zcoefsu,zcoeflo,disso_sil,silsat,disso_poc,sed_denit,disso_caco3
+  use mo_control_bgc, only: dtbgc,io_stdo_bgc,l_3Dvarsedpor
   use mo_param1_bgc,  only: ks
   use mod_xc,         only: mnproc
 
   implicit none
 
   integer, intent(in) :: kpie,kpje,kpke
   real,    intent(in) :: pddpo(kpie,kpje,kpke)
+  real,    intent(in) :: omask(kpie,kpje)
+  real,    intent(in) :: sed_por(kpie,kpje,ks)
 
   ! Local variables
   integer             :: i,j,k
@@ -79,33 +81,68 @@ subroutine bodensed(kpie,kpje,kpke,pddpo)
      write(io_stdo_bgc,*)  ' '
   endif
 
-  porwat(1) = 0.85
-  porwat(2) = 0.83
-  porwat(3) = 0.8
-  porwat(4) = 0.79
-  porwat(5) = 0.77
-  porwat(6) = 0.75
-  porwat(7) = 0.73
-  porwat(8) = 0.7
-  porwat(9) = 0.68
-  porwat(10) = 0.66
-  porwat(11) = 0.64
-  porwat(12) = 0.62
+  ! this initialization can be done later via reading a porosity map
+  if (l_3Dvarsedpor)then
+   ! lon-lat variable sediment porosity from input file
+   do k=1,ks
+   do j=1,kpje
+   do i=1,kpie
+     if(omask(i,j).gt. 0.5)then
+       porwat(i,j,k) = sed_por(i,j,k)
+     endif
+   enddo
+   enddo
+   enddo
+  else
+   porwat(:,:,1) = 0.85
+   porwat(:,:,2) = 0.83
+   porwat(:,:,3) = 0.8
+   porwat(:,:,4) = 0.79
+   porwat(:,:,5) = 0.77
+   porwat(:,:,6) = 0.75
+   porwat(:,:,7) = 0.73
+   porwat(:,:,8) = 0.7
+   porwat(:,:,9) = 0.68
+   porwat(:,:,10) = 0.66
+   porwat(:,:,11) = 0.64
+   porwat(:,:,12) = 0.62
+  endif
 
   if (mnproc == 1) then
-     write(io_stdo_bgc,*)  'Pore water in sediment: ',porwat
+     write(io_stdo_bgc,*)  'Pore water in sediment initialized'
   endif
 
   seddzi(1) = 500.
   do k = 1, ks
-     porsol(k) = 1. - porwat(k)
-     if(k >= 2) porwah(k) = 0.5 * (porwat(k) + porwat(k-1))
-     if(k == 1) porwah(k) = 0.5 * (1. + porwat(1))
      seddzi(k+1) = 1. / dzs(k+1)
      seddw(k) = 0.5 * (dzs(k) + dzs(k+1))
+     do j = 1, kpje
+     do i = 1, kpie
+        porsol(i,j,k) = 1. - porwat(i,j,k)
+        if(k >= 2) porwah(i,j,k) = 0.5 * (porwat(i,j,k) + porwat(i,j,k-1))
+        if(k == 1) porwah(i,j,k) = 0.5 * (1. + porwat(i,j,1))
+     enddo
+     enddo
   enddo
+
+  sedict = 1.e-9 * dtbgc ! Molecular diffusion coefficient
+  ! Dissolution rate constant of opal (disso) [1/(kmol Si(OH)4/m3)*1/sec]
+  ! THIS NEEDS TO BE CHANGED TO disso=3.e-8! THIS IS ONLY KEPT FOR THE MOMENT
+  ! FOR BACKWARDS COMPATIBILITY
+  !disso_sil = 3.e-8*dtbgc  ! (2011-01-04) EMR
+  !disso_sil = 1.e-6*dtbgc  ! test vom 03.03.04 half live sil ca. 20.000 yr 
+  disso_sil = 1.e-6*dtbgc
+  ! Silicate saturation concentration is 1 mol/m3
+  silsat    = 0.001
+
+  ! Degradation rate constant of POP (disso) [1/(kmol O2/m3)*1/sec]
+  disso_poc = 0.01 / 86400. * dtbgc  !  disso=3.e-5 was quite high
 
-  sedict = 1.e-9 * dtbgc
+  ! Denitrification rate constant of POP (disso) [1/sec]
+  sed_denit =  0.01/86400. * dtbgc 
+
+  ! Dissolution rate constant of CaCO3 (disso) [1/(kmol CO3--/m3)*1/sec]
+  disso_caco3 = 1.e-7 * dtbgc
 
 ! ******************************************************************
 ! densities etc. for SEDIMENT SHIFTING
@@ -131,9 +168,26 @@ subroutine bodensed(kpie,kpje,kpke,pddpo)
 
 ! determine total solid sediment volume
   solfu = 0.
+  do i = 1, kpie
+  do j = 1, kpje
   do k = 1, ks
-     solfu = solfu + seddw(k) * porsol(k)
+     solfu(i,j) = solfu(i,j) + seddw(k) * porsol(i,j,k)
+  enddo
+  enddo
   enddo
 
+! Initialize porosity-dependent diffusion coefficients of sediment
+  zcoefsu(:,:,0) = 0.0
+  do k = 1,ks
+  do j = 1, kpje
+  do i = 1, kpie
+     ! sediment diffusion coefficient * 1/dz * fraction of pore water at half depths
+     zcoefsu(i,j,k  ) = -sedict * seddzi(k) * porwah(i,j,k)
+     zcoeflo(i,j,k-1) = -sedict * seddzi(k) * porwah(i,j,k)    ! why the same ?
+  enddo
+  enddo
+  enddo
+  zcoeflo(:,:,ks) = 0.0                    ! diffusion coefficient for bottom sediment layer
 
+
 end subroutine bodensed

hamocc/dipowa.F90

@@ -56,7 +56,7 @@ subroutine dipowa(kpie,kpje,kpke,omask,lspin)
 !**********************************************************************
 
   use mo_carbch,     only: ocetra, sedfluxo 
-  use mo_sedmnt,     only: powtra,porwat,porwah,sedict,seddw,seddzi 
+  use mo_sedmnt,     only: powtra,porwat,porwah,sedict,seddw,seddzi,zcoefsu,zcoeflo 
   use mo_param1_bgc, only: ks,npowtra,map_por2octra
   use mo_vgrid,      only: kbo,bolay
 #ifdef cisonew
@@ -77,22 +77,13 @@ subroutine dipowa(kpie,kpje,kpke,omask,lspin)
   integer :: iv_oc                     ! index of ocetra in powtra loop
 
   real :: sedb1(kpie,0:ks,npowtra)     ! ????
-  real :: zcoefsu(0:ks),zcoeflo(0:ks)  ! diffusion coefficients (upper/lower)
   real :: tredsy(kpie,0:kpke,3)        ! redsy for 'reduced system'?
   real :: aprior                       ! start value of oceanic tracer in bottom layer
 
 !ik accelerated sediment
 !ik needed for boundary layer ventilation in fast sediment routine
   real :: bolven(kpie)                 ! bottom layer ventilation rate
 
-  zcoefsu(0) = 0.0
-  do k = 1,ks
-     ! sediment diffusion coefficient * 1/dz * fraction of pore water at half depths
-     zcoefsu(k  ) = -sedict * seddzi(k) * porwah(k)
-     zcoeflo(k-1) = -sedict * seddzi(k) * porwah(k)    ! why the same ?
-  enddo
-  zcoeflo(ks) = 0.0                    ! diffusion coefficient for bottom sediment layer
-
 !$OMP PARALLEL DO                            &
 !$OMP&PRIVATE(i,k,iv,l,bolven,tredsy,sedb1,aprior,iv_oc)
   j_loop: do j=1,kpje
@@ -104,8 +95,8 @@ subroutine dipowa(kpie,kpje,kpke,omask,lspin)
 
   k = 0
   do i = 1,kpie
-     tredsy(i,k,1) = zcoefsu(k)
-     tredsy(i,k,3) = zcoeflo(k)
+     tredsy(i,k,1) = zcoefsu(i,j,k)
+     tredsy(i,k,3) = zcoeflo(i,j,k)
      tredsy(i,k,2) = bolven(i)*bolay(i,j) - tredsy(i,k,1) - tredsy(i,k,3)
      !                            dz(kbo) - diff upper    - diff lower
   enddo
@@ -124,17 +115,17 @@ subroutine dipowa(kpie,kpje,kpke,omask,lspin)
 
   do k = 1,ks
      do i = 1,kpie
-        tredsy(i,k,1) = zcoefsu(k)
-        tredsy(i,k,3) = zcoeflo(k)
-        tredsy(i,k,2) = seddw(k)*porwat(k) -tredsy(i,k,1) -tredsy(i,k,3)
+        tredsy(i,k,1) = zcoefsu(i,j,k)
+        tredsy(i,k,3) = zcoeflo(i,j,k)
+        tredsy(i,k,2) = seddw(k)*porwat(i,j,k) -tredsy(i,k,1) -tredsy(i,k,3)
      enddo
   enddo
 
   do iv = 1,npowtra
      do k = 1,ks
         do i = 1,kpie
            ! tracer_concentration(k[1:ks]) * porewater fraction(k) * dz(k)
-           sedb1(i,k,iv) = powtra(i,j,k,iv) * porwat(k) * seddw(k)
+           sedb1(i,k,iv) = powtra(i,j,k,iv) * porwat(i,j,k) * seddw(k)
         enddo
      enddo
   enddo

hamocc/hamocc_init.F90

@@ -46,7 +46,7 @@ subroutine hamocc_init(read_rest,rstfnm_hamocc)
        &                    do_ndep,do_rivinpt,do_oalk,do_sedspinup,            &
        &                    sedspin_yr_s,sedspin_yr_e,sedspin_ncyc,             &
        &                    dtb,dtbgc,io_stdo_bgc,ldtbgc,                       &
-       &                    ldtrunbgc,ndtdaybgc,with_dmsph
+       &                    ldtrunbgc,ndtdaybgc,with_dmsph,l_3Dvarsedpor
   use mo_param1_bgc,  only: ks,nsedtra,npowtra,init_por2octra_mapping
   use mo_carbch,      only: alloc_mem_carbch,ocetra,atm,atm_co2
   use mo_biomod,      only: alloc_mem_biomod
@@ -58,6 +58,7 @@ subroutine hamocc_init(read_rest,rstfnm_hamocc)
   use mo_read_ndep,   only: ini_read_ndep,ndepfile
   use mo_read_oafx,   only: ini_read_oafx,oalkfile,oalkscen
   use mo_read_pi_ph,  only: ini_pi_ph,pi_ph_file
+  use mo_read_sedpor, only: read_sedpor,sedporfile
   use mo_clim_swa,    only: ini_swa_clim,swaclimfile
   use mo_Gdata_read,  only: inidic,inialk,inipo4,inioxy,inino3,                 &
        &                    inisil,inid13c,inid14c
@@ -76,13 +77,14 @@ subroutine hamocc_init(read_rest,rstfnm_hamocc)
 
   integer :: i,j,k,l,nt
   integer :: iounit
+  real    :: sed_por(idm,jdm,ks) = 0.
 
   namelist /bgcnml/ atm_co2,fedepfile,do_rivinpt,rivinfile,do_ndep,ndepfile,    &
        &   do_oalk,oalkscen,oalkfile,do_sedspinup,sedspin_yr_s,                 &
        &   sedspin_yr_e,sedspin_ncyc,                                           &
        &   inidic,inialk,inipo4,inioxy,inino3,inisil,                           &
        &   inid13c,inid14c,swaclimfile,                                         &
-       &   with_dmsph,pi_ph_file
+       &   with_dmsph,pi_ph_file,l_3Dvarsedpor,sedporfile
   !
   ! --- Set io units and some control parameters
   !
@@ -174,8 +176,9 @@ subroutine hamocc_init(read_rest,rstfnm_hamocc)
   call set_vgrid(idm,jdm,kdm,bgc_dp)
   !
   ! --- Initialize sediment layering
-  !
-  CALL BODENSED(idm,jdm,kdm,bgc_dp)
+  ! First raed the porosity, then apply it in bodensed
+  CALL read_sedpor(idm,jdm,ks,omask,sed_por)
+  CALL BODENSED(idm,jdm,kdm,bgc_dp,omask,sed_por)
   !
   ! --- Initialize parameters, sediment and ocean tracer.
   !

hamocc/inventory_bgc.F90

@@ -143,7 +143,7 @@ SUBROUTINE INVENTORY_BGC(kpie,kpje,kpke,dlxp,dlyp,ddpo,omask,iogrp)
      DO j=1,kpje
         DO i=1,kpie
            ztmp1(i,j) = ztmp1(i,j) + omask(i,j)*seddw(k)                       &
-                &       *dlxp(i,j)*dlyp(i,j)*porwat(k)
+                &       *dlxp(i,j)*dlyp(i,j)*porwat(i,j,k)
         ENDDO
      ENDDO
   ENDDO
@@ -155,7 +155,7 @@ SUBROUTINE INVENTORY_BGC(kpie,kpje,kpke,dlxp,dlyp,ddpo,omask,iogrp)
      DO k=1,ks
         DO j=1,kpje
            DO i=1,kpie
-              vol    = seddw(k)*dlxp(i,j)*dlyp(i,j)*porwat(k)
+              vol    = seddw(k)*dlxp(i,j)*dlyp(i,j)*porwat(i,j,k)
               ztmp1(i,j)= ztmp1(i,j) + omask(i,j)*powtra(i,j,k,l)*vol
            ENDDO
         ENDDO
@@ -174,7 +174,7 @@ SUBROUTINE INVENTORY_BGC(kpie,kpje,kpke,dlxp,dlyp,ddpo,omask,iogrp)
      DO k=1,ks
         DO j=1,kpje
            DO i=1,kpie
-              vol = porsol(k)*seddw(k)*dlxp(i,j)*dlyp(i,j)
+              vol = porsol(i,j,k)*seddw(k)*dlxp(i,j)*dlyp(i,j)
               ztmp1(i,j) = ztmp1(i,j) + omask(i,j)*sedlay(i,j,k,l)*vol
            ENDDO
         ENDDO
@@ -187,7 +187,7 @@ SUBROUTINE INVENTORY_BGC(kpie,kpje,kpke,dlxp,dlyp,ddpo,omask,iogrp)
   DO k=1,ks
      DO j=1,kpje
         DO i=1,kpie
-           vol = porsol(k)*seddw(k)*dlxp(i,j)*dlyp(i,j)
+           vol = porsol(i,j,k)*seddw(k)*dlxp(i,j)*dlyp(i,j)
            ztmp1(i,j) = ztmp1(i,j) + omask(i,j)*sedhpl(i,j,k)*vol
         ENDDO
      ENDDO

hamocc/meson.build

@@ -34,6 +34,7 @@ sources += files(
   'mo_read_pi_ph.F90',
   'mo_read_rivin.F90',
   'mo_read_oafx.F90',
+  'mo_read_sedpor.F90',
   'mo_sedmnt.F90',
   'mo_vgrid.F90',
   'ncout_hamocc.F90',

hamocc/mo_control_bgc.F90

@@ -57,8 +57,9 @@ MODULE mo_control_bgc
 
       REAL,    save :: rmasks = 0.0       !  value at wet cells in sediment.
       REAL,    save :: rmasko = 99999.00  !  value at wet cells in ocean.
-
+      
 ! Logical switches set via namelist
+      LOGICAL, save :: l_3Dvarsedpor = .false.  ! apply lon-lat-depth variable sediment porosity via input file
       LOGICAL, save :: do_ndep     =.true.   ! apply n-deposition
       LOGICAL, save :: do_rivinpt  =.true.   ! apply riverine input
       LOGICAL, save :: do_sedspinup=.false.  ! apply sediment spin-up