Correction of stability functions for the tke-epsilon-tpe PBL scheme (#2120)

TYPE: bug fix

KEYWORDS: tke-epsilon-tpe, PBL, stability functions

SOURCE: Andrea Zonato, Royal Netherlands Meteorological Institute (KNMI)

DESCRIPTION OF CHANGES:
Problem:
The stability functions at the surface (phim,phieps) for calculating the boundary values of tke and its dissipation rate in the tke-epsilon-tpe PBL scheme [Zonato et al., 2022](https://doi.org/10.1175/MWR-D-21-0299.1) does not include the correction by z0.

Solution:
Pass the correct stability functions, calculated in module_sf_sfclayrev.F

LIST OF MODIFIED FILES:
phys/module_bl_keps.F
phys/module_pbl_driver.F

TESTS CONDUCTED: 
The Jenkins tests are all passing.

RELEASE NOTE: 
The similarity stability functions phim and phieps, necessary for calculating the surface values of tke and dissipation rate in the tke-epsilon-tpe PBL scheme [Zonato et al., 2022](https://doi.org/10.1175/MWR-D-21-0299.1) have been updated considering the correction term accounting for the roughness length z0. No relevant differences are found in temperature, wind speed, and humidity. Regarding turbulence variables, the stable case has just negligible differences, while the unstable case shows higher values of TKE and dissipation rate and lower values of temperature variance.

Author: andreazonato

phys/module_bl_keps.F

@@ -1073,16 +1073,16 @@ subroutine surface_bl_pra_ri(kms,kme,kts,kte,dz,z,rho,g,cp,z0,sflux,raten,pi1d,p
      zz=(zl1+z0)/zl1
      zol1 = max(b_ric*psim*psim/psih,rimin)
     if(sfcflg)then
-     zol1 = min(zol1,-zfmin)
+     zol1 = zz*min(zol1,-zfmin)
     else
-     zol1 = max(zol1,zfmin)
+     zol1 = zz*max(zol1,zfmin)
     endif
      hol1 = zol1*pblh/zl1*sfcfrac
 
     if(sfcflg)then
-       phim = (1.-aphi16*zol1*zz)**(-1./4.)
-       phih = (1.-aphi16*zol1*zz)**(-1./2.)
-       phieps=1.-zol1*zz
+       phim = (1.-aphi16*zol1)**(-1./4.)
+       phih = (1.-aphi16*zol1)**(-1./2.)
+       phieps=1.-zol1
        phim_sl = (1.-aphi16*hol1)**(-1./4.)
        phih_sl = (1.-aphi16*hol1)**(-1./2.)
        bfx0 = max(sflux,0.)
@@ -1093,9 +1093,9 @@ subroutine surface_bl_pra_ri(kms,kme,kts,kte,dz,z,rho,g,cp,z0,sflux,raten,pi1d,p
      else
        phim_sl = (1.+aphi5*hol1)
        phih_sl= phim_sl
-       phim=(1.+aphi5*zol1*zz)
-       phieps=(1+2.5*(zol1*zz)**0.6)**(3./2.)
-       phih= phim_sl
+       phim=(1.+aphi5*zol1)
+       phieps=(1+2.5*(zol1)**0.6)**(3./2.)
+       phih= phim
        wstar3=0.
        wstar3_2=0.
      endif

phys/module_pbl_driver.F

@@ -1895,7 +1895,7 @@ SUBROUTINE pbl_driver(                                          &
               ,B_E_BEP=b_e_bep                                                 &
               ,SF_BEP=sf_bep,VL_BEP=vl_bep                                     &
               ,BR=br,ZNT=znt                                                   &
-              ,PSIM=psim,PSIH=psih                                             &
+              ,PSIM=fm,PSIH=fhh                                             &
               ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde                 &
               ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme                 &
               ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte  )