Merge pull request #5 from lanl/cabolt/IM1

wip - Includes partial implementation of polygonal tundra subtype land unit subtype of soil land unit types.
Includes new microtopographic parameters, calculation of inundation fraction, and calculation of subsidence associated with melting of excess ice.

components/elm/src/biogeophys/CanopyHydrologyMod.F90

@@ -794,24 +794,28 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, &
  integer , intent(in), optional :: no_update ! flag to make calculation w/o updating variables
  !
  ! !LOCAL VARIABLES:
- integer :: c,f,l ! indices
- real(r8):: d,fd,dfdd ! temporary variable for frac_h2oscs iteration
- real(r8):: sigma ! microtopography pdf sigma in mm
+ integer :: c,f,l,k ! indices
+ real(r8):: d,fd,dfdd ! temporary variable for frac_h2oscs iteration
+ real(r8):: sigma ! microtopography pdf sigma in mm
+ real(r8):: swc ! surface water content in m
  real(r8):: min_h2osfc
  !-----------------------------------------------------------------------
 
  associate( &
- micro_sigma => col_pp%micro_sigma  , & ! Input: [real(r8) (:) ] microtopography pdf sigma (m)
+ micro_sigma => col_pp%micro_sigma , & ! Input: [real(r8) (:) ] microtopography pdf sigma (m)
 
- h2osno => col_ws%h2osno , & ! Input: [real(r8) (:) ] snow water (mm H2O)
- excess_ice => col_ws%excess_ice , & ! Input: [real(r8) (:,:) ] excess ground ice [kg/m2]
+ 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)
 
- h2osoi_liq => col_ws%h2osoi_liq , & ! Output: [real(r8) (:,:) ] liquid water (col,lyr) [kg/m2]
- h2osfc => col_ws%h2osfc , & ! Output: [real(r8) (:) ] surface water (mm)
- frac_sno => col_ws%frac_sno , & ! Output: [real(r8) (:) ] fraction of ground covered by snow (0 to 1)
- frac_sno_eff => col_ws%frac_sno_eff , & ! Output: [real(r8) (:) ] eff. fraction of ground covered by snow (0 to 1)
- frac_h2osfc => col_ws%frac_h2osfc , & ! Output: [real(r8) (:) ] col fractional area with surface water greater than zero
- frac_h2osfc_act => col_ws%frac_h2osfc_act & ! Output: [real(r8) (:) ] col fractional area with surface water greater than zero
+ h2osno => col_ws%h2osno , & ! Input: [real(r8) (:) ] snow water (mm H2O)
+
+ h2osoi_liq => col_ws%h2osoi_liq , & ! Output: [real(r8) (:,:) ] liquid water (col,lyr) [kg/m2]
+ h2osfc => col_ws%h2osfc , & ! Output: [real(r8) (:) ] surface water (mm)
+ frac_sno => col_ws%frac_sno , & ! Output: [real(r8) (:) ] fraction of ground covered by snow (0 to 1)
+ frac_sno_eff => col_ws%frac_sno_eff , & ! Output: [real(r8) (:) ] eff. fraction of ground covered by snow (0 to 1)
+ frac_h2osfc => col_ws%frac_h2osfc , & ! Output: [real(r8) (:) ] col fractional area with surface water greater than zero
+ frac_h2osfc_act => col_ws%frac_h2osfc_act & ! Output: [real(r8) (:) ] col fractional area with surface water greater than zero
  )
 
  ! arbitrary lower limit on h2osfc for safer numerics...
@@ -826,25 +830,50 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, &
 
  ! Use newton-raphson method to iteratively determine frac_h20sfc
  ! based on amount of surface water storage (h2osfc) and
- ! microtopography variability (micro_sigma)
+ ! microtopography variability (micro_sigma) if nonpolygonal, or
+ ! amount of surface water (h2osfc) and microtopographic attributes
+ ! (iwp_microrel, iwp_exclvol) if polygonal tundra
 
  if (h2osfc(c) > min_h2osfc) then
  ! a cutoff is needed for numerical reasons...(nonconvergence after 5 iterations)
- d=0.0
-
- sigma=1.0e3 * micro_sigma(c) ! convert to mm
- do l=1,10
- fd = 0.5*d*(1.0_r8+erf(d/(sigma*sqrt(2.0)))) &
- +sigma/sqrt(2.0*shr_const_pi)*exp(-d**2/(2.0*sigma**2)) &
- -h2osfc(c)
- dfdd = 0.5*(1.0_r8+erf(d/(sigma*sqrt(2.0))))
-
- d = d - fd/dfdd
- enddo
- !-- update the submerged areal fraction using the new d value
- frac_h2osfc(c) = 0.5*(1.0_r8+erf(d/(sigma*sqrt(2.0))))
 
- else
+ if (lun_pp%ispolygon(l)) then
+ ! calculate water depth and inundation fraction if column is polygonal
+ swc = h2osfc(c)/1000_r8 ! convert to m
+
+ if (swc > iwp_microrel(c) - iwp_exclvol(c)) then
+ d = swc + iwp_exclvol(c)
+ else
+ d = 0.0
+ do k=1,10
+ fd = (2_r8*iwp_exclvol(c) - iwp_microrel(c)) * (d/iwp_microrel(c))**3_r8 &
+ + (2_r8*iwp_microrel(c) - 3_r8*iwp_exclvol(c)) * (d/iwp_microrel(c))**2_r8 &
+ - swc
+ dfdd = (3_r8/iwp_microrel(c)) * (2_r8*iwp_exclvol(c) - iwp_microrel(c)) * (d/iwp_microrel(c))**2_r8 &
+ + (2_r8/iwp_microrel(c)) * (2_r8*iwp_microrel(c) - 3_r8*iwp_exclvol(c)) * (d/iwp_microrel(c))
+ d = d - fd/dfdd
+ enddo
+ endif
+
+ !-- update the submerged areal fraction using the new d value
+ frac_h2osfc(c) = (3_r8/iwp_microrel(c)) * (2_r8*iwp_exclvol(c) - iwp_microrel(c)) * (d/iwp_microrel(c))**2_r8 &
+ + (2_r8/iwp_microrel(c)) * (2_r8*iwp_microrel(c) - 3_r8*iwp_exclvol(c)) * (d/iwp_microrel(c))
+
+ else ! calculate water depth and inudation fraction if column is non-polygonal
+ d=0.0
+ sigma=1.0e3 * micro_sigma(c) ! convert to mm
+ do k=1,10
+ fd = 0.5*d*(1.0_r8+erf(d/(sigma*sqrt(2.0)))) &
+ +sigma/sqrt(2.0*shr_const_pi)*exp(-d**2/(2.0*sigma**2)) &
+ -h2osfc(c)
+ dfdd = 0.5*(1.0_r8+erf(d/(sigma*sqrt(2.0))))
+
+ d = d - fd/dfdd
+ enddo
+ !-- update the submerged areal fraction using the new d value
+ frac_h2osfc(c) = 0.5*(1.0_r8+erf(d/(sigma*sqrt(2.0))))
+ endif ! end if polygonal test
+ else ! if h2osfc(c) <= min_h2osfc
  frac_h2osfc(c) = 0._r8
  h2osoi_liq(c,1) = h2osoi_liq(c,1) + h2osfc(c)
  qflx_h2osfc2topsoi(c) = h2osfc(c)/dtime
@@ -867,7 +896,6 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, &
  frac_sno_eff(c)=frac_sno(c)
 
  endif
-
  endif ! end of no_update construct
 
  else !if landunit not istsoil/istcrop, set frac_h2osfc to zero

components/elm/src/biogeophys/SoilHydrologyMod.F90

@@ -201,13 +201,18 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
 
  do fc = 1, num_hydrologyc
  c = filter_hydrologyc(fc)
-
- ! assume qinmax large relative to qflx_top_soil in control
- if (origflag == 1) then
- qflx_surf(c) =  fcov(c) * qflx_top_soil(c)
+ l = col_pp%landunit(c)
+ ! no qflx_surf in polygonal ground
+ if (lun_pp%ispolygon(l)) then
+ qflx_surf(c) = 0
  else
- ! only send fast runoff directly to streams
- qflx_surf(c) = fsat(c) * qflx_top_soil(c)
+ ! assume qinmax large relative to qflx_top_soil in control
+ if (origflag == 1) then
+ qflx_surf(c) = fcov(c) * qflx_top_soil(c)
+ else
+ ! only send fast runoff directly to streams
+ qflx_surf(c) = fsat(c) * qflx_top_soil(c)
+ endif
  endif
  end do
 
@@ -267,7 +272,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
  use elm_varcon , only : denh2o, denice, roverg, wimp, mu, tfrz
  use elm_varcon , only : pondmx, watmin
  use column_varcon , only : icol_roof, icol_road_imperv, icol_sunwall, icol_shadewall, icol_road_perv
- use landunit_varcon , only : istsoil, istcrop
+ use landunit_varcon , only : istsoil, istcrop, ilowcenpoly
  use elm_time_manager , only : get_step_size, get_nstep
  use atm2lndType , only : atm2lnd_type ! land river two way coupling
  use lnd2atmType , only : lnd2atm_type
@@ -313,6 +318,9 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
  real(r8) :: d
  real(r8) :: h2osoi_vol
  real(r8) :: basis ! temporary, variable soil moisture holding capacity
+ real(r8) :: vdep ! temporary, ice wedge polygon volumetric depression depth (m)
+ real(r8) :: phi_eff ! temporary, polygonal ground effective subsidence (maxes out at 0.4) (m)
+ real(r8) :: swc ! temporary, polygonal surface water content in m
  ! in top VIC layers for runoff calculation
  real(r8) :: rsurf_vic ! temp VIC surface runoff
  real(r8) :: top_moist(bounds%begc:bounds%endc) ! temporary, soil moisture in top VIC layers
@@ -328,7 +336,11 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
  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 
- wtgcell => col_pp%wtgcell , & ! Input: [real(r8) (:) ] weight (relative to 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)
+ 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)
 
@@ -514,15 +526,36 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
  endif
  endif
 
- ! limit runoff to value of storage above S(pc)
- if(h2osfc(c) >= h2osfc_thresh(c) .and. h2osfcflag/=0) then
- ! spatially variable k_wet
- k_wet=1.0_r8 * sin((rpi/180.) * col_pp%topo_slope(c))
- qflx_h2osfc_surf(c) = k_wet * frac_infclust * (h2osfc(c) - h2osfc_thresh(c))
-
- qflx_h2osfc_surf(c)=min(qflx_h2osfc_surf(c),(h2osfc(c) - h2osfc_thresh(c))/dtime)
+ 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
+ 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) &
+ * (710.3_r8*meangradz(c)**2 - 28.736_r8*meangradz(c) + 12.74_r8)
+ else
+ k_wet = 24.925_r8 * (710.3_r8*meangradz(c)**2 - 28.736_r8*meangradz(c) + 12.74_r8)
+ endif
+ qflx_h2osfc_surf(c) = k_wet * (swc - vdep)
+ qflx_h2osfc_surf(c) = min(qflx_h2osfc_surf(c), (swc - vdep)*1000_r8/dtime)
+ else
+ qflx_h2osfc_surf(c) = 0._r8
+ endif
+
  else
- qflx_h2osfc_surf(c)= 0._r8
+ ! limit runoff to value of storage above S(pc)
+ if(h2osfc(c) >= h2osfc_thresh(c) .and. h2osfcflag/=0) then
+ ! spatially variable k_wet
+ k_wet=1.0_r8 * sin((rpi/180.) * col_pp%topo_slope(c))
+ qflx_h2osfc_surf(c) = k_wet * frac_infclust * (h2osfc(c) - h2osfc_thresh(c))
+
+ qflx_h2osfc_surf(c)=min(qflx_h2osfc_surf(c),(h2osfc(c) - h2osfc_thresh(c))/dtime)
+ else
+ qflx_h2osfc_surf(c)= 0._r8
+ endif
  endif
 
  ! cutoff lower limit
@@ -696,7 +729,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, fil
  real(r8) :: q_perch
  real(r8) :: q_perch_max
  real(r8) :: dflag=0._r8
-  real(r8) :: qcharge_temp
+  real(r8) :: qcharge_temp
  !-----------------------------------------------------------------------
 
  associate( &
@@ -783,23 +816,23 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, fil
  ! Water table changes due to qcharge
  do fc = 1, num_hydrologyc
  c = filter_hydrologyc(fc)
-   nlevbed = nlev2bed(c)
+   nlevbed = nlev2bed(c)
 
  !scs: use analytical expression for aquifer specific yield
  rous = watsat(c,nlevbed) &
  * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevbed))**(-1./bsw(c,nlevbed)))
  rous=max(rous,0.02_r8)
 
  !-- water table is below the soil column --------------------------------------
-  g = col_pp%gridcell(c)
+  g = col_pp%gridcell(c)
  l = col_pp%landunit(c)
  qcharge_temp = qcharge(c)
 
  wa(c) = wa(c) - qflx_grnd_irrig_col(c) * dtime
  zwt(c) = zwt(c) + (qflx_grnd_irrig_col(c) * dtime)/1000._r8/rous
 
  if(jwt(c) == nlevbed) then
-  if (.not. (zengdecker_2009_with_var_soil_thick)) then
+  if (.not. (zengdecker_2009_with_var_soil_thick)) then
  wa(c) = wa(c) + qcharge(c) * dtime
  zwt(c) = zwt(c) - (qcharge(c) * dtime)/1000._r8/rous
  end if
@@ -865,7 +898,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, fil
  ! perched water table code
  do fc = 1, num_hydrologyc
  c = filter_hydrologyc(fc)
-   nlevbed = nlev2bed(c)
+   nlevbed = nlev2bed(c)
 
  ! define frost table as first frozen layer with unfrozen layer above it
  if(t_soisno(c,1) > tfrz) then
@@ -1132,7 +1165,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
 
  do fc = 1, num_hydrologyc
  c = filter_hydrologyc(fc)
-   nlevbed = nlev2bed(c)
+   nlevbed = nlev2bed(c)
  jwt(c) = nlevbed
  ! allow jwt to equal zero when zwt is in top layer
  do j = 1,nlevbed
@@ -1142,7 +1175,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
  else
  jwt(c) = j-1
  exit
-  end if
+  end if
  end if
  enddo
  end do
@@ -1153,7 +1186,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
  ! perched water table code
  do fc = 1, num_hydrologyc
  c = filter_hydrologyc(fc)
-   nlevbed = nlev2bed(c)
+   nlevbed = nlev2bed(c)
 
  ! specify maximum drainage rate
  q_perch_max = 1.e-5_r8 * sin(col_pp%topo_slope(c) * (rpi/180._r8))
@@ -1357,11 +1390,11 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
  ! make sure baseflow isn't negative
  rsub_top(c) = max(0._r8, rsub_top(c))
  else
-  if (jwt(c) == nlevbed .and. zengdecker_2009_with_var_soil_thick) then
+  if (jwt(c) == nlevbed .and. zengdecker_2009_with_var_soil_thick) then
  rsub_top(c) = 0._r8
  else
  rsub_top(c) = imped * rsub_top_max* exp(-fff(c)*zwt(c))
- end if
+ end if
  end if
 
  if (use_vsfm) rsub_top(c) = 0._r8
@@ -1373,10 +1406,10 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
 
  !-- water table is below the soil column --------------------------------------
  if(jwt(c) == nlevbed) then
-  if (zengdecker_2009_with_var_soil_thick) then
-   if (-1._r8 * smp_l(c,nlevbed) < 0.5_r8 * dzmm(c,nlevbed)) then
-   zwt(c) = z(c,nlevbed) - (smp_l(c,nlevbed) / 1000._r8)
-  end if
+  if (zengdecker_2009_with_var_soil_thick) then
+   if (-1._r8 * smp_l(c,nlevbed) < 0.5_r8 * dzmm(c,nlevbed)) then
+   zwt(c) = z(c,nlevbed) - (smp_l(c,nlevbed) / 1000._r8)
+  end if
  rsub_top(c) = imped * rsub_top_max * exp(-fff(c) * zwt(c))
  rsub_top_tot = - rsub_top(c) * dtime
  s_y = watsat(c,nlevbed) &
@@ -1391,8 +1424,8 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
  else
  zwt(c) = zi(c,nlevbed)
  end if
-  if (rsub_top_tot < 0.) then
-  rsub_top(c) = rsub_top(c) + rsub_top_tot / dtime
+  if (rsub_top_tot < 0.) then
+  rsub_top(c) = rsub_top(c) + rsub_top_tot / dtime
  rsub_top_tot = 0.
  end if
  else
@@ -1488,7 +1521,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
 
  do fc = 1, num_hydrologyc
  c = filter_hydrologyc(fc)
-   nlevbed = nlev2bed(c)
+   nlevbed = nlev2bed(c)
  do j = nlevbed,2,-1
  xsi(c) = max(h2osoi_liq(c,j)-eff_porosity(c,j)*dzmm(c,j),0._r8)
  if (use_vsfm) then
@@ -1536,8 +1569,8 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
 
  do fc = 1, num_hydrologyc
  c = filter_hydrologyc(fc)
-   nlevbed = nlev2bed(c)
-   do j = 1, nlevbed-1
+   nlevbed = nlev2bed(c)
+   do j = 1, nlevbed-1
  if (h2osoi_liq(c,j) < watmin) then
  xs(c) = watmin - h2osoi_liq(c,j)
  ! deepen water table if water is passed from below zwt layer
@@ -1555,8 +1588,8 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
  ! Get water for bottom layer from layers above if possible
  do fc = 1, num_hydrologyc
  c = filter_hydrologyc(fc)
-   nlevbed = nlev2bed(c)
-   j = nlevbed
+   nlevbed = nlev2bed(c)
+   j = nlevbed
  if (h2osoi_liq(c,j) < watmin) then
  xs(c) = watmin-h2osoi_liq(c,j)
  searchforwater: do i = nlevbed-1, 1, -1