From 5b3ec4b1eca8cc335d894e7033410fa44a2aebfc Mon Sep 17 00:00:00 2001
From: weiwangncar <weiwang@ucar.edu>
Date: Wed, 9 Aug 2017 13:52:37 -0600
Subject: [PATCH] Fix to make consistent calculation of PSIT and PSIQ when
 COARE formula is used in the revised surface layer scheme (#271)

TYPE: Bug fix

KEYWORDS: integrated similarity functions, PSIT, PSIQ

SOURCE: Internal

DESCRIPTION OF CHANGES:

When COARE 3 formula was introduced to the revised surface layer scheme (sf_sfclay_physics=1), the modified thermal roughness was not used consistently with the new formulation of the PSIT, PSIQ calculations. This has a small effect in an idealized tropical cyclone test and a real-data hurricane test in terms of track and intensity.

LIST OF MODIFIED FILES:

modified: phys/module_sf_sfclayrev.F

TESTS CONDUCTED: idealized tropical cyclone case at 3km and a real-data typhoon case. Reg tests conducted by Mike.

(cherry picked from commit 5fab1111d69277cd62a0e78315aebde862e8f341)
---
 phys/module_sf_sfclayrev.F | 53 ++++++++++++++++++++++++++++++++++----
 1 file changed, 48 insertions(+), 5 deletions(-)

diff --git a/phys/module_sf_sfclayrev.F b/phys/module_sf_sfclayrev.F
index 2a595ef3d5..9c36abb62f 100644
--- a/phys/module_sf_sfclayrev.F
+++ b/phys/module_sf_sfclayrev.F
@@ -708,11 +708,54 @@ SUBROUTINE SFCLAYREV1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d,                &
               Z0T = MAX(Z0T,2.0e-9)
               Z0Q = Z0T
 
-              PSIQ=max(ALOG((ZA(I)+Z0Q)/Z0Q)-PSIH(I), 2.)
-              PSIT=max(ALOG((ZA(I)+Z0T)/Z0T)-PSIH(I), 2.)
-              PSIQ2=max(ALOG((2.+Z0Q)/Z0Q)-PSIH2(I), 2.)
-              PSIT2=max(ALOG((2.+Z0T)/Z0T)-PSIH2(I), 2.)
-              PSIQ10=max(ALOG((10.+Z0Q)/Z0Q)-PSIH10(I), 2.)
+! following paj:
+           zolzz=zol(I)*(za(I)+z0t)/za(I)    ! (z+z0t)/L
+           zol10=zol(I)*(10.+z0t)/za(I)   ! (10+z0t)/L
+           zol2=zol(I)*(2.+z0t)/za(I)     ! (2+z0t)/L
+           zol0=zol(I)*z0t/za(I)          ! z0t/L
+!
+              if (zol(I).gt.0.) then
+              psih(I)=psih_stable(zolzz)-psih_stable(zol0)
+              psih10(I)=psih_stable(zol10)-psih_stable(zol0)
+              psih2(I)=psih_stable(zol2)-psih_stable(zol0)
+              else
+                if (zol(I).eq.0) then
+                psih(I)=0.
+                psih10(I)=0.
+                psih2(I)=0.
+                else
+                psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
+                psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
+                psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
+                endif
+              endif
+              PSIT=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I)
+              PSIT2=ALOG((2.+z0t)/Z0t)-PSIH2(I)
+
+           zolzz=zol(I)*(za(I)+z0q)/za(I)    ! (z+z0q)/L
+           zol10=zol(I)*(10.+z0q)/za(I)   ! (10+z0q)/L
+           zol2=zol(I)*(2.+z0q)/za(I)     ! (2+z0q)/L
+           zol0=zol(I)*z0q/za(I)          ! z0q/L
+!
+              if (zol(I).gt.0.) then
+              psih(I)=psih_stable(zolzz)-psih_stable(zol0)
+              psih10(I)=psih_stable(zol10)-psih_stable(zol0)
+              psih2(I)=psih_stable(zol2)-psih_stable(zol0)
+              else
+                if (zol(I).eq.0) then
+                psih(I)=0.
+                psih10(I)=0.
+                psih2(I)=0.
+                else
+                psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
+                psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
+                psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
+                endif
+              endif
+!
+              PSIQ=ALOG((ZA(I)+z0q)/Z0q)-PSIH(I)
+              PSIQ2=ALOG((2.+z0q)/Z0q)-PSIH2(I)
+              PSIQ10=ALOG((10.+z0q)/Z0q)-PSIH10(I)
         ENDIF
 
         IF ( PRESENT(ISFTCFLX) ) THEN