phys/module_bl_gbmpbl.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3 ! Grenier-Bretherton mixing PBL scheme (GBM PBL):
   4 !
   5 ! REFERENCES
   6 ! 1) Grenier, H., and C. S. Bretherton, 2001: A moist PBL parameterization
   7 !    for large-scale models and its application to subtropical cloud-topped
   8 !    marine boundary layers. Mon. Wea. Rev., 129, 357-377.
   9 !    ("GB01" in the comments)
  10 ! 2) Galperin, B., L. H. Kantha, S. Hassid, and A. Rosati, 1988:
  11 !    A quasi-equilibrium turbulent energy model for geophysical flows.
  12 !    J. Atmos. Sci., 45, 55-62.
  13 !    ("Gal88" in the comments)
  14 ! 3) Mellor, G., and T. Yamada, 1982: Development of a turbulence closure
  15 !    model for geophysical fluid problems. Rev. Astrophys.Space Phys.,
  16 !    20, 851-875.
  17 !    ("MY82" in the comments)
  18 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  19 !! Natalie Perlin, October 2012 nperlin@coas.oregonstate.edu
  20 !!
  21 !! GBM PBL Code history:
  22 !!
  23 !! Herve Grenier, Christopher Bretherton - original contribution
  24 !! Jim McCaa - late 1990s, 2000-2001
  25 !! Nicolai Thum  2004-2005(?)
  26 !! Qingtao Song - 2007-2008
  27 !! Natalie Perlin 2011-2012
  28 !
  29
  30 MODULE module_bl_gbmpbl
  31
  32    USE module_model_constants, ONLY: cp, g, rcp, r_d,            &
  33                                      r_v, svp1, svp2, svp3,      &
  34                                      svpt0, ep_1, ep_2, xlv,     &
  35                                      karman
  36 public gbmpbl
  37 public gbmpblinit
  38 private
  39
  40 CONTAINS
  41
  42   SUBROUTINE GBMPBL(U3D,V3D,TH3D,T3D,QV3D,QC3D,QI3D,P3D,PI3D,     &
  43        RUBLTEN,RVBLTEN,RTHBLTEN,                                  &
  44        RQVBLTEN,RQCBLTEN,RQIBLTEN,                                &
  45        KZM_GBM,KTH_GBM,KETHL_GBM,EL_GBM,                          &
  46        dz8w,z,PSFC,TKE_PBL,RTHRATEN,                              &
  47        ZNT,UST,ZOL,HOL,PBL,KPBL2D,REGIME,PSIM,PSIH,               &
  48        XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR,                          &
  49        DT,DTMIN,                                                  &
  50        ids,ide, jds,jde, kds,kde,                                 &
  51        ims,ime, jms,jme, kms,kme,                                 &
  52        its,ite, jts,jte, kts,kte                        )
  53     !-------------------------------------------------------------------
  54     IMPLICIT NONE
  55     !-------------------------------------------------------------------
  56     !-- U3D         3D u-velocity interpolated to theta points (m/s)
  57     !-- V3D         3D v-velocity interpolated to theta points (m/s)
  58     !-- TH3D    3D potential temperature (K)
  59     !-- T3D         temperature (K)
  60     !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
  61     !-- QC3D        3D cloud mixing ratio (Kg/Kg)
  62     !-- QI3D        3D ice mixing ratio (Kg/Kg)
  63     !-- P3D         3D pressure (Pa)
  64     !-- PI3D    3D exner function (dimensionless)
  65     !-- RUBLTEN     Rho_dU tendency due to
  66     !               PBL parameterization (kg/m^3 . m/s)
  67     !-- RVBLTEN     Rho_dV tendency due to
  68     !               PBL parameterization (kg/m^3 . m/s)
  69     !-- RTHBLTEN    Rho_dTheta_m tendency due to
  70     !               PBL parameterization (kg/m^3 . K)
  71     !-- RQVBLTEN    Rho_dQv tendency due to
  72     !               PBL parameterization (kg/m^3 . kg/kg)
  73     !-- RQCBLTEN    Rho_dQc tendency due to
  74     !               PBL parameterization (kg/m^3 . kg/kg)
  75     !-- RQIBLTEN    Rho_dQi tendency due to
  76     !               PBL parameterization (kg/m^3 . kg/kg)
  77     !
  78     !-- KZM_GBM       exchange coefficient for momentum (m^2/s)
  79     !-- KTH_GBM       exchange coefficient for heat (K m/s)
  80     !-- KETHL_GBM     exchange coeff. for TKE enhanced (m^2/s)
  81     !
  82     ! NB!! Vertical staggering of the GBM output variables
  83     !  is as follows:
  84     !  KZM_GBM, KTH_GBM are on full levels, starting from the sfc
  85     !  KETHL_GBM is on half-levels
  86     !  TKE_PBL and EL_GBM are on full levels, starting from the sfc
  87     !
  88     !-- cp          heat capacity at constant pressure for dry air (J/kg/K)
  89     !-- g           acceleration due to gravity (m/s^2)
  90     !-- rcp         R/CP
  91     !-- r_d         gas constant for dry air (J/kg/K)
  92     !-- P_QI    species index for cloud ice
  93     !-- dz8w    dz between full levels (m)
  94     !-- z               height above sea level (m)
  95     !-- XLV         latent heat of vaporization (J/kg)
  96     !-- r_v             gas constant for water vapor (J/kg/K)
  97     !-- PSFC        pressure at the surface (Pa)
  98     !-- ZNT             roughness length (m)
  99     !-- UST             u* in similarity theory (m/s)
 100     !-- ZOL             z/L height over Monin-Obukhov length
 101     !-- HOL             PBL height over Monin-Obukhov length
 102     !-- PBL             PBL height (m)
 103     !-- KPBL2D     Layer index containing PBL top within or at the base interface
 104     !-- REGIME  flag indicating PBL regime (stable, unstable, etc.)
 105     !-- PSIM        similarity stability function for momentum
 106     !-- PSIH        similarity stability function for heat
 107     !-- XLAND   land mask (1 for land, 2 for water)
 108     !-- HFX             upward heat flux at the surface (W/m^2)
 109     !-- QFX             upward moisture flux at the surface (kg/m^2/s)
 110     !-- TSK             surface temperature (K)
 111     !-- GZ1OZ0      log(z/z0) where z0 is roughness length
 112     !-- WSPD        wind speed at lowest model level (m/s)
 113     !-- BR          bulk Richardson number in surface layer
 114     !-- DT              time step (s)
 115     !-- DTMIN   time step (minute)
 116     !-- svp1        constant for saturation vapor pressure (kPa)
 117     !-- svp2        constant for saturation vapor pressure (dimensionless)
 118     !-- svp3        constant for saturation vapor pressure (K)
 119     !-- svpt0       constant for saturation vapor pressure (K)
 120     !-- ep_1        constant for virtual temperature (r_v/r_d - 1) (dimensionless)
 121     !-- ep_2        constant for specific humidity calculation
 122     !-- karman      Von Karman constant
 123     !-- ids         start index for i in domain
 124     !-- ide         end index for i in domain
 125     !-- jds         start index for j in domain
 126     !-- jde         end index for j in domain
 127     !-- kds         start index for k in domain
 128     !-- kde         end index for k in domain
 129     !-- ims         start index for i in memory
 130     !-- ime         end index for i in memory
 131     !-- jms         start index for j in memory
 132     !-- jme         end index for j in memory
 133     !-- kms         start index for k in memory
 134     !-- kme         end index for k in memory
 135     !-- its         start index for i in tile
 136     !-- ite         end index for i in tile
 137     !-- jts         start index for j in tile
 138     !-- jte         end index for j in tile
 139     !-- kts         start index for k in tile
 140     !-- kte         end index for k in tile
 141     !-------------------------------------------------------------------
 142
 143     INTEGER,  INTENT(IN   )   ::      ids,ide, jds,jde, kds,kde, &
 144          ims,ime, jms,jme, kms,kme, &
 145          its,ite, jts,jte, kts,kte
 146
 147     REAL,     INTENT(IN   )   ::      DT,DTMIN
 148
 149     REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 150          INTENT(IN   )   ::                           QV3D, &
 151          QC3D, &
 152          QI3D, &
 153          P3D, &
 154          PI3D, &
 155          TH3D, &
 156          T3D, &
 157          dz8w, &
 158          z
 159
 160     !
 161     REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 162          INTENT(INOUT)   ::                          RUBLTEN, &
 163          RVBLTEN, &
 164          RTHBLTEN, &
 165          RQVBLTEN, &
 166          RQCBLTEN, &
 167          RQIBLTEN
 168
 169     REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 170          INTENT(OUT)   ::   KZM_GBM,KTH_GBM,KETHL_GBM, EL_GBM
 171
 172     REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 173          INTENT(IN   )   ::                          XLAND, &
 174          HFX, &
 175          QFX, &
 176          REGIME
 177
 178     REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 179          INTENT(INOUT)   ::                            HOL, &
 180          UST, &
 181          PBL, &
 182          ZNT
 183     INTEGER,  DIMENSION( ims:ime, jms:jme )                    , &
 184          INTENT(INOUT)   ::                            KPBL2D
 185
 186     !
 187 !m The following 5 variables are changed to memory size from tile size--
 188 !
 189      REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  ::    &
 190                                                              PSIM, &
 191                                                              PSIH
 192
 193      REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)   ::   &
 194                                                              WSPD
 195
 196      REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(IN   )   ::   &
 197                                                            GZ1OZ0, &
 198                                                                BR
 199
 200       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 201                 INTENT(IN   )               ::               PSFC
 202
 203       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 204                 INTENT(IN   )   ::                            TSK
 205
 206       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 207                 INTENT(INOUT)   ::                            ZOL
 208
 209       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 210                 INTENT(IN   )   ::                            U3D, &
 211                                                               V3D
 212
 213     REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME)                          &
 214          ,INTENT(INOUT) ::                    TKE_PBL
 215     REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME)                          &
 216          ,INTENT(IN) ::                    RTHRATEN
 217
 218     ! LOCAL VARIABLES
 219
 220     REAL,       DIMENSION( its:ite, kts:kte )          ::   dz8w2d, &
 221          z2d,t2dten_ra
 222
 223     INTEGER ::  I,J,K,NK,pass
 224     CHARACTER(LEN=200)        :: string
 225     REAL,DIMENSION(IMS:IME,KMS:KME)   :: TKE2d_1,TKE2d_2 ,          &
 226                                          u2dblten_1,u2dblten_2,     &
 227                                          v2dblten_1,v2dblten_2
 228
 229     tke2d_1 = 0.0
 230     tke2d_2 = 0.0
 231     u2dblten_1 = 0.0
 232     u2dblten_2 = 0.0
 233     v2dblten_1 = 0.0
 234     v2dblten_2 = 0.0
 235     t2dten_ra  = 0.0
 236
 237     DO J=jts,jte
 238        DO k=kts,kte
 239           NK=kme-k
 240           DO i=its,ite
 241              dz8w2d(I,K) = dz8w(i,K,j)
 242              z2d(I,K) = z(i,NK,j)
 243              t2dten_ra(i,k) = RTHRATEN(i,k,j)*PI3D(I,K,J)
 244              !t2dten is flipped inside the code
 245 ! TKE_PBL = tke2d = 0.5 * (q**2) == e
 246 ! tke_2d_2 in gbmpbl is equivalent to "e" in GB01
 247              tke2d_2(i,k)= TKE_PBL(i,k,j)
 248           ENDDO
 249        ENDDO
 250
 251     do pass=1,2
 252        CALL GBM2D(J,U3D(ims,kms,j),V3D(ims,kms,j),T3D(ims,kms,j),&
 253             QV3D(ims,kms,j),QC3D(ims,kms,j),QI3D(ims,kms,j),     &
 254             P3D(ims,kms,j),u2dblten_2(ims,kms),v2dblten_2(ims,kms),&
 255             RTHBLTEN(ims,kms,j),RQVBLTEN(ims,kms,j),             &
 256             RQCBLTEN(ims,kms,j),RQIBLTEN(ims,kms,j),             &
 257             KZM_GBM(ims,kms,j),KTH_GBM(ims,kms,j),               &
 258             KETHL_GBM(ims,kms,j),EL_GBM(ims,kms,j),              &
 259             TKE2d_2(ims,kms),t2dten_ra,                          &
 260             dz8w2d,z2d,                                          &
 261             PSFC(ims,j),ZNT(ims,j),UST(ims,j),ZOL(ims,j),        &
 262             HOL(ims,j),PBL(ims,j),KPBL2D(ims,j),PSIM(ims,j),     &
 263             PSIH(ims,j),XLAND(ims,j),HFX(ims,j),QFX(ims,j),      &
 264             TSK(ims,j),GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),      &
 265             DT,DTMIN,                                            &
 266             ids,ide, jds,jde, kds,kde,                           &
 267             ims,ime, jms,jme, kms,kme,                           &
 268             its,ite, jts,jte, kts,kte                            )
 269       if(pass==1)then
 270         TKE2d_1=tke2d_2 !tke at n+1
 271         u2dblten_1=u2dblten_2
 272         v2dblten_1=v2dblten_2
 273       end if
 274       if(pass==2)then
 275 ! tke_2d_2 in gbmpbl is equivalent to e == 2*(q**2)
 276         TKE_PBL(:,:,j)=tke2d_1 !otherwise tke is advanced to n+2
 277         rublten(:,:,j)=0.5*(u2dblten_2+u2dblten_1)
 278         rvblten(:,:,j)=0.5*(v2dblten_2+v2dblten_1)
 279       end if
 280     end do
 281
 282
 283        DO k=kts,kte
 284           DO i=its,ite
 285              RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/PI3D(I,K,J)
 286           ENDDO
 287        ENDDO
 288     ENDDO
 289
 290   END SUBROUTINE GBMPBL
 291
 292   SUBROUTINE GBM2D(J,U2D,V2D,T2D,QV2D,QC2D,QI2D,P2D,    &
 293        U2DTEN,V2DTEN,T2DTEN,                            &
 294        QV2DTEN,QC2DTEN,QI2DTEN,                         &
 295        KZM2D,KTH2D,KETHL2D,EL2D,                        &
 296        TKE2D,T2DTEN_RA,                                 &
 297        dz8w2d,z2d,PSFCPA,                               &
 298        ZNT,UST,ZOL,HOL,PBL,kpbl1d,PSIM,PSIH,            &
 299        XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR,                &
 300        DT,DTMIN,                                        &
 301        ids,ide, jds,jde, kds,kde,                       &
 302        ims,ime, jms,jme, kms,kme,                       &
 303        its,ite, jts,jte, kts,kte                        )
 304     !-------------------------------------------------------------------
 305     implicit none
 306     !
 307     INTEGER,  INTENT(IN   )   ::      ids,ide, jds,jde, kds,kde, &
 308          ims,ime, jms,jme, kms,kme, &
 309          its,ite, jts,jte, kts,kte,J
 310     !
 311     REAL,     INTENT(IN   )   ::      DT,DTMIN
 312     REAL                      ::     SVP1PA
 313
 314     REAL,     DIMENSION( ims:ime, kms:kme )                    , &
 315          INTENT(IN   )   ::                           QV2D, &
 316          QC2D, &
 317          QI2D, &
 318          P2D, &
 319          T2D
 320
 321     !
 322     REAL,     DIMENSION( ims:ime, kms:kme )                    , &
 323          INTENT(INOUT)   ::                         U2DTEN, &
 324          V2DTEN, &
 325          T2DTEN, &
 326          QV2DTEN, &
 327          QC2DTEN, &
 328          QI2DTEN, &
 329          tke2d
 330     !
 331     REAL,     DIMENSION( ims:ime, kms:kme )                    , &
 332          INTENT(OUT)   ::    KZM2d,KTH2d,KETHL2d,EL2D
 333
 334     REAL,     DIMENSION( ims:ime )                             , &
 335          INTENT(INOUT)   ::                            HOL, &
 336          UST, &
 337          PBL, &
 338          ZNT
 339
 340     INTEGER,  DIMENSION( ims:ime )                             , &
 341          INTENT(INOUT)   ::                            kpbl1d
 342
 343     REAL,     DIMENSION( ims:ime )                             , &
 344          INTENT(IN   )   ::                          XLAND, &
 345          HFX, &
 346          QFX
 347     !
 348     REAL,     DIMENSION( ims:ime ), INTENT(IN   )   ::      PSIM, &
 349          PSIH
 350
 351     REAL,     DIMENSION( ims:ime ), INTENT(INOUT)   ::      WSPD
 352
 353     REAL,     DIMENSION( ims:ime ), INTENT(IN   )   ::    GZ1OZ0, &
 354          BR
 355
 356     REAL,     DIMENSION( ims:ime )                             , &
 357          INTENT(IN   )               ::             PSFCPA
 358
 359     REAL,     DIMENSION( ims:ime )                             , &
 360          INTENT(IN   )   ::                            TSK
 361
 362     REAL,     DIMENSION( ims:ime )                             , &
 363          INTENT(INOUT)   ::                            ZOL
 364
 365     !
 366     REAL,     DIMENSION( its:ite, kts:kte ) ,                    &
 367          INTENT(IN)      ::                         dz8w2d, &
 368          z2d,t2dten_ra
 369     !
 370     REAL,     DIMENSION( ims:ime, kms:kme )                    , &
 371          INTENT(IN   )   ::                            U2D, &
 372          V2D
 373     !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 374     !     PBL SCHEME                                                      C
 375     !     Questions? Contact mccaa@atmos.washington.edu                   C
 376     !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 377     !
 378     !     Order of events:
 379     !     0.  Preliminaries:  Some things we only need to do once
 380     !     1.  First, save initial values and calculate some derivatives.
 381     !     4.  Calculate the buoyancy jump at flux levels
 382     !     5.  Calculate the pbl height and length scale
 383     !     6.  Compute diffusivity profiles
 384     !     7.  Implicit diffusion of thetal, total water
 385     !     8.  Implicit diffusion of momentum
 386     !     9.  Implicit diffusion of cloud ice
 387     !     10. Prediction of TKE
 388     !     11. Calculation of tendencies for output to model
 389     !
 390     !     Local variables on full levels
 391     real zqx(kts:kte+2)
 392     REAL KTH(KTS:KTE+1),KZM(KTS:KTE+1),RHOXFL(KTS:KTE+1),tke(kts:kte+1),tkes(kts:kte+1), &
 393          rrhoxfl(kts:kte+1),BBLS(KTS:KTE+1),NSQUAR(KTS:KTE+1),BOUYAN(KTS:KTE+1), &
 394          DQWDZ(kts:kte+1),rdza(kts:kte+1),dza(kts:kte+1),SVS(KTS:KTE+1),presfl(kts:kte+1), &
 395          exnerfl(kts:kte+1),SHEAR(KTS:KTE+1),rexnerfl(kts:kte+1),rcldb(kts:kte+1), &
 396          epop(kts:kte+1),DTLDZ(KTS:KTE+1)
 397     !     Local variables on half levels
 398     REAL UX(KTS:KTE),VX(KTS:KTE),THX(KTS:KTE),QX(KTS:KTE),THVX(KTS:KTE),zax(kts:kte),qix(kts:kte), &
 399          KETHL(KTS:KTE),THLX(KTS:KTE),THXS(KTS:KTE),tx(kts:kte),tvx(kts:kte),rttenx(kts:kte), &
 400          PRESHL(KTS:KTE),QCX(KTS:KTE),QWX(KTS:KTE),dzq(kts:kte),rRHOXHL(KTS:KTE),UXS(KTS:KTE), &
 401          QXS(KTS:KTE),RHOXHL(KTS:KTE),exnerhl(kts:kte),rexnerhl(kts:kte),rdzq(kts:kte), &
 402          VXS(KTS:KTE),qixs(kts:kte),qcxs(kts:kte)
 403     REAL,     DIMENSION( its:ite ) :: wspd1
 404     REAL UFLXP,VFLXP,RHOXSF,Q0S, &
 405          RDT,dt2, &
 406          aone,atwo,czero,tskx, &
 407          tvcon,fracz,dudz,dvdz,rvls,thv0,dthv,xfr, &
 408          cpoxlv,r1cp,rczero, &
 409          templ,temps,gamcrt,gamcrq,cell, &
 410          rwspd,cfac, &
 411          thgb,pblx,gothv,capcd,capch,capcq, &
 412          ustx,ustxm,qfxx,hfxx,rlwp,tkeavg,wstar,xlvocp,wso,phih, &
 413          rstbl,vk,tbbls,pref
 414     integer i,k,l,iconv,ilay, &
 415          ktop(kts:kte),kpbl2dx,kmix2dx, &
 416          iteration,pass,kr,&
 417          ktop_save(kts:kte) !NT new to store original height in case layers get merged
 418     !     Arrays for tridiagonal matrix solver
 419     real aimp(kts:kte),bimp(kts:kte),cimp(kts:kte)
 420     real uimp1(kts:kte),rimp1(kts:kte)
 421     real uimp2(kts:kte),rimp2(kts:kte)
 422     !NT some temporary variables to recompute n2
 423     real  THX_t,THVX_t,DTHV_t
 424     !
 425     parameter(XFR=0.1)        !Fraction of turb layer to be considered in bbls
 426 !    parameter(aone=1.9*xfr,atwo=15.,cfac=7.8)  orig, why is atwo 15!???? atn aone 1.9
 427     parameter(aone=1.9*xfr,atwo=10.,cfac=7.8)
 428     parameter(czero=5.869)    ! b1/2**(3/2)
 429     parameter(rstbl=1.5)      ! empirical constant for l at stable interfaces
 430     PARAMETER (GAMCRT=3.,GAMCRQ=2.E-3)
 431     parameter(vk=0.4)
 432     parameter(pref = 100000.)
 433     !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 434     !
 435     !     ---0---  Preliminaries:  Some things we only need to do once
 436     !
 437     RDT=1./DT
 438     dt2=dt
 439     SVP1PA=svp1*10. ! NT actually this is svp1 in mb
 440     cpoxlv = cp/xlv
 441     xlvocp = xlv/cp
 442     r1cp = 1./cp
 443     rczero = 1./czero
 444
 445     !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 446     !     This is the beginning of big I loop
 447     !
 448     !     ---1---   First, save initial values and calculate some derivatives.
 449     !
 450     DO I=its,ite
 451        !     Copy in j-dependent variables
 452        tskx = tsk(i)
 453        pblx = pbl(i)
 454        qfxx = qfx(i)
 455        hfxx = hfx(i)
 456        zqx(kte+1)=0.0
 457        zqx(kte+1+1)=0.0
 458        tke(kte+1) = 0.0
 459        DO K = kte,kts,-1
 460           KR = kte+1-K
 461           zqx(k)= zqx(k+1) + dz8w2d(i,kr)
 462           zax(k)=0.5*(zqx(k)+zqx(k+1))
 463           tke(k) = tke2d(i,kr)
 464           !     Wind components
 465           UX(K)=U2D(I,KR)  !
 466           VX(K)=V2D(I,KR)  !
 467           TX(K)=T2D(I,KR)
 468           QX(K)=QV2D(I,KR)
 469           QCX(K)=QC2D(I,KR)
 470           qix(k)=qi2D(i,kr)
 471           rttenx(k)=t2dten_ra(i,kr)
 472        ENDDO
 473        !     Done copying j-dependent variables
 474        tkes(kte+1) = tke(kte+1)
 475        do k = kts,kte
 476           kr = kte+1-k
 477           !     Pressure at half levels
 478           PRESHL(K)=p2d(i,kr)
 479           DZQ(K)=ZQX(K)-ZQX(K+1)
 480           rdzq(k)=1./dzq(k)
 481           exnerhl(k)=(preshl(k)/pref)**rcp
 482           rexnerhl(k)=1./exnerhl(k)
 483           THX(K)=TX(K)*rexnerhl(k)
 484           QWX(K) = QX(K) + QCX(K)
 485           TVCON=(1.+ep_1*QX(K)-qcx(k))
 486           TVX(K)=TX(K)*TVCON  ! virtual temperature
 487           THVX(K)=THX(K)*TVCON
 488           THLX(K) = THX(K) - XLVOCP * rexnerhl(K) * QCX(K)
 489           !     SAVE INITIAL VALUES of winds, thx, and qx
 490           tkes(k) = tke(k)
 491           UXS(K)  = UX(K)
 492           VXS(K)  = VX(K)
 493           THXS(K) = THX(K)
 494           QXS(K) = QX(K)
 495           qixs(k) = qix(k)
 496           qcxs(k) = qcx(k)
 497           !     Density at half levels
 498           RHOXHL(K)=PRESHL(K)/(r_d*TVX(K))
 499           rrhoxhl(k)=1./rhoxhl(k)
 500        enddo
 501        !     Density and vertical derivatives at the full sigma levels.
 502        presfl(1)=0.
 503        DO K = 2, kte
 504           kr = kte+1-k
 505           !     Pressure at full levels
 506           PRESFL(K)=exp(0.5*(log(p2d(i,kr+1))+log(p2d(i,kr))))
 507           epop(k)=ep_2/presfl(k)
 508           DZA(K)=ZAX(K-1)-ZAX(K)
 509           rdza(k)=1./dza(k)
 510           exnerfl(k)=(presfl(k)/pref)**rcp
 511           rexnerfl(k)=1./exnerfl(k)
 512           FRACZ=(ZQX(K)-ZAX(K))*RDZA(K)
 513           RHOXFL(K)=RHOXHL(K)+(RHOXHL(K-1)-RHOXHL(K))*FRACZ
 514           RRHOXFL(K)=1./RHOXFL(K)
 515           DUDZ =  (UX(K-1)  - UX(K))  *RDZA(K)
 516           DVDZ =  (VX(K-1)  - VX(K))  *RDZA(K)
 517           SVS(K)= DUDZ*DUDZ + DVDZ*DVDZ
 518           DQWDZ(K) =  (QWX(K-1)  - QWX(K))  *RDZA(K)
 519           DTLDZ(K) =  (THLX(K-1)  - THLX(K))  *RDZA(K)
 520        ENDDO
 521 !write(*,*)'dzq', dzq
 522 !write(*,*)'dza', dza
 523        !     Pressure at the surface (in centibars)
 524        PRESFL(KTE+1)=psfcpa(i)
 525        rexnerfl(kte+1)=(pref/presfl(kte+1))**rcp
 526        exnerfl(kte+1)=1./rexnerfl(kte+1)
 527        !     Saturation Mixing Ratio at Surface
 528        q0s=ep_2/(PRESFL(KTE+1)/(100.*SVP1PA*EXP(svp2*(tskx-svpt0)/(tskx-svp3)))-1.)
 529        !     COMPUTE SOIL POTENTIAL TEMPERATURE
 530        THGB = TSKX * rexnerfl(kte+1)
 531        !     DENSITY AT THE SURFACE
 532        RHOXSF=(PRESFL(KTE+1))/(r_d*TVX(kte))
 533        !     More surface variables
 534        WSPD1(i)=SQRT(UX(KTE)*UX(KTE)+VX(KTE)*VX(KTE))+1.e-4
 535        THV0=THGB*(1.+ep_1*Q0S)
 536        DTHV=(THVX(KTE)-THV0)
 537        gothv = g/thvx(kte)
 538        !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 539        ustx = wspd1(i)*vk/(gz1oz0(i)-psim(i))
 540        ustxm = amax1(ustx,0.1)
 541        if ((xland(i)-1.5).lt.0)  ustx = ustxm
 542        !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 543        !
 544        !
 545        !     ---4---   Calculate the buoyancy jump at flux levels
 546        !
 547        call n2(thlx,exnerfl,epop,qwx,cpoxlv,rdza,            &
 548             rexnerfl,kts,kte,nsquar,rcldb,xlvocp,svp1pa)
 549        nsquar(kte+1)=  g/thvx(kte) * dthv / zax(kte)
 550        !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 551        !
 552        !     ---5--- Calculate the pbl height and length scale
 553        !
 554        call pblhgt( &
 555             ! input variables &
 556        zqx,kts,kte, &
 557             nsquar,tke,presfl,rhoxfl,rcldb,exnerfl, &
 558             rttenx,thvx,thlx,thx,qwx,rexnerhl,qcx, &
 559             xfr,xlvocp,aone,atwo,rstbl,            &
 560             ! output variables &
 561        bbls,pblx, &
 562             ktop,iconv,kmix2dx,kpbl2dx)
 563        !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 564        !
 565        !     ---6---   Compute diffusivity profiles
 566        !
 567        do pass = 1,2
 568           call my( &
 569                ! Input variables &
 570           bbls,nsquar,tke,iconv,ktop,thlx,thx,qwx,     &
 571                xlvocp,rcldb,rexnerhl,aone,atwo,kts,kte, &
 572                ! Output variables &
 573           kzm,kth,kethl )
 574           !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 575           !
 576           !     ---7---    Implicit diffusion of thetal, total water
 577           !
 578           !     First find the coefficients that apply to all scalars at half-levels
 579           do k = 1, kte
 580              if(k.eq.1)then
 581                 aimp(k)=0.
 582              else
 583                 aimp(k)=-(rhoxfl(k)*rrhoxhl(k))* &
 584                      kth(k) * dt2 *rdzq(k)*rdza(k)
 585              endif
 586              if(k.eq.kte)then
 587                 cimp(k)=0.
 588              else
 589                 cimp(k)=-(rhoxfl(k+1)*rrhoxhl(k))* &
 590                      kth(k+1) * dt2 *rdzq(k)*rdza(k+1)
 591              endif
 592              bimp(k) = 1 - aimp(k) - cimp(k)
 593              !     Now find right side for various scalars:
 594              !     No flux out top, so no (k.eq.1)
 595              if(k.eq.kte)then  ! at surface
 596                 !     Include surface latent heat flux
 597                 rimp2(k) = qwx(k) + dt2 * qfxx*rrhoxhl(k)*rdzq(kte)
 598                 !     Include surface sensible heat flux
 599                 rimp1(k) = thlx(k) + &
 600                      dt2 * hfxx*rrhoxhl(k)*r1cp*rdzq(kte)*rexnerhl(kte)
 601              else
 602                 rimp2(k) = qwx(k)
 603                 rimp1(k) = thlx(k)
 604              endif
 605           enddo
 606           call tridag(aimp,bimp,cimp,rimp1,uimp1,kte)
 607           call tridag(aimp,bimp,cimp,rimp2,uimp2,kte)
 608           !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 609           !
 610           !     Recompute the buoyancy jump at flux levels
 611           !
 612           call n2(uimp1,exnerfl,epop,uimp2,cpoxlv,rdza,       &
 613                rexnerfl, kts,kte,nsquar,rcldb,xlvocp,svp1pa)
 614                !NT new recompute also n2 at sfc!!
 615                ! THLX(K) = THX(K) - XLVOCP * rexnerhl(K) * QCX(K)
 616                THX_t = uimp1(KTE) + XLVOCP * rexnerhl(KTE) * QCX(KTE)
 617                THVX_t=THX_t*TVCON
 618                DTHV_t=(THVX_t-THV0)
 619                nsquar(kte+1)=  g/thvx_t * dthv_t / zax(kte)
 620                !NT end
 621           !     Update only nsquar, then go back
 622        !TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
 623        !
 624        !       Recompute pbl-height? since n2 changes; use uimp1 instead thlx
 625        !       TEST!!!
 626        !
 627 !NT       call pblhgt( &
 628 !NT            ! input variables &
 629 !NT       zqx,kts,kte, &
 630 !NT            nsquar,tke,presfl,rhoxfl,rcldb,exnerfl, &
 631 !NT            rttenx,thvx,uimp1,thx,qwx,rexnerhl,qcx, &
 632 !NT            xfr,xlvocp,aone,atwo,rstbl, &
 633 !NT            ! output variables &
 634 !NT       bbls,pblx, &
 635 !NT            ktop,iconv,kmix2dx,kpbl2dx)
 636        !TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
 637        enddo
 638        !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 639        !
 640        !     Update the scalars (only after second pass)
 641        !
 642        !     Calculate new values of thx, qx and qcx using new values of thlx and qcx.
 643        do  k = 1, kte
 644           thlx(k) = uimp1(k)
 645           qwx(k) = uimp2(k)
 646           if(k.ge.2) DQWDZ(K) = (QWX(K-1) - QWX(K)) *RDZA(K)
 647           if(k.ge.2) DTLDZ(K) = (THLX(K-1) - THLX(K)) *RDZA(K)
 648           templ = thlx(k)*exnerhl(k)
 649           temps = templ
 650           rvls = ep_2/ &
 651                (preshl(k)/(100.*svp1pa*EXP(svp2*(temps-svpt0)/(temps-svp3)))-1.)
 652           do iteration = 1, 3
 653              temps = temps  + ((templ-temps)*cp/xlv + qwx(k)-rvls)/ &
 654                   (cp/xlv+ep_2*xlv*rvls/r_d/temps/temps)
 655              rvls = ep_2/ &
 656                   (preshl(k)/(100.*svp1pa*EXP(svp2*(temps-svpt0)/(temps-svp3)))-1.)
 657           enddo
 658           qcx(k)=max(qwx(k)-rvls,0.)
 659           qx(k) = qwx(k) - qcx(k)
 660           thx(k) = (templ+qcx(k)*xlv/cp) / exnerhl(k) ! theta_l
 661           THVX(K)=THX(K)*(1.+ep_1*QX(K)-QCX(K)) !NT is this GB01:A13 then it should
 662 !NT       qx=q_v ; qwx=q_t; qcx=q_l
 663 !nt thvx is theta virtual
 664        ENDDO
 665        !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 666        !
 667        !     ---8--- Implicit diffusion of momentum
 668        !
 669        !     First find the coefficients that apply to winds at half-levels
 670        do k = 1, kte
 671           if(k.eq.1)then
 672              aimp(k)=0.
 673           else
 674              aimp(k)=-(rhoxfl(k)*rrhoxhl(k))* &
 675                   kzm(k) * dt2 *rdzq(k)*rdza(k)
 676           endif
 677           if(k.eq.kte)then
 678              cimp(k)=0.
 679           else
 680              cimp(k)=-(rhoxfl(k+1)*rrhoxhl(k))* &
 681                   kzm(k+1) * dt2 *rdzq(k)*rdza(k+1)
 682           endif
 683           bimp(k) = 1 - aimp(k) - cimp(k)
 684           !     Now find right side
 685           !     No flux out top, so no (k.eq.1)
 686           if(k.eq.kte)then
 687              !     At surface
 688              UFLXP=-USTX*USTX*UX(KTE)/WSPD1(i)
 689              VFLXP=-USTX*USTX*VX(KTE)/WSPD1(i)
 690              !     Include surface momentum fluxes
 691              rimp1(k) = ux(k) + dt2 *  &
 692                   uflxp * (rhoxsf*rrhoxhl(k)) *rdzq(kte)
 693              rimp2(k) = vx(k) + dt2 *  &
 694                   vflxp * (rhoxsf*rrhoxhl(k)) *rdzq(kte)
 695           else
 696              rimp1(k) = ux(k)
 697              rimp2(k) = vx(k)
 698           endif
 699        enddo
 700        call tridag(aimp,bimp,cimp,rimp1,uimp1,kte)
 701        call tridag(aimp,bimp,cimp,rimp2,uimp2,kte)
 702        !     Update the winds
 703        do  k = 1, kte
 704           ux(k) = uimp1(k)
 705           vx(k) = uimp2(k)
 706           if(k.ge.2)then
 707              DUDZ =  (UX(K-1)  - UX(K))  *RDZA(K)
 708              DVDZ =  (VX(K-1)  - VX(K))  *RDZA(K)
 709              SVS(K)= DUDZ*DUDZ + DVDZ*DVDZ
 710           endif
 711        enddo
 712        !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 713        !
 714        !     ---9---    Implicit diffusion of cloud ice
 715        !
 716        !     First find the coefficients that apply to all scalars at half-levels
 717        do k = 1, kte
 718           if(k.eq.1)then
 719              aimp(k)=0.
 720           else
 721              aimp(k)=-(rhoxfl(k)*rrhoxhl(k))* &
 722                   kth(k) * dt2 *rdzq(k)*rdza(k)
 723           endif
 724           if(k.eq.kte)then
 725              cimp(k)=0.
 726           else
 727              cimp(k)=-(rhoxfl(k+1)*rrhoxhl(k))* &
 728                   kth(k+1) * dt2 *rdzq(k)*rdza(k+1)
 729           endif
 730           bimp(k) = 1 - aimp(k) - cimp(k)
 731           !     Now find right side for various scalars:
 732           !     No flux out top, so no (k.eq.1)
 733           !     No flux in bottom, so no (k.eq.kte)
 734           rimp1(k) = qix(k)
 735           !            rimp2(k) = qncx(k)
 736        enddo
 737        call tridag(aimp,bimp,cimp,rimp1,qix,kte)
 738        !         call tridag(aimp,bimp,cimp,rimp2,qncx,kte)
 739        !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 740        !
 741        !     ---10---   Prediction of TKE
 742        !
 743        !     Features:
 744        !     a. Explicit calculation of buoyancy and shear terms using
 745        !     time=t+1 values of thetal, qw, and winds
 746        !     b. Surface TKE is diagnosed
 747        !     c. Semi-implicit calculation of dissipation term and
 748        !     implicit calculation of turbulent transfer term
 749        !     First, buoyancy and shear terms
 750        DO K=2,KTE
 751           !     Compute buoyancy term with new values of thetal
 752           BOUYAN(K)=-KTH(K)*nsquar(K)
 753           !     Compute shear term with new values of svs
 754           SHEAR(K)  = KZM(K) * SVS(K)
 755        ENDDO
 756        do ilay = 1,iconv
 757           k = ktop(ilay)
 758           if(qcx(k).gt.1e-8.and.k.gt.1) bouyan(k) = bouyan(k) - & !NT GB01:A14
 759                rttenx(k)*(presfl(k+1)-presfl(k)) * rrhoxfl(k) &
 760                * rexnerfl(k) / thvx(k)
 761        enddo
 762        !     TKE at top is fixed
 763        tke(1)=0.
 764        bbls(1)= 0.0
 765        !     Diagnose TKE at surface, following MY 82, B1 ** (2/3) / 2 = 3.25
 766        tke(kte+1)=3.25 * ustx * ustx ! normal
 767        !     tke(kte+1)=.33*wstar**2 ! dryrun
 768        !     Now the implicit calculations
 769        !     First find the coefficients that apply for full levels
 770        do k = 2, kte
 771           if(k.eq.2)then
 772              aimp(k-1)=0.
 773           else
 774              aimp(k-1)=-(rhoxhl(k-1)*rrhoxfl(k))* &
 775                   kethl(k-1)*dt2*rdzq(k-1)*rdza(k)
 776           endif
 777           if(k.eq.kte)then
 778              cimp(k-1)=0.
 779              !     Account for implicit part of flux between level kte and surface
 780              if(bbls(k).gt.0)then
 781                 bimp(k-1) = 1 - aimp(k-1) - cimp(k-1) + &
 782                      dt2 * ( sqrt(tke(k))*rczero/bbls(k) + &
 783                      rhoxhl(k)*rrhoxfl(k)*kethl(k)*rdzq(k)*rdza(k) )
 784              else
 785                 bimp(k-1) = 1 - aimp(k-1) - cimp(k-1) + &
 786                      dt2 * rhoxhl(k)*rrhoxfl(k)* kethl(k)*rdzq(k)*rdza(k)
 787              endif
 788
 789           else
 790              cimp(k-1)=-(rhoxhl(k)*rrhoxfl(k))* &
 791                   kethl(k)*dt2*rdzq(k)*rdza(k)
 792              tbbls = max(bbls(k),bbls(k+1))
 793              if(tbbls.gt.0)then
 794                 bimp(k-1)= 1 - aimp(k-1) - cimp(k-1) + dt2 * sqrt(tke(k))*rczero/tbbls
 795              else
 796                 bimp(k-1)= 1 - aimp(k-1) - cimp(k-1)
 797              endif
 798           endif
 799           !     Now find right side
 800           if(k.eq.kte)then
 801              !     Account for fixed part of flux between level kte and surface
 802              rimp1(k-1) = tke(k) + dt2 * ( SHEAR(K)+BOUYAN(K) + &
 803                   tke(kte+1)*(rhoxhl(k)*rrhoxfl(k))*kethl(k)*rdzq(k)*rdza(k) )
 804           else
 805              rimp1(k-1) = tke(k) + dt2 * (SHEAR(K)+BOUYAN(K))
 806           endif
 807        enddo
 808
 809        call tridag(aimp,bimp,cimp,rimp1,uimp1,kte-1)
 810        !     Update the tke
 811        do  k = 2, kte
 812           tke(k) = max(uimp1(k-1),0.001) ! background tke .001
 813        enddo
 814        !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 815        !
 816        !     ---11---   Calculation of tendencies for output to model
 817        !
 818
 819        DO K=KTS,KTE
 820           kr = kte+1-k
 821           U2DTEN(I,KR)= (UX(K)-UXS(K))*RDT
 822           V2DTEN(I,KR)=(VX(K)-VXS(K))*RDT
 823           T2DTEN(I,KR)=(THX(K)-THXS(K))*EXNERHL(K)*RDT
 824           QV2DTEN(I,KR) = (QX(K)-QXS(K))*RDT
 825           qi2Dten(i,kr) = (qix(k)-qixs(k))*rdt
 826           qc2Dten(i,kr) = (qcx(k)-qcxs(k))*rdt
 827           tke2d(i,kr)=tke(k)
 828           kzm2d(i,kr)=kzm(k)
 829           kth2d(i,kr)=kth(k)
 830           kethl2d(i,kr)=kethl(k)
 831           EL2D(i,kr) = bbls(k)
 832        ENDDO
 833        !     Copy out j-dependent variables
 834        pbl(i)=pblx
 835        ust(i)=ustx
 836        kpbl1d(i)=kte+1-kpbl2dx
 837     end do
 838     !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 839     RETURN
 840   END SUBROUTINE GBM2D
 841   !CCCCCCCCCCCCCCCCCCCCC END OF PBL SUBROUTINE CCCCCCCCCCCCCCCCCCCCCCCCC
 842   !
 843   subroutine tridag(a,b,c,r,u,n)
 844     !     Solves tridiagonal matrix
 845     !     See "Numerical Recipes in Fortran 77", p. 42
 846     implicit none
 847     integer n,nmax
 848     real a(n),b(n),c(n),r(n),u(n)
 849     parameter (nmax=100)
 850     integer j
 851     real rbet,gam(nmax)
 852     rbet=1./b(1)
 853     u(1)=r(1)*rbet
 854     do j=2,n
 855        gam(j)=c(j-1)*rbet
 856        rbet=1./(b(j)-a(j)*gam(j))
 857        u(j)=(r(j)-a(j)*u(j-1))*rbet
 858     end do
 859     do  j=n-1,1,-1
 860        u(j)=u(j)-gam(j+1)*u(j+1)
 861     end do
 862     return
 863   end subroutine tridag
 864
 865   subroutine pblhgt( &
 866        ! input variables
 867        zqx,kts,kte, &
 868        nsquar,tke,presfl,rhoxfl,rcldb,exnerfl, &
 869        rttenx,thvx,thlx,thx,qwx,rexnerhl,qcx, &
 870        xfr,xlvocp,aone,atwo,rstbl,            &
 871        ! output variables
 872        bbls,pblx, &
 873        ktop,iconv,kmix2dx,kpbl2dx &
 874        )
 875
 876     implicit none
 877     !     Input variables
 878     real zqx(kts:kte+2)
 879     real nsquar(kts:kte+1),tke(kts:kte+1),presfl(kts:kte+1),rhoxfl(kts:kte+1), &
 880          rcldb(kts:kte+1),exnerfl(kts:kte+1)
 881     real rttenx(kts:kte),thvx(kts:kte),thlx(kts:kte),thx(kts:kte),qwx(kts:kte), &
 882          rexnerhl(kts:kte),qcx(kts:kte)
 883     real xfr,xlvocp,aone,atwo,rstbl
 884     !     Output variables
 885     real bbls(kts:kte+1), pblx
 886     integer ktop(kts:kte+1), iconv,kmix2dx,kpbl2dx, &
 887             ktop_save(kts:kte+1)
 888
 889     !     Local variables
 890     integer kbot(kts:kte+1),kts,kte,kstart
 891     integer istabl,ibeg,ilay,nlev,k,itemp
 892     real blinf,rnnll,tkeavg,trnnll,radnnll,delthvl,elambda, &
 893          bige,biga,entnnll,tbbls
 894
 895     !     Find noncontiguous convectively unstable layers
 896     iconv = 0
 897     istabl = 1
 898     do k=2,kte+1   !nt nscquar is defined at kte+1 after the call to n2
 899        if(nsquar(k).le.0)then
 900           if(istabl.eq.1)then
 901              iconv = iconv + 1
 902              ktop(iconv)=k
 903           endif
 904           istabl=0
 905           kbot(iconv)=k
 906        else
 907           istabl=1
 908           BBLS(K) = MIN(rstbl*SQRT(TKE(K)/NSQUAR(K)),karman*ZQX(K))
 909        endif
 910     enddo
 911     !NT new
 912 !     if(ktop(iconv).eq.kte+1 .and. kbot(iconv).eq.kte+1)then
 913 !            iconv=iconv-1 !NT don't let the bottome layer be a convective top.
 914 !     end if
 915     !NT end
 916
 917     !     Now see if they have sufficient buoyant convection to connect
 918 !NT save ktops to be able to reset top of a merged layers
 919    ktop_save=ktop
 920 !NT put this if-statement in to make sure;
 921   IF(iconv.ge.1)THEN
 922 !NT
 923     ibeg = 1
 924 2745 do ilay = ibeg, iconv
 925        blinf = xfr*(zqx(ktop(ilay)-1) - zqx(kbot(ilay)+1))
 926        !     Find average n*n*l*l for layer
 927        rnnll = 0.
 928        tkeavg = 0.
 929        nlev = kbot(ilay)-ktop(ilay)+1
 930        do k = ktop(ilay), kbot(ilay)
 931           bbls(k) = min(blinf,karman*ZQX(K))
 932           rnnll = rnnll + nsquar(k)*bbls(k)*bbls(k)   !NT is not averaged, but correct
 933           tkeavg = tkeavg + tke(k) / nlev
 934        enddo
 935        !     First extend up
 936        kstart=ktop(ilay)-1
 937 !NT orig       do k = ktop(ilay)-1,2,-1
 938        do k = kstart,2,-1
 939           ktop(ilay)=k        ! We always go up at least one, for the entrainment interface
 940           bbls(k) = min(karman * ZQX(K),blinf)
 941           trnnll = nsquar(k)*bbls(k)*bbls(k)
 942           if(trnnll*nlev.ge.-0.5*rnnll) goto 2746 ! Is it the top?
 943 !NT test          if(trnnll*nlev.ge.-0.7*rnnll) goto 2746 ! Requires stronger jump to be top
 944 !NT doesn't workkpbl           if(trnnll*nlev.ge.-0.5*rnnll.and. &
 945 !NT n2 gets recomputed      nsquar(k).ge.1.e-6) goto 2746 ! Avoids weak layers with even weaker tops
 946 !NT after pblhgt                           ! Imagine n2 small negative over several layers
 947                                            ! then average rnnll is very small and a very weak
 948                                            ! pos n2 is enough to be 'inversion' top. This makes
 949                                            ! sure that we go up at least one more.
 950           if(ilay.gt.1) then
 951              if(ktop(ilay).eq.kbot(ilay-1))then ! did we merge with layer above?
 952                 ibeg = ilay - 1
 953    !NT orig     ktop(ibeg)=ktop(ibeg)+1 !NT not correct if one up was not inversion, the new thicker
 954                                         !NT layer might have different average properties, should
 955                                         !NT reset to original ktop
 956                 ktop(ibeg)=ktop_save(ibeg) !NT new
 957                 kbot(ibeg)=kbot(ibeg+1)
 958                 iconv = iconv - 1
 959                 do itemp = ibeg+1,iconv !NT if there's a layer below decrease layer index
 960                    ktop(itemp)=ktop(itemp+1)
 961                    kbot(itemp)=kbot(itemp+1)
 962                    ktop_save(itemp)=ktop_save(itemp+1) !NT new
 963                 enddo
 964                 goto 2745        ! recompute for the new, deeper layer
 965              endif
 966           endif
 967           rnnll = rnnll + trnnll
 968           nlev = nlev + 1 !NT moved up
 969        enddo
 970        !     Add radiative/entrainment contribution to total
 971 2746   k = ktop(ilay)
 972        radnnll = 0.
 973        if(qcx(k).gt.1e-8) radnnll =  &
 974             rttenx(k)*(presfl(k+1)-presfl(k))/ &
 975             (rhoxfl(k)*thvx(k)*exnerfl(k))
 976        entnnll = 0.
 977        if(k.ge.3)then
 978           delthvl = (thlx(k-2)+thx(k-2)*ep_1*qwx(k-2)) &
 979                - (thlx(k) + thx(k)*ep_1*qwx(k))
 980           elambda = xlvocp*rcldb(k)*rexnerhl(k)/max(delthvl,0.1)
 981           bige = 0.8 * elambda
 982           biga = aone * (1 + atwo * bige)
 983           entnnll = biga * sqrt(tkeavg**3) / bbls(k)
 984        endif
 985        if(tkeavg.gt.0.) rnnll = rnnll + min(0.,bbls(k)/sqrt(tkeavg) * (radnnll + entnnll) )
 986        !     Now extend down
 987        do k = kbot(ilay)+1,kte+1
 988           tbbls = min(karman * ZQX(K),blinf)
 989           trnnll = nsquar(k)*tbbls*tbbls
 990           if(trnnll*nlev.ge.-0.5*rnnll)goto 2747 !  Is it the bottom?
 991           kbot(ilay)=k
 992           if(ilay.lt.iconv.and.kbot(ilay).eq.ktop(ilay+1))then ! did we merge with layer below?
 993 !NT orig             ktop(ilay)=ktop(ilay)+1
 994              ktop(ilay)=ktop_save(ilay)
 995              kbot(ilay)=kbot(ilay+1)
 996              iconv = iconv - 1
 997              do itemp = ilay+1,iconv
 998                 ktop(itemp)=ktop(itemp+1)
 999                 kbot(itemp)=kbot(itemp+1)
1000                 ktop_save(itemp)=ktop_save(itemp+1)
1001              enddo
1002              goto 2745        ! recompute for the new, deeper layer
1003           endif
1004           rnnll = rnnll + trnnll
1005           bbls(k)=tbbls
1006           nlev = nlev + 1
1007        enddo
1008 2747   continue
1009     enddo !NT all sorting is done, go through layers and calk l
1010     do ilay = 1, iconv
1011        blinf = xfr*(zqx(ktop(ilay)-1) - zqx(kbot(ilay)+1))
1012        do k = ktop(ilay),kbot(ilay)
1013           bbls(k) = min(karman * ZQX(K),blinf)
1014        enddo
1015     enddo
1016 !NT added to check if iconv .ge. 1
1017  ENDIF
1018 !NT end
1019     !     We should now have tops and bottoms for iconv layers
1020     ! NT not clear how kpbl2dx should work, but doesn't matter since
1021     ! NT only kmix2dx is used no matter what kpbl2dx is.
1022     ! NT Looks like it could be used to choose between mixing pbl and
1023     ! NT convective pbl height if there are more than 1 unstable layers
1024
1025     if(iconv.gt.0)then
1026        if(kbot(iconv).eq.kte+1)then
1027           kmix2dx = ktop(iconv)
1028           if(kpbl2dx.ge.0)then
1029              if(iconv.gt.1)then
1030                 kpbl2dx = ktop(iconv-1)
1031              else
1032                 kpbl2dx = kmix2dx
1033              endif
1034           else
1035              kpbl2dx=-kpbl2dx
1036           endif
1037        else
1038           kmix2dx = kte
1039           if(kpbl2dx.ge.0)then
1040              kpbl2dx = ktop(iconv)
1041           else
1042              kpbl2dx = -kpbl2dx
1043           endif
1044        endif
1045     else
1046        kmix2dx = kte
1047        if(kpbl2dx.ge.0)then
1048           kpbl2dx = kmix2dx
1049        else
1050           kpbl2dx = -kpbl2dx
1051        endif
1052     endif
1053     pblx=ZQX(KMIX2DX)
1054     return
1055   end subroutine pblhgt
1056
1057
1058   subroutine roots(a,b,c,r1,r2)
1059     implicit none
1060     real a,b,c,r1,r2,q
1061     if(a.eq.0)then            ! form b*x + c = 0
1062        if(b.eq.0)then         ! failure: c = 0
1063           r1 = -9.99e33
1064        else                   ! b*x + c = 0
1065           r1 = -c / b
1066        endif
1067        r2 = r1
1068     else
1069        if(b.eq.0.)then        ! form a*x**2 + c = 0
1070           if(a*c.gt.0.)then   ! failure: x**2 = -c/a < 0
1071              r1 =  -9.99e33
1072           else                ! x**2 = -c/a
1073              r1 = sqrt(-c/a)
1074           endif
1075           r2 = -r1
1076        else                   ! form a*x**2 + b*x + c = 0
1077           if((b**2 - 4*a*c).lt.0.)then ! failure, no real roots
1078              r1 =  -9.99e33
1079              r2 = -r1
1080           else
1081              q = - 0.5 * ( b + sign(1.0,b) * &
1082                   sqrt(b**2 - 4*a*c) )
1083              r1 = q/a
1084              r2 = c/q
1085           endif
1086        endif
1087     endif
1088     return
1089   end subroutine roots
1090
1091   subroutine n2(thlx,exnerfl,epop,qwx,cpoxlv,rdza,        &
1092        rexnerfl, kts,kte, nsquar,rcldb,xlvocp, svp1pa)
1093     implicit none
1094     !     Input/output variables
1095     real thlx(kts:kte),exnerfl(kts:kte+1),epop(kts:kte+1),qwx(kts:kte), &
1096          rexnerfl(kts:kte+1),rdza(kts:kte+1),nsquar(kts:kte+1),rcldb(kts:kte+1)
1097     real cpoxlv,xlvocp,svp1pa
1098     !     Local variables
1099     real templ,rvls,temps,tempv,tvbl,rcld,tvab,thvxfl,dtvdz
1100     integer k,kts,kte
1101
1102     DO K = 2, KTE
1103        !     Buoyancy is jump in thetav across flux level/dza
1104        !     First, layer below, go up and see if anything condenses.
1105        templ = thlx(k)*exnerfl(k)
1106        rvls  = 100.*svp1pa*EXP(svp2*(templ-svpt0)/(templ-svp3))*epop(k)
1107        temps=templ + (qwx(k)-rvls)/(cpoxlv + &
1108             ep_2*xlv*rvls/(r_d*templ**2))
1109        rvls  = 100.*svp1pa*EXP(svp2*(temps-svpt0)/(temps-svp3))*epop(k)
1110        rcldb(k)=max(qwx(k)-rvls,0.)
1111        tempv = (templ + xlvocp*rcldb(k)) * &
1112             (1 + ep_1*(qwx(k)-rcldb(k)) - rcldb(k))
1113        tvbl = tempv*rexnerfl(k)
1114        !     Now do layer above; go down to see how much evaporates
1115        templ = thlx(k-1)*exnerfl(k)
1116        rvls  = 100.*svp1pa*EXP(svp2*(templ-svpt0)/(templ-svp3))*epop(k)
1117        temps=templ+(qwx(k-1)-rvls)/(cpoxlv+ &
1118             ep_2*xlv*rvls/(r_d*templ**2))
1119        rvls  = 100.*svp1pa*EXP(svp2*(temps-svpt0)/(temps-svp3))*epop(k)
1120        rcld=max(qwx(k-1)-rvls,0.)
1121        tempv = (templ + xlvocp*rcld) * &
1122             (1 + ep_1*(qwx(k-1)-rcld) - rcld)
1123        tvab = tempv*rexnerfl(k)
1124        !
1125        thvxfl= .5 * (tvab+tvbl)
1126        dtvdz = (tvab - tvbl) *rdza(k)
1127        nsquar(k) = g/thvxfl * dtvdz
1128     ENDDO
1129     nsquar(1)=nsquar(2)
1130     return
1131   end subroutine n2
1132
1133   subroutine my( &
1134        ! Input variables &
1135     bbls,nsquar,tke,iconv,ktop,thlx,thx,qwx,     &
1136          xlvocp,rcldb,rexnerhl,aone,atwo,kts,kte, &
1137          ! Output variables &
1138          kzm,kth,kethl)
1139     implicit none
1140     real bbls(kts:kte+1),nsquar(kts:kte+1),tke(kts:kte+1),thlx(kts:kte),thx(kts:kte), &
1141          qwx(kts:kte),kzm(kts:kte+1),kth(kts:kte+1),kethl(kts:kte)
1142     real xlvocp,rcldb(kts:kte+1),rexnerhl(kts:kte),aone,atwo
1143     integer iconv,ktop(kts:kte)
1144     !     Local variables
1145     real gh,a1,b1,c1,a2,b2,a1ob1,nuk,delthvl,elambda,bige,biga,kmax,n2min
1146     real sm(kts:kte+1),sh(kts:kte+1)
1147     integer k,ilay,kts,kte
1148     parameter(A1=0.92,B1=16.6,C1=0.08,A2=0.74,B2=10.1)
1149     parameter(nuk=5.0)        !multiplier for kethl, GB01 value
1150
1151     parameter(kmax=1000.)      ! max for kethl
1152     parameter(n2min=1.e-7)     ! min for n2
1153     a1ob1 = a1/b1
1154
1155     !     Calculate the diffusivities for momentum, thetal, qw and tke
1156     !     KTH and KZM are at full levels.
1157     !     KETHL is at half levels. NT: BMG03:6/7
1158     kzm(kte+1)=0.
1159     kth(kte+1)=0.
1160     kzm(1)=0.
1161     kth(1)=0.
1162     sm(kte+1)=1.
1163     sh(kte+1)=1.
1164     sm(1) = 1.
1165     sh(1) = 1.
1166     DO K = KTE, 2, -1
1167 ! According on how the TKE equiation is computed, the "tke" variable in gbmpbl !NP
1168 ! is equivalent to "e" in GB01 paper, or e==(q^2 / 2)  in MY82 or Gal88.       !NP
1169 ! Calculations of gh (MY82, Eq.33(b) and Gal88, Eq.26) contain q^2;            !NP
1170 ! therefore, gh in the code below should contain (2*tke)                       !NP
1171        gh =-bbls(k)*bbls(k)*nsquar(k)/(2*tke(k)+1e-9) ! MY82:33b/ Gal88:26
1172        gh = min(gh,0.0233)
1173        gh = max(gh,-0.28)  ! according to Gal88:30;  -0.28 < gh < 0.0233
1174        sm(k) = a1 * (1. - 3.*c1 - 6.*a1ob1 - 3.*a2*gh* &
1175             ((b2-3.*a2)*(1.-6.*a1ob1) - 3.*c1*(b2+6.*a1))) / &
1176             ((1. - 3.*a2*gh*(6.*a1 + b2)) * (1. - 9.*a1*a2*gh))
1177        sh(k) = a2 * (1. - 6.*a1ob1) / (1. - 3.*a2*gh*(6.*a1+b2))
1178 ! See the previous comments regarding the convention of "e", "q^2/2" and "tke".!NP
1179 ! In order to use "q" in the calculations of kzm and kzh, apply "sqrt(2*tke)"  !NP
1180        kzm(k) = min(bbls(k)*sqrt(2*tke(k))*sm(k),kmax)
1181        kth(k) = min(bbls(k)*sqrt(2*tke(k))*sh(k),kmax)
1182 !
1183        kethl(k)=nuk*sqrt(kzm(k)*kzm(k+1)) ! GB01:A7
1184        kethl(k)=min(kethl(k),kmax) ! keep in bounds
1185     ENDDO
1186     !     Special case for tops of convective layers
1187     ! NT GB01 calls this local entrainment closure : 28
1188      do ilay = 1,iconv
1189       k = ktop(ilay)
1190       IF(nsquar(k).ge.n2min)THEN !NTnew
1191       kethl(k  )=nuk*kzm(k+1)
1192       if(k.ge.3)then !NT only in combination with previous n2 test.
1193          kethl(k-1)=0.0   !NT no transport through inversion
1194          delthvl = (thlx(k-2)+thx(k-2)*ep_1*qwx(k-2)) - &   !NT GB01:18,B12
1195               (thlx(k) + thx(k) * ep_1 * qwx(k))
1196          elambda = xlvocp*rcldb(k)*rexnerhl(k)/max(delthvl,0.1)
1197          bige = 0.8 * elambda
1198          biga = aone * (1 + atwo * bige)
1199          kth(k) = min(kth(k), biga*sqrt(TKE(k)**3)/NSQUAR(k)/ &
1200               max(bbls(k),bbls(k+1)))
1201          kzm(k) = min(kzm(k), kth(k) / sh(k+1) * sm(k+1)) ! prandtl number from layer below
1202         endif
1203     ENDIF !NTnew end
1204    enddo
1205      !     Need KETHL at the top (top-down indexing)
1206      KETHL(1) = kethl(2)
1207      !     Replace KETHL at the surface with something non-zero (top-down indexing)
1208      kethl(kte) = nuk*0.5*kzm(kte)
1209     return
1210   end subroutine my
1211
1212   SUBROUTINE gbmpblinit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,        &
1213        RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR,       &
1214        TKE_PBL,KTH_GBM,                             &
1215        restart,allowed_to_read,                     &
1216        ids, ide, jds, jde, kds, kde,                &
1217        ims, ime, jms, jme, kms, kme,                &
1218        its, ite, jts, jte, kts, kte                 )
1219     !-------------------------------------------------------------------
1220     IMPLICIT NONE
1221     !-------------------------------------------------------------------
1222     LOGICAL , INTENT(IN)          :: restart,allowed_to_read
1223     INTEGER , INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
1224          ims, ime, jms, jme, kms, kme, &
1225          its, ite, jts, jte, kts, kte
1226     INTEGER , INTENT(IN)          ::  P_QI,P_FIRST_SCALAR
1227
1228     REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::         &
1229          RUBLTEN, &
1230          RVBLTEN, &
1231          RTHBLTEN, &
1232          RQVBLTEN, &
1233          RQCBLTEN, &
1234          RQIBLTEN, &
1235          TKE_PBL,  &
1236          KTH_GBM
1237     INTEGER :: i, j, k, itf, jtf, ktf
1238
1239     jtf=min0(jte,jde-1)
1240     ktf=min0(kte,kde-1)
1241     itf=min0(ite,ide-1)
1242
1243     IF(.not.restart)THEN
1244        DO j=jts,jtf
1245           DO k=kts,ktf
1246              DO i=its,itf
1247                 RUBLTEN(i,k,j)=0.
1248                 RVBLTEN(i,k,j)=0.
1249                 RTHBLTEN(i,k,j)=0.
1250                 RQVBLTEN(i,k,j)=0.
1251                 RQCBLTEN(i,k,j)=0.
1252                 TKE_PBL(i,k,j)=0.001 ! use array for TKE
1253                 KTH_GBM(i,k,j)=0.
1254              ENDDO
1255           ENDDO
1256        ENDDO
1257     ENDIF
1258
1259     IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
1260        DO j=jts,jtf
1261           DO k=kts,ktf
1262              DO i=its,itf
1263                 RQIBLTEN(i,k,j)=0.
1264              ENDDO
1265           ENDDO
1266        ENDDO
1267     ENDIF
1268
1269   END SUBROUTINE gbmpblinit
1270 !-------------------------------------------------------------------
1271
1272 END MODULE module_bl_gbmpbl
1273