Merge pull request NCAR#25 from HelinWei-NOAA/fix_tq2m_diagnoahmp

using NoahMP tiled tq2m with the bug fix
dustinswales · Dec 2, 2022 · 9d45cd1 · 9d45cd1
2 parents 18d3fda + 9a66a9e
commit 9d45cd1
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 10 deletions.
diff --git 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
@@ -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
@@ -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