Fix to make consistent calculation of PSIT and PSIQ when COARE formul…

…a 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 5fab111)

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