Scale qflx_parflow at the column instead of the gridcell level (#62)

* fix: change qflx_parflow_col unit to 1/hr, remove unit transformation for qflx_parflow_grc

* qflx_parflow unit transformation r8 compliant

* set nan values of qflx_parflow_grc to zero

* clean-up

* fix qflx_drain

* use array operation for converting qflx_parflow units

Author: s-poll

src/clm5/biogeophys/SoilHydrologyMod.F90

@@ -2362,9 +2362,13 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, &
      ! !LOCAL VARIABLES:
      character(len=32) :: subname = 'ParFlowDrainage'   ! subroutine name
      integer  :: c,j,fc                                   ! indices
+     real(r8), parameter :: m_per_mm = 1.e-3_r8          ! 0.001 meters per mm
+     real(r8), parameter :: sec_per_hr = 3600._r8        ! 3600 s in 1 hour
+
      !-----------------------------------------------------------------------
 
-     associate(                                                            & 
+     associate(                                                            &
+          dz                 =>    col%dz                                , & ! Input:  [real(r8) (:,:) ] layer depth (m)
           qflx_snwcp_liq     =>    waterflux_inst%qflx_snwcp_liq_col     , & ! excess rainfall due to snow capping (mm H2O /s) [+]
           qflx_drain         =>    waterflux_inst%qflx_drain_col         , & ! sub-surface runoff (mm H2O /s)
           qflx_drain_perched =>    waterflux_inst%qflx_drain_perched_col , & ! perched wt sub-surface runoff (mm H2O /s)         
@@ -2377,25 +2381,25 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, &
 
          ! COUP_OAS_PFL
          ! Calculate here the source/sink term for ParFlow
-         do j = 1, nlevsoi
-            do fc = 1, num_hydrologyc
-               c = filter_hydrologyc(fc)
+         do fc = 1, num_hydrologyc
+            c = filter_hydrologyc(fc)
+            do j = 1, nlevsoi
                if (j == 1) then
                   ! From SoilWaterPlantSinkMod:
                   ! qflx_rootsoi_col(c,j) = rootr_col(c,j)*qflx_tran_veg_col(c)
-                  qflx_parflow(c,j) = (qflx_infl(c) - qflx_rootsoi(c,j)) !* 3.6_r8 / dz(c,j)
+                   qflx_parflow(c,j) = (qflx_infl(c) - qflx_rootsoi(c,j)) !mm/s
                else 
-                  qflx_parflow(c,j) = -qflx_rootsoi(c,j) !* 3.6_r8 / dz(c,j)
+                   qflx_parflow(c,j) = -qflx_rootsoi(c,j) !mm/s
                end if
             end do
-         end do
-
-         do fc = 1, num_hydrologyc
-            c = filter_hydrologyc(fc)
-            qflx_drain(c)         = -sum(qflx_parflow(c,:)) !0._r8
+            ! Compute subsurface run-off (mm/s)
+            qflx_drain(c) = -sum(qflx_parflow(c,:))
             qflx_drain_perched(c) = 0._r8  
             qflx_rsub_sat(c)      = 0._r8
             qflx_qrgwl(c)         = qflx_snwcp_liq(c)   ! Set imbalance for snow capping
+            ! Convert eCLM fluxes (mm/s) to ParFlow fluxes (1/hr):
+            !         1/hr            =           [mm/s]          *   [s/hr]   *  [m/mm]  *         [1/m]
+            qflx_parflow(c,1:nlevsoi) = qflx_parflow(c,1:nlevsoi) * sec_per_hr * m_per_mm * (1._r8/dz(c,1:nlevsoi))
          end do
 
          ! No drainage for urban columns (necessary to account for water balance errors)

src/clm5/biogeophys/WaterfluxType.F90

@@ -73,7 +73,7 @@ module WaterfluxType
      real(r8), pointer :: qflx_adv_col             (:,:) ! col advective flux across different soil layer interfaces [mm H2O/s] [+ downward]
      real(r8), pointer :: qflx_rootsoi_col         (:,:) ! col root and soil water exchange [mm H2O/s] [+ into root]
 #ifdef COUP_OAS_PFL
-     real(r8), pointer :: qflx_parflow_col         (:,:) ! col source/sink flux per soil layer sent to ParFlow [mm H2O/s] [- out from root]
+     real(r8), pointer :: qflx_parflow_col         (:,:) ! col source/sink flux per soil layer sent to ParFlow [1/hr] [- out from root]
 #endif
      real(r8), pointer :: qflx_infl_col            (:)   ! col infiltration (mm H2O /s)
      real(r8), pointer :: qflx_surf_col            (:)   ! col surface runoff (mm H2O /s)

src/clm5/main/lnd2atmMod.F90

@@ -162,8 +162,6 @@ subroutine lnd2atm(bounds, &
     real(r8), parameter :: amC   = 12.0_r8                      ! Atomic mass number for Carbon
     real(r8), parameter :: amO   = 16.0_r8                      ! Atomic mass number for Oxygen
     real(r8), parameter :: amCO2 = amC + 2.0_r8*amO             ! Atomic mass number for CO2
-    real(r8), parameter :: m_per_mm = 1.e-3_r8                  ! 0.001 meters per mm
-    real(r8), parameter :: sec_per_hr = 3600                    ! 3600 s in 1 hour
     ! The following converts g of C to kg of CO2
     real(r8), parameter :: convertgC2kgCO2 = 1.0e-3_r8 * (amCO2/amC)
     !------------------------------------------------------------------------
@@ -433,18 +431,15 @@ subroutine lnd2atm(bounds, &
          lnd2atm_inst%qflx_parflow_grc   (bounds%begg:bounds%endg, :), &
          c2l_scale_type= 'unity',  l2g_scale_type='unity' )
 
-    do c = bounds%begc, bounds%endc
-     if (col%hydrologically_active(c)) then
-       if (col%itype(c) == istsoil .or. col%itype(c) == istcrop) then
-         g = col%gridcell(c)
-         do j = 1, nlevsoi
-          ! Convert eCLM fluxes (mm/s) to ParFlow fluxes (1/hr): 
-          !             1/hr                 =             [mm/s]                 *   [s/hr]   *  [m/mm]  *     [1/m] 
-          lnd2atm_inst%qflx_parflow_grc(g,j) = lnd2atm_inst%qflx_parflow_grc(g,j) * sec_per_hr * m_per_mm * (1/col%dz(c,j))
-         enddo
-       end if
-     end if
-    end do
+    ! adjust nan values after c2g, need to be rechecked
+    do g = bounds%begg, bounds%endg
+     do j = 1, nlevsoi
+      if (lnd2atm_inst%qflx_parflow_grc(g,j) == spval) then
+       lnd2atm_inst%qflx_parflow_grc(g,j) = 0._r8
+       write(iulog,*)'WARNING: qflx_parflow_grc is nan at grid point ',g,' level',j,' replaced with 0.'
+      end if
+     end do
+    enddo
 #endif
 
   end subroutine lnd2atm