Skip to content

Commit

Permalink
iwp_subsidence bug fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@rfiorella
rfiorella committed Sep 26, 2024
1 parent d5780e3 commit 2d780f9
Show file tree
Hide file tree
Showing 3 changed files with 103 additions and 60 deletions.
83 changes: 61 additions & 22 deletions components/elm/src/biogeophys/ActiveLayerMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,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, dzsoi
use elm_varcon , only : zsoi, dzsoi, zisoi
!
! !ARGUMENTS:
integer , intent(in) :: num_soilc ! number of soil columns in filter
Expand All @@ -58,17 +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 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
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, orig_excess, old_mfrac ! temporary variables
real(r8), dimension(nlevgrnd) :: melt_profile ! profile of melted excess ice
!-----------------------------------------------------------------------

! RF NOTE: use of 1989 ALT in these parameterizations is somewhat of a placeholder used to compare against
Expand All @@ -87,11 +87,12 @@ subroutine alt_calc(num_soilc, filter_soilc, &
altmax_lastyear_indx => canopystate_vars%altmax_lastyear_indx_col , & ! Output: [integer (:) ] prior year maximum annual depth of thaw
altmax_1989_indx => canopystate_vars%altmax_1989_indx_col, & ! Output: [integer (:) ] index of maximum ALT in 1989
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 (-)
excess_ice => col_ws%excess_ice , & ! Input/output:[real(r8) (:,:)] depth variable excess ice content in soil column (-)
rmax => col_ws%iwp_microrel , & ! Output: [real(r8) (:) ] ice wedge polygon microtopographic relief (m)
vexc => col_ws%iwp_exclvol , & ! Output: [real(r8) (:) ] ice wedge polygon excluded volume (m)
ddep => col_ws%iwp_ddep , & ! Output: [real(r8) (:) ] ice wedge polygon depression depth (m)
subsidence => col_ws%iwp_subsidence & ! Input/output:[real(r8) (:) ] ice wedge polygon subsidence (m)
subsidence => col_ws%iwp_subsidence , & ! Input/output:[real(r8)(:)] ice wedge polygon subsidence (m)
frac_melted => col_ws%frac_melted & ! Input/output:[real(r8)(:)] fraction of layer that has ever melted (-)
)

! on a set annual timestep, update annual maxima
Expand Down Expand Up @@ -184,15 +185,49 @@ subroutine alt_calc(num_soilc, filter_soilc, &
! update subsidence based on change in ALT
! melt_profile stores the amount of excess_ice
! melted in this timestep.
do j = 1,nlevgrnd
if (j < k_frz) then
melt_profile(j) = 0.0_r8
! note that this may cause some unexpected results
! for taliks

! initialize melt_profile as zero
melt_profile(:) = 0._r8

do j = nlevgrnd,1,-1 ! note, this will go from bottom to top
!write(iulog,*) "processing level j,",j,k_frz,excess_ice(c,j)
if (j .gt. k_frz + 1) then ! all layers below k_frz + 1 remain frozen
melt_profile(j) = 0.0_r8
else if (j .eq. k_frz + 1) then ! first layer below the 'thawed' layer
! need to check to see if the active layer thickness is is actually
! in this layer (and not between the midpoint of j_frz and bottom interface
! or else inferred melt will be negative
! also note: only have ice to melt if alt has never been this deep, otherwise
! ice will continue to be removed each time step the alt remains in this layer
if ((alt(c)-zisoi(j-1)) .ge. 0._r8 .and. (alt(c) .eq. altmax_ever(c)) .and. (frac_melted(c,j) .lt. 1._r8)) then
orig_excess = (1._r8/(1._r8-frac_melted(c,j))) * excess_ice(c,j)
old_mfrac = frac_melted(c,j)
! update frac melted
frac_melted(c,j) = min(max(frac_melted(c,j), (alt(c)-zisoi(j-1))/dzsoi(j)),1._r8)
melt_profile(j) = orig_excess*(frac_melted(c,j) - old_mfrac)
excess_ice(c,j) = excess_ice(c,j) - melt_profile(j)
! DEBUG:
write(iulog,*) "melt and excess ice here are for j and j+1", melt_profile(j), &
melt_profile(j+1), excess_ice(c,j), &
excess_ice(c,j+1),(alt(c)-zisoi(j-1))/dzsoi(j),c,j+1,alt(c),zisoi(j-1),dzsoi(j)
else
melt_profile(j) = 0._r8 ! no melt
end if
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!!
! remove melted excess ice:
excess_ice(c,j) = 0._r8
excess_ice(c,j+1) = excess_ice(c,j+1)*(1._r8 - min(1._r8,(z2-alt(c))/(z2-z1)))
else
if (alt(c) .eq. altmax_ever(c) .and. (frac_melted(c,j) .lt. 1._r8)) then
orig_excess = (1._r8/(1._r8 - frac_melted(c,j))) * excess_ice(c,j)
old_mfrac = frac_melted(c,j)
! update frac_melted:
frac_melted(c,j) = min(max(frac_melted(c,j), (alt(c)-zsoi(j-1))/dzsoi(j)),1._r8)
! remove ice, only if alt has never been this deep before:
melt_profile(j) = orig_excess*(frac_melted(c,j) - old_mfrac)
excess_ice(c,j) = excess_ice(c,j) - melt_profile(j)
else
melt_profile(j) = 0._r8
end if
else !
melt_profile(j) = excess_ice(c,j)
! remove melted excess ice
excess_ice(c,j) = 0._r8
Expand All @@ -204,6 +239,10 @@ subroutine alt_calc(num_soilc, filter_soilc, &
! subsidence is integral of melt profile:
if ((year .ge. 1989) .and. (altmax_ever(c) .ge. altmax_1989(c))) then
subsidence(c) = subsidence(c) + sum(melt_profile)
if (sum(melt_profile) .lt. 0_r8) then
write(iulog,*) "subsidence(c) is",subsidence(c),sum(melt_profile),c,k_frz
write(iulog,*) "meltprofile is:", melt_profile
end if
end if

! limit subsidence to 0.4 m
Expand Down
2 changes: 1 addition & 1 deletion components/elm/src/biogeophys/SoilHydrologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -530,7 +530,7 @@ 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(iwp_subsidence(c), 0.4) !fix this variable when available to pull from alt calculations
phi_eff = min(iwp_subsidence(c), 0.4_r8) !fix this variable when available to pull from alt calculations
swc = h2osfc(c)/1000_r8 ! convert to m

if (swc >= vdep) then
Expand Down
78 changes: 41 additions & 37 deletions components/elm/src/data_types/ColumnDataType.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ module ColumnDataType
use soilorder_varcon, only : smax, ks_sorption
use elm_time_manager, only : is_restart, get_nstep
use elm_time_manager, only : is_first_step, get_step_size, is_first_restart_step
use landunit_varcon , only : istice, istwet, istsoil, istdlak, istcrop, istice_mec
use landunit_varcon , only : istice, istwet, istsoil, istdlak, istcrop, istice_mec, istlowcenpoly, isthighcenpoly
use column_varcon , only : icol_road_perv, icol_road_imperv, icol_roof, icol_sunwall, icol_shadewall
use histFileMod , only : hist_addfld1d, hist_addfld2d, no_snow_normal
use histFileMod , only : hist_addfld_decomp
Expand Down Expand Up @@ -110,7 +110,6 @@ module ColumnDataType
real(r8), pointer :: h2osoi_liq (:,:) => null() ! liquid water (-nlevsno+1:nlevgrnd) (kg/m2)
real(r8), pointer :: h2osoi_ice (:,:) => null() ! ice lens (-nlevsno+1:nlevgrnd) (kg/m2)
real(r8), pointer :: h2osoi_vol (:,:) => null() ! volumetric soil water (0<=h2osoi_vol<=watsat) (1:nlevgrnd) (m3/m3)
real(r8), pointer :: excess_ice (:,:) => null() ! NGEE Arctic: excess ground ice in column (1:nlevgrnd) (0 to 1)
real(r8), pointer :: h2osfc (:) => null() ! surface water (kg/m2)
real(r8), pointer :: h2ocan (:) => null() ! canopy water integrated to column (kg/m2)
real(r8), pointer :: total_plant_stored_h2o(:)=> null() ! total water in plants (kg/m2)
Expand Down Expand Up @@ -170,11 +169,13 @@ module ColumnDataType
real(r8), pointer :: vsfm_soilp_col_1d (:) => null() ! 1D soil liquid pressure from VSFM [Pa]
real(r8), pointer :: h2orof (:) => null() ! floodplain inundation volume received from rof (mm)
real(r8), pointer :: frac_h2orof (:) => null() ! floodplain inundation fraction received from rof (-)
! polygonal tundra
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)
! polygonal tundra (NGEE Arctic IM1)
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 :: excess_ice (:,:) => null() ! excess ground ice in column (1:nlevgrnd) (0 to 1)
real(r8), pointer :: frac_melted (:,:) => null() ! fraction of layer that has ever thawed (for tracking excess ice removal) (0 to 1)

contains
procedure, public :: Init => col_ws_init
Expand Down Expand Up @@ -1396,7 +1397,6 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_
allocate(this%h2osoi_liq (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liq (:,:) = spval
allocate(this%h2osoi_ice (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_ice (:,:) = spval
allocate(this%h2osoi_vol (begc:endc, 1:nlevgrnd)) ; this%h2osoi_vol (:,:) = spval
allocate(this%excess_ice (begc:endc, 1:nlevgrnd)) ; this%excess_ice (:,:) = spval
allocate(this%h2osfc (begc:endc)) ; this%h2osfc (:) = spval
allocate(this%h2ocan (begc:endc)) ; this%h2ocan (:) = spval
allocate(this%wslake_col (begc:endc)) ; this%wslake_col (:) = spval
Expand Down Expand Up @@ -1454,11 +1454,15 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_
allocate(this%vsfm_soilp_col_1d (ncells)) ; this%vsfm_soilp_col_1d (:) = spval
allocate(this%h2orof (begc:endc)) ; this%h2orof (:) = spval
allocate(this%frac_h2orof (begc:endc)) ; this%frac_h2orof (:) = 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
if (use_polygonal_tundra) then
! 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%frac_melted (begc:endc,1:nlevgrnd)) ; this%frac_melted (:,:) = spval
allocate(this%excess_ice (begc:endc,1:nlevgrnd)) ; this%excess_ice (:,:) = spval
end if

!-----------------------------------------------------------------------
! initialize history fields for select members of col_ws
Expand Down Expand Up @@ -1494,19 +1498,17 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_
avgflag='A', long_name='soil ice (ice landunits only)', &
ptr_col=this%h2osoi_ice, l2g_scale_type='ice')

! RPF - polygonal tundra vars
this%excess_ice(begc:endc, :) = spval
if (use_polygonal_tundra) then
this%iwp_subsidence(begc:endc) = spval
this%iwp_ddep(begc:endc) = spval
this%iwp_exclvol(begc:endc) = spval
this%iwp_microrel(begc:endc) = spval
this%frac_melted(begc:endc,:) = spval
this%excess_ice(begc:endc,:) = spval

call hist_addfld2d (fname='EXCESS_ICE', units = '1', type2d='levgrnd', &
avgflag='A', long_name='Excess ground ice (0 to 1)', &
ptr_col=this%excess_ice, l2g_scale_type='veg') ! <- RPF: should this be natveg?
end if

this%iwp_subsidence(begc:endc) = spval
this%iwp_ddep(begc:endc) = spval
this%iwp_exclvol(begc:endc) = spval
this%iwp_microrel(begc:endc) = spval
if (use_polygonal_tundra) then
call hist_addfld1d (fname="SUBSIDENCE", units='m', avgflag='A', &
long_name='ground subsidence (m)', ptr_col=this%iwp_subsidence)
call hist_addfld1d (fname="DEPRESS_DEPTH", units='m', avgflag='A', &
Expand All @@ -1516,6 +1518,9 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_
ptr_col=this%iwp_exclvol)
call hist_addfld1d (fname="MICROREL", units='m', avgflag='A', &
long_name='microtopographic relief (m)', ptr_col=this%iwp_microrel)
call hist_addfld2d (fname="FRAC_MELTED", units='-', type2d='levgrnd', &
avgflag='A', long_name='fraction of layer that has melted (-)', &
ptr_col=this%frac_melted, l2g_scale_type='veg')
endif
!/polygonal tundra

Expand Down Expand Up @@ -1679,6 +1684,7 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_
this%h2orof(c) = 0._r8
this%frac_h2orof(c) = 0._r8
this%iwp_subsidence(c) = 0._r8
this%frac_melted(c,:) = 0._r8

if (lun_pp%urbpoi(l)) then
! From Bonan 1996 (LSM technical note)
Expand Down Expand Up @@ -1731,13 +1737,16 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_
if (j > nlevbed) then
this%h2osoi_vol(c,j) = 0.0_r8
else
if (use_fates .or. use_hydrstress) then
if (use_fates .or. use_hydrstress) then
this%h2osoi_vol(c,j) = 0.70_r8*watsat_input(c,j) !0.15_r8 to avoid very dry conditions that cause errors in FATES
else if (use_arctic_init) then
this%h2osoi_vol(c,j) = watsat_input(c,j) ! start saturated for arctic
this%h2osoi_vol(c,j) = watsat_input(c,j) ! start saturated for arctic
else
this%h2osoi_vol(c,j) = 0.15_r8
endif
if (use_polygonal_tundra) then
this%frac_melted(c,j) = 0._r8
end if
endif
end do
else if (lun_pp%urbpoi(l)) then
Expand Down Expand Up @@ -1838,23 +1847,14 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_
this%h2osoi_ice(c,j) = 0._r8
this%h2osoi_liq(c,j) = col_pp%dz(c,j)*denh2o*this%h2osoi_vol(c,j)
endif
! RPF 240713 - notes from Chuck: initialize all to 0.36_r8
this%excess_ice(c,j) = 0.36_r8
! RPF - to do: a) apply to only polygonal ground
! b) pull in ALT information to get starting point right.
! for now, assuming no excess ice in top meter, 0.5 for 1-4 m,
! 0.2 below 4 m. Also: need to check signs!
!if (col_pp%z(c,j) > -1._r8) then
! this%excess_ice(c,j) = 0._r8
!else if (col_pp%z(c,j) > -4._r8 .and. col_pp%z(c,j) < -1._r8) then
! this%excess_ice(c,j) = 0.5_r8
!else
! this%excess_ice(c,j) = 0.2_r8
!end if
end do

this%h2osoi_liq_old(c,:) = this%h2osoi_liq(c,:)
this%h2osoi_ice_old(c,:) = this%h2osoi_ice(c,:)
if (use_polygonal_tundra) then
! RPF 240713 - notes from Chuck Abolt: initialize all to 0.36_r8
this%excess_ice(c,:) = 0.36_r8
end if
end do

end subroutine col_ws_init
Expand Down Expand Up @@ -1929,6 +1929,10 @@ subroutine col_ws_restart(this, bounds, ncid, flag, watsat_input)
dim1name='column', &
long_name='microtopographic relief', units='m', &
interpinic_flag='interp', readvar=readvar, data=this%iwp_microrel)
call restartvar(ncid=ncid, flag=flag, varname='FRAC_MELTED', xtype=ncd_double, &
dim1name='column', dim2name='levgrnd', switchdim=.true., &
long_name='fraction of layer that has ever melted', units='-', &
interpinic_flag='interp', readvar=readvar, data=this%frac_melted)
end if

call restartvar(ncid=ncid, flag=flag, varname='SOILP', xtype=ncd_double, &
Expand Down

0 comments on commit 2d780f9

Please sign in to comment.