From 9a66a9eed172d0f8ee4ec77e2cf77f60cf789ff8 Mon Sep 17 00:00:00 2001
From: Helin Wei <Helin.Wei@noaa.gov>
Date: Wed, 30 Nov 2022 15:36:51 -0500
Subject: [PATCH] using NoahMP tiled tq2m with the bug fix

---
 physics/module_sf_noahmp_glacier.f90 |  5 +++++
 physics/module_sf_noahmplsm.f90      | 24 ++++++++++++++++++++++--
 physics/sfc_diag_post.F90            | 16 ++++++++--------
 3 files changed, 35 insertions(+), 10 deletions(-)

diff --git a/physics/module_sf_noahmp_glacier.f90 b/physics/module_sf_noahmp_glacier.f90
index 08fdd92ce..738a40b5c 100644
--- a/physics/module_sf_noahmp_glacier.f90
+++ b/physics/module_sf_noahmp_glacier.f90
@@ -1464,6 +1464,11 @@ subroutine glacier_flux (nsoil     ,nsnow   ,emg     ,isnow   ,df      ,dzsnso
 ! 2m air temperature
      ehb2  = fv*vkc/(log((2.+z0h)/z0h)-fh2)
      cq2b  = ehb2
+! for opt_sfc 3
+     if (opt_sfc ==3) then
+       ehb2  = fv*vkc/fh2
+       cq2b  = ehb2
+     endif
 
      if (opt_sfc == 4) then
        ehb2 = ch2 * wspd1i ! need conductance,z0h from sfcdif4 
diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90
index 46b572b65..1da30f156 100644
--- a/physics/module_sf_noahmplsm.f90
+++ b/physics/module_sf_noahmplsm.f90
@@ -4289,13 +4289,19 @@ subroutine vege_flux(parameters,nsnow   ,nsoil   ,isnow   ,vegtyp  ,veg     , &
 !     qfx = (qsfc-qair)*rhoair*caw !*cpair/gammag
 
 ! 2m temperature over vegetation ( corrected for low cq2v values )
-   if (opt_sfc == 1 .or. opt_sfc == 2 .or. opt_sfc ==3 ) then
+   if (opt_sfc == 1 .or. opt_sfc == 2 ) then
 !      cah2 = fv*1./vkc*log((2.+z0h)/z0h)
       cah2 = fv*vkc/log((2.+z0h)/z0h)
       cah2 = fv*vkc/(log((2.+z0h)/z0h)-fh2)
       cq2v = cah2
    endif
 
+! opt_sfc 3: fh2 is the stability
+    if (opt_sfc ==3) then
+        cah2 = fv*vkc/fh2
+        cq2v = cah2
+    endif
+
     if (opt_sfc == 4 ) then
        rahc2 = max(1.,1./(ch2v*wspdv))
        rawc2 = rahc2
@@ -4742,7 +4748,7 @@ subroutine bare_flux (parameters,nsnow   ,nsoil   ,isnow   ,dt      ,sag     , &
 !jref:start; errors in original equation corrected.
 ! 2m air temperature
 
-     if(opt_sfc == 1 .or. opt_sfc ==2 .or. opt_sfc == 3) then
+     if(opt_sfc == 1 .or. opt_sfc ==2 ) then
        ehb2  = fv*vkc/log((2.+z0h)/z0h)
        ehb2  = fv*vkc/(log((2.+z0h)/z0h)-fh2)
        cq2b  = ehb2
@@ -4756,6 +4762,20 @@ subroutine bare_flux (parameters,nsnow   ,nsoil   ,isnow   ,dt      ,sag     , &
        if (parameters%urban_flag) q2b = qsfc
      end if
 
+! opt_sfc 3: fh2 is the stability
+     if(opt_sfc == 3 ) then
+       ehb2  = fv*vkc/fh2
+       cq2b  = ehb2
+       if (ehb2.lt.1.e-5 ) then
+         t2mb  = tgb
+         q2b   = qsfc
+       else
+         t2mb  = tgb - shb/(rhoair*cpair) * 1./ehb2
+         q2b   = qsfc - evb/(lathea*rhoair)*(1./cq2b + rsurf)
+       endif
+       if (parameters%urban_flag) q2b = qsfc
+     end if
+
     if(opt_sfc == 4) then ! consistent with veg
 
          rahb2 = max(1.,1./(ch2b*wspdb))
diff --git a/physics/sfc_diag_post.F90 b/physics/sfc_diag_post.F90
index 18936cc8b..139f30137 100644
--- a/physics/sfc_diag_post.F90
+++ b/physics/sfc_diag_post.F90
@@ -40,14 +40,14 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, dry, lssav, dtf, con_eps, con
         errmsg = ''
         errflg = 0
 
-!       if (lsm == lsm_noahmp) then
-!         do i=1,im
-!           if(dry(i)) then
-!             t2m(i) = t2mmp(i)
-!             q2m(i) = q2mp(i)
-!           endif
-!         enddo
-!       endif
+        if (lsm == lsm_noahmp) then
+          do i=1,im
+            if(dry(i)) then
+              t2m(i) = t2mmp(i)
+              q2m(i) = q2mp(i)
+            endif
+          enddo
+        endif
 
         if (lssav) then
           do i=1,im