Cabolt/im1

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
+     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) (:)   ] microtopography pdf sigma (m)
+          iwp_ddep     => col_pp%iwp_ddep     , & ! Input:  [real(r8) (:)   ] microtopography pdf sigma (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/data_types/ColumnType.F90

@@ -57,6 +57,9 @@ module ColumnType
      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)
 
      ! vertical levels
      integer , pointer :: snl          (:)   => null() ! number of snow layers

components/elm/src/data_types/LandunitType.F90

@@ -47,6 +47,8 @@ module LandunitType
      logical , pointer :: urbpoi       (:) => null() ! true=>urban point
      logical , pointer :: glcmecpoi    (:) => null() ! true=>glacier_mec point
      logical , pointer :: active       (:) => null() ! true=>do computations on this landunit
+     logical , pointer :: ispolygon    (:) => null() ! true=>is polygonal tundra
+     integer , pointer :: polygontype  (:) => null() ! ice wedge polygon initial condition (LCP, FCP, or HCP)
 
      ! urban properties
      real(r8), pointer :: canyon_hwr   (:) => null() ! urban landunit canyon height to width ratio (-)
@@ -93,6 +95,8 @@ subroutine lun_pp_Init(this, begl, endl)
     allocate(this%lakpoi       (begl:endl)); this%lakpoi    (:) = .false.
     allocate(this%urbpoi       (begl:endl)); this%urbpoi    (:) = .false.
     allocate(this%glcmecpoi    (begl:endl)); this%glcmecpoi (:) = .false.
+    allocate(this%ispolygon    (begl:endl)); this%ispolygon (:) = .false.
+    allocate(this%polygontype  (begl:endl)); this%polygontype(:) = ispval
 
     ! The following is initialized in routine setActive in module reweightMod
     allocate(this%active       (begl:endl))
@@ -129,6 +133,8 @@ subroutine lun_pp_clean(this)
     deallocate(this%lakpoi       )
     deallocate(this%urbpoi       )
     deallocate(this%glcmecpoi    )
+    deallocate(this%ispolygon    )
+    deallocate(this%polygontype  )
     deallocate(this%active       )
     deallocate(this%canyon_hwr   )
     deallocate(this%wtroad_perv  )

components/elm/src/main/landunit_varcon.F90

@@ -35,7 +35,20 @@ module landunit_varcon
 
   integer, parameter, public                   :: landunit_name_length = 40  ! max length of landunit names
   character(len=landunit_name_length), public  :: landunit_names(max_lunit)  ! name of each landunit type
+
+  ! land unit polygonal ground types
+  integer, parameter, public :: ilowcenpoly     = 1     ! low-centered polygons
+  integer, parameter, public :: iflatcenpoly    = 2	    ! flat-centered polygons
+  integer, parameter, public :: ihighcenpoly    = 3	    ! high-centered polygons
+
+  integer, parameter, public :: max_poly  = 3  !maximum value that lun_pp%polygontype can have
+                                        !(i.e., largest value in the above list)
 
+  integer, parameter, public                   :: polygon_name_length = 40  ! max length of landunit names
+  character(len=polygon_name_length), public   :: polygon_names(max_poly)  ! name of each landunit type
+
+
+
   ! parameters that depend on the above constants
 
   integer, parameter, public :: numurbl = isturb_MAX - isturb_MIN + 1   ! number of urban landunits
@@ -48,6 +61,7 @@ module landunit_varcon
   !
   ! !PRIVATE MEMBER FUNCTIONS:
   private :: set_landunit_names   ! set the landunit_names vector
+  private :: set_polygon_names    ! set the polygon_names vector
 !-----------------------------------------------------------------------
 
 contains
@@ -98,6 +112,30 @@ function landunit_is_special(ltype) result(is_special)
 
   end function landunit_is_special
 
+  subroutine set_polygon_names
+    !
+    ! !DESCRIPTION:
+    ! Set the polygon_names vector
+    !
+    ! !USES:
+    use shr_sys_mod, only : shr_sys_abort
+    !
+    character(len=*), parameter :: not_set = 'NOT_SET'
+    character(len=*), parameter :: subname = 'set_polygon_names'
+    !-----------------------------------------------------------------------
+
+    polygon_names(:) = not_set
+
+    polygon_names(ilowcenpoly) = 'low_centered_polygons'
+    polygon_names(iflatcenpoly) = 'flat_centered_polygons'
+    polygon_names(ihighcenpoly) = 'high_centered_polygons'
+
+    if (any(polygon_names == not_set)) then
+       call shr_sys_abort(trim(subname)//': Not all polygon names set')
+    end if
+
+  end subroutine set_polygon_names
+
   !-----------------------------------------------------------------------
   subroutine set_landunit_names
     !