eCLM-ICON coupling (ESCOMP#18)

Model developments by @s-poll :

* 9f00031 introduce precompiler COUP_OAS_ICON and COUP_OAS_PFL
* 59702f3 Oasis define for ICON ()
* a98d2cd Oasis snd/rcv for ICON (CLM3 vars)
* 1641475 consistent coupler precompiler naming
* f3c0cad include oas_receive_icon in eclm time stepping
* f7be360 renaming of coupled ICON-ECLM vars in oas_defineMod
* 27d984d Update snd/rcv fields from ICON
* ed1af99 Implementation of coupled variable t_sf.
* 47c2068 Couple rain_rate and snow_rate in seperate variables.
* 60850cf Restructure oasis_def_var.
* a42a3eb Inclusion of urban landunit.

Merges and bugfixes by @kvrigor :

* d79df54 Merge clm5nl-gen v0.6 (ESCOMP#16)
* 8ffe78f Merge eCLM-ParFlow coupling (ESCOMP#17)
* dfb6e69 Added compile definition COUP_OAS_PFL for ParFlow
* 0263ae8 Removed `use_parflow_soilwater()` in `BalanceCheckMod.F90` and `TotalWaterAndHeatMod.F90`
* 9ea5952 Fixed missing `psit` calculation when not coupled with ParFlow

Co-authored-by: Stefan Poll <spoll@users.noreply.github.com>

src/clm5/biogeophys/BalanceCheckMod.F90

@@ -184,7 +184,6 @@ subroutine BalanceCheck( bounds, &
      use clm_time_manager     , only : get_nstep_since_startup_or_lastDA_restart_or_pause
      use CanopyStateType      , only : canopystate_type
      use SurfaceAlbedoType    , only : surfalb_type
-     use SoilWaterMovementMod , only : use_parflow_soilwater
      use subgridAveMod
      !
      ! !ARGUMENTS:
@@ -228,7 +227,6 @@ subroutine BalanceCheck( bounds, &
           errh2osno               =>    waterstate_inst%errh2osno_col           , & ! Output: [real(r8) (:)   ]  error in h2osno (kg m-2)                
           endwb                   =>    waterstate_inst%endwb_col               , & ! Output: [real(r8) (:)   ]  water mass end of the time step         
           total_plant_stored_h2o_col => waterstate_inst%total_plant_stored_h2o_col, & ! Input: [real(r8) (:)   ]  water mass in plant tissues (kg m-2)
-          pfl_psi                 =>    waterstate_inst%pfl_psi_col             , & ! Input:  [real(r8) (:,:) ]  COUP_OAS_PFL
 
           qflx_rootsoi_col        =>    waterflux_inst%qflx_rootsoi_col         , & ! Input   [real(r8) (:)   ]  water loss in soil layers to root uptake (mm H2O/s)
                                                                                     !                            (ie transpiration demand, often = transpiration)
@@ -404,14 +402,13 @@ subroutine BalanceCheck( bounds, &
              write(iulog,*)'deltawb/dtime    = ',(endwb(indexc)-begwb(indexc))/dtime
              write(iulog,*)'deltaflux        = ',forc_rain_col(indexc)+forc_snow_col(indexc) - (qflx_evap_tot(indexc) + &
                   qflx_surf(indexc)+qflx_h2osfc_surf(indexc)+qflx_drain(indexc))
-
-             if (use_parflow_soilwater()) then
-               ! TODO: Balance errors must be fixed for fully coupled model (ICON-eCLM-ParFlow)
-               write(iulog,*)'Ignoring water balance error...'
-             else
-               write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
-               call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
-             end if
+#ifdef COUP_OAS_PFL
+             ! TODO: Balance errors must be fixed for fully coupled model (ICON-eCLM-ParFlow)
+             write(iulog,*)'Ignoring water balance error...'
+#else
+             write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
+             call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
+#endif
           else if (abs(errh2o(indexc)) > 1.e-5_r8 .and. (DAnstep > skip_steps) ) then
 
 
@@ -440,13 +437,13 @@ subroutine BalanceCheck( bounds, &
              write(iulog,*)'qflx_drain_perched         = ',qflx_drain_perched(indexc)*dtime
              write(iulog,*)'qflx_glcice_dyn_water_flux = ',qflx_glcice_dyn_water_flux(indexc)*dtime
              write(iulog,*)'qflx_rootsoi_col(1:nlevsoil)  = ',qflx_rootsoi_col(indexc,:)*dtime
-             if (use_parflow_soilwater()) then
-               ! TODO: Balance errors must be fixed for fully coupled model (ICON-eCLM-ParFlow)
-               write(iulog,*)'Ignoring water balance error...'
-             else
-               write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
-               call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
-             end if
+#ifdef COUP_OAS_PFL
+             ! TODO: Balance errors must be fixed for fully coupled model (ICON-eCLM-ParFlow)
+             write(iulog,*)'Ignoring water balance error...'
+#else
+             write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
+             call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
+#endif
           end if
        end if
 
@@ -693,7 +690,6 @@ subroutine BalanceCheck( bounds, &
        end if
 
        ! Soil energy balance check
-      ! COUP_OAS_PFL
        found = .false.
        do c = bounds%begc,bounds%endc
           if (col%active(c)) then
@@ -707,13 +703,13 @@ subroutine BalanceCheck( bounds, &
           write(iulog,*)'WARNING: BalanceCheck: soil balance error (W/m2)'
           write(iulog,*)'nstep         = ',nstep
           write(iulog,*)'errsoi_col    = ',errsoi_col(indexc)
-          if (use_parflow_soilwater()) then
-            ! TODO: Balance errors must be fixed for fully coupled model (ICON-eCLM-ParFlow)
-            write(iulog,*)'Ignoring soil balance error...'
-          else
-            write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
-            call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
-          end if
+#ifdef COUP_OAS_PFL
+          ! TODO: Balance errors must be fixed for fully coupled model (ICON-eCLM-ParFlow)
+          write(iulog,*)'Ignoring soil balance error...'
+#else
+          write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
+          call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
+#endif
        end if
 
      end associate

src/clm5/biogeophys/BareGroundFluxesMod.F90

@@ -195,7 +195,10 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, &
          t_ref2m                => temperature_inst%t_ref2m_patch               , & ! Output: [real(r8) (:)   ]  2 m height surface air temperature (Kelvin)                           
          t_ref2m_r              => temperature_inst%t_ref2m_r_patch             , & ! Output: [real(r8) (:)   ]  Rural 2 m height surface air temperature (Kelvin)                     
          t_veg                  => temperature_inst%t_veg_patch                 , & ! Output: [real(r8) (:)   ]  vegetation temperature (Kelvin)                                       
-
+#ifdef COUP_OAS_ICON
+         t_sf_patch             => temperature_inst%t_sf_patch                  , & ! Output: [real(r8) (:)   ]  patch surface temperature (K)
+!         q_sf_patch             => waterstate_inst%q_sf_patch                   , & ! Output: [real(r8) (:)   ]  patch surface humidity (kg/kg)
+#endif
          q_ref2m                => waterstate_inst%q_ref2m_patch                , & ! Output: [real(r8) (:)   ]  2 m height surface specific humidity (kg/kg)                          
          rh_ref2m_r             => waterstate_inst%rh_ref2m_r_patch             , & ! Output: [real(r8) (:)   ]  Rural 2 m height surface relative humidity (%)                        
          rh_ref2m               => waterstate_inst%rh_ref2m_patch               , & ! Output: [real(r8) (:)   ]  2 m height surface relative humidity (%)                              
@@ -374,6 +377,11 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, &
          qflx_ev_soil(p)   = -raiw*(forc_q(c) - qg_soil(c))
          qflx_ev_h2osfc(p) = -raiw*(forc_q(c) - qg_h2osfc(c))
 
+#ifdef COUP_OAS_ICON
+         t_sf_patch(p)     = t_grnd(c)
+!         q_sf_patch(p)     = qg(c)
+#endif
+
          ! 2 m height air temperature
          t_ref2m(p) = thm(p) + temp1(p)*dth(p)*(1._r8/temp12m(p) - 1._r8/temp1(p))
 

src/clm5/biogeophys/CanopyFluxesMod.F90

@@ -453,6 +453,10 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
          t_ref2m                => temperature_inst%t_ref2m_patch               , & ! Output: [real(r8) (:)   ]  2 m height surface air temperature (Kelvin)                           
          t_ref2m_r              => temperature_inst%t_ref2m_r_patch             , & ! Output: [real(r8) (:)   ]  Rural 2 m height surface air temperature (Kelvin)                     
          t_skin_patch           => temperature_inst%t_skin_patch                , & ! Output: [real(r8) (:)   ]  patch skin temperature (K)  
+#ifdef COUP_OAS_ICON
+         t_sf_patch             => temperature_inst%t_sf_patch                  , & ! Output: [real(r8) (:)   ]  patch surface temperature (K)
+!         q_sf_patch             => waterstate_inst%q_sf_patch                   , & ! Output: [real(r8) (:)   ]  patch surface humidity (kg/kg)
+#endif
 
          frac_h2osfc            => waterstate_inst%frac_h2osfc_col              , & ! Input:  [real(r8) (:)   ]  fraction of surface water                                             
          fwet                   => waterstate_inst%fwet_patch                   , & ! Input:  [real(r8) (:)   ]  fraction of canopy that is wet (0 to 1)                               
@@ -1177,6 +1181,10 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
          delq_h2osfc = wtalq(p)*qg_h2osfc(c)-wtlq0(p)*qsatl(p)-wtaq0(p)*forc_q(c)
          qflx_ev_h2osfc(p) = forc_rho(c)*wtgq(p)*delq_h2osfc
 
+#ifdef COUP_OAS_ICON
+         t_sf_patch(p)  = taf(p)
+!         q_sf(p)  = qaf(p)
+#endif
          ! 2 m height air temperature
 
          t_ref2m(p) = thm(p) + temp1(p)*dth(p)*(1._r8/temp12m(p) - 1._r8/temp1(p))

src/clm5/biogeophys/CanopyTemperatureMod.F90

@@ -151,9 +151,10 @@ subroutine CanopyTemperature(bounds, &
          qg_soil          =>    waterstate_inst%qg_soil_col           , & ! Output: [real(r8) (:)   ] specific humidity at soil surface [kg/kg]
          qg               =>    waterstate_inst%qg_col                , & ! Output: [real(r8) (:)   ] ground specific humidity [kg/kg]         
          qg_h2osfc        =>    waterstate_inst%qg_h2osfc_col         , & ! Output: [real(r8) (:)   ]  specific humidity at h2osfc surface [kg/kg]
-         dqgdT            =>    waterstate_inst%dqgdT_col             , & ! Output: [real(r8) (:)   ] d(qg)/dT                                 
+         dqgdT            =>    waterstate_inst%dqgdT_col             , & ! Output: [real(r8) (:)   ] d(qg)/dT
+#ifdef COUP_OAS_PFL
          pfl_psi          =>    waterstate_inst%pfl_psi_col           , & ! Input:  [real(r8) (:,:) ] COUP_OAS_PFL
-
+#endif
          qflx_evap_tot    =>    waterflux_inst%qflx_evap_tot_patch    , & ! Output: [real(r8) (:)   ] qflx_evap_soi + qflx_evap_can + qflx_tran_veg
          qflx_evap_veg    =>    waterflux_inst%qflx_evap_veg_patch    , & ! Output: [real(r8) (:)   ] vegetation evaporation (mm H2O/s) (+ = to atm)
          qflx_tran_veg    =>    waterflux_inst%qflx_tran_veg_patch    , & ! Output: [real(r8) (:)   ] vegetation transpiration (mm H2O/s) (+ = to atm)
@@ -257,15 +258,14 @@ subroutine CanopyTemperature(bounds, &
                wx   = (h2osoi_liq(c,1)/denh2o+h2osoi_ice(c,1)/denice)/dz(c,1)
                fac  = min(1._r8, wx/watsat(c,1))
                fac  = max( fac, 0.01_r8 )
+#ifdef COUP_OAS_PFL
                ! clm3.5/bld/usr.src/Biogeophysics1Mod.F90
-               ! if COUP_OAS_PFL
                if (pfl_psi(c,1)>= 0.0_r8)  psit = 0._r8
                if (pfl_psi(c,1) < 0.0_r8)  psit = pfl_psi(c,1)
-               !else
-               !psit = -sucsat(c,1) * fac ** (-bsw(c,1))
-               !psit = max(smpmin(c), psit)
-               !end if
-
+#else
+               psit = -sucsat(c,1) * fac ** (-bsw(c,1))
+               psit = max(smpmin(c), psit)
+#endif
                ! modify qred to account for h2osfc
                hr   = exp(psit/roverg/t_soisno(c,1))
                qred = (1._r8 - frac_sno(c) - frac_h2osfc(c))*hr &

src/clm5/biogeophys/HydrologyDrainageMod.F90

@@ -52,8 +52,13 @@ subroutine HydrologyDrainage(bounds,               &
     use clm_varctl       , only : use_vichydro
     use clm_varpar       , only : nlevgrnd, nlevurb, nlevsoi
     use clm_time_manager , only : get_step_size, get_nstep
+#if defined(COUP_OAS_PFL)
     use SoilHydrologyMod , only : CLMVICMap, Drainage, PerchedLateralFlow, LateralFlowPowerLaw, ParFlowDrainage
     use SoilWaterMovementMod , only : use_aquifer_layer, use_parflow_soilwater
+#else
+    use SoilHydrologyMod , only : CLMVICMap, Drainage, PerchedLateralFlow, LateralFlowPowerLaw
+    use SoilWaterMovementMod , only : use_aquifer_layer
+#endif
     !
     ! !ARGUMENTS:
     type(bounds_type)        , intent(in)    :: bounds               
@@ -125,7 +130,7 @@ subroutine HydrologyDrainage(bounds,               &
          call CLMVICMap(bounds, num_hydrologyc, filter_hydrologyc, &
               soilhydrology_inst, waterstate_inst)
       endif
-
+#if defined(COUP_OAS_PFL)
       if (use_parflow_soilwater()) then
          ! clm3.5/bld/usr.src/SoilHydrologyMod.F90
          ! ignore drainage calculations in eCLM and instead pass these fluxes to ParFlow 
@@ -148,7 +153,24 @@ subroutine HydrologyDrainage(bounds,               &
                 soilhydrology_inst, soilstate_inst, &
                 waterstate_inst, waterflux_inst)
 
-         endif
+      endif
+#else
+      if (use_aquifer_layer()) then 
+         call Drainage(bounds, num_hydrologyc, filter_hydrologyc, &
+            num_urbanc, filter_urbanc,&
+            temperature_inst, soilhydrology_inst, soilstate_inst, &
+            waterstate_inst, waterflux_inst)
+      else
+         call PerchedLateralFlow(bounds, num_hydrologyc, filter_hydrologyc, &
+            num_urbanc, filter_urbanc,&
+            soilhydrology_inst, soilstate_inst, &
+            waterstate_inst, waterflux_inst)
+
+         call LateralFlowPowerLaw(bounds, num_hydrologyc, filter_hydrologyc, &
+            num_urbanc, filter_urbanc,&
+            soilhydrology_inst, soilstate_inst, &
+            waterstate_inst, waterflux_inst)
+#endif
       endif
 
       do j = 1, nlevgrnd

src/clm5/biogeophys/HydrologyNoDrainageMod.F90

@@ -72,7 +72,11 @@ subroutine HydrologyNoDrainage(bounds, &
     use SnowHydrologyMod     , only : SnowWater, BuildSnowFilter 
     use SoilHydrologyMod     , only : CLMVICMap, SurfaceRunoff, Infiltration, WaterTable, PerchedWaterTable
     use SoilHydrologyMod     , only : ThetaBasedWaterTable, RenewCondensation
+#ifdef COUP_OAS_PFL
     use SoilWaterMovementMod , only : SoilWater, use_parflow_soilwater
+#else
+    use SoilWaterMovementMod , only : SoilWater
+#endif
     use SoilWaterRetentionCurveMod, only : soil_water_retention_curve_type
     use SoilWaterMovementMod , only : use_aquifer_layer
     use SoilWaterPlantSinkMod , only : Compute_EffecRootFrac_And_VertTranSink
@@ -155,9 +159,7 @@ subroutine HydrologyNoDrainage(bounds, &
          h2osoi_vol         => waterstate_inst%h2osoi_vol_col         , & ! Output: [real(r8) (:,:) ]  volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3]
          h2osno_top         => waterstate_inst%h2osno_top_col         , & ! Output: [real(r8) (:)   ]  mass of snow in top layer (col) [kg]    
          wf                 => waterstate_inst%wf_col                 , & ! Output: [real(r8) (:)   ]  soil water as frac. of whc for top 0.05 m 
-         wf2                => waterstate_inst%wf2_col                , & ! Output: [real(r8) (:)   ]  soil water as frac. of whc for top 0.17 m 
-         pfl_psi            => waterstate_inst%pfl_psi_col            , & ! Input:  [real(r8) (:,:) ]  COUP_OAS_PFL
-
+         wf2                => waterstate_inst%wf2_col                , & ! Output: [real(r8) (:)   ]  soil water as frac. of whc for top 0.17 m
          watsat             => soilstate_inst%watsat_col              , & ! Input:  [real(r8) (:,:) ]  volumetric soil water at saturation (porosity)
          sucsat             => soilstate_inst%sucsat_col              , & ! Input:  [real(r8) (:,:) ]  minimum soil suction (mm)             
          bsw                => soilstate_inst%bsw_col                 , & ! Input:  [real(r8) (:,:) ]  Clapp and Hornberger "b"              
@@ -476,6 +478,7 @@ subroutine HydrologyNoDrainage(bounds, &
       ! ZMS: Note, this form, which seems to be the same as used in SoilWater, DOES NOT distinguish between
       ! ice and water volume, in contrast to the soilpsi calculation above. It won't be used in ch4Mod if
       ! t_soisno <= tfrz, though.
+#ifdef COUP_OAS_PFL
       if (.not. use_parflow_soilwater()) then
          do j = 1, nlevgrnd
             do fc = 1, num_hydrologyc
@@ -489,7 +492,19 @@ subroutine HydrologyNoDrainage(bounds, &
             end do
          end do
       endif
+#else
+      do j = 1, nlevgrnd
+         do fc = 1, num_hydrologyc
+            c = filter_hydrologyc(fc)
 
+            s_node = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8)
+            s_node = min(1.0_r8, s_node)
+
+            smp_l(c,j) = -sucsat(c,j)*s_node**(-bsw(c,j))
+            smp_l(c,j) = max(smpmin(c), smp_l(c,j))
+         end do
+      end do
+#endif
 
  !     if (use_cn) then
          ! Available soil water up to a depth of 0.05 m.

src/clm5/biogeophys/SoilFluxesMod.F90

@@ -331,11 +331,12 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
             ! for evaporation partitioning between liquid evap and ice sublimation, 
             ! use the ratio of liquid to (liquid+ice) in the top layer to determine split
             if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0.) then
+#ifdef COUP_OAS_PFL
                ! clm3.5/bld/usr.src/Biogeophysics2Mod.F90
-               ! if COUP_OAS_PFL
                qflx_evap_grnd(p) = qflx_ev_snow(p)*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j)))
-               !else
-               !qflx_evap_grnd(p) = max(qflx_ev_snow(p)*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
+#else
+               qflx_evap_grnd(p) = max(qflx_ev_snow(p)*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
+#endif
             else
                qflx_evap_grnd(p) = 0.
             end if

src/clm5/biogeophys/SoilHydrologyMod.F90

@@ -257,18 +257,19 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
       do fc = 1, num_hydrologyc
          c = filter_hydrologyc(fc)
 
+
+#ifdef COUP_OAS_PFL
          ! clm3.5/bld/usr.src/SoilHydrologyMod.F90
-         ! if COUP_OAS_PFL
          qflx_surf(c) = 0._r8
-
-         !else
+#else
          ! assume qinmax large relative to qflx_top_soil in control
-         !if (origflag == 1) then
-         !   qflx_surf(c) =  fcov(c) * qflx_top_soil(c)
-         !else
-         !   ! only send fast runoff directly to streams
-         !   qflx_surf(c) =   fsat(c) * qflx_top_soil(c)
-         !endif
+         if (origflag == 1) then
+           qflx_surf(c) =  fcov(c) * qflx_top_soil(c)
+         else
+           ! only send fast runoff directly to streams
+           qflx_surf(c) =   fsat(c) * qflx_top_soil(c)
+         endif
+#endif
       end do
 
       ! Determine water in excess of ponding limit for urban roof and impervious road.
@@ -2339,6 +2340,7 @@ subroutine RenewCondensation(bounds, num_hydrologyc, filter_hydrologyc, &
    end subroutine RenewCondensation
 
 !#8
+#ifdef COUP_OAS_PFL
    !-----------------------------------------------------------------------
    subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, &
         num_urbanc, filter_urbanc, waterflux_inst)
@@ -2407,5 +2409,6 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, &
          end do
      end associate
    end subroutine ParFlowDrainage
+#endif
 !#0
 end module SoilHydrologyMod