From 7524a496be0cdaf09a7f6d6312da0f047190b728 Mon Sep 17 00:00:00 2001 From: Rich Fiorella Date: Wed, 10 Jul 2024 17:27:03 -0600 Subject: [PATCH] Read gradz from surfdata file when polygonal tundra is turned on. 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. --- components/elm/src/main/initVerticalMod.F90 | 86 ++++++++++++--------- 1 file changed, 50 insertions(+), 36 deletions(-) diff --git a/components/elm/src/main/initVerticalMod.F90 b/components/elm/src/main/initVerticalMod.F90 index 96e064b267a4..223314842cc0 100755 --- a/components/elm/src/main/initVerticalMod.F90 +++ b/components/elm/src/main/initVerticalMod.F90 @@ -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 @@ -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 @@ -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 !------------------------------------------------------------------------ @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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) @@ -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 @@ -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) @@ -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 !----------------------------------------------- @@ -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) @@ -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) @@ -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