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.

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

Adds polygonal tundra subsidence and micro topographic parameter calculations #6

Merged
merged 3 commits into from
May 17, 2024
+105 −43
Merged

Adds polygonal tundra subsidence and micro topographic parameter calculations #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
de55d88
wip - subsidence calculations in polygonal tundra
rfiorella May 15, 2024
f14ed6f
Add ground subsidence variable to SoilHydrologyMod
rfiorella May 15, 2024
8e93bbb
Add variable calculation of microtopographic parameters
rfiorella May 15, 2024
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
79 changes: 64 additions & 15 deletions components/elm/src/biogeophys/ActiveLayerMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,9 @@ module ActiveLayerMod
use CanopyStateType , only : canopystate_type
use GridcellType , only : grc_pp
use ColumnType , only : col_pp
use ColumnDataType , only : col_es
use ColumnDataType , only : col_es, col_ws
use LandunitType , only : lun_pp
use landunit_varcon , only : ilowcenpoly, iflatcenpoly, ihighcenpoly
!
implicit none
save
Expand Down Expand Up @@ -47,7 +49,7 @@ subroutine alt_calc(num_soilc, filter_soilc, &
use elm_varpar , only : nlevgrnd
use elm_time_manager , only : get_curr_date, get_step_size
use elm_varctl , only : iulog
use elm_varcon , only : zsoi
use elm_varcon , only : zsoi, dzsoi
!
! !ARGUMENTS:
integer , intent(in) :: num_soilc ! number of soil columns in filter
Expand All @@ -56,16 +58,17 @@ subroutine alt_calc(num_soilc, filter_soilc, &
type(canopystate_type) , intent(inout) :: canopystate_vars
!
! !LOCAL VARIABLES:
integer :: c, j, fc, g ! counters
integer :: alt_ind ! index of base of activel layer
integer :: year ! year (0, ...) for nstep+1
integer :: mon ! month (1, ..., 12) for nstep+1
integer :: day ! day of month (1, ..., 31) for nstep+1
integer :: sec ! seconds into current date for nstep+1
integer :: dtime ! time step length in seconds
integer :: k_frz ! index of first nonfrozen soil layer
logical :: found_thawlayer ! used to break loop when first unfrozen layer reached
real(r8) :: t1, t2, z1, z2 ! temporary variables
integer :: c, j, fc, g ! counters
integer :: alt_ind ! index of base of active layer
integer :: year ! year (0, ...) for nstep+1
integer :: mon ! month (1, ..., 12) for nstep+1
integer :: day ! day of month (1, ..., 31) for nstep+1
integer :: sec ! seconds into current date for nstep+1
integer :: dtime ! time step length in seconds
integer :: k_frz ! index of first nonfrozen soil layer
logical :: found_thawlayer ! used to break loop when first unfrozen layer reached
real(r8) :: t1, t2, z1, z2 ! temporary variables
real(r8), dimension(nlevgrnd) :: melt_profile ! profile of melted excess ice
!-----------------------------------------------------------------------

associate( &
Expand All @@ -78,7 +81,12 @@ subroutine alt_calc(num_soilc, filter_soilc, &
alt_indx => canopystate_vars%alt_indx_col , & ! Output: [integer (:) ] current depth of thaw
altmax_indx => canopystate_vars%altmax_indx_col , & ! Output: [integer (:) ] maximum annual depth of thaw
altmax_lastyear_indx => canopystate_vars%altmax_lastyear_indx_col , & ! Output: [integer (:) ] prior year maximum annual depth of thaw
altmax_ever_indx => canopystate_vars%altmax_ever_indx_col & ! Output: [integer (:) ] maximum thaw depth since initialization
altmax_ever_indx => canopystate_vars%altmax_ever_indx_col, & ! Output: [integer (:) ] maximum thaw depth since initialization
excess_ice => col_ws%excess_ice , & ! Input: [real(r8) (:,:) ] depth variable excess ice content in soil column (-)
rmax => col_pp%iwp_microrel , & ! Output: [real(r8) (:) ] ice wedge polygon microtopographic relief (m)
vexc => col_pp%iwp_exclvol , & ! Output: [real(r8) (:) ] ice wedge polygon excluded volume (m)
ddep => col_pp%iwp_ddep , & ! Output: [real(r8) (:) ] ice wedge polygon depression depth (m)
subsidence => col_pp%iwp_subsidence & ! Input/output:[real(r8) (:) ] ice wedge polygon subsidence (m)
)

! on a set annual timestep, update annual maxima
Expand Down Expand Up @@ -153,15 +161,56 @@ subroutine alt_calc(num_soilc, filter_soilc, &
altmax_indx(c) = alt_indx(c)
endif
if (alt(c) > altmax_ever(c)) then
if (spinup_state .eq. 0) then !overwrite if in spinup
if (spinup_state .eq. 0) then
altmax_ever(c) = alt(c)
altmax_ever_indx(c) = alt_indx(c)
else
else !overwrite if in spinup
altmax_ever(c) = 0._r8
altmax_ever_indx(c) = 0
endif
endif

! update subsidence based on change in ALT
! melt_profile stores the amount of excess_ice
! melted in this timestep.
do j = nlevgrnd-1,1,-1
if (j < k_frz) then
melt_profile(j) = 0.0_r8
else if (j .eq. k_frz) then
melt_profile(j) = excess_ice(c,j) + ((z2-alt(c))/(z2-z1))*excess_ice(c,j+1) ! TODO: check indices here!!
else
melt_profile(j) = excess_ice(c,j)
end if
! calculate subsidence at this layer:
melt_profile(j) = melt_profile(j) * dzsoi(j)
end do

! subsidence is integral of melt profile:
subsidence(c) = subsidence(c) + sum(melt_profile)

! limit subsidence to 0.4 m
subsidence(c) = min(0.4_r8, subsidence(c))

! update ice wedge polygon microtopographic parameters if in polygonal ground
!rf - min/max logic may be redunant w/ subsidence limiter above
if (lun_pp%ispolygon(col_pp%landunit(c))) then
if (lun_pp%polygontype(col_pp%landunit(c)) .eq. ilowcenpoly) then
rmax(c) = 0.4_r8
vexc(c) = 0.2_r8
ddep(c) = min(0.05_r8, max(0.15_r8 - 0.25_r8*subsidence(c), 0.15_r8))
elseif (lun_pp%polygontype(col_pp%landunit(c)) .eq. iflatcenpoly) then
rmax(c) = min(0.1_r8, max(0.4_r8, 0.1_r8 + 0.75_r8*subsidence(c)))
vexc(c) = min(0.05_r8, max(0.2_r8, 0.05_r8 + 0.375_r8*subsidence(c)))
ddep(c) = min(0.01_r8, max(0.05_r8, 0.01_r8 + 0.1_r8*subsidence(c)))
elseif (lun_pp%polygontype(col_pp%landunit(c)) .eq. ihighcenpoly) then
rmax(c) = 0.4_r8
vexc(c) = 0.2_r8
ddep(c) = 0.05_r8
else
!call endrun !<- TODO: needed? Potential way to prevent unintended updating of microtopography
! if polygonal ground is misspecified on surface file.
endif
endif
end do

end associate
Expand Down
29 changes: 15 additions & 14 deletions components/elm/src/biogeophys/SoilHydrologyMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
integer , intent(in) :: num_urbanc ! number of column urban points in column filter
integer , intent(in) :: filter_urbanc(:) ! column filter for urban points
type(atm2lnd_type) , intent(in) :: atm2lnd_vars ! land river two way coupling
type(lnd2atm_type) , intent(in) :: lnd2atm_vars
type(lnd2atm_type) , intent(in) :: lnd2atm_vars
type(energyflux_type) , intent(in) :: energyflux_vars
type(soilhydrology_type) , intent(inout) :: soilhydrology_vars
type(soilstate_type) , intent(inout) :: soilstate_vars
Expand Down Expand Up @@ -335,11 +335,12 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
snl => col_pp%snl , & ! Input: [integer (:) ] minus number of snow layers
dz => col_pp%dz , & ! Input: [real(r8) (:,:) ] layer depth (m)
nlev2bed => col_pp%nlevbed , & ! Input: [integer (:) ] number of layers to bedrock
cgridcell => col_pp%gridcell , & ! Input: [integer (:) ] column's gridcell
cgridcell => col_pp%gridcell , & ! Input: [integer (:) ] column's gridcell
wtgcell => col_pp%wtgcell , & ! Input: [real(r8) (:) ] weight (relative to gridcell)
iwp_microrel => col_pp%iwp_microrel , & ! Input: [real(r8) (:) ] ice wedge polygon microtopographic relief (m)
iwp_exclvol => col_pp%iwp_exclvol , & ! Input: [real(r8) (:) ] ice wedge polygon excluded volume (m)
iwp_ddep => col_pp%iwp_ddep , & ! Input: [real(r8) (:) ] ice wedge polygon depression depth (m)
iwp_subsidence => col_pp%iwp_subsidence , & ! Input: [real(r8) (:) ] ice wedge polygon ground subsidence (m)
meangradz => col_pp%meangradz , & ! Input: [real(r8) (:) ] mean topographic gradient at the column level (unitless)

t_soisno => col_es%t_soisno , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin)
Expand All @@ -365,7 +366,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
qflx_infl => col_wf%qflx_infl , & ! Output: [real(r8) (:) ] infiltration (mm H2O /s)
qflx_gross_infl_soil => col_wf%qflx_gross_infl_soil , & ! Output: [real(r8) (:)] gross infiltration (mm H2O/s)
qflx_gross_evap_soil => col_wf%qflx_gross_evap_soil , & ! Output: [real(r8) (:)] gross evaporation (mm H2O/s)
qflx_h2orof_drain => col_wf%qflx_h2orof_drain , & ! Output: [real(r8) (:)] drainange from floodplain inundation volume (mm H2O/s)
qflx_h2orof_drain => col_wf%qflx_h2orof_drain , & ! Output: [real(r8) (:)] drainange from floodplain inundation volume (mm H2O/s)

smpmin => soilstate_vars%smpmin_col , & ! Input: [real(r8) (:) ] restriction for min of soil potential (mm)
sucsat => soilstate_vars%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm)
Expand Down Expand Up @@ -408,7 +409,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
c = filter_hydrologyc(fc)
g = cgridcell(c)
pc = pc_grid(g)

! partition moisture fluxes between soil and h2osfc
if (lun_pp%itype(col_pp%landunit(c)) == istsoil .or. lun_pp%itype(col_pp%landunit(c))==istcrop) then

Expand All @@ -428,7 +429,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
h2orof(c) = atm2lnd_vars%h2orof_grc(g) * wtgcell(c)
frac_h2orof(c) = atm2lnd_vars%frac_h2orof_grc(g) * wtgcell(c)
endif
! TODO: add inundfrac from ocean
! TODO: add inundfrac from ocean
if ( frac_h2orof(c) > 1.0_r8 - fsno - frac_h2osfc(c) ) then
frac_h2orof(c) = 1.0_r8 - fsno - frac_h2osfc(c)
endif
Expand Down Expand Up @@ -481,7 +482,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
qinmax=(1._r8 - fsat(c)) * minval(10._r8**(-e_ice*(icefrac(c,1:3)))*hksat(c,1:3))
end if
end if

if ( use_modified_infil ) then
! Assume frac_h2osfc occurs on fsat
if ( frac_h2osfc(c) >= fsat(c) ) then
Expand All @@ -492,7 +493,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
else
qflx_infl_excess(c) = max(0._r8,qflx_in_soil(c) - (1.0_r8 - frac_h2osfc(c))*qinmax)
end if

if (use_lnd_rof_two_way) then
qflx_infl_excess(c) = max(0._r8,qflx_in_soil(c) - (1.0_r8 - frac_h2osfc(c) - frac_h2orof(c))*qinmax)
else
Expand All @@ -508,7 +509,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
if (h2osfcflag==1) then
! calculate runoff from h2osfc -------------------------------------
if (use_modified_infil) then
if (frac_h2osfc_act(c) <= pc .and. frac_h2osfc(c) <= pc) then
if (frac_h2osfc_act(c) <= pc .and. frac_h2osfc(c) <= pc) then
frac_infclust=0.0_r8
else
if (frac_h2osfc(c) <= pc) then
Expand All @@ -529,9 +530,9 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
if (lun_pp%ispolygon(col_pp%landunit(c))) then
vdep = (2_r8*iwp_exclvol(c) - iwp_microrel(c)) * (iwp_ddep(c)/iwp_microrel(c))**3_r8 &
+ (2_r8*iwp_microrel(c) - 3_r8*iwp_exclvol(c)) * (iwp_ddep(c)/iwp_microrel(c))**2_r8
phi_eff = min(subsidence, 0.4) !fix this variable when available to pull from alt calculations
phi_eff = min(iwp_subsidence(c), 0.4) !fix this variable when available to pull from alt calculations
swc = h2osfc(c)/1000_r8 ! convert to m

if (swc >= vdep) then
if (lun_pp%polygontype(col_pp%landunit(c)) == ilowcenpoly) then
k_wet = (2890_r8*phi_eff**4 - 1171.1_r8*phi_eff**3 + 144.94_r8*phi_eff**2 + 1.682_r8*phi_eff + 2.028) &
Expand All @@ -544,7 +545,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
else
qflx_h2osfc_surf(c) = 0._r8
endif

else
! limit runoff to value of storage above S(pc)
if(h2osfc(c) >= h2osfc_thresh(c) .and. h2osfcflag/=0) then
Expand Down Expand Up @@ -602,7 +603,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
!7. remove drainage from h2osfc and add to qflx_infl
h2osfc(c) = h2osfc(c) - qflx_h2osfc_drain(c) * dtime
qflx_infl(c) = qflx_infl(c) + qflx_h2osfc_drain(c)

!8. add drainage from river inundation to qflx_infl (land river two way coupling)
if (use_lnd_rof_two_way) then

Expand All @@ -625,8 +626,8 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
! remove drainage from inundation volume
h2orof(c) = h2orof(c) - qflx_h2orof_drain(c) * dtime

qflx_infl(c) = qflx_infl(c) + qflx_h2orof_drain(c)
qflx_gross_infl_soil(c) = qflx_gross_infl_soil(c) + qflx_h2osfc_drain(c) + qflx_h2orof_drain(c)
qflx_infl(c) = qflx_infl(c) + qflx_h2orof_drain(c)
qflx_gross_infl_soil(c) = qflx_gross_infl_soil(c) + qflx_h2osfc_drain(c) + qflx_h2orof_drain(c)

else
qflx_gross_infl_soil(c) = qflx_gross_infl_soil(c) + qflx_h2osfc_drain(c)
Expand Down
40 changes: 26 additions & 14 deletions components/elm/src/data_types/ColumnType.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,18 +49,19 @@ module ColumnType
logical , pointer :: active (:) => null() ! true=>do computations on this column

! topography
real(r8), pointer :: glc_topo (:) => null() ! surface elevation (m)
real(r8), pointer :: micro_sigma (:) => null() ! microtopography pdf sigma (m)
real(r8), pointer :: n_melt (:) => null() ! SCA shape parameter
real(r8), pointer :: topo_slope (:) => null() ! gridcell topographic slope
real(r8), pointer :: topo_std (:) => null() ! gridcell elevation standard deviation
real(r8), pointer :: hslp_p10 (:,:) => null() ! hillslope slope percentiles (unitless)
integer, pointer :: nlevbed (:) => null() ! number of layers to bedrock
real(r8), pointer :: zibed (:) => null() ! bedrock depth in model (interface level at nlevbed)
real(r8), pointer :: iwp_microrel (:) => null() ! ice wedge polygon microtopographic relief (m)
real(r8), pointer :: iwp_exclvol (:) => null() ! ice wedge polygon excluded volume (m)
real(r8), pointer :: iwp_ddep (:) => null() ! ice wedge polygon depression depth (m)
real(r8), pointer :: meangradz (:) => null() ! mean topographic gradient at the column level
real(r8), pointer :: glc_topo (:) => null() ! surface elevation (m)
real(r8), pointer :: micro_sigma (:) => null() ! microtopography pdf sigma (m)
real(r8), pointer :: n_melt (:) => null() ! SCA shape parameter
real(r8), pointer :: topo_slope (:) => null() ! gridcell topographic slope
real(r8), pointer :: topo_std (:) => null() ! gridcell elevation standard deviation
real(r8), pointer :: hslp_p10 (:,:) => null() ! hillslope slope percentiles (unitless)
integer, pointer :: nlevbed (:) => null() ! number of layers to bedrock
real(r8), pointer :: zibed (:) => null() ! bedrock depth in model (interface level at nlevbed)
real(r8), pointer :: iwp_microrel (:) => null() ! ice wedge polygon microtopographic relief (m)
real(r8), pointer :: iwp_exclvol (:) => null() ! ice wedge polygon excluded volume (m)
real(r8), pointer :: iwp_ddep (:) => null() ! ice wedge polygon depression depth (m)
real(r8), pointer :: iwp_subsidence(:) => null() ! ice wedge polygon ground subsidence (m)
real(r8), pointer :: meangradz (:) => null() ! mean topographic gradient at the column level

! vertical levels
integer , pointer :: snl (:) => null() ! number of snow layers
Expand Down Expand Up @@ -135,12 +136,18 @@ subroutine col_pp_init(this, begc, endc)
allocate(this%hslp_p10 (begc:endc,nlevslp)) ; this%hslp_p10 (:,:) = spval
allocate(this%nlevbed (begc:endc)) ; this%nlevbed (:) = ispval
allocate(this%zibed (begc:endc)) ; this%zibed (:) = spval
! polygonal tundra/ice wedge polygons:
allocate(this%iwp_microrel (begc:endc)) ; this%iwp_microrel (:) = spval
allocate(this%iwp_exclvol (begc:endc)) ; this%iwp_exclvol (:) = spval
allocate(this%iwp_ddep (begc:endc)) ; this%iwp_ddep (:) = spval
allocate(this%iwp_subsidence(begc:endc)) ; this%iwp_subsidence(:) = spval
allocate(this%meangradz (begc:endc)) ; this%meangradz (:) = spval

allocate(this%hydrologically_active(begc:endc)) ; this%hydrologically_active(:) = .false.

! Assume that columns are not fates columns until fates initialization begins
allocate(this%is_fates(begc:endc)); this%is_fates(:) = .false.

end subroutine col_pp_init

!------------------------------------------------------------------------
Expand Down Expand Up @@ -177,9 +184,14 @@ subroutine col_pp_clean(this)
deallocate(this%hslp_p10 )
deallocate(this%nlevbed )
deallocate(this%zibed )
deallocate(this%iwp_microrel)
deallocate(this%iwp_exclvol )
deallocate(this%iwp_ddep )
deallocate(this%iwp_subsidence)
deallocate(this%meangradz )
deallocate(this%hydrologically_active)
deallocate(this%is_fates)

end subroutine col_pp_clean

end module ColumnType