Updates for component land model

CODEOWNERS

@@ -196,6 +196,8 @@ physics/ysuvdif.*                              @Qingfu-Liu @WeiguoWang-NOAA
 physics/zhaocarr_gscond.*                      @RuiyuSun                            @grantfirl @Qingfu-Liu @dustinswales
 physics/zhaocarr_precpd.*                      @RuiyuSun                            @grantfirl @Qingfu-Liu @dustinswales
 
+physics/sfc_land.*                             @uturuncoglu @barlage
+
 ########################################################################
 
 # Lines starting with '#' are comments.

physics/CONV/Grell_Freitas/cu_gf_deep.F90

@@ -142,13 +142,13 @@ subroutine cu_gf_deep_run(        &
                                 !! betwee -1 and +1
               ,do_capsuppress,cap_suppress_j    &    !         
               ,k22                              &    !
-              ,jmin,tropics)                         !
+              ,jmin,kdt,tropics)                         !
 
    implicit none
 
      integer                                                &
         ,intent (in   )                   ::                &
-        nranflag,itf,ktf,its,ite, kts,kte,ipr,imid
+        nranflag,itf,ktf,its,ite, kts,kte,ipr,imid,kdt
      integer, intent (in   )              ::                &
         ichoice,nchem
      real(kind=kind_phys),  dimension (its:ite,4)                 &
@@ -591,6 +591,7 @@ subroutine cu_gf_deep_run(        &
          sig(i)=(1.-frh)**2
          !frh_out(i) = frh
          if(forcing(i,7).eq.0.)sig(i)=1.
+         if(kdt.le.(3600./dtime))sig(i)=1.
          frh_out(i) = frh*sig(i)
       enddo
 !$acc end kernels

physics/CONV/Grell_Freitas/cu_gf_driver.F90

@@ -68,7 +68,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
                dfi_radar_max_intervals,ldiag3d,qci_conv,do_cap_suppress,        &
                maxupmf,maxMF,do_mynnedmf,ichoice_in,ichoicem_in,ichoice_s_in,   &
                spp_cu_deep,spp_wts_cu_deep,nchem,chem3d,fscav,wetdpc_deep,      &
-               do_smoke_transport,errmsg,errflg)
+               do_smoke_transport,kdt,errmsg,errflg)
 !-------------------------------------------------------------
       implicit none
       integer, parameter :: maxiens=1
@@ -95,7 +95,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
       integer            :: ishallow_g3 ! depend on imfshalcnv
 !-------------------------------------------------------------
    integer      :: its,ite, jts,jte, kts,kte
-   integer, intent(in   ) :: im,km,ntracer, nchem
+   integer, intent(in   ) :: im,km,ntracer,nchem,kdt
    integer, intent(in   ) :: ichoice_in,ichoicem_in,ichoice_s_in
    logical, intent(in   ) :: flag_init, flag_restart, do_mynnedmf
    logical, intent(in   ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend
@@ -766,7 +766,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
                                ! betwee -1 and +1
               ,do_cap_suppress_here,cap_suppress_j &
               ,k22m          &
-              ,jminm,tropics)
+              ,jminm,kdt,tropics)
 !$acc kernels
             do i=its,itf
              do k=kts,ktf
@@ -853,7 +853,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
                                ! betwee -1 and +1
               ,do_cap_suppress_here,cap_suppress_j &
               ,k22          &
-              ,jmin,tropics)
+              ,jmin,kdt,tropics)
           jpr=0
           ipr=0
 !$acc kernels

physics/CONV/Grell_Freitas/cu_gf_driver.meta

@@ -651,6 +651,13 @@
   type = real
   kind = kind_phys
   intent = inout
+[kdt]
+  standard_name = index_of_timestep
+  long_name = current forecast iteration
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP

physics/PBL/MYNN_EDMF/module_bl_mynn.F90

@@ -2001,9 +2001,9 @@ SUBROUTINE  mym_length (                     &
         uonset= 15. 
         wt_u = (1.0 - min(max(ugrid - uonset, 0.0)/30.0, 0.5)) 
         cns  = 2.7 !was 3.5
-        alp1 = 0.22
+        alp1 = 0.23
         alp2 = 0.3
-        alp3 = 2.0 * wt_u !taper off bouyancy enhancement in shear-driven pbls
+        alp3 = 2.5 * wt_u !taper off bouyancy enhancement in shear-driven pbls
         alp4 = 5.0
         alp5 = 0.3
         alp6 = 50.
@@ -2059,12 +2059,12 @@ SUBROUTINE  mym_length (                     &
 
            !   **  Length scale limited by the buoyancy effect  **
            IF ( dtv(k) .GT. 0.0 ) THEN
-              bv  = max( sqrt( gtr*dtv(k) ), 0.001)
+              bv  = max( sqrt( gtr*dtv(k) ), 0.0001)
               elb = MAX(alp2*qkw(k),                          &
                   &    alp6*edmf_a1(k-1)*edmf_w1(k-1)) / bv   &
                   &  *( 1.0 + alp3*SQRT( vsc/(bv*elt) ) )
               elb = MIN(elb, zwk)
-              elf = 0.80 * qkw(k)/bv
+              elf = 1.0 * qkw(k)/bv
               elBLavg(k) = MAX(elBLavg(k), alp6*edmf_a1(k-1)*edmf_w1(k-1)/bv)
            ELSE
               elb = 1.0e10
@@ -2084,8 +2084,10 @@ SUBROUTINE  mym_length (                     &
            !add blending to use BouLac mixing length in free atmos;
            !defined relative to the PBLH (zi) + transition layer (h1)
            !el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf)
-           !try squared-blending
-           el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2)))
+           !try squared-blending - but take out elb (makes it underdiffusive)
+           !el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2)))
+           el(k) = sqrt( els**2/(1. + (els**2/elt**2)))
+           el(k) = min(el(k), elb)
            el(k) = MIN (el(k), elf)
            el(k) = el(k)*(1.-wt) + alp5*elBLavg(k)*wt
 
@@ -3633,13 +3635,13 @@ SUBROUTINE  mym_condensation (kts,kte,   &
 
     real(kind_phys):: qsl,esat,qsat,dqsl,cld0,q1k,qlk,eq1,qll,           &
          &q2p,pt,rac,qt,t,xl,rsl,cpm,Fng,qww,alpha,beta,bb,              &
-         &ls,wt,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water,       &
+         &ls,wt,wt2,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water,   &
          &qmq,qsat_tk,q1_rh,rh_hack,dzm1,zsl,maxqc
     real(kind_phys), parameter :: qpct_sfc=0.025
     real(kind_phys), parameter :: qpct_pbl=0.030
     real(kind_phys), parameter :: qpct_trp=0.040
     real(kind_phys), parameter :: rhcrit  =0.83 !for cloudpdf = 2
-    real(kind_phys), parameter :: rhmax   =1.01 !for cloudpdf = 2
+    real(kind_phys), parameter :: rhmax   =1.02 !for cloudpdf = 2
     integer :: i,j,k
 
     real(kind_phys):: erf
@@ -3864,25 +3866,18 @@ SUBROUTINE  mym_condensation (kts,kte,   &
            !Add condition for falling/settling into low-RH layers, so at least
            !some cloud fraction is applied for all qc, qs, and qi.
            rh_hack= rh(k)
+           wt2    = min(max( zagl - pblh2, 0.0 )/300., 1.0)
            !ensure adequate RH & q1 when qi is at least 1e-9 (above the PBLH)
-           if (qi(k)>1.e-9 .and. zagl .gt. pblh2) then
-              rh_hack =min(rhmax, rhcrit + 0.07*(9.0 + log10(qi(k))))
+           if ((qi(k)+qs(k))>1.e-9 .and. (zagl .gt. pblh2)) then
+              rh_hack =min(rhmax, rhcrit + wt2*0.045*(9.0 + log10(qi(k)+qs(k))))
               rh(k)   =max(rh(k), rh_hack)
               !add rh-based q1
               q1_rh   =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
               q1(k)   =max(q1_rh, q1(k) )
            endif
-           !ensure adequate RH & q1 when qc is at least 1e-6
-           if (qc(k)>1.e-6) then
-              rh_hack =min(rhmax, rhcrit + 0.09*(6.0 + log10(qc(k))))
-              rh(k)   =max(rh(k), rh_hack)
-              !add rh-based q1
-              q1_rh   =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
-              q1(k)   =max(q1_rh, q1(k) )
-           endif
-           !ensure adequate RH & q1 when qs is at least 1e-8 (above the PBLH)
-           if (qs(k)>1.e-8 .and. zagl .gt. pblh2) then
-              rh_hack =min(rhmax, rhcrit + 0.07*(8.0 + log10(qs(k))))
+           !ensure adequate rh & q1 when qc is at least 1e-6 (above the PBLH)
+           if (qc(k)>1.e-6 .and. (zagl .gt. pblh2)) then
+              rh_hack =min(rhmax, rhcrit + wt2*0.08*(6.0 + log10(qc(k))))
               rh(k)   =max(rh(k), rh_hack)
               !add rh-based q1
               q1_rh   =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
@@ -3994,7 +3989,7 @@ SUBROUTINE  mym_condensation (kts,kte,   &
            fac_damp = min(zagl * 0.0025, 1.0)
            !cld_factor = 1.0 + fac_damp*MAX(0.0, ( RH(k) - 0.75 ) / 0.26 )**1.9 !HRRRv4
            !cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.25 )**2, 0.3)
-           cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.35)
+           cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.37)
            cldfra_bl1D(K) = min( 1., cld_factor*cldfra_bl1D(K) )
         enddo
 
@@ -4181,38 +4176,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i,       &
 
     k=kts
 
-!original approach (drag in b-vector):
-!    a(1)=0.
-!    b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff
-!    c(1)=-dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff
-!    d(1)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff + &
-!         sub_u(k)*delt + det_u(k)*delt
-
 !rho-weighted (drag in b-vector):
     a(k)=  -dtz(k)*kmdz(k)*rhoinv(k)
-    b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) &
-           & - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
+    b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k)     &
+           & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff             &
+           & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
     c(k)=  -dtz(k)*kmdz(k+1)*rhoinv(k) &
-           & - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
-    d(k)=u(k)  + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff - &
-       & dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt
-
-!rho-weighted with drag term moved out of b-array
-!    a(k)=  -dtz(k)*kmdz(k)*rhoinv(k)
-!    b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
-!    c(k)=  -dtz(k)*kmdz(k+1)*rhoinv(k)   - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
-!    d(k)=u(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*uoce*ust**2/wspd - &
-!    !!!d(k)=u(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*uoce*ust**2/wspd - &
-!      &  dtz(k)*rhoinv(k)*s_awu(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt
-
-    DO k=kts+1,kte-1
-       a(k)=  -dtz(k)*kmdz(k)*rhoinv(k)     + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff 
-       b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + &
-           &    0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
-       c(k)=  -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
-       d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff + &
-           &    sub_u(k)*delt + det_u(k)*delt
-    ENDDO
+           & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff             &
+           & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
+    d(k)=u(k)  + dtz(k)*uoce*ust**2/wspd                        &
+           & - dtz(k)*rhoinv(k)*s_awu(k+1)*onoff                &
+           & + dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff               &
+           & + sub_u(k)*delt + det_u(k)*delt
+
+    do k=kts+1,kte-1
+       a(k)=  -dtz(k)*kmdz(k)*rhoinv(k)                         &
+           &  + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff              &
+           &  + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
+       b(k)=1.+ dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k)            &
+           &  + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff  &
+           &  + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
+       c(k)=  - dtz(k)*kmdz(k+1)*rhoinv(k)                      &
+           &  - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff            &
+           &  - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
+       d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff &
+           &  - dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff  &
+           &  + sub_u(k)*delt + det_u(k)*delt
+    enddo
 
 !! no flux at the top
 !    a(kte)=-1.
@@ -4247,37 +4237,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i,       &
 
     k=kts
 
-!original approach (drag in b-vector):
-!    a(1)=0.
-!    b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff
-!    c(1)=   - dtz(k)*dfm(k+1)               - 0.5*dtz(k)*s_aw(k+1)*onoff
-!    d(1)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff + &
-!          sub_v(k)*delt + det_v(k)*delt
-
 !rho-weighted (drag in b-vector):
     a(k)=  -dtz(k)*kmdz(k)*rhoinv(k)
-    b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) &
-        &  - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
-    c(k)=  -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
-    d(k)=v(k)  + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + &
-       & sub_v(k)*delt + det_v(k)*delt
-
-!rho-weighted with drag	term moved out of b-array
-!    a(k)=  -dtz(k)*kmdz(k)*rhoinv(k)
-!    b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k)  - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
-!    c(k)=  -dtz(k)*kmdz(k+1)*rhoinv(k)    - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
-!    d(k)=v(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*voce*ust**2/wspd - &
-!    !!!d(k)=v(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*voce*ust**2/wspd - &
-!      &  dtz(k)*rhoinv(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + sub_v(k)*delt + det_v(k)*delt
-
-    DO k=kts+1,kte-1
-       a(k)=  -dtz(k)*kmdz(k)*rhoinv(k)   + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
-       b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + &
-           &    0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
-       c(k)=  -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff 
-       d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff + &
-           &    sub_v(k)*delt + det_v(k)*delt
-    ENDDO
+    b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k)   &
+        &  - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff               &
+        &  - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
+    c(k)=  -dtz(k)*kmdz(k+1)*rhoinv(k)                          &
+        &  - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff               &
+        &  - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
+    d(k)=v(k)  + dtz(k)*voce*ust**2/wspd                        &
+        &  - dtz(k)*rhoinv(k)*s_awv(k+1)*onoff                  &
+        &  + dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff                 &
+        &  + sub_v(k)*delt + det_v(k)*delt
+
+    do k=kts+1,kte-1
+       a(k)=  -dtz(k)*kmdz(k)*rhoinv(k)                         &
+         & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff                 &
+         & + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
+       b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k)             &
+         & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff     &
+         & + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
+       c(k)=  -dtz(k)*kmdz(k+1)*rhoinv(k)                       &
+         & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff               &
+         & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
+       d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff &
+         & - dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff     &
+         & + sub_v(k)*delt + det_v(k)*delt
+    enddo
 
 !! no flux at the top
 !    a(kte)=-1.

physics/SFC_Models/Land/Noahmp/noahmpdrv.F90

@@ -136,7 +136,7 @@ subroutine noahmpdrv_run                                       &
       iopt_trs,iopt_diag,xlatin, xcoszin, iyrlen, julian, garea, &
       rainn_mp, rainc_mp, snow_mp, graupel_mp, ice_mp, rhonewsn1,&
       con_hvap, con_cp, con_jcal, rhoh2o, con_eps, con_epsm1,    &
-      con_fvirt, con_rd, con_hfus, thsfc_loc,                    &
+      con_fvirt, con_rd, con_hfus, thsfc_loc, cpllnd, cpllnd2atm,&
 
 !  ---  in/outs:
       weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc,        &
@@ -310,6 +310,9 @@ subroutine noahmpdrv_run                                       &
 
   logical                                , intent(in)    :: thsfc_loc  ! Flag for reference pressure in theta calculation
 
+  logical                                , intent(in)    :: cpllnd     ! Flag for land coupling (atm->lnd)
+  logical                                , intent(in)    :: cpllnd2atm ! Flag for land coupling (lnd->atm)
+
   real(kind=kind_phys), dimension(:)     , intent(inout) :: weasd      ! water equivalent accumulated snow depth [mm]
   real(kind=kind_phys), dimension(:)     , intent(inout) :: snwdph     ! snow depth [mm]
   real(kind=kind_phys), dimension(:)     , intent(inout) :: tskin      ! ground surface skin temperature [K]
@@ -684,7 +687,12 @@ subroutine noahmpdrv_run                                       &
   errmsg = ''
   errflg = 0
 
-do i = 1, im
+!
+!  --- Just return if external land component is activated for two-way interaction
+!
+  if (cpllnd .and. cpllnd2atm) return
+
+  do i = 1, im
 
     if (flag_iter(i) .and. dry(i)) then
 

physics/SFC_Models/Land/Noahmp/noahmpdrv.meta

@@ -635,6 +635,20 @@
   dimensions = ()
   type = logical
   intent = in
+[cpllnd]
+  standard_name = flag_for_land_coupling
+  long_name = flag controlling cpllnd collection (default off)
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[cpllnd2atm]
+  standard_name = flag_for_one_way_land_coupling_to_atmosphere
+  long_name = flag controlling land coupling to the atmosphere (default off)
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
 [weasd]
   standard_name = water_equivalent_accumulated_snow_depth_over_land
   long_name = water equiv of acc snow depth over land

physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90

@@ -1687,7 +1687,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
       endif
 
        if(newsn > zero ) then
-         SNOWFRACnewsn=MIN(one,SNHEI/SNHEI_CRIT_newsn)
+         SNOWFRACnewsn=MIN(one,snowfallac*1.e-3_kind_phys/SNHEI_CRIT_newsn)
        endif
 
        !-- due to steep slopes and blown snow, limit snow fraction in the
@@ -1700,7 +1700,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
        if(snowfrac < 0.75_kind_phys) snow_mosaic = one
 
        KEEP_SNOW_ALBEDO = zero
-       IF (NEWSN > zero .and. snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys) THEN
+       IF (snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys) THEN
        ! new snow
              KEEP_SNOW_ALBEDO = one
              ! turn off separate treatment of snow covered and snow-free portions of the grid cell
@@ -1735,7 +1735,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
 ! hwlps with these biases.. 
      if( snow_mosaic == one) then
          ALBsn=alb_snow
-         if(newsn > zero .and. KEEP_SNOW_ALBEDO > 0.9_kind_phys .and. albsn < 0.4_kind_phys) then
+         if(KEEP_SNOW_ALBEDO > 0.9_kind_phys .and. albsn < 0.4_kind_phys) then
          !-- Albedo correction with fresh snow and deep snow pack
          !-- will reduce warm bias in western Canada
          !-- and US West coast, where max snow albedo is low (0.3-0.5).