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 all commits
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
87 changes: 66 additions & 21 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_SEDSPINUP = .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 @@ -526,6 +541,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 @@ -1341,33 +1363,41 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
! 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)
! 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
! 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.
&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
WITH_DMSPH = $WITH_DMSPH
PI_PH_FILE = $PI_PH_FILE
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
WITH_DMSPH = $WITH_DMSPH
PI_PH_FILE = $PI_PH_FILE
/

! IO-NAMELIST FOR DIAGNOSTIC iHAMOCC OUTPUT
Expand Down Expand Up @@ -1494,11 +1524,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 @@ -1688,6 +1726,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
36 changes: 31 additions & 5 deletions hamocc/accfields.F90
Original file line number Diff line number Diff line change
Expand Up @@ -44,23 +44,25 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
! *REAL* *omask* - land/ocean mask
!
!**********************************************************************
use mo_carbch, only: atm,atmflx,co2fxd,co2fxu,co3,hi,kwco2sol,ndepflx,ocetra,omegaa,omegac,pco2d,satoxy
use mo_carbch, only: atm,atmflx,co2fxd,co2fxu,co3,hi,kwco2sol,ndepflx,ocetra,omegaa,omegac,pco2d,satoxy,sedfluxo
use mo_biomod, only: bsiflx_bot,bsiflx0100,bsiflx0500,bsiflx1000,bsiflx2000,bsiflx4000,calflx_bot,calflx0100,calflx0500,&
& calflx1000,calflx2000,calflx4000,carflx_bot,carflx0100,carflx0500,carflx1000,carflx2000,carflx4000,&
& expoca,expoor,exposi,intdms_bac,intdms_uv,intdmsprod,intdnit,intnfix,intphosy,phosy3d
use mod_dia, only: ddm
use mod_xc, only: mnproc
use mo_bgcmean, only: domassfluxes,jalkali,jano3,jasize,jatmco2,jbsiflx0100,jbsiflx0500,jbsiflx1000,jbsiflx2000, &
& jbsiflx4000,jbsiflx_bot,jcalc,jcalflx0100,jcalflx0500,jcalflx1000,jcalflx2000,jcalflx4000, &
& jcalflx_bot,jcarflx0100,jcarflx0500,jcarflx1000,jcarflx2000,jcarflx4000,jcarflx_bot,jco2flux, &
& jcalflx_bot,jcarflx0100,jcarflx0500,jcarflx1000,jcarflx2000,jcarflx4000,jcarflx_bot, &
& jsediffic,jsediffal,jsediffph,jsediffox,jsediffn2,jsediffno3,jsediffsi,jco2flux, &
& jco2fxd,jco2fxu,jco3,jdic,jdicsat,jdms,jdms_bac,jdms_uv,jdmsflux,jdmsprod,jdoc,jdp,jeps,jexpoca, &
& jexport,jexposi,jgrazer,jintdnit,jintnfix,jintphosy,jiralk,jirdet,jirdin,jirdip,jirdoc,jiriron, &
& jiron,jirsi,jkwco2,jlvlalkali,jlvlano3,jlvlasize,jlvlbigd14c,jlvlbromo,jlvlcalc,jlvlcalc13, &
& jlvlcfc11,jlvlcfc12,jlvlco3,jlvld13c,jlvld14c,jlvldic,jlvldic13,jlvldic14,jlvldicsat,jlvldoc, &
& jlvldoc13,jlvleps,jlvlgrazer,jlvlgrazer13,jlvliron,jlvln2o,jlvlnatalkali,jlvlnatcalc,jlvlnatco3, &
& jlvlnatdic,jlvlnatomegaa,jlvlnatomegac,jlvlnos,jlvlo2sat,jlvlomegaa,jlvlomegac,jlvlopal,jlvloxygen,&
& jlvlph,jlvlphosph,jlvlphosy,jlvlphyto,jlvlphyto13,jlvlpoc,jlvlpoc13,jlvlprefalk,jlvlprefdic, &
& jlvlprefo2,jlvlprefpo4,jlvlsf6,jlvlsilica,jlvlwnos,jlvlwphy,jn2flux,jn2o,jn2oflux,jn2ofx,jndep, &
& jlvlprefo2,jlvlprefpo4,jlvlsf6,jlvlsilica,jlvlwnos,jlvlwphy,jn2flux,jn2o,jn2oflux,jn2ofx, &
& jprorca,jprcaca,jsilpro,jpodiic,jpodial,jpodiph,jpodiox,jpodin2,jpodino3,jpodisi,jndep, &
& jniflux,jnos,jo2flux,jo2sat,jomegaa,jomegac,jopal,joxflux,joxygen,jpco2,jph,jphosph,jphosy,jphyto, &
& jpoc,jprefalk,jprefdic,jprefo2,jprefpo4,jsilica,jsrfalkali,jsrfano3,jsrfdic,jsrfiron,jsrfoxygen, &
& jsrfphosph,jsrfphyto,jsrfsilica,jwnos,jwphy,nbgc,nacc_bgc,bgcwrt,glb_inventory,bgct2d,acclvl, &
Expand Down Expand Up @@ -152,11 +154,25 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
do j=1,kpje
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
! N-deposition and riverine input fluxes
bgct2d(i,j,jndep) = bgct2d(i,j,jndep) + ndepflx(i,j)/2.0
bgct2d(i,j,jirdin) = bgct2d(i,j,jirdin) + rivinflx(i,j,irdin)/2.0
bgct2d(i,j,jirdip) = bgct2d(i,j,jirdip) + rivinflx(i,j,irdip)/2.0
Expand All @@ -165,7 +181,7 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
bgct2d(i,j,jiriron) = bgct2d(i,j,jiriron) + rivinflx(i,j,iriron)/2.0
bgct2d(i,j,jirdoc) = bgct2d(i,j,jirdoc) + rivinflx(i,j,irdoc)/2.0
bgct2d(i,j,jirdet) = bgct2d(i,j,jirdet) + rivinflx(i,j,irdet)/2.0

endif
enddo
enddo
Expand Down Expand Up @@ -262,6 +278,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
36 changes: 2 additions & 34 deletions hamocc/bodensed.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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,sedac,sedict,sedifti,solfu,orgwei
use mo_control_bgc, only: dtbgc,io_stdo_bgc,isac
& porsol,dzs,seddw,sedict,solfu,orgwei
use mo_control_bgc, only: dtbgc,io_stdo_bgc
use mo_param1_bgc, only: ks
use mod_xc, only: mnproc

Expand Down Expand Up @@ -105,39 +105,7 @@ subroutine bodensed(kpie,kpje,kpke,pddpo)
seddw(k) = 0.5 * (dzs(k) + dzs(k+1))
enddo

! ******************************************************************
! the following section is to include the SEDIMENT ACCELERATION
! mechanism to accelerate the sediment:

sedac = 1. / real(isac)

! determine total solid sediment thickness sumsed
! and reduced sediment volume
sumsed = 0.
do k = 1, ks
porwat(k) = porwat(k) * sedac
porsol(k) = porsol(k) * sedac
sumsed = sumsed + seddw(k)
enddo

! determine reduced diffusion rate sedict,
! and scaling for bottom layer ventilation, sedifti

sedict = 1.e-9 * dtbgc
sedifti = sedict / (sumsed**2)
sedict = sedict * sedac

if (mnproc == 1) then
write(io_stdo_bgc,*) 'sediment acceleration factor: ',isac
write(io_stdo_bgc,*) 'sediment area reduction: ',sedac
write(io_stdo_bgc,*) 'new diffusion is: ',sedict
write(io_stdo_bgc,*) 'total sediment thickness: ',sumsed
write(io_stdo_bgc,*) 'sedict / sumsed**2: ',sedifti
write(io_stdo_bgc,*) 'new pore water fraction: ',porwat
endif


! ******************************************************************

! ******************************************************************
! densities etc. for SEDIMENT SHIFTING
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 @@ -69,7 +69,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 @@ -174,7 +175,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 @@ -187,8 +187,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 @@ -208,9 +208,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 @@ -222,8 +222,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
Loading