ICON-eCLM: exchange-coefficient coupling scheme

src/clm5/biogeophys/BareGroundFluxesMod.F90

@@ -197,7 +197,9 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, &
          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)
+         q_sf_patch             => waterstate_inst%q_sf_patch                   , & ! Output: [real(r8) (:)   ]  patch surface humidity (kg/kg)
+         rah1                   => frictionvel_inst%rah1_patch                  , & ! Output: [real(r8) (:)   ]  heat resistance (s/m)
+         raw1                   => frictionvel_inst%raw1_patch                  , & ! Output: [real(r8) (:)   ]  moisture resistance (s/m)
 #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 (%)                        
@@ -378,8 +380,16 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, &
          qflx_ev_h2osfc(p) = -raiw*(forc_q(c) - qg_h2osfc(c))
 
 #ifdef COUP_OAS_ICON
+         !pass value to global variable
          t_sf_patch(p)     = t_grnd(c)
-!         q_sf_patch(p)     = qg(c)
+         q_sf_patch(p)     = qg(c) ! sfc specific humidity for cplscheme_exchcoef
+         rah1(p) = rah  ! heat resistance for cplscheme_exchcoef
+         ! effective moisture resistance
+         if (raiw > 0._r8) then
+            raw1(p) = forc_rho(c)/raiw
+         else
+            raw1(p) = 1.e6_r8 ! soilresis imprevious ground
+         end if
 #endif
 
          ! 2 m height air temperature

src/clm5/biogeophys/CanopyFluxesMod.F90

@@ -455,7 +455,9 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
          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)
+         q_sf_patch             => waterstate_inst%q_sf_patch                   , & ! Output: [real(r8) (:)   ]  patch surface humidity (kg/kg)
+         rah1                   => frictionvel_inst%rah1_patch                  , & ! Output: [real(r8) (:)   ]  heat resistance (s/m)
+         raw1                   => frictionvel_inst%raw1_patch                  , & ! Output: [real(r8) (:)   ]  moisture resistance (s/m)
 #endif
 
          frac_h2osfc            => waterstate_inst%frac_h2osfc_col              , & ! Input:  [real(r8) (:)   ]  fraction of surface water                                             
@@ -784,7 +786,6 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
             ram1(p)  = 1._r8/(ustar(p)*ustar(p)/um(p))
             rah(p,1) = 1._r8/(temp1(p)*ustar(p))
             raw(p,1) = 1._r8/(temp2(p)*ustar(p))
-
             ! Bulk boundary layer resistance of leaves
 
             uaf(p) = um(p)*sqrt( 1._r8/(ram1(p)*um(p)) )
@@ -1183,7 +1184,9 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
 
 #ifdef COUP_OAS_ICON
          t_sf_patch(p)  = taf(p)
-!         q_sf(p)  = qaf(p)
+         q_sf_patch(p)  = qaf(p)
+         rah1(p) = rah(p,1)
+         raw1(p) = raw(p,1)
 #endif
          ! 2 m height air temperature
 

src/clm5/biogeophys/FrictionVelocityMod.F90

@@ -45,6 +45,10 @@ module FrictionVelocityMod
      real(r8), pointer, public :: rb1_patch        (:)   ! patch aerodynamical resistance (s/m) (for dry deposition of chemical tracers)
      real(r8), pointer, public :: rb10_patch       (:)   ! 10-day mean patch aerodynamical resistance (s/m) (for LUNA model)
      real(r8), pointer, public :: ram1_patch       (:)   ! patch aerodynamical resistance (s/m)
+#ifdef COUP_OAS_ICON
+     real(r8), pointer, public :: rah1_patch       (:)   ! patch heat resistance (s/m)
+     real(r8), pointer, public :: raw1_patch       (:)   ! patch moisture resistance (s/m)
+#endif
      real(r8), pointer, public :: z0m_patch        (:)   ! patch momentum roughness length (m)
      real(r8), pointer, public :: z0mv_patch       (:)   ! patch roughness length over vegetation, momentum [m]
      real(r8), pointer, public :: z0hv_patch       (:)   ! patch roughness length over vegetation, sensible heat [m]
@@ -122,6 +126,10 @@ subroutine InitAllocate(this, bounds)
     allocate(this%rb1_patch        (begp:endp)) ; this%rb1_patch        (:)   = nan
     allocate(this%rb10_patch       (begp:endp)) ; this%rb10_patch       (:)   = spval
     allocate(this%ram1_patch       (begp:endp)) ; this%ram1_patch       (:)   = nan
+#ifdef COUP_OAS_ICON
+    allocate(this%rah1_patch       (begp:endp)) ; this%rah1_patch       (:)   = nan
+    allocate(this%raw1_patch       (begp:endp)) ; this%raw1_patch       (:)   = nan
+#endif
     allocate(this%z0m_patch        (begp:endp)) ; this%z0m_patch        (:)   = nan
     allocate(this%z0mv_patch       (begp:endp)) ; this%z0mv_patch       (:)   = nan
     allocate(this%z0hv_patch       (begp:endp)) ; this%z0hv_patch       (:)   = nan
@@ -187,13 +195,31 @@ subroutine InitHistory(this, bounds)
          avgflag='A', long_name='10-m wind for dust model', &
          ptr_patch=this%u10_patch)
 
+#ifndef COUP_OAS_ICON
     if (use_cn) then
        this%ram1_patch(begp:endp) = spval
        call hist_addfld1d (fname='RAM1', units='s/m', &
             avgflag='A', long_name='aerodynamical resistance ', &
             ptr_patch=this%ram1_patch, default='inactive')
     end if
 
+#else
+    this%ram1_patch(begp:endp) = spval
+    call hist_addfld1d (fname='RAM1', units='s/m', &
+         avgflag='A', long_name='aerodynamical resistance ', &
+         ptr_patch=this%ram1_patch, default='inactive')
+
+    this%rah1_patch(begp:endp) = spval
+    call hist_addfld1d (fname='RAH1', units='s/m', &
+         avgflag='A', long_name='aerodynamical resistance for heat', &
+         ptr_patch=this%rah1_patch, default='inactive')
+
+    this%raw1_patch(begp:endp) = spval
+    call hist_addfld1d (fname='RAW1', units='s/m', &
+         avgflag='A', long_name='aerodynamical resistance for moisture', &
+         ptr_patch=this%raw1_patch, default='inactive')
+#endif
+
     if (use_cn) then
        this%fv_patch(begp:endp) = spval
        call hist_addfld1d (fname='FV', units='m/s', &

src/clm5/biogeophys/LakeFluxesMod.F90

@@ -201,7 +201,14 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
          zetamax          =>    frictionvel_parms_inst%zetamaxstable   , & ! Input:  [real(r8)       ]  max zeta value under stable conditions
          ram1             =>    frictionvel_inst%ram1_patch            , & ! Output: [real(r8) (:)   ]  aerodynamical resistance (s/m)                    
 
-         q_ref2m          =>    waterstate_inst%q_ref2m_patch          , & ! Output: [real(r8) (:)   ]  2 m height surface specific humidity (kg/kg)      
+#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)
+         rah1             =>    frictionvel_inst%rah1_patch            , & ! Output: [real(r8) (:)   ]  heat resistance (s/m)
+         raw1             =>    frictionvel_inst%raw1_patch            , & ! Output: [real(r8) (:)   ]  moisture resistance (s/m)
+#endif
+
+         q_ref2m          =>    waterstate_inst%q_ref2m_patch          , & ! Output: [real(r8) (:)   ]  2 m height surface specific humidity (kg/kg)
          rh_ref2m         =>    waterstate_inst%rh_ref2m_patch         , & ! Output: [real(r8) (:)   ]  2 m height surface relative humidity (%)          
 
          tc_ref2m         =>    humanindex_inst%tc_ref2m_patch         , & ! Output: [real(r8) (:)]  2 m height surface air temperature (C)
@@ -615,6 +622,14 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
          eflx_lh_tot(p)   = htvp(c)*qflx_evap_soi(p)
          eflx_lh_grnd(p)  = htvp(c)*qflx_evap_soi(p)
 
+#ifdef COUP_OAS_ICON
+         t_sf_patch(p) = t_grnd(c)
+         ! invert the actual flux formula used above
+         q_sf_patch(p) = forc_q(c) + qflx_evap_soi(p)*raw(p)/forc_rho(c)
+         rah1(p) = rah(p)
+         raw1(p) = raw(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/UrbanFluxesMod.F90

@@ -221,6 +221,9 @@ subroutine UrbanFluxes (bounds, num_nourbanl, filter_nourbanl,
          taf                 =>   temperature_inst%taf_lun                  , & ! Output: [real(r8) (:)   ]  urban canopy air 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)
+         rah1                =>   frictionvel_inst%rah1_patch               , & ! Output: [real(r8) (:)   ]  heat resistance (s/m)
+         raw1                =>   frictionvel_inst%raw1_patch               , & ! Output: [real(r8) (:)   ]  moisture resistance (s/m)
 #endif
 
          tc_ref2m            => humanindex_inst%tc_ref2m_patch              , & ! Output: [real(r8) (:)   ]  2 m height surface air temperature (C)
@@ -769,6 +772,9 @@ subroutine UrbanFluxes (bounds, num_nourbanl, filter_nourbanl,
 
 #ifdef COUP_OAS_ICON
          t_sf_patch(p)  = taf(l)
+         q_sf_patch(p)  = qaf(l)
+         rah1(p) = rahu(l) !forc_rho(g) * cpair / cgrnds(p)
+         raw1(p) = rawu(l) !forc_rho(g) / qflx_evap_soi(p)
 #endif
       end do
 

src/clm5/biogeophys/WaterStateType.F90

@@ -121,6 +121,9 @@ module WaterstateType
      real(r8), pointer :: sno_liq_top_col        (:)   ! col snow liquid water fraction (mass), top layer  [fraction]
 
      real(r8), pointer :: q_ref2m_patch          (:)   ! patch 2 m height surface specific humidity (kg/kg)
+#ifdef COUP_OAS_ICON
+     real(r8), pointer :: q_sf_patch             (:)   ! patch surface humidity (kg/kg)
+#endif
      real(r8), pointer :: rh_ref2m_patch         (:)   ! patch 2 m height surface relative humidity (%)
      real(r8), pointer :: rh_ref2m_r_patch       (:)   ! patch 2 m height surface relative humidity - rural (%)
      real(r8), pointer :: rh_ref2m_u_patch       (:)   ! patch 2 m height surface relative humidity - urban (%)
@@ -312,6 +315,9 @@ subroutine InitAllocate(this, bounds)
     allocate(this%dqgdT_col              (begc:endc))                     ; this%dqgdT_col              (:)   = nan   
     allocate(this%qaf_lun                (begl:endl))                     ; this%qaf_lun                (:)   = nan
     allocate(this%q_ref2m_patch          (begp:endp))                     ; this%q_ref2m_patch          (:)   = nan
+#ifdef COUP_OAS_ICON
+    allocate(this%q_sf_patch             (begp:endp))                     ; this%q_sf_patch             (:)   = nan
+#endif
     allocate(this%rh_ref2m_patch         (begp:endp))                     ; this%rh_ref2m_patch         (:)   = nan
     allocate(this%rh_ref2m_u_patch       (begp:endp))                     ; this%rh_ref2m_u_patch       (:)   = nan
     allocate(this%rh_ref2m_r_patch       (begp:endp))                     ; this%rh_ref2m_r_patch       (:)   = nan
@@ -511,6 +517,13 @@ subroutine InitHistory(this, bounds)
          avgflag='A', long_name='2m specific humidity', &
          ptr_patch=this%q_ref2m_patch)
 
+#ifdef COUP_OAS_ICON
+    this%q_sf_patch(begp:endp) = spval
+    call hist_addfld1d (fname='QSF', units='kg/kg',  &
+         avgflag='A', long_name='surface specific humidity', &
+         ptr_patch=this%q_sf_patch)
+#endif
+
     this%rh_ref2m_patch(begp:endp) = spval
     call hist_addfld1d (fname='RH2M', units='%',  &
          avgflag='A', long_name='2m relative humidity', &

src/clm5/main/clm_varctl.F90

@@ -99,6 +99,12 @@ module clm_varctl
 
   logical, public :: ndep_from_cpl = .false.
 
+#ifdef COUP_OAS_ICON
+  !----------------------------------------------------------
+  ! Logical flag to enable exchange-coefficient coupling scheme
+  logical, public :: use_exchcoef_oascoupling = .false.
+#endif
+
   !----------------------------------------------------------
   ! Interpolation of finidat if requested
   !----------------------------------------------------------

src/clm5/main/controlMod.F90

@@ -259,6 +259,11 @@ subroutine control_init( )
     ! max number of plant functional types in naturally vegetated landunit
     namelist /clm_inparm/ maxpatch_pft
 
+#ifdef COUP_OAS_ICON
+    ! Enable exchange-coefficient coupling scheme when using OASIS ECLM coupling
+    namelist /clm_inparm/ use_exchcoef_oascoupling
+#endif
+
     ! flag for SSRE diagnostic
     namelist /clm_inparm/ use_SSRE
 
@@ -600,6 +605,9 @@ subroutine control_spmd()
     call mpi_bcast (use_mexicocity, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_noio, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_SSRE, 1, MPI_LOGICAL, 0, mpicom, ier)
+#ifdef COUP_OAS_ICON
+    call mpi_bcast (use_exchcoef_oascoupling, 1, MPI_LOGICAL, 0, mpicom, ier)
+#endif
 
     ! initial file variables
     call mpi_bcast (nrevsn, len(nrevsn), MPI_CHARACTER, 0, mpicom, ier)

src/clm5/main/lnd2atmMod.F90

@@ -199,6 +199,16 @@ subroutine lnd2atm(bounds, &
     REAL, allocatable :: h2ocan_grc(:)
     logical, allocatable :: found(:)
     logical, allocatable :: found_patch(:)
+#endif
+#ifdef COUP_OAS_ICON
+    ! patch-level aerodynamic conductances (1/resistance) for rah1/raw1/ram1: 
+    ! combine parallel resistances by area-weighting their conductances
+    integer  :: pp
+    real(r8) :: rah1_cond_patch(bounds%begp:bounds%endp)
+    real(r8) :: raw1_cond_patch(bounds%begp:bounds%endp)
+    real(r8) :: ram1_cond_patch(bounds%begp:bounds%endp)
+    ! resistance where no land patch contribute
+    real(r8), parameter :: rmax_oascoupling = 1.e6_r8
 #endif
     !------------------------------------------------------------------------
 
@@ -357,15 +367,55 @@ subroutine lnd2atm(bounds, &
          lnd2atm_inst%t_sf_grc       (bounds%begg:bounds%endg), &
          p2c_scale_type='unity', c2l_scale_type= 'unity', l2g_scale_type='unity')
 
-!    call p2g(bounds, &
-!         waterstate_inst%q_sf_patch (bounds%begp:bounds%endp), &
-!         lnd2atm_inst%q_sf_grc      (bounds%begg:bounds%endg), &
-!         p2c_scale_type='unity', c2l_scale_type= 'unity', l2g_scale_type='unity')
+    call p2g(bounds, &
+         waterstate_inst%q_sf_patch (bounds%begp:bounds%endp), &
+         lnd2atm_inst%q_sf_grc      (bounds%begg:bounds%endg), &
+         p2c_scale_type='unity', c2l_scale_type= 'unity', l2g_scale_type='unity')
+
+    ! combine parallel resistances by area-weighting conductances (1/R)
+    ! unset / nan or non-positive resistances stay spval and are skipped by p2g
+    rah1_cond_patch(bounds%begp:bounds%endp) = spval
+    raw1_cond_patch(bounds%begp:bounds%endp) = spval
+    ram1_cond_patch(bounds%begp:bounds%endp) = spval
+    do pp = bounds%begp, bounds%endp
+       if (frictionvel_inst%rah1_patch(pp) > 0._r8) rah1_cond_patch(pp) = 1._r8 / frictionvel_inst%rah1_patch(pp)
+       if (frictionvel_inst%raw1_patch(pp) > 0._r8) raw1_cond_patch(pp) = 1._r8 / frictionvel_inst%raw1_patch(pp)
+       if (frictionvel_inst%ram1_patch(pp) > 0._r8) ram1_cond_patch(pp) = 1._r8 / frictionvel_inst%ram1_patch(pp)
+    end do
+
+    call p2g(bounds, &
+         rah1_cond_patch       (bounds%begp:bounds%endp), &
+         lnd2atm_inst%rah1_grc (bounds%begg:bounds%endg), &
+         p2c_scale_type='unity', c2l_scale_type= 'unity', l2g_scale_type='unity')
+
+    call p2g(bounds, &
+         raw1_cond_patch       (bounds%begp:bounds%endp), &
+         lnd2atm_inst%raw1_grc (bounds%begg:bounds%endg), &
+         p2c_scale_type='unity', c2l_scale_type= 'unity', l2g_scale_type='unity')
+
+    call p2g(bounds, &
+         ram1_cond_patch           (bounds%begp:bounds%endp), &
+         lnd2atm_inst%ram1_cpl_grc (bounds%begg:bounds%endg), &
+         p2c_scale_type='unity', c2l_scale_type= 'unity', l2g_scale_type='unity')
 
-!    call p2g(bounds, &
-!         frictionvel_inst%rah1_patch (bounds%begp:bounds%endp), &
-!         lnd2atm_inst%rah1_grc       (bounds%begg:bounds%endg), &
-!         p2c_scale_type='unity', c2l_scale_type= 'unity', l2g_scale_type='unity')
+    ! back to resistancy after aggregation. Limit p2g for non-contributing patch
+    do g = bounds%begg, bounds%endg
+       if (lnd2atm_inst%rah1_grc(g) /= spval .and. lnd2atm_inst%rah1_grc(g) > 0._r8) then
+          lnd2atm_inst%rah1_grc(g)     = 1._r8 / lnd2atm_inst%rah1_grc(g)
+       else
+          lnd2atm_inst%rah1_grc(g)     = rmax_oascoupling
+       end if
+       if (lnd2atm_inst%raw1_grc(g) /= spval .and. lnd2atm_inst%raw1_grc(g) > 0._r8) then
+          lnd2atm_inst%raw1_grc(g)     = 1._r8 / lnd2atm_inst%raw1_grc(g)
+       else
+          lnd2atm_inst%raw1_grc(g)     = rmax_oascoupling
+       end if
+       if (lnd2atm_inst%ram1_cpl_grc(g) /= spval .and. lnd2atm_inst%ram1_cpl_grc(g) > 0._r8) then
+          lnd2atm_inst%ram1_cpl_grc(g) = 1._r8 / lnd2atm_inst%ram1_cpl_grc(g)
+       else
+          lnd2atm_inst%ram1_cpl_grc(g) = rmax_oascoupling
+       end if
+    end do
 
 #endif
 

src/clm5/main/lnd2atmType.F90

@@ -63,8 +63,10 @@ module lnd2atmType
      real(r8), pointer :: flux_ch4_grc       (:)   => null() ! net CH4 flux (kg C/m**2/s) [+ to atm]
 #ifdef COUP_OAS_ICON
      real(r8), pointer :: t_sf_grc           (:)   => null() ! surface temperature (Kelvin)
-!     real(r8), pointer :: q_sf_grc           (:)   => null() ! surface humidity (kg/kg)
-!     real(r8), pointer :: rah1_grc           (:)   => null() ! aerodynamical resistance for heat (s/m)
+     real(r8), pointer :: q_sf_grc           (:)   => null() ! surface humidity (kg/kg)
+     real(r8), pointer :: rah1_grc           (:)   => null() ! aerodynamical resistance for heat (s/m)
+     real(r8), pointer :: raw1_grc           (:)   => null() ! aerodynamical resistance for moisture (s/m)
+     real(r8), pointer :: ram1_cpl_grc       (:)   => null() ! aerodynamical resistance for momentum, conductance-weighted (s/m)
 #endif
      ! lnd->rof
      real(r8), pointer :: qflx_rofliq_grc         (:)   => null() ! rof liq forcing
@@ -178,8 +180,10 @@ subroutine InitAllocate(this, bounds)
     allocate(this%flux_ch4_grc       (begg:endg))            ; this%flux_ch4_grc       (:)   =ival
 #ifdef COUP_OAS_ICON
     allocate(this%t_sf_grc           (begg:endg))            ; this%t_sf_grc           (:)   =ival
-!    allocate(this%q_sf_grc           (begg:endg))            ; this%q_sf_grc           (:)   =ival
-!    allocate(this%rah1_grc           (begg:endg))            ; this%rah1_grc           (:)   =ival
+    allocate(this%q_sf_grc           (begg:endg))            ; this%q_sf_grc           (:)   =ival
+    allocate(this%rah1_grc           (begg:endg))            ; this%rah1_grc           (:)   =ival
+    allocate(this%raw1_grc           (begg:endg))            ; this%raw1_grc           (:)   =ival
+    allocate(this%ram1_cpl_grc       (begg:endg))            ; this%ram1_cpl_grc       (:)   =ival
 #endif
     allocate(this%qflx_rofliq_grc    (begg:endg))            ; this%qflx_rofliq_grc    (:)   =ival
     allocate(this%qflx_rofliq_qsur_grc    (begg:endg))       ; this%qflx_rofliq_qsur_grc    (:)   =ival