Skip to content

Commit 0b0a6b9

Browse files
rfiorellarfiorella
committed
wip - subsidence calculations in polygonal tundra
Adds infrastructure for calculating subsidence and related microtopographic parameters based on changes in MAXALT over time.
1 parent cb668de commit 0b0a6b9

File tree

1 file changed

+59
-14
lines changed

1 file changed

+59
-14
lines changed

‎components/elm/src/biogeophys/ActiveLayerMod.F90

Lines changed: 59 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@ module ActiveLayerMod
1313
use GridcellType , only : grc_pp
1414
use ColumnType , only : col_pp
1515
use ColumnDataType , only : col_es
16+
use LandunitType , only : lun_pp
1617
!
1718
implicit none
1819
save
@@ -47,7 +48,7 @@ subroutine alt_calc(num_soilc, filter_soilc, &
4748
use elm_varpar , only : nlevgrnd
4849
use elm_time_manager , only : get_curr_date, get_step_size
4950
use elm_varctl , only : iulog
50-
use elm_varcon , only : zsoi
51+
use elm_varcon , only : zsoi, dzsoi
5152
!
5253
! !ARGUMENTS:
5354
integer , intent(in) :: num_soilc ! number of soil columns in filter
@@ -56,16 +57,17 @@ subroutine alt_calc(num_soilc, filter_soilc, &
5657
type(canopystate_type) , intent(inout) :: canopystate_vars
5758
!
5859
! !LOCAL VARIABLES:
59-
integer :: c, j, fc, g ! counters
60-
integer :: alt_ind ! index of base of activel layer
61-
integer :: year ! year (0, ...) for nstep+1
62-
integer :: mon ! month (1, ..., 12) for nstep+1
63-
integer :: day ! day of month (1, ..., 31) for nstep+1
64-
integer :: sec ! seconds into current date for nstep+1
65-
integer :: dtime ! time step length in seconds
66-
integer :: k_frz ! index of first nonfrozen soil layer
67-
logical :: found_thawlayer ! used to break loop when first unfrozen layer reached
68-
real(r8) :: t1, t2, z1, z2 ! temporary variables
60+
integer :: c, j, fc, g ! counters
61+
integer :: alt_ind ! index of base of active layer
62+
integer :: year ! year (0, ...) for nstep+1
63+
integer :: mon ! month (1, ..., 12) for nstep+1
64+
integer :: day ! day of month (1, ..., 31) for nstep+1
65+
integer :: sec ! seconds into current date for nstep+1
66+
integer :: dtime ! time step length in seconds
67+
integer :: k_frz ! index of first nonfrozen soil layer
68+
logical :: found_thawlayer ! used to break loop when first unfrozen layer reached
69+
real(r8) :: t1, t2, z1, z2 ! temporary variables
70+
real(r8), dimension(nlevgrnd) :: melt_profile ! profile of melted excess ice
6971
!-----------------------------------------------------------------------
7072

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

8491
! on a set annual timestep, update annual maxima
@@ -153,15 +160,53 @@ subroutine alt_calc(num_soilc, filter_soilc, &
153160
altmax_indx(c) = alt_indx(c)
154161
endif
155162
if (alt(c) > altmax_ever(c)) then
156-
if (spinup_state .eq. 0) then !overwrite if in spinup
163+
if (spinup_state .eq. 0) then
157164
altmax_ever(c) = alt(c)
158165
altmax_ever_indx(c) = alt_indx(c)
159-
else
166+
else !overwrite if in spinup
160167
altmax_ever(c) = 0._r8
161168
altmax_ever_indx(c) = 0
162169
endif
163170
endif
164171

172+
! update subsidence based on change in ALT
173+
! melt_profile stores the amount of excess_ice
174+
! melted in this timestep.
175+
do j = nlevgrnd-1,1,-1
176+
if (j < k_frz) then
177+
melt_profile(j) = 0.0_r8
178+
else if (j .eq. k_frz) then
179+
melt_profile(j) = excess_ice(c,j) + ((z2-alt(c))/(z2-z1))*excess_ice(c,j+1) ! TODO: check indices here!!
180+
else
181+
melt_profile(j) = excess_ice(c,j)
182+
end if
183+
! calculate subsidence at this layer:
184+
melt_profile(j) = melt_profile(j) * dzsoi(j)
185+
end do
186+
187+
! subsidence is integral of melt profile:
188+
subsidence(c) = subsidence(c) + sum(melt_profile)
189+
190+
! update ice wedge polygon microtopographic parameters if in polygonal ground
191+
! TODO: need to retrieve landunit this column is on.
192+
if (lun_pp%ispolygon) then
193+
if (lun_pp%polygontype(c) .eq. ilowcenpoly) then
194+
rmax(c) = 0.4_r8
195+
vexc(c) = 0.2_r8
196+
ddep(c) = 0.15_r8 ! TODO - update based on subsidence calcs.
197+
elseif (lun_pp%polygontype(c) .eq. iflatcenpoly) then
198+
rmax(c) = 0.1_r8 ! TODO - update based on subsidence calcs.
199+
vexc(c) = 0.05_r8 ! TODO - update based on subsidence calcs.
200+
ddep(c) = 0.01_r8 ! TODO - update based on subsidence calcs.
201+
elseif (lun_pp%polygontype(c) .eq. ihighcenpoly) then
202+
rmax(c) = 0.4_r8
203+
vexc(c) = 0.2_r8
204+
ddep(c) = 0.05_r8
205+
else
206+
!call endrun !<- TODO: needed? Potential way to prevent unintended updating of microtopography
207+
! if polygonal ground is misspecified on surface file.
208+
endif
209+
endif
165210
end do
166211

167212
end associate

0 commit comments

Comments
 (0)