diff --git a/components/elm/src/biogeophys/ActiveLayerMod.F90 b/components/elm/src/biogeophys/ActiveLayerMod.F90 index 17480014398b..f88292864dc2 100644 --- a/components/elm/src/biogeophys/ActiveLayerMod.F90 +++ b/components/elm/src/biogeophys/ActiveLayerMod.F90 @@ -49,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, dzsoi + use elm_varcon , only : zsoi, dzsoi, zisoi ! ! !ARGUMENTS: integer , intent(in) :: num_soilc ! number of soil columns in filter @@ -58,17 +58,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 - real(r8), dimension(nlevgrnd) :: melt_profile ! profile of melted excess ice + 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 !----------------------------------------------------------------------- ! RF NOTE: use of 1989 ALT in these parameterizations is somewhat of a placeholder used to compare against @@ -87,11 +87,12 @@ subroutine alt_calc(num_soilc, filter_soilc, & 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: [real(r8) (:,:) ] depth variable excess ice content in soil column (-) + 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) + 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 @@ -184,15 +185,49 @@ subroutine alt_calc(num_soilc, filter_soilc, & ! update subsidence based on change in ALT ! melt_profile stores the amount of excess_ice ! melted in this timestep. - do j = 1,nlevgrnd - if (j < k_frz) then - melt_profile(j) = 0.0_r8 + ! 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 + !write(iulog,*) "processing level j,",j,k_frz,excess_ice(c,j) + 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) + ! DEBUG: + write(iulog,*) "melt and excess ice here are for j and j+1", melt_profile(j), & + melt_profile(j+1), excess_ice(c,j), & + excess_ice(c,j+1),(alt(c)-zisoi(j-1))/dzsoi(j),c,j+1,alt(c),zisoi(j-1),dzsoi(j) + else + melt_profile(j) = 0._r8 ! no melt + end if else if (j .eq. k_frz) then - melt_profile(j) = excess_ice(c,j) + ((z2-alt(c))/(z2-z1))*excess_ice(c,j+1) ! TODO: check indices here!! - ! remove melted excess ice: - excess_ice(c,j) = 0._r8 - excess_ice(c,j+1) = excess_ice(c,j+1)*(1._r8 - min(1._r8,(z2-alt(c))/(z2-z1))) - else + 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 @@ -204,6 +239,10 @@ subroutine alt_calc(num_soilc, filter_soilc, & ! 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) + if (sum(melt_profile) .lt. 0_r8) then + write(iulog,*) "subsidence(c) is",subsidence(c),sum(melt_profile),c,k_frz + write(iulog,*) "meltprofile is:", melt_profile + end if end if ! limit subsidence to 0.4 m diff --git a/components/elm/src/biogeophys/SoilHydrologyMod.F90 b/components/elm/src/biogeophys/SoilHydrologyMod.F90 index de8aaf8ebade..548117c7faa3 100644 --- a/components/elm/src/biogeophys/SoilHydrologyMod.F90 +++ b/components/elm/src/biogeophys/SoilHydrologyMod.F90 @@ -530,7 +530,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f 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(iwp_subsidence(c), 0.4) !fix this variable when available to pull from alt calculations + phi_eff = min(iwp_subsidence(c), 0.4_r8) !fix this variable when available to pull from alt calculations swc = h2osfc(c)/1000_r8 ! convert to m if (swc >= vdep) then diff --git a/components/elm/src/data_types/ColumnDataType.F90 b/components/elm/src/data_types/ColumnDataType.F90 index eff24c80a5e0..144aff38515f 100644 --- a/components/elm/src/data_types/ColumnDataType.F90 +++ b/components/elm/src/data_types/ColumnDataType.F90 @@ -35,7 +35,7 @@ module ColumnDataType use soilorder_varcon, only : smax, ks_sorption use elm_time_manager, only : is_restart, get_nstep use elm_time_manager, only : is_first_step, get_step_size, is_first_restart_step - use landunit_varcon , only : istice, istwet, istsoil, istdlak, istcrop, istice_mec + use landunit_varcon , only : istice, istwet, istsoil, istdlak, istcrop, istice_mec, istlowcenpoly, isthighcenpoly use column_varcon , only : icol_road_perv, icol_road_imperv, icol_roof, icol_sunwall, icol_shadewall use histFileMod , only : hist_addfld1d, hist_addfld2d, no_snow_normal use histFileMod , only : hist_addfld_decomp @@ -110,7 +110,6 @@ module ColumnDataType real(r8), pointer :: h2osoi_liq (:,:) => null() ! liquid water (-nlevsno+1:nlevgrnd) (kg/m2) real(r8), pointer :: h2osoi_ice (:,:) => null() ! ice lens (-nlevsno+1:nlevgrnd) (kg/m2) real(r8), pointer :: h2osoi_vol (:,:) => null() ! volumetric soil water (0<=h2osoi_vol<=watsat) (1:nlevgrnd) (m3/m3) - real(r8), pointer :: excess_ice (:,:) => null() ! NGEE Arctic: excess ground ice in column (1:nlevgrnd) (0 to 1) real(r8), pointer :: h2osfc (:) => null() ! surface water (kg/m2) real(r8), pointer :: h2ocan (:) => null() ! canopy water integrated to column (kg/m2) real(r8), pointer :: total_plant_stored_h2o(:)=> null() ! total water in plants (kg/m2) @@ -170,11 +169,13 @@ module ColumnDataType real(r8), pointer :: vsfm_soilp_col_1d (:) => null() ! 1D soil liquid pressure from VSFM [Pa] real(r8), pointer :: h2orof (:) => null() ! floodplain inundation volume received from rof (mm) real(r8), pointer :: frac_h2orof (:) => null() ! floodplain inundation fraction received from rof (-) - ! polygonal tundra - 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) - real(r8), pointer :: iwp_subsidence(:) => null() ! ice wedge polygon ground subsidence (m) + ! polygonal tundra (NGEE Arctic IM1) + 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) + real(r8), pointer :: iwp_subsidence (:) => null() ! ice wedge polygon ground subsidence (m) + real(r8), pointer :: excess_ice (:,:) => null() ! excess ground ice in column (1:nlevgrnd) (0 to 1) + real(r8), pointer :: frac_melted (:,:) => null() ! fraction of layer that has ever thawed (for tracking excess ice removal) (0 to 1) contains procedure, public :: Init => col_ws_init @@ -1396,7 +1397,6 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ allocate(this%h2osoi_liq (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liq (:,:) = spval allocate(this%h2osoi_ice (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_ice (:,:) = spval allocate(this%h2osoi_vol (begc:endc, 1:nlevgrnd)) ; this%h2osoi_vol (:,:) = spval - allocate(this%excess_ice (begc:endc, 1:nlevgrnd)) ; this%excess_ice (:,:) = spval allocate(this%h2osfc (begc:endc)) ; this%h2osfc (:) = spval allocate(this%h2ocan (begc:endc)) ; this%h2ocan (:) = spval allocate(this%wslake_col (begc:endc)) ; this%wslake_col (:) = spval @@ -1454,11 +1454,15 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ allocate(this%vsfm_soilp_col_1d (ncells)) ; this%vsfm_soilp_col_1d (:) = spval allocate(this%h2orof (begc:endc)) ; this%h2orof (:) = spval allocate(this%frac_h2orof (begc:endc)) ; this%frac_h2orof (:) = spval - ! polygonal tundra/ice wedge polygons: - allocate(this%iwp_microrel (begc:endc)) ; this%iwp_microrel (:) = spval - allocate(this%iwp_exclvol (begc:endc)) ; this%iwp_exclvol (:) = spval - allocate(this%iwp_ddep (begc:endc)) ; this%iwp_ddep (:) = spval - allocate(this%iwp_subsidence (begc:endc)) ; this%iwp_subsidence(:) = spval + if (use_polygonal_tundra) then + ! polygonal tundra/ice wedge polygons: + allocate(this%iwp_microrel (begc:endc)) ; this%iwp_microrel (:) = spval + allocate(this%iwp_exclvol (begc:endc)) ; this%iwp_exclvol (:) = spval + allocate(this%iwp_ddep (begc:endc)) ; this%iwp_ddep (:) = spval + allocate(this%iwp_subsidence (begc:endc)) ; this%iwp_subsidence (:) = spval + allocate(this%frac_melted (begc:endc,1:nlevgrnd)) ; this%frac_melted (:,:) = spval + allocate(this%excess_ice (begc:endc,1:nlevgrnd)) ; this%excess_ice (:,:) = spval + end if !----------------------------------------------------------------------- ! initialize history fields for select members of col_ws @@ -1494,19 +1498,17 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ avgflag='A', long_name='soil ice (ice landunits only)', & ptr_col=this%h2osoi_ice, l2g_scale_type='ice') - ! RPF - polygonal tundra vars - this%excess_ice(begc:endc, :) = spval if (use_polygonal_tundra) then + this%iwp_subsidence(begc:endc) = spval + this%iwp_ddep(begc:endc) = spval + this%iwp_exclvol(begc:endc) = spval + this%iwp_microrel(begc:endc) = spval + this%frac_melted(begc:endc,:) = spval + this%excess_ice(begc:endc,:) = spval + call hist_addfld2d (fname='EXCESS_ICE', units = '1', type2d='levgrnd', & avgflag='A', long_name='Excess ground ice (0 to 1)', & ptr_col=this%excess_ice, l2g_scale_type='veg') ! <- RPF: should this be natveg? - end if - - this%iwp_subsidence(begc:endc) = spval - this%iwp_ddep(begc:endc) = spval - this%iwp_exclvol(begc:endc) = spval - this%iwp_microrel(begc:endc) = spval - if (use_polygonal_tundra) then call hist_addfld1d (fname="SUBSIDENCE", units='m', avgflag='A', & long_name='ground subsidence (m)', ptr_col=this%iwp_subsidence) call hist_addfld1d (fname="DEPRESS_DEPTH", units='m', avgflag='A', & @@ -1516,6 +1518,9 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ ptr_col=this%iwp_exclvol) call hist_addfld1d (fname="MICROREL", units='m', avgflag='A', & long_name='microtopographic relief (m)', ptr_col=this%iwp_microrel) + call hist_addfld2d (fname="FRAC_MELTED", units='-', type2d='levgrnd', & + avgflag='A', long_name='fraction of layer that has melted (-)', & + ptr_col=this%frac_melted, l2g_scale_type='veg') endif !/polygonal tundra @@ -1679,6 +1684,7 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ this%h2orof(c) = 0._r8 this%frac_h2orof(c) = 0._r8 this%iwp_subsidence(c) = 0._r8 + this%frac_melted(c,:) = 0._r8 if (lun_pp%urbpoi(l)) then ! From Bonan 1996 (LSM technical note) @@ -1731,13 +1737,16 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ if (j > nlevbed) then this%h2osoi_vol(c,j) = 0.0_r8 else - if (use_fates .or. use_hydrstress) then + if (use_fates .or. use_hydrstress) then this%h2osoi_vol(c,j) = 0.70_r8*watsat_input(c,j) !0.15_r8 to avoid very dry conditions that cause errors in FATES else if (use_arctic_init) then - this%h2osoi_vol(c,j) = watsat_input(c,j) ! start saturated for arctic + this%h2osoi_vol(c,j) = watsat_input(c,j) ! start saturated for arctic else this%h2osoi_vol(c,j) = 0.15_r8 endif + if (use_polygonal_tundra) then + this%frac_melted(c,j) = 0._r8 + end if endif end do else if (lun_pp%urbpoi(l)) then @@ -1838,23 +1847,14 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ this%h2osoi_ice(c,j) = 0._r8 this%h2osoi_liq(c,j) = col_pp%dz(c,j)*denh2o*this%h2osoi_vol(c,j) endif - ! RPF 240713 - notes from Chuck: initialize all to 0.36_r8 - this%excess_ice(c,j) = 0.36_r8 - ! RPF - to do: a) apply to only polygonal ground - ! b) pull in ALT information to get starting point right. - ! for now, assuming no excess ice in top meter, 0.5 for 1-4 m, - ! 0.2 below 4 m. Also: need to check signs! - !if (col_pp%z(c,j) > -1._r8) then - ! this%excess_ice(c,j) = 0._r8 - !else if (col_pp%z(c,j) > -4._r8 .and. col_pp%z(c,j) < -1._r8) then - ! this%excess_ice(c,j) = 0.5_r8 - !else - ! this%excess_ice(c,j) = 0.2_r8 - !end if end do this%h2osoi_liq_old(c,:) = this%h2osoi_liq(c,:) this%h2osoi_ice_old(c,:) = this%h2osoi_ice(c,:) + if (use_polygonal_tundra) then + ! RPF 240713 - notes from Chuck Abolt: initialize all to 0.36_r8 + this%excess_ice(c,:) = 0.36_r8 + end if end do end subroutine col_ws_init @@ -1929,6 +1929,10 @@ subroutine col_ws_restart(this, bounds, ncid, flag, watsat_input) dim1name='column', & long_name='microtopographic relief', units='m', & interpinic_flag='interp', readvar=readvar, data=this%iwp_microrel) + call restartvar(ncid=ncid, flag=flag, varname='FRAC_MELTED', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='fraction of layer that has ever melted', units='-', & + interpinic_flag='interp', readvar=readvar, data=this%frac_melted) end if call restartvar(ncid=ncid, flag=flag, varname='SOILP', xtype=ncd_double, &