wip - first attempt at subsidence calc

components/elm/src/biogeophys/ActiveLayerMod.F90

@@ -48,7 +48,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
@@ -57,16 +57,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
+ 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( &
@@ -83,7 +84,8 @@ subroutine alt_calc(num_soilc, filter_soilc, &
  excess_ice => col_pp%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)
+ 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
@@ -158,15 +160,33 @@ 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 = k_frz) then
+ melt_profile(j) = excess_ice(j) + (z2-alt(c))/(z2-z1)*excess_ice(j+1) ! TODO: check indices here!!
+ else
+ melt_profile(j) = excess_ice(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)
+
  ! update ice wedge polygon microtopographic parameters if in polygonal ground
  ! TODO: need to retrieve landunit this column is on.
  if (lun_pp%ispolygon) then
@@ -181,7 +201,7 @@ subroutine alt_calc(num_soilc, filter_soilc, &
  elseif (lun_pp%polygontype(c) .eq. ihighcenpoly) then
  rmax(c) = 0.4_r8
  vexc(c) = 0.2_r8
- ddep(c) = 0.05_m
+ 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.