Skip to content

Commit

Permalink
Read gradz from surfdata file when polygonal tundra is turned on.
Browse files Browse the repository at this point in the history 
Adds logic to read gradz from surfdata file when NGEE Arctic
polygonal tundra is turned on. There are already topo_slope variables
so should make sure that this additional variable is actually needed.
  • Loading branch information
@rfiorella
rfiorella committed Jul 10, 2024
1 parent 496931b commit 7524a49
Showing 1 changed file with 50 additions and 36 deletions.
86 changes: 50 additions & 36 deletions components/elm/src/main/initVerticalMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,16 +14,16 @@ module initVerticalMod
use spmdMod , only : masterproc
use elm_varpar , only : more_vertlayers, nlevsno, nlevgrnd, nlevlak
use elm_varpar , only : toplev_equalspace, nlev_equalspace
use elm_varpar , only : nlevsoi, nlevsoifl, nlevurb, nlevslp
use elm_varpar , only : nlevsoi, nlevsoifl, nlevurb, nlevslp
use elm_varctl , only : fsurdat, iulog, use_var_soil_thick
use elm_varctl , only : use_vancouver, use_mexicocity, use_vertsoilc, use_extralakelayers, use_extrasnowlayers
use elm_varctl , only : use_erosion
use elm_varcon , only : zlak, dzlak, zsoi, dzsoi, zisoi, dzsoi_decomp, spval, grlnd
use elm_varctl , only : use_erosion, use_polygonal_tundra
use elm_varcon , only : zlak, dzlak, zsoi, dzsoi, zisoi, dzsoi_decomp, spval, grlnd
use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv
use landunit_varcon, only : istdlak, istice_mec
use fileutils , only : getfil
use LandunitType , only : lun_pp
use ColumnType , only : col_pp
use LandunitType , only : lun_pp
use ColumnType , only : col_pp
use SnowHydrologyMod, only : InitSnowLayers
use ncdio_pio
use topounit_varcon , only : max_topounits
Expand All @@ -50,13 +50,14 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
real(r8) , intent(in) :: thick_roof(bounds%begl:)
!
! LOCAL VARAIBLES:
integer :: c,l,t,ti,topi,g,i,j,lev ! indices
integer :: c,l,t,ti,topi,g,i,j,lev ! indices
type(file_desc_t) :: ncid ! netcdf id
logical :: readvar
logical :: readvar
integer :: dimid ! dimension id
character(len=256) :: locfn ! local filename
real(r8) ,pointer :: std (:) ! read in - topo_std
real(r8) ,pointer :: tslope (:) ! read in - topo_slope
real(r8) ,pointer :: std (:) ! read in - topo_std
real(r8) ,pointer :: tslope (:) ! read in - topo_slope
read(r8) ,pointer :: gradz(:) ! read in - gradz (polygonal tundra only)
real(r8) ,pointer :: hslp_p10 (:,:,:) ! read in - hillslope slope percentiles
real(r8) ,pointer :: dtb (:,:) ! read in - DTB
real(r8) :: beddep ! temporary
Expand All @@ -68,14 +69,14 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
integer :: ier ! error status
real(r8) :: scalez = 0.025_r8 ! Soil layer thickness discretization (m)
real(r8) :: thick_equal = 0.2
real(r8) ,pointer :: lakedepth_in(:,:) ! read in - lakedepth
real(r8) ,pointer :: lakedepth_in(:,:) ! read in - lakedepth
real(r8), allocatable :: zurb_wall(:,:) ! wall (layer node depth)
real(r8), allocatable :: zurb_roof(:,:) ! roof (layer node depth)
real(r8), allocatable :: dzurb_wall(:,:) ! wall (layer thickness)
real(r8), allocatable :: dzurb_roof(:,:) ! roof (layer thickness)
real(r8), allocatable :: ziurb_wall(:,:) ! wall (layer interface)
real(r8), allocatable :: ziurb_roof(:,:) ! roof (layer interface)
real(r8) :: depthratio ! ratio of lake depth to standard deep lake depth
real(r8) :: depthratio ! ratio of lake depth to standard deep lake depth
integer :: begc, endc
integer :: begl, endl
!------------------------------------------------------------------------
Expand All @@ -87,7 +88,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
SHR_ASSERT_ALL((ubound(thick_wall) == (/endl/)), errMsg(__FILE__, __LINE__))
SHR_ASSERT_ALL((ubound(thick_roof) == (/endl/)), errMsg(__FILE__, __LINE__))

! Open surface dataset to read in data below
! Open surface dataset to read in data below

call getfil (fsurdat, locfn, 0)
call ncd_pio_openfile (ncid, locfn, 0)
Expand All @@ -103,19 +104,19 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
end if

! --------------------------------------------------------------------
! Define layer structure for soil, lakes, urban walls and roof
! Define layer structure for soil, lakes, urban walls and roof
! Vertical profile of snow is not initialized here - but below
! --------------------------------------------------------------------

! Soil layers and interfaces (assumed same for all non-lake patches)
! "0" refers to soil surface and "nlevsoi" refers to the bottom of model soil

if ( more_vertlayers )then
! replace standard exponential grid with a grid that starts out exponential,
! then has several evenly spaced layers, then finishes off exponential.
! this allows the upper soil to behave as standard, but then continues
! replace standard exponential grid with a grid that starts out exponential,
! then has several evenly spaced layers, then finishes off exponential.
! this allows the upper soil to behave as standard, but then continues
! with higher resolution to a deeper depth, so that, for example, permafrost
! dynamics are not lost due to an inability to resolve temperature, moisture,
! dynamics are not lost due to an inability to resolve temperature, moisture,
! and biogeochemical dynamics at the base of the active layer
do j = 1, toplev_equalspace
zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8) !node depths
Expand Down Expand Up @@ -148,7 +149,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
zisoi(nlevgrnd) = zsoi(nlevgrnd) + 0.5_r8*dzsoi(nlevgrnd)

if (masterproc) then
write(iulog, *) 'zsoi', zsoi(:)
write(iulog, *) 'zsoi', zsoi(:)
write(iulog, *) 'zisoi: ', zisoi(:)
write(iulog, *) 'dzsoi: ', dzsoi(:)
end if
Expand All @@ -164,7 +165,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
dzsoi_decomp(1) = 1.
end if
if (masterproc) then
write(iulog, *) 'dzsoi_decomp', dzsoi_decomp(:)
write(iulog, *) 'dzsoi_decomp', dzsoi_decomp(:)
end if

if (nlevurb > 0) then
Expand All @@ -187,7 +188,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)

! "0" refers to urban wall/roof surface and "nlevsoi" refers to urban wall/roof bottom
if (lun_pp%urbpoi(l)) then
if (use_vancouver) then
if (use_vancouver) then
zurb_wall(l,1) = 0.010_r8/2._r8
zurb_wall(l,2) = zurb_wall(l,1) + 0.010_r8/2._r8 + 0.020_r8/2._r8
zurb_wall(l,3) = zurb_wall(l,2) + 0.020_r8/2._r8 + 0.070_r8/2._r8
Expand Down Expand Up @@ -281,13 +282,13 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)

dzurb_roof(l,1) = 0.5*(zurb_roof(l,1)+zurb_roof(l,2)) !thickness b/n two interfaces
do j = 2,nlevurb-1
dzurb_roof(l,j)= 0.5*(zurb_roof(l,j+1)-zurb_roof(l,j-1))
dzurb_roof(l,j)= 0.5*(zurb_roof(l,j+1)-zurb_roof(l,j-1))
enddo
dzurb_roof(l,nlevurb) = zurb_roof(l,nlevurb)-zurb_roof(l,nlevurb-1)

dzurb_wall(l,1) = 0.5*(zurb_wall(l,1)+zurb_wall(l,2)) !thickness b/n two interfaces
do j = 2,nlevurb-1
dzurb_wall(l,j)= 0.5*(zurb_wall(l,j+1)-zurb_wall(l,j-1))
dzurb_wall(l,j)= 0.5*(zurb_wall(l,j+1)-zurb_wall(l,j-1))
enddo
dzurb_wall(l,nlevurb) = zurb_wall(l,nlevurb)-zurb_wall(l,nlevurb-1)

Expand Down Expand Up @@ -362,7 +363,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
t = col_pp%topounit(c)
topi = grc_pp%topi(g)
ti = t - topi + 1

col_pp%lakedepth(c) = lakedepth_in(g,ti)
end do
deallocate(lakedepth_in)
Expand Down Expand Up @@ -435,7 +436,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)

else if (col_pp%lakedepth(c) > 1._r8 .and. col_pp%lakedepth(c) < 5000._r8) then

depthratio = col_pp%lakedepth(c) / (zlak(nlevlak) + 0.5_r8*dzlak(nlevlak))
depthratio = col_pp%lakedepth(c) / (zlak(nlevlak) + 0.5_r8*dzlak(nlevlak))
col_pp%z_lake(c,1) = zlak(1)
col_pp%dz_lake(c,1) = dzlak(1)
col_pp%dz_lake(c,2:nlevlak-1) = dzlak(2:nlevlak-1)*depthratio
Expand Down Expand Up @@ -466,12 +467,12 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
end do

!-----------------------------------------------
! Set cold-start values for snow levels, snow layers and snow interfaces
! Set cold-start values for snow levels, snow layers and snow interfaces
!-----------------------------------------------

if (.not. use_extrasnowlayers) then
associate(snl => col_pp%snl) ! Output: [integer (:) ] number of snow layers

associate(snl => col_pp%snl) ! Output: [integer (:) ] number of snow layers

do c = bounds%begc,bounds%endc
l = col_pp%landunit(c)
Expand Down Expand Up @@ -547,12 +548,12 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
col_pp%zi(c,-nlevsno+0:0) = 0._r8
end if
end do

end associate

else ! use extra (16) snow layers, for firn model
call InitSnowLayers(bounds, snow_depth(bounds%begc:bounds%endc))

end if

!-----------------------------------------------
Expand All @@ -563,7 +564,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
call ncd_io(ncid=ncid, varname='SLOPE', flag='read', data=tslope, dim1name=grlnd, readvar=readvar)
if (.not. readvar) then
call shr_sys_abort(' ERROR: TOPOGRAPHIC SLOPE NOT on surfdata file'//&
errMsg(__FILE__, __LINE__))
errMsg(__FILE__, __LINE__))
end if
do c = begc,endc
g = col_pp%gridcell(c)
Expand All @@ -572,11 +573,24 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
end do
deallocate(tslope)

if (use_polygonal_tundra) then
allocate(gradz(bounds%begg:bounds%endg))
call ncd_io(ncid=ncid, varname='GRADZ', flag='read', data=gradz, dim1name=grlnd, readvar=readvar)
if (.not. readvar) then
call shr_sys_abort(' ERROR: polygonal tundra turned on, but GRADZ not on surfdata file'//&
errMsg(__FILE__, __LINE__))
end if
do c = begc,endc
g = col_pp%gridcell(c)
col_pp%meangradz(c) = gradz(g)
end do
end if

allocate(std(bounds%begg:bounds%endg))
call ncd_io(ncid=ncid, varname='STD_ELEV', flag='read', data=std, dim1name=grlnd, readvar=readvar)
if (.not. readvar) then
call shr_sys_abort(' ERROR: TOPOGRAPHIC STDdev (STD_ELEV) NOT on surfdata file'//&
errMsg(__FILE__, __LINE__))
errMsg(__FILE__, __LINE__))
end if
do c = begc,endc
g = col_pp%gridcell(c)
Expand Down Expand Up @@ -623,11 +637,11 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
else
do c = begc,endc
g = col_pp%gridcell(c)
l = col_pp%landunit(c)
l = col_pp%landunit(c)
t = col_pp%topounit(c)
topi = grc_pp%topi(g)
ti = t - topi + 1

if (lun_pp%urbpoi(l) .and. col_pp%itype(c) /= icol_road_imperv .and. col_pp%itype(c) /= icol_road_perv) then
col_pp%nlevbed(c) = nlevurb
else if (lun_pp%itype(l) == istdlak) then
Expand Down

0 comments on commit 7524a49

Please sign in to comment.