Merge branch 'cwhicker/elm/use_extrasnowlayers_bugfix_snowdiags' (PR #6859)

Minor mods to fix segmentation fault errors in use_extrasnowlayers errors
that occur when the snow depth is very shallow.

[B4B]

Conflicts:
	components/elm/src/biogeophys/BalanceCheckMod.F90
	components/elm/src/biogeophys/HydrologyDrainageMod.F90
	components/elm/src/biogeophys/SoilTemperatureMod.F90
	components/elm/src/data_types/ColumnDataType.F90

Author: bishtgautam

components/elm/src/biogeophys/BalanceCheckMod.F90

@@ -164,10 +164,10 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
      ! error = abs(precipitation - change of water storage - evaporation - runoff)
      !
      ! !USES:
-     use elm_varcon        , only : spval
+     use elm_varcon        , only : spval, h2osno_max
      use column_varcon     , only : icol_roof, icol_sunwall, icol_shadewall
      use column_varcon     , only : icol_road_perv, icol_road_imperv
-     use landunit_varcon   , only : istice_mec, istdlak, istsoil,istcrop,istwet
+     use landunit_varcon   , only : istice_mec, istice, istdlak, istsoil,istcrop,istwet
      use elm_varctl        , only : create_glacier_mec_landunit, use_IM2_hillslope_hydrology
      use elm_initializeMod , only : surfalb_vars  
      use CanopyStateType   , only : canopystate_type
@@ -221,6 +221,7 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
           qflx_surf_irrig_col        =>    col_wf%qflx_surf_irrig         , & ! Input:  [real(r8) (:)   ]  real surface irrigation flux (mm H2O /s)
           qflx_over_supply_col       =>    col_wf%qflx_over_supply        , & ! Input:  [real(r8) (:)   ]  over supply irrigation flux (mm H2O /s)
           qflx_snwcp_ice             =>    col_wf%qflx_snwcp_ice          , & ! Input:  [real(r8) (:)   ]  excess snowfall due to snow capping (mm H2O /s) [+]`
+          qflx_ice_runoff_xs         =>    col_wf%qflx_ice_runoff_xs      , & ! Input:  [real(r8) (:)   ] 
           qflx_evap_tot              =>    col_wf%qflx_evap_tot           , & ! Input:  [real(r8) (:)   ]  qflx_evap_soi + qflx_evap_can + qflx_tran_veg
           qflx_dew_snow              =>    col_wf%qflx_dew_snow           , & ! Input:  [real(r8) (:)   ]  surface dew added to snow pack (mm H2O /s) [+]
           qflx_sub_snow              =>    col_wf%qflx_sub_snow           , & ! Input:  [real(r8) (:)   ]  sublimation rate from snow pack (mm H2O /s) [+]
@@ -334,7 +335,7 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
                   - (forc_rain_col(c) + forc_snow_col(c)  + qflx_floodc(c) + qflx_from_uphill(c) &
                   + qflx_surf_irrig_col(c) + qflx_over_supply_col(c) &
                   - qflx_evap_tot(c) - qflx_surf(c)  - qflx_h2osfc_surf(c) - qflx_to_downhill(c) &
-                  - qflx_qrgwl(c) - qflx_drain(c) - qflx_drain_perched(c) - qflx_snwcp_ice(c) &
+                  - qflx_qrgwl(c) - qflx_drain(c) - qflx_drain_perched(c) - qflx_snwcp_ice(c) - qflx_ice_runoff_xs(c) &
                   - qflx_lateral(c) + qflx_h2orof_drain(c)) * dtime
              dwb(c) = (endwb(c)-begwb(c))/dtime
 
@@ -406,6 +407,7 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
              write(iulog,*)'qflx_lateral               = ',qflx_lateral(indexc)
              write(iulog,*)'total_plant_stored_h2o_col = ',total_plant_stored_h2o_col(indexc)
              write(iulog,*)'qflx_h2orof_drain          = ',qflx_h2orof_drain(indexc)
+             write(iulog,*)'qflx_ice_runoff_xs          = ',qflx_ice_runoff_xs(indexc)
              write(iulog,*)'elm model is stopping'
              call endrun(decomp_index=indexc, elmlevel=namec, msg=errmsg(__FILE__, __LINE__))
 
@@ -436,6 +438,7 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
              write(iulog,*)'qflx_lateral               = ',qflx_lateral(indexc)
              write(iulog,*)'total_plant_stored_h2o_col = ',total_plant_stored_h2o_col(indexc)
              write(iulog,*)'qflx_h2orof_drain          = ',qflx_h2orof_drain(indexc)
+             write(iulog,*)'qflx_ice_runoff_xs          = ',qflx_ice_runoff_xs(indexc)
              write(iulog,*)'elm model is stopping'
              call endrun(decomp_index=indexc, elmlevel=namec, msg=errmsg(__FILE__, __LINE__))
           end if
@@ -559,15 +562,23 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
              write(iulog,*)'h2osno_old       = ',h2osno_old(indexc)
              write(iulog,*)'snow_sources     = ',snow_sources(indexc)
              write(iulog,*)'snow_sinks       = ',snow_sinks(indexc)
-             write(iulog,*)'qflx_prec_grnd   = ',qflx_prec_grnd(indexc)*dtime
-             write(iulog,*)'qflx_sub_snow    = ',qflx_sub_snow(indexc)*dtime
-             write(iulog,*)'qflx_evap_grnd   = ',qflx_evap_grnd(indexc)*dtime
-             write(iulog,*)'qflx_top_soil    = ',qflx_top_soil(indexc)*dtime
-             write(iulog,*)'qflx_dew_snow    = ',qflx_dew_snow(indexc)*dtime
-             write(iulog,*)'qflx_dew_grnd    = ',qflx_dew_grnd(indexc)*dtime
-             write(iulog,*)'qflx_snwcp_ice   = ',qflx_snwcp_ice(indexc)*dtime
-             write(iulog,*)'qflx_snwcp_liq   = ',qflx_snwcp_liq(indexc)*dtime
-             write(iulog,*)'qflx_sl_top_soil = ',qflx_sl_top_soil(indexc)*dtime
+             write(iulog,*)'qflx_prec_grnd   = ',qflx_prec_grnd(indexc)
+             write(iulog,*)'qflx_sub_snow    = ',qflx_sub_snow(indexc)
+             write(iulog,*)'qflx_evap_grnd   = ',qflx_evap_grnd(indexc)
+             write(iulog,*)'qflx_top_soil    = ',qflx_top_soil(indexc)
+             write(iulog,*)'qflx_dew_snow    = ',qflx_dew_snow(indexc)
+             write(iulog,*)'qflx_dew_grnd    = ',qflx_dew_grnd(indexc)
+             write(iulog,*)'qflx_snwcp_ice   = ',qflx_snwcp_ice(indexc)
+             write(iulog,*)'qflx_snwcp_liq   = ',qflx_snwcp_liq(indexc)
+             write(iulog,*)'qflx_sl_top_soil = ',qflx_sl_top_soil(indexc)
+             write(iulog,*)'qflx_snow_melt   = ',qflx_snow_melt(indexc)
+             write(iulog,*)'frac_sno_eff     = ',frac_sno_eff(indexc)
+             write(iulog,*)'qflx_rain_grnd_col= ',qflx_rain_grnd_col(indexc)
+             write(iulog,*)'qflx_snow_grnd_col= ',qflx_snow_grnd_col(indexc)
+             write(iulog,*)'qflx_h2osfc_to_ice= ',qflx_h2osfc_to_ice(indexc)
+             write(iulog,*)'qflx_snow_h2osfc= ',qflx_snow_h2osfc(indexc)*dtime
+             write(iulog,*)'(glc_dyn_runoff_routing(g))',(glc_dyn_runoff_routing(g))
+
              if (create_glacier_mec_landunit) then
                 write(iulog,*)'qflx_glcice_frz  = ',qflx_glcice_frz(indexc)*dtime
              end if

components/elm/src/biogeophys/HydrologyDrainageMod.F90

@@ -7,7 +7,7 @@ module HydrologyDrainageMod
   use shr_kind_mod      , only : r8 => shr_kind_r8
   use shr_log_mod       , only : errMsg => shr_log_errMsg
   use decompMod         , only : bounds_type
-  use elm_varctl        , only : iulog, use_vichydro
+  use elm_varctl        , only : iulog, use_vichydro, use_firn_percolation_and_compaction
   use elm_varcon        , only : e_ice, denh2o, denice, rpi, spval
   use atm2lndType       , only : atm2lnd_type
   use glc2lndMod        , only : glc2lnd_type
@@ -238,6 +238,12 @@ subroutine HydrologyDrainage(bounds,    &
            qflx_glcice(c) = qflx_glcice(c) + qflx_glcice_frz(c)
            if (glc_dyn_runoff_routing(g)) qflx_snwcp_ice(c) = 0._r8
          end if
+
+         !if (lun_pp%itype(l)==istice) 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
 
       ! Determine wetland and land ice hydrology (must be placed here
@@ -270,8 +276,9 @@ subroutine HydrologyDrainage(bounds,    &
             ! glc_dyn_runoff_routing = true: in this case, melting ice runs off, and excess
             ! snow is sent to CISM, where it is converted to ice. These corrections are
             ! done here:
-
             if (glc_dyn_runoff_routing(g) .and. lun_pp%itype(l)==istice_mec) then
+            ! this allows GLC melt to runoff to qflx_qrgwl! 
+
                ! If glc_dyn_runoff_routing=T, add meltwater from istice_mec ice columns to the runoff.
                !    Note: The meltwater contribution is computed in PhaseChanges (part of Biogeophysics2)
                qflx_qrgwl(c) = qflx_qrgwl(c) + qflx_glcice_melt(c)

components/elm/src/biogeophys/LakeTemperatureMod.F90

@@ -1291,6 +1291,7 @@ subroutine PhaseChange_Lake (bounds, num_lakec, filter_lakec, cv, cv_lake, lhabs
           qflx_snofrz_lyr => col_wf%qflx_snofrz_lyr , & ! Input:  [real(r8)  (:,:) ] snow freezing rate (positive definite) (col,lyr) [kg m-2 s-1]
           qflx_snow_melt  => col_wf%qflx_snow_melt  , & ! Output: [real(r8)  (:)   ] net snow melt
           qflx_snomelt    => col_wf%qflx_snomelt    , & ! Output: [real(r8)  (:)   ] snow melt (mm H2O /s)
+          qflx_snomelt_lyr=> col_wf%qflx_snomelt_lyr, & ! Output: [real(r8)  (:)   ] snow melt per layer (mm H2O /s)
           qflx_snofrz_col => col_wf%qflx_snofrz     , & ! Output: [real(r8)  (:)   ] column-integrated snow freezing rate (kg m-2 s-1) [+]
 
           t_soisno        => col_es%t_soisno      , & ! Input:  [real(r8)  (:,:) ] soil temperature (Kelvin)
@@ -1306,10 +1307,11 @@ subroutine PhaseChange_Lake (bounds, num_lakec, filter_lakec, cv, cv_lake, lhabs
        do fc = 1,num_lakec
           c = filter_lakec(fc)
 
-          qflx_snomelt(c)   = 0._r8
-          eflx_snomelt(c)   = 0._r8
-          lhabs(c)          = 0._r8
-          qflx_snow_melt(c) = 0._r8
+          qflx_snomelt(c)      = 0._r8
+          qflx_snomelt_lyr(c,:)= 0._r8
+          eflx_snomelt(c)      = 0._r8
+          lhabs(c)             = 0._r8
+          qflx_snow_melt(c)    = 0._r8
        end do
 
        do j = -nlevsno+1,0
@@ -1399,6 +1401,7 @@ subroutine PhaseChange_Lake (bounds, num_lakec, filter_lakec, cv, cv_lake, lhabs
                    if (j <= 0) then !snow
                       imelt(c,j) = 1
                       qflx_snomelt(c) = qflx_snomelt(c) + melt/dtime
+                      qflx_snomelt_lyr(c,j) = melt/dtime
                    end if
                 else if (t_soisno(c,j) < tfrz .and. h2osoi_liq(c,j) > 0._r8) then !freezing
                    dophasechangeflag = .true.

components/elm/src/biogeophys/SnowHydrologyMod.F90

@@ -2273,6 +2273,7 @@ subroutine SnowCapping(bounds, num_nolakec, filter_initc, num_snowc, filter_snow
      associate( &
          qflx_snwcp_ice     => col_wf%qflx_snwcp_ice               , & ! Output: [real(r8) (:)   ]  excess solid h2o due to snow capping (outgoing) (mm H2O /s) [+]
          qflx_snwcp_liq     => col_wf%qflx_snwcp_liq               , & ! Output: [real(r8) (:)   ]  excess liquid h2o due to snow capping (outgoing) (mm H2O /s) [+]
+         qflx_ice_runoff_xs =>    col_wf%qflx_ice_runoff_xs        , & ! Input:  [real(r8) (:)   ] 
          h2osoi_ice         => col_ws%h2osoi_ice                   , & ! In/Out: [real(r8) (:,:) ] ice lens (kg/m2)                       
          h2osoi_liq         => col_ws%h2osoi_liq                   , & ! In/Out: [real(r8) (:,:) ] liquid water (kg/m2)                   
          h2osno             => col_ws%h2osno                       , & ! Input:  [real(r8) (:)   ] snow water (mm H2O)
@@ -2294,6 +2295,7 @@ subroutine SnowCapping(bounds, num_nolakec, filter_initc, num_snowc, filter_snow
         c = filter_initc(fc)
         qflx_snwcp_ice(c) = 0.0_r8
         qflx_snwcp_liq(c) = 0.0_r8
+        qflx_ice_runoff_xs(c) = 0.0_r8
      end do
 
      loop_columns: do fc = 1, num_snowc

components/elm/src/biogeophys/SoilHydrologyMod.F90

@@ -1016,9 +1016,11 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
      use elm_varpar       , only : nlevsoi, nlevgrnd, nlayer, nlayert
      use elm_varcon       , only : pondmx, tfrz, watmin,rpi, secspday, nlvic
      use column_varcon    , only : icol_roof, icol_road_imperv, icol_road_perv
-     use elm_varctl       , only : use_vsfm, use_var_soil_thick
+     use elm_varctl       , only : use_vsfm, use_var_soil_thick, use_firn_percolation_and_compaction
      use SoilWaterMovementMod, only : zengdecker_2009_with_var_soil_thick
      use pftvarcon        , only : rsub_top_globalmax
+     use LandunitType      , only : lun_pp
+     use landunit_varcon  , only : istice_mec, istice
      !
      ! !ARGUMENTS:
      type(bounds_type)        , intent(in)    :: bounds
@@ -1121,6 +1123,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
 
           qflx_snwcp_liq     =>    col_wf%qflx_snwcp_liq     , & ! Output: [real(r8) (:)   ] excess rainfall due to snow capping (mm H2O /s) [+]
           qflx_snwcp_ice     =>    col_wf%qflx_snwcp_ice     , & ! Output: [real(r8) (:)   ] excess snowfall due to snow capping (mm H2O /s) [+]
+          qflx_ice_runoff_xs =>    col_wf%qflx_ice_runoff_xs , & ! Output: solid runoff from excess ice in soil (mm H2O /s) [+]
           !qflx_dew_grnd      =>    col_wf%qflx_dew_grnd      , & ! Output: [real(r8) (:)   ] ground surface dew formation (mm H2O /s) [+]
           !qflx_dew_snow      =>    col_wf%qflx_dew_snow      , & ! Output: [real(r8) (:)   ] surface dew added to snow pack (mm H2O /s) [+]
           !qflx_sub_snow      =>    col_wf%qflx_sub_snow      , & ! Output: [real(r8) (:)   ] sublimation rate from snow pack (mm H2O /s) [+]
@@ -1559,7 +1562,11 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
           ! add in ice check
           xs1(c)          = max(max(h2osoi_ice(c,1),0._r8)-max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1))),0._r8)
           h2osoi_ice(c,1) = min(max(0._r8,pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1)), h2osoi_ice(c,1))
-          qflx_snwcp_ice(c) = qflx_snwcp_ice(c) + xs1(c) / dtime
+          if ( (lun_pp%itype(col_pp%landunit(c)) == istice .or. lun_pp%itype(col_pp%landunit(c)) == istice_mec) .or. (.not. use_firn_percolation_and_compaction)) then      
+                qflx_snwcp_ice(c) = qflx_snwcp_ice(c) + xs1(c) / dtime
+          else
+                qflx_ice_runoff_xs(c) = qflx_ice_runoff_xs(c) + xs1(c) / dtime
+          endif
        end do
 
        ! Limit h2osoi_liq to be greater than or equal to watmin.

components/elm/src/biogeophys/SoilTemperatureMod.F90

@@ -1372,6 +1372,7 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
          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)
+         qflx_snomelt_lyr     =>    col_wf%qflx_snomelt_lyr     , & ! Output: [real(r8) (:)   ] snow melt (mm H2O /s)
 
          eflx_snomelt     =>    col_ef%eflx_snomelt    , & ! Output: [real(r8) (:)   ] snow melt heat flux (W/m**2)
          eflx_snomelt_r   =>    col_ef%eflx_snomelt_r  , & ! Output: [real(r8) (:)   ] rural snow melt heat flux (W/m**2)
@@ -1391,6 +1392,7 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
          l = col_pp%landunit(c)
 
          qflx_snomelt(c) = 0._r8
+         qflx_snomelt_lyr(c,-nlevsno+1:0) = 0._r8
          xmf(c) = 0._r8
          qflx_snofrz_lyr(c,-nlevsno+1:0) = 0._r8
          qflx_snofrz_col(c) = 0._r8
@@ -1626,7 +1628,7 @@ subroutine Phasechange_beta (bounds, num_nolakec, filter_nolakec, dhsdT, &
 
                      if (imelt(c,j) == 1 .AND. j < 1) then
                         qflx_snomelt(c) = qflx_snomelt(c) + max(0._r8,(wice0(c,j)-h2osoi_ice(c,j)))/dtime
-
+                        qflx_snomelt_lyr(c,j) = max(0._r8,(wice0(c,j)-h2osoi_ice(c,j)))/dtime 
 
                      endif
 
@@ -1646,7 +1648,6 @@ 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 (j>=1 .and. h2osoi_liq(c,j) > 0._r8) then   ! ice layer with meltwater
                   ! melting corresponds to a negative ice flux