Skip to content

Added flexibility on soil layer definitions in MPAS #33

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.

Sign up for GitHub

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

Draft
mpaiao wants to merge 2 commits into
base: feature/monan_2.0.0
Choose a base branch
from
Draft

Added flexibility on soil layer definitions in MPAS #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8fcda2b
This fixes a cryptic bug that only occurs when the code is modified t…
mpaiao May 26, 2026
b592e0c
First commit implementing a flexible soil grid.
mpaiao May 26, 2026
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
12 changes: 12 additions & 0 deletions src/core_atmosphere/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -591,6 +591,8 @@
<var name="precipw"/>
<var name="precipci"/>
<var name="ctt"/>
<var name="dzstop"/>
<var name="kzs"/>
<var name="dzs"/>
<var name="sh2o"/>
<var name="smois"/>
Expand Down Expand Up @@ -883,6 +885,8 @@
<var name="seaice"/>
<var name="xice"/>
<var name="xland"/>
<var name="dzstop"/>
<var name="kzs"/>
<var name="dzs"/>
<var name="smcrel"/>
<var name="sh2o"/>
Expand Down Expand Up @@ -1240,6 +1244,8 @@
<var name="q2"/>
<var name="t2m"/>
<var name="precipw"/>
<var name="dzstop"/>
<var name="kzs"/>
<var name="dzs"/>
<var name="zs"/>
<var name="ter"/>
Expand Down Expand Up @@ -3898,6 +3904,12 @@
<var name="greenfrac" type="real" dimensions="nMonths nCells" units="percent"
description="monthly-mean climatological greenness fraction"/>

<var name="dzstop" type="real" dimensions="nCells Time" units="m"
description="soil layer thickness of the top most soil layer "/>

<var name="kzs" type="real" dimensions="nCells Time" units=""
description="exponential factor for thickening deeper soil layers"/>

<var name="dzs" type="real" dimensions="nSoilLevels nCells Time" units="m"
description="soil layer thickness"/>

Expand Down
2 changes: 1 addition & 1 deletion src/core_atmosphere/mpas_atm_core.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -572,7 +572,7 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
call mpas_pool_get_subpool(block % structs, 'atm_input', atm_input)
call mpas_pool_get_subpool(block % structs, 'output_noahmp', output_noahmp)
call physics_tables_init(dminfo, block % configs)
call physics_registry_init(mesh, block % configs, sfc_input)
call physics_registry_init(mesh, block % dimensions, block % configs, sfc_input)
call physics_run_init(block % configs, mesh, state, clock, stream_manager)

!initialization of all physics:
Expand Down
2 changes: 1 addition & 1 deletion src/core_atmosphere/physics/Registry_noahmp.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
<dim name="nSnowLevels" definition="3"
description="The number of snow layers used by the NOAHMP land-surface scheme"/>

<dim name="nzSnowLevels" definition="7"
<dim name="nzSnowLevels" definition="nSoilLevels+nSnowLevels"
description="The number of snow layer depths used by the NOAHMP land-surface scheme"/>

</dims>
Expand Down
41 changes: 27 additions & 14 deletions src/core_atmosphere/physics/mpas_atmphys_control.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -358,11 +358,12 @@ subroutine physics_namelist_check(configs)
end subroutine physics_namelist_check

!=================================================================================================================
subroutine physics_registry_init(mesh,configs,sfc_input)
subroutine physics_registry_init(mesh,dims,configs,sfc_input)
!=================================================================================================================

!input arguments:
type(mpas_pool_type),intent(in):: mesh
type(mpas_pool_type),intent(in):: dims
type(mpas_pool_type),intent(in):: configs

!inout arguments:
Expand All @@ -372,40 +373,52 @@ subroutine physics_registry_init(mesh,configs,sfc_input)
logical,pointer:: config_do_restart
character(len=StrKIND),pointer:: config_lsm_scheme
integer,pointer:: nCells
integer, pointer :: nSoilLevels
integer,dimension(:),pointer:: landmask

real(kind=RKIND),dimension(:,:),pointer:: dzs
real(kind=RKIND),dimension(:,:), pointer :: dzs
real(kind=RKIND),dimension(:) , pointer :: dzstop, kzs

!local variables:
integer:: iCell
integer:: iSoil

!-----------------------------------------------------------------------------------------------------------------

call mpas_pool_get_config(configs,'config_do_restart',config_do_restart)
call mpas_pool_get_config(configs,'config_lsm_scheme',config_lsm_scheme)
call mpas_pool_get_dimension(dims, 'nSoilLevels', nSoilLevels)

call mpas_pool_get_dimension(mesh,'nCells',nCells)

call mpas_pool_get_array(sfc_input,'landmask',landmask)
call mpas_pool_get_array(sfc_input,'dzs' , dzs )
call mpas_pool_get_array(sfc_input,'dzstop' , dzstop )
call mpas_pool_get_array(sfc_input,'kzs' , kzs )

!initialization of input variables, if needed:



!initialization of input variables, if needed:
if(.not. config_do_restart) then

lsm_select: select case(trim(config_lsm_scheme))
case("sf_noah","sf_noahmp")
! Report soil layers as defined during the initialisation step.
call mpas_log_write(' ')
call mpas_log_write('---~---')
call mpas_log_write(' Soil depth layers (physics_registry_init)')
call mpas_log_write('---~---')
iCell = min(2,nCells)
do iSoil=1, nSoilLevels
call mpas_log_write(' Layer $i - Thickness = $r' &
,intArgs=(/iSoil/),realArgs=(/dzs(iSoil,iCell)/))
end do
call mpas_log_write('---~---')
call mpas_log_write(' ')

case default

case("sf_noah","sf_noahmp")
!initialize the thickness of the soil layers for the Noah scheme:
do iCell = 1, nCells
dzs(1,iCell) = 0.10_RKIND
dzs(2,iCell) = 0.30_RKIND
dzs(3,iCell) = 0.60_RKIND
dzs(4,iCell) = 1.00_RKIND
enddo

case default

end select lsm_select

endif
Expand Down
6 changes: 1 addition & 5 deletions src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -842,14 +842,10 @@ subroutine lsm_noahmp_toMPAS(diag_physics,diag_physics_noahmp,output_noahmp,sfc_
snicexy(ns,i) = mpas_noahmp%snicexy(i,n)
snliqxy(ns,i) = mpas_noahmp%snliqxy(i,n)
enddo
do ns = 1,nsnow
do ns = 1,nzsnow
n = ns - nsnow
zsnsoxy(ns,i) = mpas_noahmp%zsnsoxy(i,n)
enddo
do ns = nsnow+1,nzsnow
n = ns - nsoil + 1
zsnsoxy(ns,i) = mpas_noahmp%zsnsoxy(i,n)
enddo
enddo


Expand Down
79 changes: 65 additions & 14 deletions src/core_atmosphere/physics/mpas_atmphys_initialize_real.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -285,9 +285,13 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)
integer :: iCell,iSoil
integer, pointer :: nCellsSolve,nSoilLevels,nFGSoilLevels
integer, pointer :: config_nsoillevels
real(kind=RKIND), pointer :: config_dzstop
real(kind=RKIND), pointer :: config_kzs

real(kind=RKIND),dimension(:,:),pointer:: dzs_fg,zs_fg
real(kind=RKIND),dimension(:,:),pointer:: dzs,zs
real(kind=RKIND),dimension(:),pointer:: dzstop,kzs
logical :: config_fine

!-----------------------------------------------------------------------------------------------------------------
!call mpas_log_write('')
Expand All @@ -301,11 +305,32 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)
call mpas_pool_get_array(fg, 'dzs_fg', dzs_fg)
call mpas_pool_get_array(fg, 'zs', zs)
call mpas_pool_get_array(fg, 'dzs', dzs)
call mpas_pool_get_array(fg, 'dzstop', dzstop)
call mpas_pool_get_array(fg, 'kzs', kzs)

call mpas_pool_get_config(configs, 'config_nsoillevels', config_nsoillevels)

if(config_nsoillevels .ne. 4) &
call physics_error_fatal('NOAH lsm uses 4 soil layers. Correct config_nsoillevels.')

call mpas_pool_get_config(configs, 'config_dzstop', config_dzstop)
call mpas_pool_get_config(configs, 'config_kzs' , config_kzs)

config_fine = ( config_nsoillevels >= 4 ) .and. ( config_dzstop >= 0.001 ) &
.and. ( config_kzs > 0.0 )

if (.not. config_fine) then
call mpas_log_write( ' ' )
call mpas_log_write( '---~---' )
call mpas_log_write( ' Invalid soil depth settings.' )
call mpas_log_write( '---~---' )
call mpas_log_write( ' config_nsoillevels = $i (It must be at least 4).' &
, intArgs=(/config_nsoillevels/))
call mpas_log_write( ' config_dztop = $r (It must be >= 0.001).' &
, realArgs=(/config_dzstop/))
call mpas_log_write( ' config_kzs = $r (It must be >= 0.0).' &
, realArgs=(/ config_kzs /))
call mpas_log_write(' ')
call physics_error_fatal(' At least one of the settings for soil depth is incorrect.')
end if

do iCell = 1, nCellsSolve
iSoil = 1
Expand All @@ -318,10 +343,11 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)
enddo

do iCell = 1, nCellsSolve
dzs(1,iCell) = 0.10_RKIND
dzs(2,iCell) = 0.30_RKIND
dzs(3,iCell) = 0.60_RKIND
dzs(4,iCell) = 1.00_RKIND
dzstop(iCell) = config_dzstop
kzs(iCell) = config_kzs
do iSoil = 1, nSoilLevels
dzs(iSoil,iCell) = dzstop(iCell) * exp( kzs(iCell) * real(iSoil-1,kind=RKIND) )
end do

iSoil = 1
zs(iSoil,iCell) = 0.5_RKIND * dzs(iSoil,iCell)
Expand All @@ -333,6 +359,18 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)

enddo

call mpas_log_write(' ')
call mpas_log_write('---~---')
call mpas_log_write(' Soil depth layers (init_soil_layers_depth)')
call mpas_log_write('---~---')
iCell = min(2,nCellsSolve)
do iSoil=1, nSoilLevels
call mpas_log_write(' Layer $i -- Thickness = $r -- MidDepth = $r ' &
,intArgs=(/iSoil/),realArgs=(/dzs(iSoil,iCell),zs(iSoil,iCell)/))
end do
call mpas_log_write('---~---')
call mpas_log_write(' ')

end subroutine init_soil_layers_depth

!=================================================================================================================
Expand All @@ -351,7 +389,7 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
type (mpas_pool_type), intent(inout) :: fg

!local variables:
integer:: iCell,ifgSoil,iSoil
integer:: iCell,ifgSoil,iSoil,iCellShow
integer, pointer:: nCellsSolve,nFGSoilLevels,nSoilLevels
integer:: num_sm,num_st
integer,dimension(:),pointer:: landmask
Expand Down Expand Up @@ -409,8 +447,8 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
call mpas_log_write('Error in interpolation of sm_fg to MPAS grid: num_sm = $i', messageType=MPAS_LOG_CRIT, intArgs=(/num_sm/))
endif

if(config_nsoillevels .ne. 4) &
call physics_error_fatal('NOAH lsm uses 4 soil layers. Correct config_nsoillevels.')
if(config_nsoillevels < 4) &
call physics_error_fatal('NOAH and Noah-MP require at least 4 soil layers. Correct config_nsoillevels.')

if(.not.allocated(zhave) ) allocate(zhave(nFGSoilLevels+2,nCellsSolve) )
if(.not.allocated(st_input)) allocate(st_input(nFGSoilLevels+2,nCellsSolve))
Expand All @@ -429,28 +467,29 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
sm_input(ifgSoil+1,iCell) = sm_fg(ifgSoil,iCell)
enddo

zhave(nFGSoilLevels+2,iCell) = 300._RKIND/100._RKIND
zhave(nFGSoilLevels+2,iCell) = max( zs(nSoilLevels,iCell) + dzs(nSoilLevels,iCell), 300._RKIND/100._RKIND)
st_input(nFGSoilLevels+2,iCell) = tmn(iCell)
sm_input(nFGSoilLevels+2,iCell) = sm_input(nFGSoilLevels,iCell)

if(iCell .eq. 1) then
do ifgSoil = 1,nFGSoilLevels+2
call mpas_log_write('$i $r', intArgs=(/ifgSoil/), realArgs=(/zhave(ifgSoil,iCell)/))
call mpas_log_write(' ifgSoil = $i -- zhave = $r', intArgs=(/ifgSoil/), realArgs=(/zhave(ifgSoil,iCell)/))
enddo
endif

enddo

!... interpolate the soil temperature, soil moisture, and soil liquid temperature to the four
! layers used in the NOAH land surface scheme:

iCellShow = 0
do iCell = 1, nCellsSolve

if(landmask(iCell) .eq. 1) then
if (iCellShow == 0) iCellShow = iCell

noah: do iSoil = 1 , nSoilLevels
input: do ifgSoil = 1 , nFGSoilLevels+2-1
if(iCell .eq. 1) call mpas_log_write('$i $i $r $r $r', &
if(iCell .eq. 1) call mpas_log_write(' iSoil = $i -- ifgSoil = $i -- zs = $r -- zhave = $r zhaveBelow = $r', &
intArgs=(/iSoil,ifgSoil/), &
realArgs=(/zs(iSoil,iCell), zhave(ifgSoil,iCell),zhave(ifgSoil+1,iCell)/))

Expand All @@ -461,7 +500,7 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
(st_input(ifgSoil,iCell) * (zhave(ifgSoil+1,iCell)-zs(iSoil,iCell)) &
+ st_input(ifgSoil+1,iCell) * (zs(iSoil,iCell)-zhave(ifgSoil,iCell))) &
/ (zhave(ifgSoil+1,iCell)-zhave(ifgSoil,iCell))
if(iCell .eq. 1) call mpas_log_write('$i $i $r $r $r', &
if(iCell .eq. 1) call mpas_log_write(' iSoil = $i -- ifgSoil = $i -- zs = $r -- zhave = $r zhaveBelow = $r', &
intArgs=(/iSoil,ifgSoil/), &
realArgs=(/zs(iSoil,iCell), zhave(ifgSoil,iCell),zhave(ifgSoil+1,iCell)/))

Expand Down Expand Up @@ -506,6 +545,18 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)

enddo

call mpas_log_write(' ')
call mpas_log_write('---~---')
call mpas_log_write(' Soil initial conditions (init_soil_layers_properties)')
call mpas_log_write('---~---')
iCell = max(1,iCellShow)
do iSoil=1, nSoilLevels
call mpas_log_write(' Layer $i -- Thickness = $r -- midDepth = $r -- Temperature = $r -- Moisture = $r' &
,intArgs=(/iSoil/),realArgs=(/dzs(iSoil,iCell),zs(iSoil,iCell),tslb(iSoil,iCell),smois(iSoil,iCell)/))
end do
call mpas_log_write('---~---')
call mpas_log_write(' ')

if(allocated(zhave) ) deallocate(zhave )
if(allocated(st_input)) deallocate(st_input)
if(allocated(sm_input)) deallocate(sm_input)
Expand Down
8 changes: 4 additions & 4 deletions src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,10 +114,10 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
if ( noahmp%config%nmlist%OptSoilProperty == 1 ) then
noahmp%config%domain%SoilType(1:NumSoilLayer) = NoahmpIO%ISLTYP(I) ! soil type same in all layers
elseif ( noahmp%config%nmlist%OptSoilProperty == 2 ) then
noahmp%config%domain%SoilType(1) = nint(NoahmpIO%SOILCL1(I)) ! soil type in layer1
noahmp%config%domain%SoilType(2) = nint(NoahmpIO%SOILCL2(I)) ! soil type in layer2
noahmp%config%domain%SoilType(3) = nint(NoahmpIO%SOILCL3(I)) ! soil type in layer3
noahmp%config%domain%SoilType(4) = nint(NoahmpIO%SOILCL4(I)) ! soil type in layer4
noahmp%config%domain%SoilType(1) = nint(NoahmpIO%SOILCL1(I)) ! soil type in layer1
noahmp%config%domain%SoilType(2) = nint(NoahmpIO%SOILCL2(I)) ! soil type in layer2
noahmp%config%domain%SoilType(3) = nint(NoahmpIO%SOILCL3(I)) ! soil type in layer3
noahmp%config%domain%SoilType(4:NumSoilLayer) = nint(NoahmpIO%SOILCL4(I)) ! soil type in layers 4 and below
elseif ( noahmp%config%nmlist%OptSoilProperty == 3 ) then
noahmp%config%domain%SoilType(1:NumSoilLayer) = NoahmpIO%ISLTYP(I) ! to initialize with default
endif
Expand Down
Loading