eCLM-ParFlow: Couple effective porosity for frozen soil

src/clm5/biogeophys/SoilHydrologyType.F90

@@ -29,6 +29,9 @@ Module SoilHydrologyType
      real(r8), pointer :: qcharge_col       (:)     ! col aquifer recharge rate (mm/s) 
      real(r8), pointer :: fracice_col       (:,:)   ! col fractional impermeability (-)
      real(r8), pointer :: icefrac_col       (:,:)   ! col fraction of ice       
+#ifdef COUP_OAS_PFL
+     real(r8), pointer :: pfl_eff_porosity_col (:,:) ! col effective porosity sent to ParFlow = watsat - vol_ice (m3/m3), unfloored
+#endif
      real(r8), pointer :: fcov_col          (:)     ! col fractional impermeable area
      real(r8), pointer :: fsat_col          (:)     ! col fractional area with water table at surface
      real(r8), pointer :: h2osfc_thresh_col (:)     ! col level at which h2osfc "percolates"   (time constant)
@@ -130,6 +133,9 @@ subroutine InitAllocate(this, bounds)
     allocate(this%qcharge_col       (begc:endc))                 ; this%qcharge_col       (:)     = nan
     allocate(this%fracice_col       (begc:endc,nlevgrnd))        ; this%fracice_col       (:,:)   = nan
     allocate(this%icefrac_col       (begc:endc,nlevgrnd))        ; this%icefrac_col       (:,:)   = nan
+#ifdef COUP_OAS_PFL
+    allocate(this%pfl_eff_porosity_col (begc:endc,nlevgrnd))     ; this%pfl_eff_porosity_col(:,:) = 0._r8
+#endif
     allocate(this%fcov_col          (begc:endc))                 ; this%fcov_col          (:)     = nan   
     allocate(this%fsat_col          (begc:endc))                 ; this%fsat_col          (:)     = nan
     allocate(this%h2osfc_thresh_col (begc:endc))                 ; this%h2osfc_thresh_col (:)     = nan
@@ -235,6 +241,7 @@ subroutine InitCold(this, bounds)
     ! averaging for the accum field
     do c = bounds%begc, bounds%endc
        this%num_substeps_col(c) = spval
+       this%icefrac_col(c,:) = spval
     end do
 
   end subroutine InitCold

src/clm5/biogeophys/SoilStateType.F90

@@ -293,8 +293,12 @@ subroutine InitHistory(this, bounds)
          avgflag='A', long_name='urban factor limiting ground evap', &
          ptr_col=this%soilalpha_u_col, set_nourb=spval, default='inactive')
 
+#ifdef COUP_OAS_PFL
+    if (.true.) then
+#else
     if (use_cn) then
-       this%watsat_col(begc:endc,:) = spval 
+#endif
+       this%watsat_col(begc:endc,:) = spval
        call hist_addfld2d (fname='watsat', units='m^3/m^3', type2d='levgrnd', &
             avgflag='A', long_name='water saturated', &
             ptr_col=this%watsat_col, default='inactive')

src/clm5/biogeophys/SoilWaterMovementMod.F90

@@ -1503,20 +1503,26 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, &
          h2osoi_ice        =>    waterstate_inst%h2osoi_ice_col     , & ! Output: [real(r8) (:,:) ]  ice lens (kg/m2)
          smpmin            =>    soilstate_inst%smpmin_col          , & ! Input:  [real(r8) (:)   ]  restriction for min of soil potential (mm)
          pfl_h2osoi_liq    =>    waterstate_inst%pfl_h2osoi_liq_col , & ! Input:  [real(r8) (:,:) ]  ParFlow soil water (mm)
-         pfl_psi           =>    waterstate_inst%pfl_psi_col          & ! Input:  [real(r8) (:,:) ]  ParFlow pressure head (mm)
+         pfl_psi           =>    waterstate_inst%pfl_psi_col        , & ! Input:  [real(r8) (:,:) ]  ParFlow pressure head (mm)
+         h2osoi_ice_prev   =>    waterstate_inst%h2osoi_ice_prev_col, & ! Input:  [real(r8) (:,:) ]  ice lens before PhaseChange (kg/m2)
+         pfl_eff_porosity  =>    soilhydrology_inst%pfl_eff_porosity_col & ! Output: [real(r8) (:,:) ]  effective porosity sent to ParFlow (m3/m3)
          )  ! end associate statement
 
 
          ! Exchange of soil water and pressure head between ParFlow and eCLM
-         ! ParFlow has no ice phase, so that the received water it the total water content.
-         ! Keep the ice what PhaseChanged diagnosed earlier in the time step and
-         ! attribute the remainder to liquid, so that liq and ice matches the received total.
+         ! ParFlow carries the liquid phase only, so the received water is the liquid
+         ! water content and h2osoi_ice stays as PhaseChange left it. ParFlow applies the
+         ! freeze one coupling interval later, so subtract it here to keep eCLM in sync.
          do fc = 1, num_hydrologyc
             c = filter_hydrologyc(fc)
 
             do j = 1, nlevgrnd
-               h2osoi_ice(c,j) = min(h2osoi_ice(c,j), pfl_h2osoi_liq(c,j))
-               h2osoi_liq(c,j) = max(0._r8, pfl_h2osoi_liq(c,j) - h2osoi_ice(c,j))
+               ! PhaseChange bounds the ice by the available water. Capping the state
+               ! keeps the frozen mass ParFlow derives from it exact.
+               h2osoi_ice(c,j) = min(h2osoi_ice(c,j), watsat(c,j)*dz(c,j)*denice)
+               h2osoi_liq(c,j) = max(0._r8, pfl_h2osoi_liq(c,j) &
+                                            - (h2osoi_ice(c,j) - h2osoi_ice_prev(c,j)))
+               pfl_eff_porosity(c,j) = watsat(c,j) - h2osoi_ice(c,j)/(dz(c,j)*denice)
                if (pfl_psi(c,j) <= 0) then
                   smp_l(c,j) = max(smpmin(c), pfl_psi(c,j))
                end if

src/clm5/biogeophys/WaterStateType.F90

@@ -108,6 +108,7 @@ module WaterstateType
 #ifdef COUP_OAS_PFL
      real(r8), pointer :: pfl_psi_col            (:,:) ! ParFlow pressure head   COUP_OAS_PFL
      real(r8), pointer :: pfl_h2osoi_liq_col     (:,:) ! ParFlow soil liquid     COUP_OAS_PFL
+     real(r8), pointer :: h2osoi_ice_prev_col    (:,:) ! ice lens at the start of the time step, before PhaseChange (kg/m2)
 #endif
      real(r8), pointer :: total_plant_stored_h2o_col(:)! col water that is bound in plants, including roots, sapwood, leaves, etc
                                                        ! in most cases, the vegetation scheme does not have a dynamic
@@ -247,6 +248,7 @@ subroutine InitAllocate(this, bounds)
 #ifdef COUP_OAS_PFL
     allocate(this%pfl_psi_col            (begc:endc,1:nlevgrnd))          ; this%pfl_psi_col           (:,:) = nan
     allocate(this%pfl_h2osoi_liq_col     (begc:endc,1:nlevgrnd))          ; this%pfl_h2osoi_liq_col    (:,:) = nan
+    allocate(this%h2osoi_ice_prev_col    (begc:endc,1:nlevgrnd))          ; this%h2osoi_ice_prev_col   (:,:) = nan
 #endif
     allocate(this%h2osoi_ice_tot_col     (begc:endc))                     ; this%h2osoi_ice_tot_col     (:)   = nan
     allocate(this%h2osoi_liq_tot_col     (begc:endc))                     ; this%h2osoi_liq_tot_col     (:)   = nan
@@ -861,6 +863,7 @@ subroutine InitCold(this, bounds, &
 #ifdef COUP_OAS_PFL
       this%pfl_psi_col(bounds%begc:bounds%endc,            1:) = -1000._r8
       this%pfl_h2osoi_liq_col(bounds%begc:bounds%endc,     1:) = spval
+      this%h2osoi_ice_prev_col(bounds%begc:bounds%endc,    1:) = 0._r8
 #endif
       this%h2osoi_vol_prs_grc(bounds%begg:bounds%endg,     1:) = spval
       this%h2osoi_liq_col(bounds%begc:bounds%endc,-nlevsno+1:) = spval

src/clm5/cpl/lnd_comp_mct.F90

@@ -18,7 +18,7 @@ module lnd_comp_mct
   use oas_sendReceiveMod, only : oas_receive_icon, oas_send_icon
 #endif                   
 #ifdef COUP_OAS_PFL
-   use oas_sendReceiveMod, only : oas_receive, oas_send
+   use oas_sendReceiveMod, only : oas_receive_parflow, oas_send_parflow
 #endif
 #endif
   !
@@ -474,7 +474,7 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l)
 #endif
 
 #ifdef COUP_OAS_PFL
-       call oas_receive(bounds, time_elapsed, atm2lnd_inst)
+       call oas_receive_parflow(bounds, time_elapsed, atm2lnd_inst)
 #endif
        ! Run clm 
        call t_barrierf('sync_clm_run1', mpicom)
@@ -493,7 +493,7 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l)
 #endif
 
 #if defined(COUP_OAS_PFL)
-       call oas_send(bounds, time_elapsed, lnd2atm_inst)
+       call oas_send_parflow(bounds, time_elapsed, lnd2atm_inst)
 #endif
        ! Create l2x_l export state - add river runoff input to l2x_l if appropriate
        call t_startf ('lc_lnd_export')

src/clm5/main/clm_driver.F90

@@ -1068,7 +1068,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
          solarabs_inst, drydepvel_inst,       &
          vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, &
          lnd2atm_inst, &
-#ifdef USE_PDAF
+#if defined(USE_PDAF) || defined(COUP_OAS_PFL)
          soilhydrology_inst, soilstate_inst, &
 #endif
          net_carbon_exchange_grc = net_carbon_exchange_grc(bounds_proc%begg:bounds_proc%endg))
@@ -1219,7 +1219,7 @@ subroutine clm_drv_init(bounds, &
     ! !USES:
     use shr_kind_mod       , only : r8 => shr_kind_r8
     use shr_infnan_mod     , only : nan => shr_infnan_nan, assignment(=)
-    use clm_varpar         , only : nlevsno, nlevsoi
+    use clm_varpar         , only : nlevsno, nlevsoi, nlevgrnd
     use CanopyStateType    , only : canopystate_type
     use WaterStateType     , only : waterstate_type
     use WaterFluxType      , only : waterflux_type
@@ -1257,6 +1257,7 @@ subroutine clm_drv_init(bounds, &
 #ifdef COUP_OAS_PFL
          pfl_psi            => waterstate_inst%pfl_psi_col               , & ! Input:  [real(r8) (:,:) ]  COUP_OAS_PFL
          pfl_h2osoi_liq     => waterstate_inst%pfl_h2osoi_liq_col        , & ! Input:  [real(r8) (:,:) ]  COUP_OAS_PFL
+         h2osoi_ice_prev    => waterstate_inst%h2osoi_ice_prev_col       , & ! Output: [real(r8) (:,:) ]  COUP_OAS_PFL
 #endif
          elai               => canopystate_inst%elai_patch               , & ! Input:  [real(r8) (:)   ]  one-sided leaf area index with burying by snow    
          esai               => canopystate_inst%esai_patch               , & ! Input:  [real(r8) (:)   ]  one-sided stem area index with burying by snow    
@@ -1314,6 +1315,8 @@ subroutine clm_drv_init(bounds, &
           pfl_psi(c,:) = atm2lnd_inst%pfl_psi_grc(g,:)
           pfl_h2osoi_liq(c,:) = atm2lnd_inst%pfl_h2osoi_liq_grc(g,:)
         end if
+        ! Get ice content before PhaseChange, to make it consistent with ParFlow
+        h2osoi_ice_prev(c,1:nlevgrnd) = h2osoi_ice(c,1:nlevgrnd)
       end do
 #endif
     end associate

src/clm5/main/lnd2atmMod.F90

@@ -43,9 +43,11 @@ module lnd2atmMod
   use landunit_varcon      , only : istice_mec, istsoil, istcrop
 #ifdef USE_PDAF
   use clm_time_manager     , only : get_nstep
+  use PatchType            , only : patch
+#endif
+#if defined(USE_PDAF) || defined(COUP_OAS_PFL)
   use SoilHydrologyType    , only : soilhydrology_type
   use SoilStateType        , only : soilstate_type
-  use PatchType            , only : patch
 #endif
   !
   ! !PUBLIC TYPES:
@@ -60,6 +62,10 @@ module lnd2atmMod
   ! !PRIVATE MEMBER FUNCTIONS:
   private :: handle_ice_runoff
 
+#ifdef COUP_OAS_PFL
+  real(r8), parameter, private :: pfl_uncoupled = -9999._r8
+#endif
+
   character(len=*), parameter, private :: sourcefile = &
        __FILE__
   !------------------------------------------------------------------------
@@ -136,7 +142,7 @@ subroutine lnd2atm(bounds, &
        solarabs_inst, drydepvel_inst,  &
        vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, &
        lnd2atm_inst, &
-#ifdef USE_PDAF
+#if defined(USE_PDAF) || defined(COUP_OAS_PFL)
        soilhydrology_inst, soilstate_inst, &
 #endif
        net_carbon_exchange_grc) 
@@ -165,16 +171,17 @@ subroutine lnd2atm(bounds, &
     type(ch4_type)              , intent(in)    :: ch4_inst
     type(glc_behavior_type)     , intent(in)    :: glc_behavior
     type(lnd2atm_type)          , intent(inout) :: lnd2atm_inst 
-#ifdef USE_PDAF
-    ! Yorck
+#if defined(USE_PDAF) || defined(COUP_OAS_PFL)
     type(soilhydrology_type)    , intent(inout) :: soilhydrology_inst
     type(soilstate_type)        , intent(inout) :: soilstate_inst
-    ! end Yorck
 #endif
     real(r8)                    , intent(in)    :: net_carbon_exchange_grc( bounds%begg: )  ! net carbon exchange between land and atmosphere, positive for source (gC/m2/s)
     !
     ! !LOCAL VARIABLES:
     integer  :: c, g, j  ! indices
+#ifdef COUP_OAS_PFL
+    real(r8) :: pfl_eff_porosity_wt(bounds%begg:bounds%endg) ! summed weight of the hydrologically active columns
+#endif
 #ifdef USE_PDAF
     integer  :: p, l, index, counter  ! indices
 #endif
@@ -466,6 +473,38 @@ subroutine lnd2atm(bounds, &
          waterflux_inst%qflx_parflow_col (bounds%begc:bounds%endc, :), &
          lnd2atm_inst%qflx_parflow_grc   (bounds%begg:bounds%endg, :), &
          c2l_scale_type= 'unity',  l2g_scale_type='unity' )
+
+    ! Porosity does not scale with size, so it cannot use c2g. Average over the hydrologically
+    ! active columns only. Not coupled cells get out of range value, must be <0.
+    pfl_eff_porosity_wt(bounds%begg:bounds%endg) = 0._r8
+    do c = bounds%begc, bounds%endc
+       if (col%hydrologically_active(c) .and. col%wtgcell(c) > 0._r8) then
+          g = col%gridcell(c)
+          pfl_eff_porosity_wt(g) = pfl_eff_porosity_wt(g) + col%wtgcell(c)
+       end if
+    end do
+    do g = bounds%begg, bounds%endg
+       if (pfl_eff_porosity_wt(g) > 0._r8) then
+          lnd2atm_inst%pfl_eff_porosity_grc(g,:) = 0._r8
+       else
+          lnd2atm_inst%pfl_eff_porosity_grc(g,:) = pfl_uncoupled
+       end if
+    end do
+    do c = bounds%begc, bounds%endc
+       if (col%hydrologically_active(c) .and. col%wtgcell(c) > 0._r8) then
+          g = col%gridcell(c)
+          do j = 1, nlevgrnd
+             lnd2atm_inst%pfl_eff_porosity_grc(g,j) = lnd2atm_inst%pfl_eff_porosity_grc(g,j) &
+                  + col%wtgcell(c) * soilhydrology_inst%pfl_eff_porosity_col(c,j)
+          end do
+       end if
+    end do
+    do g = bounds%begg, bounds%endg
+       if (pfl_eff_porosity_wt(g) > 0._r8) then
+          lnd2atm_inst%pfl_eff_porosity_grc(g,:) = lnd2atm_inst%pfl_eff_porosity_grc(g,:) &
+               / pfl_eff_porosity_wt(g)
+       end if
+    end do
 #endif
 
 #ifdef USE_PDAF

src/clm5/main/lnd2atmType.F90

@@ -77,6 +77,7 @@ module lnd2atmType
      real(r8), pointer :: qflx_liq_from_ice_col   (:)   => null() ! liquid runoff from converted ice runoff
 #ifdef COUP_OAS_PFL
      real(r8), pointer :: qflx_parflow_grc        (:,:) => null() ! source/sink flux per soil layer sent to ParFlow [1/hr] [- out from root]
+     real(r8), pointer :: pfl_eff_porosity_grc    (:,:) => null() ! effective porosity per soil layer sent to ParFlow [m3/m3]
 #endif
      real(r8), pointer :: qirrig_grc              (:)   => null() ! irrigation flux
 
@@ -191,6 +192,7 @@ subroutine InitAllocate(this, bounds)
     allocate(this%qflx_liq_from_ice_col(begc:endc))          ; this%qflx_liq_from_ice_col(:)   =ival
 #ifdef COUP_OAS_PFL
     allocate(this%qflx_parflow_grc     (begg:endg,1:nlevsoi)); this%qflx_parflow_grc     (:,:) =ival
+    allocate(this%pfl_eff_porosity_grc (begg:endg,1:nlevgrnd)); this%pfl_eff_porosity_grc(:,:) =-9999._r8
 #endif
     allocate(this%qirrig_grc           (begg:endg))          ; this%qirrig_grc           (:)   =ival
 
@@ -344,6 +346,10 @@ subroutine InitHistory(this, bounds)
     call hist_addfld2d (fname='QPARFLOW_TO_OASIS', units='1/hr', type2d='levsoi', &
          avgflag='A', long_name='source/sink flux per soil layer sent to ParFlow', &
          ptr_lnd=this%qflx_parflow_grc, default='inactive')
+
+    call hist_addfld2d (fname='EFF_POROSITY_TO_OASIS', units='m3/m3', type2d='levgrnd', &
+         avgflag='A', long_name='effective porosity per soil layer sent to ParFlow', &
+         ptr_lnd=this%pfl_eff_porosity_grc, default='inactive')
 #endif
   end subroutine InitHistory
 

src/clm5/oasis3/oas_defineMod.F90

@@ -102,12 +102,13 @@ subroutine oas_definitions_init(bounds)
     var_nodims(1) = 1         ! unused
     var_nodims(2) = nlevsoi   ! number of fields in a bundle
 
-    call oasis_def_var(oas_et_loss_id, "ECLM_ET", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) 
-    
+    call oasis_def_var(oas_et_loss_id, "ECLM_ET", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror)
+
     var_nodims(2) = nlevgrnd         ! number of fields in a bundle
     call oasis_def_var(oas_sat_id, "ECLM_SOILLIQ", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror)
     call oasis_def_var(oas_psi_id, "ECLM_PSI", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror)
-#endif 
+    call oasis_def_var(oas_eff_porosity_id, "ECLM_EFFPOROSITY", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror)
+#endif
 
 #ifdef COUP_OAS_ICON
 

src/clm5/oasis3/oas_sendReceiveMod.F90

@@ -11,8 +11,8 @@ module oas_sendReceiveMod
   private
 
 #ifdef COUP_OAS_PFL
-  public  :: oas_send
-  public  :: oas_receive
+  public  :: oas_send_parflow
+  public  :: oas_receive_parflow
 #endif
 
 #ifdef COUP_OAS_ICON
@@ -23,7 +23,7 @@ module oas_sendReceiveMod
 contains
 
 #ifdef COUP_OAS_PFL
-  subroutine oas_receive(bounds, seconds_elapsed, atm2lnd_inst)
+  subroutine oas_receive_parflow(bounds, seconds_elapsed, atm2lnd_inst)
     use atm2lndType, only: atm2lnd_type
 
     type(bounds_type),  intent(in)    :: bounds
@@ -34,9 +34,9 @@ subroutine oas_receive(bounds, seconds_elapsed, atm2lnd_inst)
     call oasis_get(oas_psi_id, seconds_elapsed, atm2lnd_inst%pfl_psi_grc, info)
     call oasis_get(oas_sat_id, seconds_elapsed, atm2lnd_inst%pfl_h2osoi_liq_grc, info)
 
-  end subroutine oas_receive
+  end subroutine oas_receive_parflow
 
-  subroutine oas_send(bounds, seconds_elapsed, lnd2atm_inst)
+  subroutine oas_send_parflow(bounds, seconds_elapsed, lnd2atm_inst)
     use lnd2atmType, only : lnd2atm_type
     use spmdMod,     only : mpicom
     use shr_mpi_mod, only: shr_mpi_barrier
@@ -48,8 +48,9 @@ subroutine oas_send(bounds, seconds_elapsed, lnd2atm_inst)
     integer                           :: info
 
     call oasis_put(oas_et_loss_id, seconds_elapsed, lnd2atm_inst%qflx_parflow_grc, info)
-
-  end subroutine oas_send
+    call oasis_put(oas_eff_porosity_id, seconds_elapsed, lnd2atm_inst%pfl_eff_porosity_grc, info)
+
+  end subroutine oas_send_parflow
 #endif
 
 #ifdef COUP_OAS_ICON

src/clm5/oasis3/oas_vardefMod.F90

@@ -3,7 +3,7 @@ module oas_vardefMod
   save
 
 #ifdef COUP_OAS_PFL
-  integer :: oas_psi_id, oas_et_loss_id, oas_sat_id
+  integer :: oas_psi_id, oas_et_loss_id, oas_sat_id, oas_eff_porosity_id
 #endif
 
 #ifdef COUP_OAS_ICON