Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature hamocc sediment spinup #162

Merged
merged 5 commits into from
May 18, 2022
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
99 changes: 72 additions & 27 deletions cime_config/buildnml
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,10 @@ set BLOM_N_DEPOSITION = `./xmlquery BLOM_N_DEPOSITION --value`
set BLOM_NDEP_SCENARIO = `./xmlquery BLOM_NDEP_SCENARIO --value`
set HAMOCC_VSLS = `./xmlquery HAMOCC_VSLS --value`
set HAMOCC_CISO = `./xmlquery HAMOCC_CISO --value`
set HAMOCC_SEDSPINUP = `./xmlquery HAMOCC_SEDSPINUP --value`
set HAMOCC_SEDSPINUP_YR_START = `./xmlquery HAMOCC_SEDSPINUP_YR_START --value`
set HAMOCC_SEDSPINUP_YR_END = `./xmlquery HAMOCC_SEDSPINUP_YR_END --value`
set HAMOCC_SEDSPINUP_NCYCLE = `./xmlquery HAMOCC_SEDSPINUP_NCYCLE --value`
set RUN_STARTDATE = `./xmlquery RUN_STARTDATE --value`
set PIO_TYPENAME_OCN = `./xmlquery PIO_TYPENAME_OCN --value`
set PIO_NETCDF_FORMAT_OCN = `./xmlquery PIO_NETCDF_FORMAT_OCN --value`
Expand Down Expand Up @@ -185,6 +189,17 @@ else
set DO_NDEP = .false.
set NDEPFNAME = ""
endif
if ($HAMOCC_SEDSPINUP == TRUE) then
set DO_SEDSPINUP = .true.
set SEDSPIN_YR_S = $HAMOCC_SEDSPINUP_YR_START
set SEDSPIN_YR_E = $HAMOCC_SEDSPINUP_YR_END
set SEDSPIN_NCYC = $HAMOCC_SEDSPINUP_NCYCLE
else
set DO_SETSPINUP = .false.
set SEDSPIN_YR_S = -1
set SEDSPIN_YR_E = -1
set SEDSPIN_NCYC = -1
endif
# VSLS is currently only available for the tnx1v4 grid, since a climatlogy file for
# short wave radiation is needed.
if ($HAMOCC_VSLS == TRUE && $OCN_GRID != tnx1v4) then
Expand Down Expand Up @@ -524,6 +539,13 @@ set LVL_WNOS = '0, 2, 2'
set LVL_EPS = '0, 0, 0'
set LVL_ASIZE = '0, 0, 0'
set LVL_BROMO = '0, 2, 2'
set FLX_SEDIFFIC = '0, 0, 2'
set FLX_SEDIFFAL = '0, 0, 2'
set FLX_SEDIFFPH = '0, 0, 2'
set FLX_SEDIFFOX = '0, 0, 2'
set FLX_SEDIFFN2 = '0, 0, 2'
set FLX_SEDIFFNO3 = '0, 0, 2'
set FLX_SEDIFFSI = '0, 0, 2'
set SDM_POWAIC = '0, 0, 2'
set SDM_POWAAL = '0, 0, 2'
set SDM_POWAPH = '0, 0, 2'
Expand Down Expand Up @@ -1336,32 +1358,40 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
!
! CONTENTS:
!
! ATM_CO2 : Atmospheric CO2 concentration [ppmv]
! FEDEPFILE : File name (incl. full path) for iron (dust) deposition data
! SWACLIMFILE : File name (incl. full path) for swa climatology field (needed
! if bromoform scheme is activated)
! DO_RIVINPT : Logical switch to activate riverine input
! RIVINFILE : File name (incl. full path) for riverine input data
! DO_NDEP : Logical switch to activate N-deposition
! NDEPFILE : File name (incl. full path) for atmopheric N-deposition data
! INIXXX : Initial condition file for iHAMOCC, where XXX=DIC, ALK, PO4,
! OXY, NO3, SIL, D13C, and D14C
! ATM_CO2 : Atmospheric CO2 concentration [ppmv]
! FEDEPFILE : File name (incl. full path) for iron (dust) deposition data
! SWACLIMFILE : File name (incl. full path) for swa climatology field (needed
! if bromoform scheme is activated)
! DO_RIVINPT : Logical switch to activate riverine input
! RIVINFILE : File name (incl. full path) for riverine input data
! DO_NDEP : Logical switch to activate N-deposition
! NDEPFILE : File name (incl. full path) for atmopheric N-deposition data
! DO_SEDSPINUP : Logical switch to activate sediment spin-up
! SEDSPIN_YR_S : Start year for sediment spinup
! SEDSPIN_YR_E : End year for sediment spinup
! SEDSPIN_NCYC : Number of subcyles per time-step for sediment spinup
! INIXXX : Initial condition file for iHAMOCC, where XXX=DIC, ALK, PO4,
! OXY, NO3, SIL, D13C, and D14C
&BGCNML
ATM_CO2 = $CCSM_CO2_PPMV
FEDEPFILE = $FEDEPFILE
SWACLIMFILE= $SWACLIMFILE
DO_RIVINPT = $DO_RIVINPT
RIVINFILE = $RIVINFILE
DO_NDEP = $DO_NDEP
NDEPFILE = $NDEPFILE
INIDIC = $INIDIC
INIALK = $INIALK
INIPO4 = $INIPO4
INIOXY = $INIOXY
ININO3 = $ININO3
INISIL = $INISIL
INID13C = $INID13C
INID14C = $INID14C
ATM_CO2 = $CCSM_CO2_PPMV
FEDEPFILE = $FEDEPFILE
SWACLIMFILE = $SWACLIMFILE
DO_RIVINPT = $DO_RIVINPT
RIVINFILE = $RIVINFILE
DO_NDEP = $DO_NDEP
NDEPFILE = $NDEPFILE
DO_SEDSPINUP = $DO_SEDSPINUP
SEDSPIN_YR_S = $SEDSPIN_YR_S
SEDSPIN_YR_E = $SEDSPIN_YR_E
SEDSPIN_NCYC = $SEDSPIN_NCYC
INIDIC = $INIDIC
INIALK = $INIALK
INIPO4 = $INIPO4
INIOXY = $INIOXY
ININO3 = $ININO3
INISIL = $INISIL
INID13C = $INID13C
INID14C = $INID14C
/

! IO-NAMELIST FOR DIAGNOSTIC iHAMOCC OUTPUT
Expand Down Expand Up @@ -1485,11 +1515,19 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
! NFIX - Vertically integrated nitrogen fixation
! DNIT - Vertically integrated denitrification
!
! Particle fluxes (FLX, where ****=0100,0500,1000,2000,4000, or _BOT)
! Particle fluxes (FLX, e.g CARFLX****, where ****=0100,0500,1000,2000,4000, or _BOT)
! and diffusive fluxes at the sediment - water-column interface (SEDIFF*)
! CARFLX**** - POC flux at **** metres depth [mol C m-2 s-1]
! BSIFLX**** - Biogenic silica flux at **** metres depth [mol Si m-2 s-1]
! CALFLX**** - Calcium carbonate flux at **** metres depth [mol C m-2 s-1]
!
! SEDIFFIC - sediment - water-column diffusive flux of DIC [mol C m-2 s-1]
! SEDIFFAL - sediment - water-column diffusive flux of alkalinity [mol m-2 s-1]
! SEDIFFPH - sediment - water-column diffusive flux of phosphate [mol PO3 m-2 s-1]
! SEDIFFOX - sediment - water-column diffusive flux of oxygen [mol O2 m-2 s-1]
! SEDIFFN2 - sediment - water-column diffusive flux of N2 [mol N2 m-2 s-1]
! SEDIFFNO3 - sediment - water-column diffusive flux of nitrate [mol NO3 m-2 s-1]
! SEDIFFSI - sediment - water-column diffusive flux of silica [mol Si m-2 s-1]
!
! Sediment fields (SDM)
! POWAIC - (powdic) [mol C m-3]
! POWAAL - (powalk) [eq m-3]
Expand Down Expand Up @@ -1679,6 +1717,13 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
LVL_EPS = $LVL_EPS
LVL_ASIZE = $LVL_ASIZE
LVL_BROMO = $LVL_BROMO
FLX_SEDIFFIC = $FLX_SEDIFFIC
FLX_SEDIFFAL = $FLX_SEDIFFAL
FLX_SEDIFFPH = $FLX_SEDIFFPH
FLX_SEDIFFOX = $FLX_SEDIFFOX
FLX_SEDIFFN2 = $FLX_SEDIFFN2
FLX_SEDIFFNO3 = $FLX_SEDIFFNO3
FLX_SEDIFFSI = $FLX_SEDIFFSI
SDM_POWAIC = $SDM_POWAIC
SDM_POWAAL = $SDM_POWAAL
SDM_POWAPH = $SDM_POWAPH
Expand Down
40 changes: 40 additions & 0 deletions cime_config/config_component.xml
Original file line number Diff line number Diff line change
Expand Up @@ -147,7 +147,47 @@
<file>env_build.xml</file>
<desc>Set preprocessor option to activate the carbon isotope code. Requires module ecosys</desc>
</entry>

<entry id="HAMOCC_SEDSPINUP">
<type>logical</type>
<valid_values>TRUE,FALSE</valid_values>
<default_value>FALSE</default_value>
<group>run_component_blom</group>
<file>env_run.xml</file>
<desc>Activate sediment spinup. HAMOCC_SEDSPINUP_YR_START and HAMOCC_SEDSPINUP_YR_END
need to be set to valid values for this option to take effect. Requires module ecosys</desc>
</entry>

<entry id="HAMOCC_SEDSPINUP_YR_START">
<type>integer</type>
<valid_values/>
<default_value>-1</default_value>
<group>run_component_blom</group>
<file>env_run.xml</file>
<desc>Set start year for HAMOCC sediment spin-up if HAMOCC_SEDSPINUP == TRUE.
Requires module ecosys</desc>
</entry>

<entry id="HAMOCC_SEDSPINUP_YR_END">
<type>integer</type>
<valid_values/>
<default_value>-1</default_value>
<group>run_component_blom</group>
<file>env_run.xml</file>
<desc>Set end year for HAMOCC sediment spin-up if HAMOCC_SEDSPINUP == TRUE.
Requires module ecosys</desc>
</entry>

<entry id="HAMOCC_SEDSPINUP_NCYCLE">
<type>integer</type>
<valid_values/>
<default_value>-1</default_value>
<group>run_component_blom</group>
<file>env_run.xml</file>
<desc>Set the number of sub-cycles for HAMOCC sediment spin-up if HAMOCC_SEDSPINUP == TRUE.
Requires module ecosys</desc>
</entry>

<entry id="HAMOCC_VSLS">
<type>logical</type>
<valid_values>TRUE,FALSE</valid_values>
Expand Down
23 changes: 23 additions & 0 deletions hamocc/accfields.F90
Original file line number Diff line number Diff line change
Expand Up @@ -99,10 +99,23 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
do i=1,kpie
if(omask(i,j).gt.0.5) then

! Atmosphere-ocean fluxes
bgct2d(i,j,jco2flux) = bgct2d(i,j,jco2flux) + atmflx(i,j,iatmco2)/2.0
bgct2d(i,j,jo2flux) = bgct2d(i,j,jo2flux) + atmflx(i,j,iatmo2)/2.0
bgct2d(i,j,jn2flux) = bgct2d(i,j,jn2flux) + atmflx(i,j,iatmn2)/2.0
bgct2d(i,j,jn2oflux) = bgct2d(i,j,jn2oflux) + atmflx(i,j,iatmn2o)/2.0
! Particle fluxes between water-column and sediment
bgct2d(i,j,jprorca) = bgct2d(i,j,jprorca) + carflx_bot(i,j)/2.0
bgct2d(i,j,jprcaca) = bgct2d(i,j,jprcaca) + calflx_bot(i,j)/2.0
bgct2d(i,j,jsilpro) = bgct2d(i,j,jsilpro) + bsiflx_bot(i,j)/2.0
! Diffusive fluxes between water-column and sediment
bgct2d(i,j,jpodiic) = bgct2d(i,j,jpodiic) + sedfluxo(i,j,ipowaic)/2.0
bgct2d(i,j,jpodial) = bgct2d(i,j,jpodial) + sedfluxo(i,j,ipowaal)/2.0
bgct2d(i,j,jpodiph) = bgct2d(i,j,jpodiph) + sedfluxo(i,j,ipowaph)/2.0
bgct2d(i,j,jpodiox) = bgct2d(i,j,jpodiox) + sedfluxo(i,j,ipowaox)/2.0
bgct2d(i,j,jpodin2) = bgct2d(i,j,jpodin2) + sedfluxo(i,j,ipown2)/2.0
bgct2d(i,j,jpodino3) = bgct2d(i,j,jpodino3) + sedfluxo(i,j,ipowno3)/2.0
bgct2d(i,j,jpodisi) = bgct2d(i,j,jpodisi) + sedfluxo(i,j,ipowasi)/2.0

endif
enddo
Expand Down Expand Up @@ -200,6 +213,16 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
call accsrf(jcalflx_bot,calflx_bot,omask,0)
ENDIF

! Accumulate diffusive fluxes between water column and sediment
call accsrf(jsediffic,sedfluxo(1,1,ipowaic),omask,0)
call accsrf(jsediffal,sedfluxo(1,1,ipowaal),omask,0)
call accsrf(jsediffph,sedfluxo(1,1,ipowaph),omask,0)
call accsrf(jsediffox,sedfluxo(1,1,ipowaox),omask,0)
call accsrf(jsediffn2,sedfluxo(1,1,ipown2),omask,0)
call accsrf(jsediffno3,sedfluxo(1,1,ipowno3),omask,0)
call accsrf(jsediffsi,sedfluxo(1,1,ipowasi),omask,0)


! Accumulate layer diagnostics
call acclyr(jdp,pddpo,pddpo,0)
call acclyr(jphyto,ocetra(1,1,1,iphy),pddpo,1)
Expand Down
15 changes: 8 additions & 7 deletions hamocc/dipowa.F90
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
! along with BLOM. If not, see https://www.gnu.org/licenses/.


subroutine dipowa(kpie,kpje,kpke,omask)
subroutine dipowa(kpie,kpje,kpke,omask,lspin)
!**********************************************************************
!
!**** *DIPOWA* - 'diffusion of pore water'
Expand Down Expand Up @@ -65,7 +65,8 @@ subroutine dipowa(kpie,kpje,kpke,omask)
implicit none

integer, intent(in) :: kpie, kpje, kpke
real, dimension(kpie,kpje), intent(in) :: omask
real, intent(in) :: omask(kpie,kpje)
logical, intent(in) :: lspin

! Local variables
integer :: i,j,k,l,iv
Expand Down Expand Up @@ -170,7 +171,6 @@ subroutine dipowa(kpie,kpje,kpke,omask)
enddo
enddo

! call maschk(kpie,kpje,kpke,23)
! sediment column
do iv = 1,npowtra
do k = 1,ks-1
Expand All @@ -183,8 +183,8 @@ subroutine dipowa(kpie,kpje,kpke,omask)
enddo
enddo
enddo
! call maschk(kpie,kpje,kpke,24)

if(.not. lspin) THEN
! sediment ocean interface
!
! CAUTION - the following assumes same indecees for ocetra and powtra
Expand All @@ -204,9 +204,9 @@ subroutine dipowa(kpie,kpje,kpke,omask)
& ( sedb1(i,l,iv) - tredsy(i,l,3) * powtra(i,j,l+1,iv) ) &
& / tredsy(i,l,2)

! used in inventory_bgc/maschk (diagnostics)
! diffusive fluxes (positive downward)
sedfluxo(i,j,iv) = sedfluxo(i,j,iv) &
& + ocetra(i,j,kbo(i,j),iv) - aprior
& -(ocetra(i,j,kbo(i,j),iv) - aprior)* bolay(i,j)
Comment on lines 212 to +213
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be without accumulation in sedfluxo?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is correct as it is: In powach sedfluxo is set to zero when entering the sediment. Then it is receiving a value for silicate and oxygen (but not the others) in powach. After that dipowa is called, where sedfluxo is updated in a loop. In order to avoid overwriting the values for oxygen and silicate, it is written as above.

#ifdef natDIC
if (iv==isco212) ocetra(i,j,kbo(i,j),inatsco212) = &
& ocetra(i,j,kbo(i,j),inatsco212) + &
Expand All @@ -218,8 +218,9 @@ subroutine dipowa(kpie,kpje,kpke,omask)
endif
enddo
enddo
! call maschk(kpie,kpje,kpke,25)

endif ! .not. lspin

enddo j_loop

end subroutine dipowa
30 changes: 27 additions & 3 deletions hamocc/hamocc4bcm.F90
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,8 @@ SUBROUTINE HAMOCC4BCM(kpie,kpje,kpke,kbnd,kplyear,kplmon,kplday,kldtday,&
REAL, intent(out) :: pflxbromo(1-kbnd:kpie+kbnd,1-kbnd:kpje+kbnd)

INTEGER :: i,j,k,l
INTEGER :: nspin,it
LOGICAL :: lspin

IF (mnproc.eq.1) THEN
write(io_stdo_bgc,*) 'iHAMOCC',KLDTDAY,LDTRUNBGC,NDTDAYBGC
Expand Down Expand Up @@ -173,7 +175,7 @@ SUBROUTINE HAMOCC4BCM(kpie,kpje,kpke,kbnd,kplyear,kplmon,kplday,kldtday,&
#endif

#ifdef BROMO
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(i)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Separate openMP bugfix?

DO j=1,kpje
DO i=1,kpie
IF (patmbromo(i,j).gt.0.) THEN
Expand Down Expand Up @@ -293,8 +295,30 @@ SUBROUTINE HAMOCC4BCM(kpie,kpje,kpke,kbnd,kplyear,kplmon,kplday,kldtday,&
#ifndef sedbypass
! jump over sediment if sedbypass is defined

CALL POWACH(kpie,kpje,kpke,kbnd,prho,omask,psao)
if(do_sedspinup .and. kplyear>=sedspin_yr_s &
.and. kplyear<=sedspin_yr_e) then
nspin = sedspin_ncyc
if(mnproc == 1) then
write(io_stdo_bgc,*)
write(io_stdo_bgc,*) 'iHAMOCC: sediment spinup activated with ', &
nspin, ' subcycles'
endif
else
nspin = 1
endif

! Loop for sediment spinup. If deactivated then nspin=1 and lspin=.false.
do it=1,nspin

if( it<nspin ) then
lspin=.true.
else
lspin=.false.
endif

call POWACH(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin)

enddo

#ifdef PBGC_CK_TIMESTEP
IF (mnproc.eq.1) THEN
Expand Down Expand Up @@ -353,7 +377,7 @@ SUBROUTINE HAMOCC4BCM(kpie,kpje,kpke,kbnd,kplyear,kplmon,kplday,kldtday,&
!--------------------------------------------------------------------
! Pass bromoform flux. Convert unit from kmol CHBr3/m^2 to kg/m^2/s.
! Negative values to the atmosphere
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(i)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Separate openMP bugfix?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, I should have done this in a separate commit, but now I don't have the time to disentangle this so I'll leave it as it is.

DO j=1,kpje
DO i=1,kpie
if(omask(i,j) .gt. 0.5) pflxbromo(i,j)=-252.7*atmflx(i,j,iatmbromo)/dtbgc
Expand Down
Loading