Merge branch 'cwhicker/elm/smb_h0_fields' (PR #6682)

Introduce land ice surface mass balance fields in default runs without MECs/GLC

Introduced new diagnostic Surface Mass Balance fields, they are
calculated as if there is a coupled land ice component but are only
diagnostic - it does not modify the energy/water budgets in default
or coupled GLC simulations.

[BFB]

components/elm/src/biogeophys/HydrologyDrainageMod.F90

@@ -47,7 +47,7 @@ subroutine HydrologyDrainage(bounds,    &
     !
     ! !USES:
       !$acc routine seq
-    use landunit_varcon  , only : istice, istwet, istsoil, istice_mec, istcrop
+    use landunit_varcon  , only : istice, istwet, istsoil, istice_mec, istcrop, istice
     use column_varcon    , only : icol_roof, icol_road_imperv, icol_road_perv, icol_sunwall, icol_shadewall
     use elm_varcon       , only : denh2o, denice, secspday
     use elm_varctl       , only : glc_snow_persistence_max_days, use_vichydro, use_betr
@@ -120,7 +120,9 @@ subroutine HydrologyDrainage(bounds,    &
          qflx_runoff_r          => col_wf%qflx_runoff_r           , & ! Output: [real(r8) (:)   ]  Rural total runoff (qflx_drain+qflx_surf+qflx_qrgwl) (mm H2O /s)
          qflx_snwcp_ice         => col_wf%qflx_snwcp_ice          , & ! Output: [real(r8) (:)   ]  excess snowfall due to snow capping (mm H2O /s) [+]`
          qflx_glcice            => col_wf%qflx_glcice             , & ! Output: [real(r8) (:)   ]  flux of new glacier ice (mm H2O /s)
-         qflx_glcice_frz        => col_wf%qflx_glcice_frz           & ! Output: [real(r8) (:)   ]  ice growth (positive definite) (mm H2O/s)
+         qflx_glcice_frz        => col_wf%qflx_glcice_frz         , & ! Output: [real(r8) (:)   ]  ice growth (positive definite) (mm H2O/s)
+         qflx_glcice_diag       => col_wf%qflx_glcice_diag        , & ! Output: [real(r8) (:)   ]  flux of new glacier ice (mm H2O/s) - diagnostic, no MECs or GLC
+         qflx_glcice_frz_diag   => col_wf%qflx_glcice_frz_diag      & ! Output: [real(r8) (:)   ]  ice growth (positive definite) (mm H2O/s)) - diagnostic, no MECs or GLC
          )
 
       ! Determine time step and step size
@@ -215,17 +217,23 @@ subroutine HydrologyDrainage(bounds,    &
 
       do c = bounds%begc,bounds%endc
          qflx_glcice_frz(c) = 0._r8
+         qflx_glcice_frz_diag(c) = 0._r8
+
+         if (lun_pp%itype(l)==istice .and. qflx_snwcp_ice(c) > 0.0_r8) then
+               qflx_glcice_frz_diag(c) = qflx_snwcp_ice(c)
+               qflx_glcice_diag(c) = qflx_glcice_diag(c) + qflx_glcice_frz_diag(c)
+         endif
       end do
       do fc = 1,num_do_smb_c
          c = filter_do_smb_c(fc)
          l = col_pp%landunit(c)
          g = col_pp%gridcell(c)
          ! In the following, we convert glc_snow_persistence_max_days to r8 to avoid overflow
          if ( (snow_persistence(c) >= (real(glc_snow_persistence_max_days, r8) * secspday)) &
-              .or. lun_pp%itype(l) == istice_mec) then
-            qflx_glcice_frz(c) = qflx_snwcp_ice(c)
-            qflx_glcice(c) = qflx_glcice(c) + qflx_glcice_frz(c)
-            if (glc_dyn_runoff_routing(g)) qflx_snwcp_ice(c) = 0._r8
+              .or. lun_pp%itype(l) == istice_mec ) then
+           qflx_glcice_frz(c) = qflx_snwcp_ice(c)
+           qflx_glcice(c) = qflx_glcice(c) + qflx_glcice_frz(c)
+           if (glc_dyn_runoff_routing(g)) qflx_snwcp_ice(c) = 0._r8
          end if
       end do
 

components/elm/src/biogeophys/SnowHydrologyMod.F90

@@ -670,7 +670,7 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, &
                       if (bi > dm) ddz1 = ddz1*exp(-46.0e-3_r8*(bi-dm))
                    else
                       ddz1_fresh = (-grav * (burden(c) + wx/2._r8)) / &
-                                   (0.007_r8 * bi**(4.75_r8 + td/40._r8))
+                                   (0.007_r8 * min(max(bi,dm),denice)**(4.75_r8 + min(td,0._r8)/40._r8))
                       snw_ssa = 3.e6_r8 / (denice * snw_rds(c,j))
                       if (snw_ssa < 50._r8) then
                           ddz1_fresh = ddz1_fresh * exp(-46.e-2_r8 * (50._r8 - snw_ssa))

components/elm/src/biogeophys/SoilTemperatureMod.F90

@@ -1316,7 +1316,7 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
     use elm_varctl       , only : iulog
     use elm_varcon       , only : tfrz, hfus, grav
     use column_varcon    , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv
-    use landunit_varcon  , only : istsoil, istcrop, istice_mec
+    use landunit_varcon  , only : istsoil, istcrop, istice_mec,istice
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds
@@ -1369,6 +1369,8 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
          qflx_snofrz_col  =>    col_wf%qflx_snofrz      , & ! Output: [real(r8) (:)   ] column-integrated snow freezing rate (positive definite) [kg m-2 s-1]
          qflx_glcice      =>    col_wf%qflx_glcice      , & ! Output: [real(r8) (:)   ] flux of new glacier ice (mm H2O/s) [+ = ice grows]
          qflx_glcice_melt =>    col_wf%qflx_glcice_melt , & ! Output: [real(r8) (:)   ] ice melt (positive definite) (mm H2O/s)
+         qflx_glcice_diag      =>    col_wf%qflx_glcice_diag      , & ! Output: [real(r8) (:)   ] flux of new glacier ice (mm H2O/s) [+ = ice grows]
+         qflx_glcice_melt_diag =>    col_wf%qflx_glcice_melt_diag , & ! Output: [real(r8) (:)   ] ice melt (positive definite) (mm H2O/s)
          qflx_snomelt     =>    col_wf%qflx_snomelt     , & ! Output: [real(r8) (:)   ] snow melt (mm H2O /s)
 
          eflx_snomelt     =>    col_ef%eflx_snomelt    , & ! Output: [real(r8) (:)   ] snow melt heat flux (W/m**2)
@@ -1393,6 +1395,7 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
          qflx_snofrz_lyr(c,-nlevsno+1:0) = 0._r8
          qflx_snofrz_col(c) = 0._r8
          qflx_glcice_melt(c) = 0._r8
+         qflx_glcice_melt_diag(c) = 0._r8
          qflx_snow_melt(c) = 0._r8
       end do
 
@@ -1643,8 +1646,8 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
             ! as computed in HydrologyDrainageMod.F90.
 
             l = col_pp%landunit(c)
-            if (lun_pp%itype(l)==istice_mec) then
 
+            if ( lun_pp%itype(l)==istice_mec) then
                if (j>=1 .and. h2osoi_liq(c,j) > 0._r8) then   ! ice layer with meltwater
                   ! melting corresponds to a negative ice flux
                   qflx_glcice_melt(c) = qflx_glcice_melt(c) + h2osoi_liq(c,j)/dtime
@@ -1656,6 +1659,16 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
 
                endif  ! liquid water is present
             endif     ! istice_mec
+            ! for diagnostic QICE SMB output only - 
+            ! these are to calculate SMB even without MECs 
+            if ( lun_pp%itype(l)==istice) then
+               if (j>=1 .and. h2osoi_liq(c,j) > 0._r8) then   ! ice layer with meltwater
+                  ! melting corresponds to a negative ice flux
+                  qflx_glcice_melt_diag(c) = qflx_glcice_melt_diag(c) + h2osoi_liq(c,j)/dtime
+                  qflx_glcice_diag(c) = qflx_glcice_diag(c) - h2osoi_liq(c,j)/dtime
+               endif  ! liquid water is present
+            endif     ! istice_mec
+
 
          end do   ! end of column-loop
       enddo   ! end of level-loop

components/elm/src/data_types/ColumnDataType.F90

@@ -502,6 +502,9 @@ module ColumnDataType
     real(r8), pointer :: qflx_glcice          (:)   => null() ! net flux of new glacial ice (growth - melt) (mm H2O/s), passed to GLC
     real(r8), pointer :: qflx_glcice_frz      (:)   => null() ! ice growth (positive definite) (mm H2O/s)
     real(r8), pointer :: qflx_glcice_melt     (:)   => null() ! ice melt (positive definite) (mm H2O/s)
+    real(r8), pointer :: qflx_glcice_diag     (:)   => null() ! net flux of new glacial ice (growth - melt) (mm H2O/s), passed to GLC
+    real(r8), pointer :: qflx_glcice_frz_diag (:)   => null() ! ice growth (positive definite) (mm H2O/s)
+    real(r8), pointer :: qflx_glcice_melt_diag(:)   => null() ! ice melt (positive definite) (mm H2O/s)
     real(r8), pointer :: qflx_drain_vr        (:,:) => null() ! liquid water lost as drainage (m /time step)
     real(r8), pointer :: qflx_h2osfc2topsoi   (:)   => null() ! liquid water coming from surface standing water top soil (mm H2O/s)
     real(r8), pointer :: qflx_snow2topsoi     (:)   => null() ! liquid water coming from residual snow to topsoil (mm H2O/s)
@@ -5725,6 +5728,9 @@ subroutine col_wf_init(this, begc, endc)
     allocate(this%qflx_glcice            (begc:endc))             ; this%qflx_glcice          (:)   = spval
     allocate(this%qflx_glcice_frz        (begc:endc))             ; this%qflx_glcice_frz      (:)   = spval
     allocate(this%qflx_glcice_melt       (begc:endc))             ; this%qflx_glcice_melt     (:)   = spval
+    allocate(this%qflx_glcice_diag       (begc:endc))             ; this%qflx_glcice_diag     (:)   = spval
+    allocate(this%qflx_glcice_frz_diag   (begc:endc))             ; this%qflx_glcice_frz_diag (:)   = spval
+    allocate(this%qflx_glcice_melt_diag  (begc:endc))             ; this%qflx_glcice_melt_diag(:)   = spval
     allocate(this%qflx_drain_vr          (begc:endc,1:nlevgrnd))  ; this%qflx_drain_vr        (:,:) = spval
     allocate(this%qflx_h2osfc2topsoi     (begc:endc))             ; this%qflx_h2osfc2topsoi   (:)   = spval
     allocate(this%qflx_snow2topsoi       (begc:endc))             ; this%qflx_snow2topsoi     (:)   = spval
@@ -5842,23 +5848,39 @@ subroutine col_wf_init(this, begc, endc)
      call hist_addfld1d (fname='QSNOFRZ', units='kg/m2/s', &
           avgflag='A', long_name='column-integrated snow freezing rate', &
            ptr_col=this%qflx_snofrz, set_lake=spval, c2l_scale_type='urbanf', default='inactive')
-
+   
     if (create_glacier_mec_landunit) then
-       this%qflx_glcice(begc:endc) = spval
-       call hist_addfld1d (fname='QICE',  units='mm/s',  &
-            avgflag='A', long_name='ice growth/melt', &
-            ptr_col=this%qflx_glcice, l2g_scale_type='ice')
-
-       this%qflx_glcice_frz(begc:endc) = spval
-       call hist_addfld1d (fname='QICE_FRZ',  units='mm/s',  &
-            avgflag='A', long_name='ice growth', &
-            ptr_col=this%qflx_glcice_frz, l2g_scale_type='ice')
-
-       this%qflx_glcice_melt(begc:endc) = spval
-       call hist_addfld1d (fname='QICE_MELT',  units='mm/s',  &
-            avgflag='A', long_name='ice melt', &
-            ptr_col=this%qflx_glcice_melt, l2g_scale_type='ice')
-    endif
+            this%qflx_glcice(begc:endc) = spval
+             call hist_addfld1d (fname='QICE',  units='mm/s',  &
+                  avgflag='A', long_name='ice growth/melt (with active GLC/MECs)', &
+                   ptr_col=this%qflx_glcice, l2g_scale_type='ice')
+
+            this%qflx_glcice_frz(begc:endc) = spval
+             call hist_addfld1d (fname='QICE_FRZ',  units='mm/s',  &
+                  avgflag='A', long_name='ice growth (with active GLC/MECs)', &
+                   ptr_col=this%qflx_glcice_frz, l2g_scale_type='ice')
+
+            this%qflx_glcice_melt(begc:endc) = spval
+             call hist_addfld1d (fname='QICE_MELT',  units='mm/s',  &
+                  avgflag='A', long_name='ice melt (with active GLC/MECs)', &
+                   ptr_col=this%qflx_glcice_melt, l2g_scale_type='ice')
+    else 
+             this%qflx_glcice_diag(begc:endc) = spval
+             call hist_addfld1d (fname='QICE',  units='mm/s',  &
+                  avgflag='A', long_name='diagnostic ice growth/melt (no active GLC/MECs)', &
+                   ptr_col=this%qflx_glcice_diag, l2g_scale_type='ice')
+
+            this%qflx_glcice_frz_diag(begc:endc) = spval
+             call hist_addfld1d (fname='QICE_FRZ',  units='mm/s',  &
+                  avgflag='A', long_name='diagnostic ice growth (no active GLC/MECs)', &
+                   ptr_col=this%qflx_glcice_frz_diag, l2g_scale_type='ice')
+
+            this%qflx_glcice_melt_diag(begc:endc) = spval
+             call hist_addfld1d (fname='QICE_MELT',  units='mm/s',  &
+                  avgflag='A', long_name='diagnostic ice melt (no active GLC/MECs)', &
+                   ptr_col=this%qflx_glcice_melt_diag, l2g_scale_type='ice')
+   end if
+
 
     ! As defined here, snow_sources - snow_sinks will equal the change in h2osno at any
     ! given time step but only if there is at least one snow layer (for all landunits

components/elm/src/main/elm_driver.F90

@@ -1604,6 +1604,8 @@ subroutine elm_drv_init(bounds, &
 
          qflx_glcice        => col_wf%qflx_glcice            , & ! Output: [real(r8) (:)   ]  flux of new glacier ice (mm H2O/s) [+ = ice grows]
 
+         qflx_glcice_diag   => col_wf%qflx_glcice_diag       , & ! Output: [real(r8) (:)   ]  flux of new glacier ice (mm H2O/s) [+ = ice grows]
+
          eflx_bot           => col_ef%eflx_bot              , & ! Output: [real(r8) (:)   ]  heat flux from beneath soil/ice column (W/m**2)
 
          cisun_z            => photosyns_vars%cisun_z_patch              , & ! Output: [real(r8) (:)   ]  intracellular sunlit leaf CO2 (Pa)
@@ -1637,6 +1639,7 @@ subroutine elm_drv_init(bounds, &
 
          ! Initialize qflx_glcice everywhere, to zero.
          qflx_glcice(c) = 0._r8
+         qflx_glcice_diag(c) = 0._r8
 
       end do