NGEE Arctic IM1 - Polygonal tundra

cime_config/tests.py

@@ -49,6 +49,9 @@
             "ERS.f19_g16.I20TRGSWCNPECACNTBC.elm-eca_f19_g16_I20TRGSWCNPECACNTBC",
             "ERS.f19_g16.I20TRGSWCNPRDCTCBC.elm-ctc_f19_g16_I20TRGSWCNPRDCTCBC",
             "ERS.r05_r05.ICNPRDCTCBC.elm-cbudget",
+            "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_default_I1850CNPRDCTCBC",
+            "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_codetest_I1850CNPRDCTCBC",
+            "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-polygonal_tundra"
             )
         },
 
@@ -94,8 +97,6 @@
             "SMS.r05_r05.IELM.elm-topounit",
             "ERS.ELM_USRDAT.I1850ELM.elm-usrdat",
             "ERS.r05_r05.IELM.elm-lnd_rof_2way",
-            "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_default_I1850CNPRDCTCBC",
-            "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_codetest_I1850CNPRDCTCBC",
             "ERS.r05_r05.IELM.elm-V2_ELM_MOSART_features",
             "ERS.ELM_USRDAT.IELM.elm-surface_water_dynamics"
             )

components/elm/bld/namelist_files/namelist_defaults.xml

@@ -2228,4 +2228,8 @@ this mask will have smb calculated over the entire global land surface
 <use_fan use_fates=".true." >.false.</use_fan>
 <use_fan fan_mode='disable'>.false.</use_fan>
 
+<!-- NGEE Arctic Options -->
+<use_polygonal_tundra>.false.</use_polygonal_tundra>;
+<use_arctic_init>.false.</use_arctic_init>;
+
 </namelist_defaults>

components/elm/bld/namelist_files/namelist_definition.xml

@@ -2160,6 +2160,16 @@ written only if do_budgets variable is .true.,
 default: 0
 </entry>
 
+<!-- ========================================================================================  -->
+<!-- NGEE Arctic Options                                                                       -->
+<!-- ========================================================================================  -->
+
+<entry id="use_arctic_init" type="logical" category="default_settings"
+	group="elm_inparm" value=".false.">
+Cold-start initialization conditions in tundra start saturated and frozen.
+<default>Default: .false.</default>
+</entry>
+
 <!-- ========================================================================================  -->
 <!-- Namelist options for soil erosion model                                                   -->
 <!-- ========================================================================================  -->
@@ -2222,4 +2232,14 @@ not be called.
 <default>Default: 0</default>
 </entry>
 
+<!-- ========================================================================================  -->
+<!-- NGEE Arctic Options                                                                       -->
+<!-- ========================================================================================  -->
+
+<entry id="use_polygonal_tundra" type="logical" category="default_settings"
+	group="elm_inparm" value=".false.">
+Use polygonal tundra parameterizations for ice wedge polygon microtopography and inindation fraction
+<default>Default: .false.</default>
+</entry>
+
 </namelist_definition>

components/elm/cime_config/testdefs/testmods_dirs/elm/polygonal_tundra/shell_commands

@@ -0,0 +1,9 @@
+./xmlchange LND_DOMAIN_FILE=domain.lnd.1x1pt_icycape_c240909.nc
+./xmlchange ATM_DOMAIN_FILE=domain.lnd.1x1pt_icycape_c240909.nc
+./xmlchange LND_DOMAIN_PATH="\$DIN_LOC_ROOT/share/domains/domain.clm"
+./xmlchange ATM_DOMAIN_PATH="\$DIN_LOC_ROOT/share/domains/domain.clm"
+./xmlchange DATM_MODE=CLMGSWP3v1
+./xmlchange DATM_CLMNCEP_YR_START=2000
+./xmlchange DATM_CLMNCEP_YR_END=2000
+./xmlchange NTASKS=1
+./xmlchange NTHRDS=1

components/elm/cime_config/testdefs/testmods_dirs/elm/polygonal_tundra/user_nl_elm

@@ -0,0 +1,3 @@
+use_polygonal_tundra = .true.
+use_arctic_init = .true.
+fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_1x1pt_IcyCape_polygons_c240909.nc'

components/elm/src/biogeophys/ActiveLayerMod.F90

@@ -7,12 +7,14 @@ module ActiveLayerMod
   ! !USES:
   use shr_kind_mod    , only : r8 => shr_kind_r8
   use shr_const_mod   , only : SHR_CONST_TKFRZ
-  use elm_varctl      , only : iulog
+  use elm_varctl      , only : iulog, spinup_state, use_polygonal_tundra
   use TemperatureType , only : temperature_type
   use CanopyStateType , only : canopystate_type
   use GridcellType    , only : grc_pp       
   use ColumnType      , only : col_pp
-  use ColumnDataType  , only : col_es  
+  use ColumnDataType  , only : col_es, col_ws
+  use LandunitType    , only : lun_pp
+  use landunit_varcon , only : ilowcenpoly, iflatcenpoly, ihighcenpoly
   !
   implicit none
   save
@@ -47,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
+    use elm_varcon       , only : zsoi, dzsoi, zisoi
     !
     ! !ARGUMENTS:
     integer                , intent(in)    :: num_soilc       ! number of soil columns in filter
@@ -56,27 +58,41 @@ 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 activel 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, orig_excess, old_mfrac  ! temporary variables
+    real(r8), dimension(nlevgrnd) :: melt_profile       ! profile of melted excess ice
     !-----------------------------------------------------------------------
 
-    associate(                                                                & 
-         t_soisno             =>    col_es%t_soisno        ,    & ! Input:   [real(r8) (:,:) ]  soil temperature (Kelvin)  (-nlevsno+1:nlevgrnd)                    
-
-         alt                  =>    canopystate_vars%alt_col             ,    & ! Output:  [real(r8) (:)   ]  current depth of thaw                                                 
-         altmax               =>    canopystate_vars%altmax_col          ,    & ! Output:  [real(r8) (:)   ]  maximum annual depth of thaw                                          
-         altmax_lastyear      =>    canopystate_vars%altmax_lastyear_col ,    & ! Output:  [real(r8) (:)   ]  prior year maximum annual depth of thaw                               
-         alt_indx             =>    canopystate_vars%alt_indx_col        ,    & ! Output:  [integer  (:)   ]  current depth of thaw                                                  
-         altmax_indx          =>    canopystate_vars%altmax_indx_col     ,    & ! Output:  [integer  (:)   ]  maximum annual depth of thaw                                           
-         altmax_lastyear_indx =>    canopystate_vars%altmax_lastyear_indx_col & ! Output:  [integer  (:)   ]  prior year maximum annual depth of thaw                                
+    ! RF NOTE: use of 1989 ALT in these parameterizations is somewhat of a placeholder used to compare against
+    ! ATS runs for NGEE Arctic Phase 3. It should ultimately be replaced by a more physically based threshold
+    ! later on.
+    associate(                                                                &
+         t_soisno             =>    col_es%t_soisno        ,    & ! Input:   [real(r8) (:,:) ]  soil temperature (Kelvin)  (-nlevsno+1:nlevgrnd)
+
+         alt                  =>    canopystate_vars%alt_col             ,      & ! Output:  [real(r8) (:)   ]  current depth of thaw
+         altmax               =>    canopystate_vars%altmax_col          ,      & ! Output:  [real(r8) (:)   ]  maximum annual depth of thaw
+         altmax_lastyear      =>    canopystate_vars%altmax_lastyear_col ,      & ! Output:  [real(r8) (:)   ]  prior year maximum annual depth of thaw
+         altmax_1989          =>    canopystate_vars%altmax_1989_col     ,      & ! Output:  [real(r8) (:)   ]  maximum ALT in 1989
+         altmax_ever          =>    canopystate_vars%altmax_ever_col     ,      & ! Output:  [real(r8) (:)   ]  maximum thaw depth since initialization
+         alt_indx             =>    canopystate_vars%alt_indx_col        ,      & ! Output:  [integer  (:)   ]  current depth of thaw
+         altmax_indx          =>    canopystate_vars%altmax_indx_col     ,      & ! Output:  [integer  (:)   ]  maximum annual depth of thaw
+         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/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)
+         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
@@ -144,14 +160,111 @@ subroutine alt_calc(num_soilc, filter_soilc, &
             endif
          endif
 
-
          ! if appropriate, update maximum annual active layer thickness
          if (alt(c) > altmax(c)) then
             altmax(c) = alt(c)
             altmax_indx(c) = alt_indx(c)
          endif
+         if (alt(c) > altmax_ever(c)) then
+            if (spinup_state .eq. 0) then
+                altmax_ever(c) = alt(c)
+                altmax_ever_indx(c) = alt_indx(c)
+            else !overwrite if in spinup
+                altmax_ever(c) = 0._r8
+                altmax_ever_indx(c) = 0
+            endif
+         endif
+
+         ! special loop for if year = 1989, see above note regarding
+         ! replacing with more physically based mechanism.
+         if (year .eq. 1989) then
+            altmax_1989(c) = altmax(c)
+            altmax_1989_indx(c) = altmax_indx(c)
+         endif
 
-      end do
+         if (use_polygonal_tundra) then
+           ! update subsidence based on change in ALT
+           ! melt_profile stores the amount of excess_ice
+           ! melted in this timestep.
+           ! 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
+              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)
+                else
+                  melt_profile(j) = 0._r8 ! no melt
+                end if
+              else if (j .eq. k_frz) then
+                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
+              end if
+              ! calculate subsidence at this layer:
+              melt_profile(j) = melt_profile(j) * dzsoi(j)
+           end do
+
+           ! 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)
+           end if
+
+           ! limit subsidence to 0.4 m
+           subsidence(c) = min(0.4_r8, subsidence(c))
+
+           ! update ice wedge polygon microtopographic parameters if in polygonal ground
+           !rpf - min/max logic may be redunant w/ subsidence limiter above
+           !rpf - TODO: many of these (particularly the invariant ones) should be
+           ! moved to data_types/ColumnDataType.F90 to avoid the highcenpoly elseif condition
+           ! each time step.
+           if (lun_pp%ispolygon(col_pp%landunit(c))) then
+             if (lun_pp%polygontype(col_pp%landunit(c)) .eq. ilowcenpoly) then
+               rmax(c) = 0.4_r8
+               vexc(c) = 0.2_r8
+               ddep(c) = max(0.05_r8, 0.15_r8 - 0.25_r8*subsidence(c))
+             elseif (lun_pp%polygontype(col_pp%landunit(c)) .eq. iflatcenpoly) then
+               rmax(c) = min(0.4_r8, 0.1_r8 + 0.75_r8*subsidence(c))
+               vexc(c) = min(0.2_r8, 0.05_r8 + 0.375_r8*subsidence(c))
+               ddep(c) = min(0.05_r8, 0.01_r8 + 0.1_r8*subsidence(c))
+             elseif (lun_pp%polygontype(col_pp%landunit(c)) .eq. ihighcenpoly) then
+               rmax(c) = 0.4_r8
+               vexc(c) = 0.2_r8
+               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.
+             endif
+           endif
+         endif
+       end do
 
     end associate