From f60e840d1275dc23031dd9a1d672a0c4b661e302 Mon Sep 17 00:00:00 2001 From: jmaerz <92309038+jmaerz@users.noreply.github.com> Date: Thu, 15 Sep 2022 16:28:54 +0200 Subject: [PATCH] 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 | 16 ++++++ hamocc/bodensed.F90 | 96 +++++++++++++++++++++++++-------- hamocc/dipowa.F90 | 23 +++----- hamocc/hamocc_init.F90 | 11 ++-- hamocc/inventory_bgc.F90 | 8 +-- hamocc/meson.build | 1 + hamocc/mo_control_bgc.F90 | 3 +- hamocc/mo_read_sedpor.F90 | 108 ++++++++++++++++++++++++++++++++++++++ hamocc/mo_sedmnt.F90 | 77 +++++++++++++++++++++++++-- hamocc/powach.F90 | 88 +++++++++++++------------------ hamocc/powadi.F90 | 12 ++--- hamocc/sedshi.F90 | 12 ++--- 12 files changed, 342 insertions(+), 113 deletions(-) create mode 100644 hamocc/mo_read_sedpor.F90 diff --git a/cime_config/buildnml b/cime_config/buildnml index 6f3c1e30..bd3d33a3 100755 --- a/cime_config/buildnml +++ b/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 diff --git a/hamocc/bodensed.F90 b/hamocc/bodensed.F90 index 99bc782a..74cb9335 100644 --- a/hamocc/bodensed.F90 +++ b/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,8 +44,8 @@ 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 @@ -53,6 +53,8 @@ subroutine bodensed(kpie,kpje,kpke,pddpo) 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 diff --git a/hamocc/dipowa.F90 b/hamocc/dipowa.F90 index be53c421..e6fb22a2 100644 --- a/hamocc/dipowa.F90 +++ b/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,7 +77,6 @@ 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 @@ -85,14 +84,6 @@ subroutine dipowa(kpie,kpje,kpke,omask,lspin) !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,9 +115,9 @@ 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 @@ -134,7 +125,7 @@ subroutine dipowa(kpie,kpje,kpke,omask,lspin) 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 diff --git a/hamocc/hamocc_init.F90 b/hamocc/hamocc_init.F90 index 0f753f30..f7f2dcf3 100644 --- a/hamocc/hamocc_init.F90 +++ b/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. ! diff --git a/hamocc/inventory_bgc.F90 b/hamocc/inventory_bgc.F90 index 8ca17e61..c9f104f5 100644 --- a/hamocc/inventory_bgc.F90 +++ b/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 diff --git a/hamocc/meson.build b/hamocc/meson.build index ed3461de..acc6319d 100644 --- a/hamocc/meson.build +++ b/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', diff --git a/hamocc/mo_control_bgc.F90 b/hamocc/mo_control_bgc.F90 index b59c19ee..c7058aa5 100644 --- a/hamocc/mo_control_bgc.F90 +++ b/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 diff --git a/hamocc/mo_read_sedpor.F90 b/hamocc/mo_read_sedpor.F90 new file mode 100644 index 00000000..6ea984c6 --- /dev/null +++ b/hamocc/mo_read_sedpor.F90 @@ -0,0 +1,108 @@ +! Copyright (C) 2020 S. Gao, I. Bethke, J. Tjiputra, J. Schwinger +! +! This file is part of BLOM/iHAMOCC. +! +! BLOM is free software: you can redistribute it and/or modify it under the +! terms of the GNU Lesser General Public License as published by the Free +! Software Foundation, either version 3 of the License, or (at your option) +! any later version. +! +! BLOM is distributed in the hope that it will be useful, but WITHOUT ANY +! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +! FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for +! more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with BLOM. If not, see https://www.gnu.org/licenses/. + +module mo_read_sedpor +!***************************************************************************** +! Purpose +! ------- +! - Routine for reading sediment porosity from netcdf file +! +! Description +! ----------- +! Public routines and variable of this module: +! +! - subroutine ini_read_sedpor +! read sediment porosity file +! +! L_SED_POR must be set to true in nml to activate +! lon-lat variable sediment porosity. +! +! The model attempts to read lon-lat-sediment depth variable porosity +! from the input file 'SEDPORFILE' (incl. full path) +! +! sed_por holds then the porosity that can be applied later +! via mo_apply_sedpor +! +!***************************************************************************** + +implicit none + +private + +public :: read_sedpor,sedporfile + +character(len=512),save :: sedporfile = '' + +contains + +subroutine read_sedpor(kpie,kpje,ks,omask,sed_por) + use mod_xc, only: mnproc,xchalt + use mod_dia, only: iotype + use mo_control_bgc, only: io_stdo_bgc,l_3Dvarsedpor + use mod_nctools, only: ncfopn,ncread,ncfcls + + implicit none + + integer, intent(in) :: kpie,kpje,ks + real, intent(in) :: omask(kpie,kpje) + real, intent(inout) :: sed_por(kpie,kpje,ks) + + !local variables + integer :: i,j,k,errstat,dummymask(2) + real :: sed_por_in(kpie,kpje,ks) + logical :: file_exists = .false. + + ! Return if l_3Dvarsedpor is turned off + if (.not. l_3Dvarsedpor) then + if (mnproc.eq.1) then + write(io_stdo_bgc,*) '' + write(io_stdo_bgc,*) 'read_sedpor: spatially variable sediment porosity is not activated.' + endif + return + endif + + ! Check if sediment porosity file exists. If not, abort. + inquire(file=sedporfile,exist=file_exists) + if (.not. file_exists .and. mnproc.eq.1) then + write(io_stdo_bgc,*) '' + write(io_stdo_bgc,*) 'read_sedpor: Cannot find sediment porosity file... ' + call xchalt('(read_sedpor)') + stop '(read_sedpor)' + endif + + ! read sediment porosity from file + if (mnproc.eq.1) then + write(io_stdo_bgc,*) '' + write(io_stdo_bgc,*) 'read_sedpor: read sediment porosity from ', & + trim(sedporfile) + endif + call ncfopn(trim(sedporfile),'r',' ',1,iotype) + call ncread('sedpor',sed_por_in,dummymask,0,0.) + call ncfcls + + do k=1,ks + do j=1,kpje + do i=1,kpie + if(omask(i,j).gt. 0.5)then + sed_por(i,j,k)=sed_por_in(i,j,k) + endif + enddo + enddo + enddo + +end subroutine read_sedpor +end module mo_read_sedpor diff --git a/hamocc/mo_sedmnt.F90 b/hamocc/mo_sedmnt.F90 index 1eb66ade..a1286f22 100644 --- a/hamocc/mo_sedmnt.F90 +++ b/hamocc/mo_sedmnt.F90 @@ -70,13 +70,16 @@ MODULE mo_sedmnt REAL, save :: dzs(ksp) = 0.0 REAL, save :: seddzi(ksp) = 0.0 REAL, save :: seddw(ks) = 0.0 - REAL, save :: porsol(ks) = 0.0 - REAL, save :: porwah(ks) = 0.0 - REAL, save :: porwat(ks) = 0.0 REAL, DIMENSION (:,:,:,:), ALLOCATABLE :: sedlay REAL, DIMENSION (:,:,:,:), ALLOCATABLE :: powtra REAL, DIMENSION (:,:,:), ALLOCATABLE :: sedhpl + REAL, DIMENSION (:,:,:), ALLOCATABLE :: porsol + REAL, DIMENSION (:,:,:), ALLOCATABLE :: porwah + REAL, DIMENSION (:,:,:), ALLOCATABLE :: porwat + REAL, DIMENSION (:,:), ALLOCATABLE :: solfu + REAL, DIMENSION (:,:,:), ALLOCATABLE :: zcoefsu + REAL, DIMENSION (:,:,:), ALLOCATABLE :: zcoeflo REAL, DIMENSION (:,:), ALLOCATABLE :: silpro REAL, DIMENSION (:,:), ALLOCATABLE :: prorca @@ -91,7 +94,8 @@ MODULE mo_sedmnt REAL :: sedict,rno3,o2ut,ansed REAL :: calcwei, opalwei, orgwei REAL :: calcdens, opaldens, orgdens, claydens - REAL :: calfa, oplfa, orgfa, clafa, solfu + REAL :: calfa, oplfa, orgfa, clafa + REAL :: disso_sil,silsat,disso_poc,sed_denit,disso_caco3 CONTAINS @@ -195,6 +199,71 @@ SUBROUTINE ALLOC_MEM_SEDMNT(kpie,kpje) if(errstat.ne.0) stop 'not enough memory sedhpl' sedhpl(:,:,:) = 0.0 + IF (mnproc.eq.1) THEN + WRITE(io_stdo_bgc,*)'Memory allocation for variable porsol ...' + WRITE(io_stdo_bgc,*)'First dimension : ',kpie + WRITE(io_stdo_bgc,*)'Second dimension : ',kpje + WRITE(io_stdo_bgc,*)'Third dimension : ',ks + ENDIF + + ALLOCATE (porsol(kpie,kpje,ks),stat=errstat) + if(errstat.ne.0) stop 'not enough memory porsol' + porsol(:,:,:) = 0.0 + + IF (mnproc.eq.1) THEN + WRITE(io_stdo_bgc,*)'Memory allocation for variable porwah ...' + WRITE(io_stdo_bgc,*)'First dimension : ',kpie + WRITE(io_stdo_bgc,*)'Second dimension : ',kpje + WRITE(io_stdo_bgc,*)'Third dimension : ',ks + ENDIF + + ALLOCATE (porwah(kpie,kpje,ks),stat=errstat) + if(errstat.ne.0) stop 'not enough memory porwah' + porwah(:,:,:) = 0.0 + + IF (mnproc.eq.1) THEN + WRITE(io_stdo_bgc,*)'Memory allocation for variable porwat ...' + WRITE(io_stdo_bgc,*)'First dimension : ',kpie + WRITE(io_stdo_bgc,*)'Second dimension : ',kpje + WRITE(io_stdo_bgc,*)'Third dimension : ',ks + ENDIF + + ALLOCATE (porwat(kpie,kpje,ks),stat=errstat) + if(errstat.ne.0) stop 'not enough memory porwat' + porwat(:,:,:) = 0.0 + + IF (mnproc.eq.1) THEN + WRITE(io_stdo_bgc,*)'Memory allocation for variable solfu ...' + WRITE(io_stdo_bgc,*)'First dimension : ',kpie + WRITE(io_stdo_bgc,*)'Second dimension : ',kpje + ENDIF + + ALLOCATE (solfu(kpie,kpje),stat=errstat) + if(errstat.ne.0) stop 'not enough memory solfu' + solfu(:,:) = 0.0 + + IF (mnproc.eq.1) THEN + WRITE(io_stdo_bgc,*)'Memory allocation for variable zcoefsu ...' + WRITE(io_stdo_bgc,*)'First dimension : ',kpie + WRITE(io_stdo_bgc,*)'Second dimension : ',kpje + WRITE(io_stdo_bgc,*)'Third dimension : ',ks + ENDIF + + ALLOCATE (zcoefsu(kpie,kpje,0:ks),stat=errstat) + if(errstat.ne.0) stop 'not enough memory zcoefsu' + zcoefsu(:,:,:) = 0.0 + + IF (mnproc.eq.1) THEN + WRITE(io_stdo_bgc,*)'Memory allocation for variable zcoeflo ...' + WRITE(io_stdo_bgc,*)'First dimension : ',kpie + WRITE(io_stdo_bgc,*)'Second dimension : ',kpje + WRITE(io_stdo_bgc,*)'Third dimension : ',ks + ENDIF + + ALLOCATE (zcoeflo(kpie,kpje,0:ks),stat=errstat) + if(errstat.ne.0) stop 'not enough memory zcoeflo' + zcoeflo(:,:,:) = 0.0 + IF (mnproc.eq.1) THEN WRITE(io_stdo_bgc,*)'Memory allocation for variable burial ...' WRITE(io_stdo_bgc,*)'First dimension : ',kpie diff --git a/hamocc/powach.F90 b/hamocc/powach.F90 index 5c7f6fbc..540d4c94 100644 --- a/hamocc/powach.F90 +++ b/hamocc/powach.F90 @@ -61,7 +61,8 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) !****************************************************************************** use mo_carbch, only: co3,keqb,ocetra,sedfluxo use mo_chemcon, only: calcon - use mo_sedmnt, only: porwat,porsol,powtra,produs,prcaca,prorca,rno3,seddw,sedhpl,sedlay,silpro + use mo_sedmnt, only: porwat,porsol,powtra,produs,prcaca,prorca,rno3,seddw,sedhpl,sedlay,silpro,disso_sil,silsat,disso_poc, & + & sed_denit,disso_caco3 use mo_biomod, only: rnit,ro2ut use mo_control_bgc, only: dtbgc use mo_param1_bgc, only: ioxygen,ipowaal,ipowaic,ipowaox,ipowaph,ipowasi,ipown2,ipowno3,isilica,isssc12,issso12,issssil, & @@ -92,7 +93,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) real :: aerob13(kpie,ks),anaerob13(kpie,ks) real :: aerob14(kpie,ks),anaerob14(kpie,ks) #endif - real :: disso, dissot, undsa, silsat, posol + real :: dissot, undsa, posol real :: umfa, denit, saln, rrho, alk, c, sit, pt real :: K1, K2, Kb, Kw, Ks1, Kf, Ksi, K1p, K2p, K3p real :: ah1, ac, cu, cb, cc, satlev @@ -117,7 +118,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) !$OMP PARALLEL DO & !$OMP&PRIVATE(sedb1,sediso,solrat,powcar,aerob,anaerob, & -!$OMP& disso,dissot,undsa,silsat,posol, & +!$OMP& dissot,undsa,posol, & !$OMP& umfa,denit,saln,rrho,alk,c,sit,pt, & !$OMP& K1,K2,Kb,Kw,Ks1,Kf,Ksi,K1p,K2p,K3p, & !$OMP& ah1,ac,cu,cb,cc,satlev,bolven, & @@ -158,17 +159,8 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) ! Calculate silicate-opal cycle and simultaneous silicate diffusion !****************************************************************** -! 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=3.e-8 ! (2011-01-04) EMR - disso=1.e-6 ! test vom 03.03.04 half live sil ca. 20.000 yr - dissot=disso*dtbgc - -! Silicate saturation concentration is 1 mol/m3 - - silsat = 0.001 +! Dissolution rate constant of opal (disso) [1/(kmol Si(OH)4/m3)*1/sec]*dtbgc + dissot=disso_sil ! Evaluate boundary conditions for sediment-water column exchange. ! Current undersaturation of bottom water: sedb(i,0) and @@ -180,8 +172,8 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) sedb1(i,0) = bolay(i,j) * (silsat - ocetra(i,j,kbo(i,j),isilica)) & & * bolven(i) solrat(i,1) = ( sedlay(i,j,1,issssil) & - & + silpro(i,j) / (porsol(1) * seddw(1)) ) & - & * dissot / (1. + dissot * undsa) * porsol(1) / porwat(1) + & + silpro(i,j) / (porsol(i,j,1) * seddw(1)) ) & + & * dissot / (1. + dissot * undsa) * porsol(i,j,1) / porwat(i,j,1) endif enddo @@ -194,9 +186,9 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do i = 1, kpie if(omask(i,j) > 0.5) then undsa = silsat - powtra(i,j,k,ipowasi) - sedb1(i,k) = seddw(k) * porwat(k) * (silsat - powtra(i,j,k,ipowasi)) + sedb1(i,k) = seddw(k) * porwat(i,j,k) * (silsat - powtra(i,j,k,ipowasi)) if ( k > 1 ) solrat(i,k) = sedlay(i,j,k,issssil) & - & * dissot / (1. + dissot * undsa) * porsol(k) / porwat(k) + & * dissot / (1. + dissot * undsa) * porsol(i,j,k) / porwat(i,j,k) endif enddo enddo @@ -218,7 +210,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) ocetra(i,j,kbo(i,j),isilica) = silsat - sediso(i,0) endif sedlay(i,j,1,issssil) = & - & sedlay(i,j,1,issssil) + silpro(i,j) / (porsol(1) * seddw(1)) + & sedlay(i,j,1,issssil) + silpro(i,j) / (porsol(i,j,1) * seddw(1)) endif enddo @@ -230,7 +222,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do k = 1, ks do i = 1, kpie if(omask(i,j) > 0.5) then - umfa = porsol(k)/porwat(k) + umfa = porsol(i,j,k)/porwat(i,j,k) solrat(i,k) = sedlay(i,j,k,issssil) * dissot & & / (1. + dissot * sediso(i,k)) posol = sediso(i,k) * solrat(i,k) @@ -243,10 +235,8 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) ! Calculate oxygen-POC cycle and simultaneous oxygen diffusion !************************************************************* -! Degradation rate constant of POP (disso) [1/(kmol O2/m3)*1/sec] - - disso = 0.01 / 86400. ! disso=3.e-5 was quite high - dissot = disso * dtbgc +! Degradation rate constant of POP (disso) [1/(kmol O2/m3)*1/sec]*dtbgc + dissot = disso_poc ! This scheme is not based on undersaturation, but on O2 itself @@ -259,9 +249,9 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) undsa = powtra(i,j,1,ipowaox) sedb1(i,0) = bolay(i,j) * ocetra(i,j,kbo(i,j),ioxygen) * bolven(i) solrat(i,1) = ( sedlay(i,j,1,issso12) + prorca(i,j) & - & / (porsol(1) * seddw(1)) ) & + & / (porsol(i,j,1) * seddw(1)) ) & & * ro2ut * dissot / (1. + dissot * undsa) & - & * porsol(1) / porwat(1) + & * porsol(i,j,1) / porwat(i,j,1) endif enddo @@ -273,9 +263,9 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do i = 1, kpie if(omask(i,j) > 0.5) then undsa = powtra(i,j,k,ipowaox) - sedb1(i,k) = seddw(k) * porwat(k) * powtra(i,j,k,ipowaox) + sedb1(i,k) = seddw(k) * porwat(i,j,k) * powtra(i,j,k,ipowaox) if (k > 1) solrat(i,k) = sedlay(i,j,k,issso12) * ro2ut * dissot & - & / (1. + dissot*undsa) * porsol(k) / porwat(k) + & / (1. + dissot*undsa) * porsol(i,j,k) / porwat(i,j,k) endif enddo enddo @@ -297,12 +287,12 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) ocetra(i,j,kbo(i,j),ioxygen) = sediso(i,0) endif sedlay(i,j,1,issso12) = & - & sedlay(i,j,1,issso12) + prorca(i,j) / (porsol(1)*seddw(1)) + & sedlay(i,j,1,issso12) + prorca(i,j) / (porsol(i,j,1)*seddw(1)) #ifdef cisonew sedlay(i,j,1,issso13) = & - & sedlay(i,j,1,issso13) + pror13(i,j) / (porsol(1)*seddw(1)) + & sedlay(i,j,1,issso13) + pror13(i,j) / (porsol(i,j,1)*seddw(1)) sedlay(i,j,1,issso14) = & - & sedlay(i,j,1,issso14) + pror14(i,j) / (porsol(1)*seddw(1)) + & sedlay(i,j,1,issso14) + pror14(i,j) / (porsol(i,j,1)*seddw(1)) #endif endif enddo @@ -315,7 +305,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do k = 1, ks do i = 1, kpie if(omask(i,j) > 0.5) then - umfa = porsol(k) / porwat(k) + umfa = porsol(i,j,k) / porwat(i,j,k) solrat(i,k) = sedlay(i,j,k,issso12) * dissot/(1. + dissot*sediso(i,k)) posol = sediso(i,k)*solrat(i,k) aerob(i,k) = posol*umfa !this has P units: kmol P/m3 of pore water @@ -345,18 +335,15 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) ! Calculate nitrate reduction under anaerobic conditions explicitely !******************************************************************* -! Denitrification rate constant of POP (disso) [1/sec] -! Store flux in array anaerob, for later computation of DIC and alkalinity. - -!ik denit = 1.e-6*dtbgc - denit = 0.01/86400. *dtbgc + ! Denitrification rate constant of POP (disso) [1/sec]*dtbgc + denit = sed_denit do k = 1, ks do i = 1, kpie if(omask(i,j) > 0.5) then if(powtra(i,j,k,ipowaox) < 1.e-6) then posol = denit * MIN(0.5*powtra(i,j,k,ipowno3)/114., & & sedlay(i,j,k,issso12)) - umfa = porsol(k)/porwat(k) + umfa = porsol(i,j,k)/porwat(i,j,k) anaerob(i,k) = posol*umfa !this has P units: kmol P/m3 of pore water #ifdef cisonew rato13 = sedlay(i,j,k,issso13) / (sedlay(i,j,k,issso12) + safediv) @@ -389,7 +376,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) if(omask(i,j) > 0.5) then if(powtra(i,j,k,ipowaox) < 3.e-6 .and. powtra(i,j,k,ipowno3) < 3.e-6) then posol = denit * sedlay(i,j,k,issso12) ! remineralization of poc - umfa = porsol(k) / porwat(k) + umfa = porsol(i,j,k) / porwat(i,j,k) !this overwrites anaerob from denitrification. added =anaerob+..., works anaerob(i,k) = anaerob(i,k) + posol*umfa !this has P units: kmol P/m3 of pore water #ifdef cisonew @@ -456,9 +443,8 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) enddo -! Dissolution rate constant of CaCO3 (disso) [1/(kmol CO3--/m3)*1/sec] - disso = 1.e-7 - dissot = disso * dtbgc +! Dissolution rate constant of CaCO3 (disso) [1/(kmol CO3--/m3)*1/sec]*dtbgc + dissot = disso_caco3 ! Evaluate boundary conditions for sediment-water column exchange. ! Current undersaturation of bottom water: sedb(i,0) and @@ -473,8 +459,8 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) undsa = MAX( satlev-powcar(i,1), 0. ) sedb1(i,0) = bolay(i,j) * (satlev-co3(i,j,kbo(i,j))) * bolven(i) solrat(i,1) = (sedlay(i,j,1,isssc12) & - & + prcaca(i,j) / (porsol(1)*seddw(1))) & - & * dissot / (1.+dissot*undsa) * porsol(1) / porwat(1) + & + prcaca(i,j) / (porsol(i,j,1)*seddw(1))) & + & * dissot / (1.+dissot*undsa) * porsol(i,j,1) / porwat(i,j,1) endif enddo @@ -486,9 +472,9 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do i = 1, kpie if(omask(i,j) > 0.5) then undsa = MAX( keqb(11,i,j) / calcon - powcar(i,k), 0. ) - sedb1(i,k) = seddw(k) * porwat(k) * undsa + sedb1(i,k) = seddw(k) * porwat(i,j,k) * undsa if (k > 1) solrat(i,k) = sedlay(i,j,k,isssc12) & - & * dissot/(1.+dissot*undsa) * porsol(k)/porwat(k) + & * dissot/(1.+dissot*undsa) * porsol(i,j,k)/porwat(i,j,k) if (undsa <= 0.) solrat(i,k) = 0. endif enddo @@ -504,12 +490,12 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do i = 1, kpie if(omask(i,j) > 0.5) then sedlay(i,j,1,isssc12) = & - & sedlay(i,j,1,isssc12) + prcaca(i,j) / (porsol(1)*seddw(1)) + & sedlay(i,j,1,isssc12) + prcaca(i,j) / (porsol(i,j,1)*seddw(1)) #ifdef cisonew sedlay(i,j,1,isssc13) = & - & sedlay(i,j,1,isssc13) + prca13(i,j) / (porsol(1)*seddw(1)) + & sedlay(i,j,1,isssc13) + prca13(i,j) / (porsol(i,j,1)*seddw(1)) sedlay(i,j,1,isssc14) = & - & sedlay(i,j,1,isssc14) + prca14(i,j) / (porsol(1)*seddw(1)) + & sedlay(i,j,1,isssc14) + prca14(i,j) / (porsol(i,j,1)*seddw(1)) #endif endif enddo @@ -523,7 +509,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do k = 1, ks do i = 1, kpie if(omask(i,j) > 0.5) then - umfa = porsol(k) / porwat(k) + umfa = porsol(i,j,k) / porwat(i,j,k) solrat(i,k) = sedlay(i,j,k,isssc12) & & * dissot / (1. + dissot * sediso(i,k)) posol = sediso(i,k) * solrat(i,k) @@ -565,7 +551,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do j = 1, kpje do i = 1, kpie sedlay(i,j,1,issster) = sedlay(i,j,1,issster) & - & + produs(i,j) / (porsol(1) * seddw(1)) + & + produs(i,j) / (porsol(i,j,1) * seddw(1)) enddo enddo !$OMP END PARALLEL DO diff --git a/hamocc/powadi.F90 b/hamocc/powadi.F90 index a5828d91..413c3046 100644 --- a/hamocc/powadi.F90 +++ b/hamocc/powadi.F90 @@ -76,21 +76,21 @@ subroutine powadi(j,kpie,kpje,solrat,sedb1,sediso,bolven,omask) !********************************************************************** do k = 1, ks - asu = sedict * seddzi(k) * porwah(k) - alo = 0. - if(k < ks) alo = sedict * seddzi(k+1) * porwah(k+1) do i = 1, kpie + asu = sedict * seddzi(k) * porwah(i,j,k) + alo = 0. + if(k < ks) alo = sedict * seddzi(k+1) * porwah(i,j,k+1) tredsy(i,k,1) = -asu tredsy(i,k,3) = -alo - tredsy(i,k,2) = seddw(k) * porwat(k) - tredsy(i,k,1) & - & - tredsy(i,k,3) + solrat(i,k) * porwat(k) * seddw(k) + tredsy(i,k,2) = seddw(k) * porwat(i,j,k) - tredsy(i,k,1) & + & - tredsy(i,k,3) + solrat(i,k) * porwat(i,j,k) * seddw(k) enddo enddo k = 0 asu = 0. - alo = sedict * seddzi(1) * porwah(1) do i = 1, kpie + alo = sedict * seddzi(1) * porwah(i,j,1) if(omask(i,j) > 0.5) then tredsy(i,k,1) = -asu tredsy(i,k,3) = -alo diff --git a/hamocc/sedshi.F90 b/hamocc/sedshi.F90 index 51eea483..44058447 100644 --- a/hamocc/sedshi.F90 +++ b/hamocc/sedshi.F90 @@ -100,7 +100,7 @@ SUBROUTINE SEDSHI(kpie,kpje,omask) uebers=wsed(i,j)*sedlay(i,j,k,iv) sedlay(i,j,k ,iv)=sedlay(i,j,k ,iv)-uebers sedlay(i,j,k+1,iv)=sedlay(i,j,k+1,iv)+uebers & - & *(seddw(k)*porsol(k))/(seddw(k+1)*porsol(k+1)) + & *(seddw(k)*porsol(i,j,k))/(seddw(k+1)*porsol(i,j,k+1)) endif enddo !end i-loop enddo !end j-loop @@ -140,7 +140,7 @@ SUBROUTINE SEDSHI(kpie,kpje,omask) !ka if(bolay(i,j).gt.0.) then uebers=wsed(i,j)*sedlay(i,j,k,iv) sedlay(i,j,ks ,iv)=sedlay(i,j,ks ,iv)-uebers - burial(i,j,iv)=burial(i,j,iv)+uebers*seddw(k)*porsol(k) + burial(i,j,iv)=burial(i,j,iv)+uebers*seddw(k)*porsol(i,j,k) endif enddo !end i-loop enddo !end j-loop @@ -178,7 +178,7 @@ SUBROUTINE SEDSHI(kpie,kpje,omask) & +calfa*sedlay(i,j,k,isssc12) & & +oplfa*sedlay(i,j,k,issssil) & & +clafa*sedlay(i,j,k,issster) - fulsed(i,j)=fulsed(i,j)+porsol(k)*seddw(k)*sedlo + fulsed(i,j)=fulsed(i,j)+porsol(i,j,k)*seddw(k)*sedlo endif enddo !end i-loop enddo !end j-loop @@ -197,7 +197,7 @@ SUBROUTINE SEDSHI(kpie,kpje,omask) ! deficiency to fully loaded sediment packed in sedlay(i,j,ks) ! this is the volume required from the buried layer - seddef=solfu-fulsed(i,j) + seddef=solfu(i,j)-fulsed(i,j) ! total volume of solid constituents in buried layer spresent=orgfa*rcar*burial(i,j,issso12) & @@ -219,7 +219,7 @@ SUBROUTINE SEDSHI(kpie,kpje,omask) ! fill the last active layer refill=seddef/(buried+1.e-10) - frac = porsol(ks)*seddw(ks) !changed k to ks, ik + frac = porsol(i,j,ks)*seddw(ks) sedlay(i,j,ks,issso12)=sedlay(i,j,ks,issso12) & & +refill*burial(i,j,issso12)/frac @@ -269,7 +269,7 @@ SUBROUTINE SEDSHI(kpie,kpje,omask) if(omask(i,j).gt.0.5) then !ka if(bolay(i,j).gt.0.) then uebers=sedlay(i,j,k,iv)*wsed(i,j) - frac=porsol(k)*seddw(k)/(porsol(k-1)*seddw(k-1)) + frac=porsol(i,j,k)*seddw(k)/(porsol(i,j,k-1)*seddw(k-1)) sedlay(i,j,k,iv)=sedlay(i,j,k,iv)-uebers sedlay(i,j,k-1,iv)=sedlay(i,j,k-1,iv)+uebers*frac #ifdef cisonew