phys/module_bl_ysu.F

   1 !=================================================================================================================
   2 !module_bl_ysu.F was modified to accomodate both the WRF and MPAS models / 2018-12-7
   3 !=================================================================================================================
   4 !WRF:model_layer:physics
   5 !
   6 !
   7 !
   8 !
   9 !
  10 !
  11 !
  12 module module_bl_ysu
  13 contains
  14 !
  15 !
  16 !-------------------------------------------------------------------------------
  17 !
  18    subroutine ysu(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d,               &
  19                   rublten,rvblten,rthblten,                                    &
  20                   rqvblten,rqcblten,rqiblten,flag_qi,                          &
  21                   cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv,                    &
  22                   dz8w,psfc,                                                   &
  23                   znt,ust,hpbl,psim,psih,                                      &
  24                   xland,hfx,qfx,wspd,br,                                       &
  25                   dt,kpbl2d,                                                   &
  26                   exch_h,exch_m,                                               &
  27                   wstar,delta,                                                 &
  28                   u10,v10,                                                     &
  29                   uoce,voce,                                                   &
  30                   rthraten,ysu_topdown_pblmix,                                 &
  31                   ctopo,ctopo2,                                                &
  32                   idiff,flag_bep,frc_urb2d,                                    &
  33                   a_u_bep,a_v_bep,a_t_bep,                                     &
  34                   a_q_bep,                                                     &
  35                   a_e_bep,b_u_bep,b_v_bep,                                     &
  36                   b_t_bep,b_q_bep,                                             &
  37                   b_e_bep,dlg_bep,                                             &
  38                   dl_u_bep,sf_bep,vl_bep,                                      &
  39                   ids,ide, jds,jde, kds,kde,                                   &
  40                   ims,ime, jms,jme, kms,kme,                                   &
  41                   its,ite, jts,jte, kts,kte,                                   &
  42                 !optional
  43                   regime                                                       &
  44                  )
  45 !-------------------------------------------------------------------------------
  46    implicit none
  47 !-------------------------------------------------------------------------------
  48 !-- u3d         3d u-velocity interpolated to theta points (m/s)
  49 !-- v3d         3d v-velocity interpolated to theta points (m/s)
  50 !-- th3d        3d potential temperature (k)
  51 !-- t3d         temperature (k)
  52 !-- qv3d        3d water vapor mixing ratio (kg/kg)
  53 !-- qc3d        3d cloud mixing ratio (kg/kg)
  54 !-- qi3d        3d ice mixing ratio (kg/kg)
  55 !               (note: if P_QI<PARAM_FIRST_SCALAR this should be zero filled)
  56 !-- p3d         3d pressure (pa)
  57 !-- p3di        3d pressure (pa) at interface level
  58 !-- pi3d        3d exner function (dimensionless)
  59 !-- rr3d        3d dry air density (kg/m^3)
  60 !-- rublten     u tendency due to
  61 !               pbl parameterization (m/s/s)
  62 !-- rvblten     v tendency due to
  63 !               pbl parameterization (m/s/s)
  64 !-- rthblten    theta tendency due to
  65 !               pbl parameterization (K/s)
  66 !-- rqvblten    qv tendency due to
  67 !               pbl parameterization (kg/kg/s)
  68 !-- rqcblten    qc tendency due to
  69 !               pbl parameterization (kg/kg/s)
  70 !-- rqiblten    qi tendency due to
  71 !               pbl parameterization (kg/kg/s)
  72 !-- cp          heat capacity at constant pressure for dry air (j/kg/k)
  73 !-- g           acceleration due to gravity (m/s^2)
  74 !-- rovcp       r/cp
  75 !-- rd          gas constant for dry air (j/kg/k)
  76 !-- rovg        r/g
  77 !-- dz8w        dz between full levels (m)
  78 !-- xlv         latent heat of vaporization (j/kg)
  79 !-- rv          gas constant for water vapor (j/kg/k)
  80 !-- psfc        pressure at the surface (pa)
  81 !-- znt         roughness length (m)
  82 !-- ust         u* in similarity theory (m/s)
  83 !-- hpbl        pbl height (m)
  84 !-- psim        similarity stability function for momentum
  85 !-- psih        similarity stability function for heat
  86 !-- xland       land mask (1 for land, 2 for water)
  87 !-- hfx         upward heat flux at the surface (w/m^2)
  88 !-- qfx         upward moisture flux at the surface (kg/m^2/s)
  89 !-- wspd        wind speed at lowest model level (m/s)
  90 !-- u10         u-wind speed at 10 m (m/s)
  91 !-- v10         v-wind speed at 10 m (m/s)
  92 !-- uoce        sea surface zonal currents (m s-1)
  93 !-- voce        sea surface meridional currents (m s-1)
  94 !-- br          bulk richardson number in surface layer
  95 !-- dt          time step (s)
  96 !-- rvovrd      r_v divided by r_d (dimensionless)
  97 !-- ep1         constant for virtual temperature (r_v/r_d - 1)
  98 !-- ep2         constant for specific humidity calculation
  99 !-- karman      von karman constant
 100 !-- ids         start index for i in domain
 101 !-- ide         end index for i in domain
 102 !-- jds         start index for j in domain
 103 !-- jde         end index for j in domain
 104 !-- kds         start index for k in domain
 105 !-- kde         end index for k in domain
 106 !-- ims         start index for i in memory
 107 !-- ime         end index for i in memory
 108 !-- jms         start index for j in memory
 109 !-- jme         end index for j in memory
 110 !-- kms         start index for k in memory
 111 !-- kme         end index for k in memory
 112 !-- its         start index for i in tile
 113 !-- ite         end index for i in tile
 114 !-- jts         start index for j in tile
 115 !-- jte         end index for j in tile
 116 !-- kts         start index for k in tile
 117 !-- kte         end index for k in tile
 118 !-- idiff       diff3d BEP/BEM+BEM diffusion flag
 119 !-- flag_bep    flag to use BEP/BEP+BEM
 120 !-- frc_urb2d   urban fraction
 121 !-- a_u_bep     BEP/BEP+BEM implicit component u-mom
 122 !-- a_v_bep     BEP/BEP+BEM implicit component v-mom
 123 !-- a_t_bep     BEP/BEP+BEM implicit component pot. temp.
 124 !-- a_q_bep     BEP/BEP+BEM implicit component vapor mixing ratio
 125 !-- a_e_bep     BEP/BEP+BEM implicit component TKE
 126 !-- b_u_bep     BEP/BEP+BEM explicit component u-mom
 127 !-- b_v_bep     BEP/BEP+BEM explicit component v-mom
 128 !-- b_t_bep     BEP/BEP+BEM explicit component pot.temp.
 129 !-- b_q_bep     BEP/BEP+BEM explicit component vapor mixing ratio
 130 !-- b_e_bep     BEP/BEP+BEM explicit component TKE
 131 !-- dlg_bep     Height above ground Martilli et al. (2002) Eq. 24
 132 !-- dl_u_bep    modified length scale Martilli et al. (2002) Eq. 22
 133 !-- sf_bep      fraction of vertical surface not occupied by buildings
 134 !-- vl_bep      volume fraction of grid cell not occupied by buildings
 135 !-------------------------------------------------------------------------------
 136 !
 137    integer,parameter ::  ndiff = 3
 138    real,parameter    ::  rcl = 1.0
 139 !
 140    integer,  intent(in   )   ::      ids,ide, jds,jde, kds,kde,                &
 141                                      ims,ime, jms,jme, kms,kme,                &
 142                                      its,ite, jts,jte, kts,kte
 143
 144    integer,  intent(in)      ::      ysu_topdown_pblmix
 145 !
 146    real,     intent(in   )   ::      dt,cp,g,rovcp,rovg,rd,xlv,rv
 147 !
 148    real,     intent(in )     ::      ep1,ep2,karman
 149 !
 150    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &
 151              intent(in   )   ::                                          qv3d, &
 152                                                                          qc3d, &
 153                                                                          qi3d, &
 154                                                                           p3d, &
 155                                                                          pi3d, &
 156                                                                          th3d, &
 157                                                                           t3d, &
 158                                                                          dz8w, &
 159                                                                      rthraten
 160    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &
 161              intent(in   )   ::                                          p3di
 162 !
 163    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &
 164              intent(inout)   ::                                       rublten, &
 165                                                                       rvblten, &
 166                                                                      rthblten, &
 167                                                                      rqvblten, &
 168                                                                      rqcblten
 169 !
 170    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &
 171              intent(inout)   ::                                        exch_h, &
 172                                                                        exch_m
 173    real,     dimension( ims:ime, jms:jme )                                   , &
 174              intent(inout)   ::                                         wstar
 175    real,     dimension( ims:ime, jms:jme )                                   , &
 176              intent(inout)   ::                                         delta
 177    real,     dimension( ims:ime, jms:jme )                                   , &
 178              intent(inout)   ::                                           u10, &
 179                                                                           v10
 180    real,     dimension( ims:ime, jms:jme )                                   , &
 181              intent(in   )   ::                                          uoce, &
 182                                                                          voce
 183 !
 184    real,     dimension( ims:ime, jms:jme )                                   , &
 185              intent(in   )   ::                                         xland, &
 186                                                                           hfx, &
 187                                                                           qfx, &
 188                                                                            br, &
 189                                                                          psfc
 190    real,     dimension( ims:ime, jms:jme )                                   , &
 191              intent(in   )   ::                                                &
 192                                                                          psim, &
 193                                                                          psih
 194    real,     dimension( ims:ime, jms:jme )                                   , &
 195              intent(inout)   ::                                           znt, &
 196                                                                           ust, &
 197                                                                          hpbl, &
 198                                                                           wspd
 199 !
 200    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &
 201              intent(in   )   ::                                           u3d, &
 202                                                                           v3d
 203 !
 204    integer,  dimension( ims:ime, jms:jme )                                   , &
 205              intent(out  )   ::                                        kpbl2d
 206    logical,  intent(in)      ::                                       flag_qi
 207    integer, intent(in) ::                                               idiff
 208    logical, intent(in) ::                                            flag_bep
 209    real,dimension(ims:ime,kms:kme,jms:jme),intent(in) ::              a_u_bep, &
 210                                                               a_v_bep,a_t_bep, &
 211                                                               a_e_bep,b_u_bep, &
 212                                                               a_q_bep,b_q_bep, &
 213                                                               b_v_bep,b_t_bep, &
 214                                                               b_e_bep,dlg_bep, &
 215                                                                      dl_u_bep, &
 216                                                                 vl_bep,sf_bep
 217    real, dimension(ims:ime,jms:jme),intent(in) :: frc_urb2d
 218 !
 219 !optional
 220 !
 221    real,     dimension( ims:ime, jms:jme )                                   , &
 222              optional                                                        , &
 223              intent(inout)   ::                                        regime
 224 !
 225    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &
 226              optional                                                        , &
 227              intent(inout)   ::                                       rqiblten
 228 !
 229    real,     dimension( ims:ime, jms:jme )                                   , &
 230              optional                                                        , &
 231              intent(in   )   ::                                         ctopo, &
 232                                                                        ctopo2
 233 !local
 234    integer ::  i,j,k
 235    real,     dimension( its:ite, kts:kte*ndiff )  ::                 rqvbl2dt, &
 236                                                                          qv2d
 237    real,     dimension( its:ite, kts:kte )  ::                            pdh
 238    real,     dimension( its:ite, kts:kte+1 )  ::                         pdhi
 239    real,     dimension( its:ite )  ::                                          &
 240                                                                         dusfc, &
 241                                                                         dvsfc, &
 242                                                                         dtsfc, &
 243                                                                         dqsfc
 244    real,dimension(its:ite,kts:kte,jts:jte) :: a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e, &
 245                                                      a_q,b_q,dlg,dl_u,sfk,vlk
 246    real,dimension(its:ite,jts:jte)         :: frcurb
 247    real :: bepswitch ! 0 if not using bep or bep+bem, 1 if using
 248 !
 249    qv2d(its:ite,:) = 0.0
 250 !
 251    bepswitch = 0.0
 252    a_u(:,:,:)=0.0
 253    a_v(:,:,:)=0.0
 254    a_t(:,:,:)=0.0
 255    a_q(:,:,:)=0.0
 256    a_e(:,:,:)=0.0
 257    b_u(:,:,:)=0.0
 258    b_v(:,:,:)=0.0
 259    b_t(:,:,:)=0.0
 260    b_q(:,:,:)=0.0
 261    b_e(:,:,:)=0.0
 262    sfk(:,:,:)=1.0
 263    vlk(:,:,:)=1.0
 264    dl_u(:,:,:)=0.0
 265    dlg(:,:,:)=0.0
 266    frcurb(:,:)=0.0
 267
 268    do j = jts,jte
 269       do k = kts,kte+1
 270         do i = its,ite
 271           if(k.le.kte)pdh(i,k) = p3d(i,k,j)
 272           pdhi(i,k) = p3di(i,k,j)
 273         enddo
 274       enddo
 275
 276       do k = kts,kte
 277         do i = its,ite
 278           qv2d(i,k) = qv3d(i,k,j)
 279           qv2d(i,k+kte) = qc3d(i,k,j)
 280           if(flag_qi) qv2d(i,k+kte+kte) = qi3d(i,k,j)
 281         enddo
 282       enddo
 283
 284       if(flag_bep) then
 285         bepswitch=1.0
 286         do k=kts,kte
 287           do i=its,ite
 288             a_u(i,k,j)=a_u_bep(i,k,j)
 289             a_v(i,k,j)=a_v_bep(i,k,j)
 290             a_t(i,k,j)=a_t_bep(i,k,j)
 291             a_q(i,k,j)=a_q_bep(i,k,j)
 292             a_e(i,k,j)=a_e_bep(i,k,j)
 293             b_u(i,k,j)=b_u_bep(i,k,j)
 294             b_v(i,k,j)=b_v_bep(i,k,j)
 295             b_t(i,k,j)=b_t_bep(i,k,j)
 296             b_q(i,k,j)=b_q_bep(i,k,j)
 297             b_e(i,k,j)=b_e_bep(i,k,j)
 298             sfk(i,k,j)=sf_bep(i,k,j)
 299             vlk(i,k,j)=vl_bep(i,k,j)
 300             dl_u(i,k,j)=dl_u_bep(i,k,j)
 301             dlg(i,k,j)=dlg_bep(i,k,j)
 302             frcurb(i,j)=frc_urb2d(i,j)
 303           enddo
 304         enddo
 305       endif
 306 !
 307       call ysu2d(J=j,ux=u3d(ims,kms,j),vx=v3d(ims,kms,j)                       &
 308               ,tx=t3d(ims,kms,j)                                               &
 309               ,qx=qv2d(its,kts)                                                &
 310               ,p2d=pdh(its,kts),p2di=pdhi(its,kts)                             &
 311               ,pi2d=pi3d(ims,kms,j)                                            &
 312               ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j)                 &
 313               ,ttnp=rthblten(ims,kms,j),qtnp=rqvbl2dt(its,kts),ndiff=ndiff     &
 314               ,cp=cp,g=g,rovcp=rovcp,rd=rd,rovg=rovg                           &
 315               ,xlv=xlv,rv=rv                                                   &
 316               ,ep1=ep1,ep2=ep2,karman=karman                                   &
 317               ,dz8w2d=dz8w(ims,kms,j)                                          &
 318               ,psfcpa=psfc(ims,j),znt=znt(ims,j),ust=ust(ims,j)                &
 319               ,hpbl=hpbl(ims,j)                                                &
 320               ,regime=regime(ims,j),psim=psim(ims,j)                           &
 321               ,psih=psih(ims,j),xland=xland(ims,j)                             &
 322               ,hfx=hfx(ims,j),qfx=qfx(ims,j)                                   &
 323               ,wspd=wspd(ims,j),br=br(ims,j)                                   &
 324               ,dusfc=dusfc,dvsfc=dvsfc,dtsfc=dtsfc,dqsfc=dqsfc                 &
 325               ,dt=dt,rcl=1.0,kpbl1d=kpbl2d(ims,j)                              &
 326               ,exch_hx=exch_h(ims,kms,j)                                       &
 327               ,exch_mx=exch_m(ims,kms,j)                                       &
 328               ,wstar=wstar(ims,j)                                              &
 329               ,delta=delta(ims,j)                                              &
 330               ,u10=u10(ims,j),v10=v10(ims,j)                                   &
 331               ,uox=uoce(ims,j),vox=voce(ims,j)                                 &
 332               ,rthraten=rthraten(ims,kms,j),p2diORG=p3di(ims,kms,j)            &
 333               ,ysu_topdown_pblmix=ysu_topdown_pblmix                           &
 334               ,ctopo=ctopo(ims,j),ctopo2=ctopo2(ims,j)                         &
 335               ,a_u2d=a_u(its,kts,j), a_v2d=a_v(its,kts,j)                      &
 336               ,a_t2d=a_t(its,kts,j), a_q2d=a_q(its,kts,j)                      &
 337               ,b_u2d=b_u(its,kts,j), b_v2d=b_v(its,kts,j)                      &
 338               ,b_t2d=b_t(its,kts,j), b_q2d=b_q(its,kts,j)                      &
 339               ,b_e2d=b_e(its,kts,j), a_e2d=a_e(its,kts,j)                      &
 340               ,sfk2d=sfk(its,kts,j), vlk2d=vlk(its,kts,j)                      &
 341               ,dlu2d=dl_u(its,kts,j), dlg2d=dlg(its,kts,j)                     &
 342               ,frc_urb1d=frcurb(its,j), bepswitch=bepswitch                    &
 343               ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde               &
 344               ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme               &
 345               ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte   )
 346 !
 347      do k = kts,kte
 348        do i = its,ite
 349          rthblten(i,k,j) = rthblten(i,k,j)/pi3d(i,k,j)
 350          rqvblten(i,k,j) = rqvbl2dt(i,k)
 351          rqcblten(i,k,j) = rqvbl2dt(i,k+kte)
 352          if(flag_qi) rqiblten(i,k,j) = rqvbl2dt(i,k+kte+kte)
 353        enddo
 354      enddo
 355 !
 356    enddo
 357 !
 358    end subroutine ysu
 359 !
 360 !-------------------------------------------------------------------------------
 361 !
 362    subroutine ysu2d(j,ux,vx,tx,qx,p2d,p2di,pi2d,                               &
 363                   utnp,vtnp,ttnp,qtnp,ndiff,                                   &
 364                   cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv,                    &
 365                   dz8w2d,psfcpa,                                               &
 366                   znt,ust,hpbl,psim,psih,                                      &
 367                   xland,hfx,qfx,wspd,br,                                       &
 368                   dusfc,dvsfc,dtsfc,dqsfc,                                     &
 369                   dt,rcl,kpbl1d,                                               &
 370                   exch_hx,exch_mx,                                             &
 371                   wstar,delta,                                                 &
 372                   u10,v10,                                                     &
 373                   uox,vox,                                                     &
 374                   rthraten,p2diORG,                                            &
 375                   ysu_topdown_pblmix,                                          &
 376                   ctopo,ctopo2,                                                &
 377                   a_u2d, a_v2d, a_t2d, a_q2d,                                  &
 378                   b_u2d, b_v2d, b_t2d, b_q2d,                                  &
 379                   b_e2d, a_e2d, sfk2d, vlk2d,                                  &
 380                   dlu2d, dlg2d,                                                &
 381                   frc_urb1d, bepswitch,                                        &
 382                   ids,ide, jds,jde, kds,kde,                                   &
 383                   ims,ime, jms,jme, kms,kme,                                   &
 384                   its,ite, jts,jte, kts,kte,                                   &
 385                 !optional
 386                   regime                                                       &
 387                    )
 388 !-------------------------------------------------------------------------------
 389    implicit none
 390 !-------------------------------------------------------------------------------
 391 !
 392 !     this code is a revised vertical diffusion package ("ysupbl")
 393 !     with a nonlocal turbulent mixing in the pbl after "mrfpbl".
 394 !     the ysupbl (hong et al. 2006) is based on the study of noh
 395 !     et al.(2003) and accumulated realism of the behavior of the
 396 !     troen and mahrt (1986) concept implemented by hong and pan(1996).
 397 !     the major ingredient of the ysupbl is the inclusion of an explicit
 398 !     treatment of the entrainment processes at the entrainment layer.
 399 !     this routine uses an implicit approach for vertical flux
 400 !     divergence and does not require "miter" timesteps.
 401 !     it includes vertical diffusion in the stable atmosphere
 402 !     and moist vertical diffusion in clouds.
 403 !
 404 !     mrfpbl:
 405 !     coded by song-you hong (ncep), implemented by jimy dudhia (ncar)
 406 !              fall 1996
 407 !
 408 !     ysupbl:
 409 !     coded by song-you hong (yonsei university) and implemented by
 410 !              song-you hong (yonsei university) and jimy dudhia (ncar)
 411 !              summer 2002
 412 !
 413 !     further modifications :
 414 !              an enhanced stable layer mixing, april 2008
 415 !              ==> increase pbl height when sfc is stable (hong 2010)
 416 !              pressure-level diffusion, april 2009
 417 !              ==> negligible differences
 418 !              implicit forcing for momentum with clean up, july 2009
 419 !              ==> prevents model blowup when sfc layer is too low
 420 !              incresea of lamda, maximum (30, 0.1 x del z) feb 2010
 421 !              ==> prevents model blowup when delz is extremely large
 422 !              revised prandtl number at surface, peggy lemone, feb 2010
 423 !              ==> increase kh, decrease mixing due to counter-gradient term
 424 !              revised thermal, shin et al. mon. wea. rev. , songyou hong, aug 2011
 425 !              ==> reduce the thermal strength when z1 < 0.1 h
 426 !              revised prandtl number for free convection, dudhia, mar 2012
 427 !              ==> pr0 = 1 + bke (=0.272) when neutral, kh is reduced
 428 !              minimum kzo = 0.01, lo = min (30m,delz), hong, mar 2012
 429 !              ==> weaker mixing when stable, and les resolution in vertical
 430 !              gz1oz0 is removed, and psim psih are ln(z1/z0)-psim,h, hong, mar 2012
 431 !              ==> consider thermal z0 when differs from mechanical z0
 432 !              a bug fix in wscale computation in stable bl, sukanta basu, jun 2012
 433 !              ==> wscale becomes small with height, and less mixing in stable bl
 434 !              revision in background diffusion (kzo), jan 2016
 435 !              ==> kzo = 0.1 for momentum and = 0.01 for mass to account for
 436 !                  internal wave mixing of large et al. (1994), songyou hong, feb 2016
 437 !              ==> alleviate superious excessive mixing when delz is large
 438 !              add multilayer urban canopy models of BEP and BEP+BEM, jan 2021
 439 !
 440 !     references:
 441 !
 442 !        hendricks, knievel, and wang (2020), j. appl. meteor. clim.
 443 !        hong (2010) quart. j. roy. met. soc
 444 !        hong, noh, and dudhia (2006), mon. wea. rev.
 445 !        hong and pan (1996), mon. wea. rev.
 446 !        noh, chun, hong, and raasch (2003), boundary layer met.
 447 !        troen and mahrt (1986), boundary layer met.
 448 !
 449 !-------------------------------------------------------------------------------
 450 !
 451    real,parameter    ::  xkzminm = 0.1,xkzminh = 0.01
 452    real,parameter    ::  xkzmin = 0.01,xkzmax = 1000.,rimin = -100.
 453    real,parameter    ::  rlam = 30.,prmin = 0.25,prmax = 4.
 454    real,parameter    ::  brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4
 455    real,parameter    ::  afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0
 456    real,parameter    ::  phifac = 8.,sfcfrac = 0.1
 457    real,parameter    ::  d1 = 0.02, d2 = 0.05, d3 = 0.001
 458    real,parameter    ::  h1 = 0.33333333, h2 = 0.6666667
 459    real,parameter    ::  zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16.
 460    real,parameter    ::  tmin=1.e-2
 461    real,parameter    ::  gamcrt = 3.,gamcrq = 2.e-3
 462    real,parameter    ::  xka = 2.4e-5
 463    integer,parameter ::  imvdif = 1
 464 !
 465    integer,  intent(in   )   ::     ids,ide, jds,jde, kds,kde,                 &
 466                                     ims,ime, jms,jme, kms,kme,                 &
 467                                     its,ite, jts,jte, kts,kte,                 &
 468                                     j,ndiff
 469
 470    integer,  intent(in)      ::     ysu_topdown_pblmix
 471 !
 472    real,     intent(in   )   ::     dt,rcl,cp,g,rovcp,rovg,rd,xlv,rv
 473 !
 474    real,     intent(in )     ::     ep1,ep2,karman
 475 !
 476    real,     dimension( ims:ime, kms:kme ),                                    &
 477              intent(in)      ::                                        dz8w2d, &
 478                                                                          pi2d, &
 479                                                                       p2diorg
 480 !
 481    real,     dimension( ims:ime, kms:kme )                                   , &
 482              intent(in   )   ::                                            tx
 483    real,     dimension( its:ite, kts:kte*ndiff )                             , &
 484              intent(in   )   ::                                            qx
 485 !
 486    real,     dimension( ims:ime, kms:kme )                                   , &
 487              intent(inout)   ::                                          utnp, &
 488                                                                          vtnp, &
 489                                                                          ttnp
 490    real,     dimension( its:ite, kts:kte*ndiff )                             , &
 491              intent(inout)   ::                                          qtnp
 492 !
 493    real,     dimension( its:ite, kts:kte+1 )                                 , &
 494              intent(in   )   ::                                          p2di
 495 !
 496    real,     dimension( its:ite, kts:kte )                                   , &
 497              intent(in   )   ::                                           p2d
 498 !
 499    real,     dimension( ims:ime )                                            , &
 500              intent(inout)   ::                                           ust, &
 501                                                                          hpbl, &
 502                                                                           znt
 503    real,     dimension( ims:ime )                                            , &
 504              intent(in   )   ::                                         xland, &
 505                                                                           hfx, &
 506                                                                           qfx
 507 !
 508    real,     dimension( ims:ime ), intent(inout)   ::                    wspd
 509    real,     dimension( ims:ime ), intent(in  )    ::                      br
 510 !
 511    real,     dimension( ims:ime ), intent(in   )   ::                    psim, &
 512                                                                          psih
 513 !
 514    real,     dimension( ims:ime ), intent(in   )   ::                  psfcpa
 515    integer,  dimension( ims:ime ), intent(out  )   ::                  kpbl1d
 516 !
 517    real,     dimension( ims:ime, kms:kme )                                   , &
 518              intent(in   )   ::                                            ux, &
 519                                                                            vx, &
 520                                                                       rthraten
 521    real,     dimension( ims:ime )                                            , &
 522              optional                                                        , &
 523              intent(in   )   ::                                         ctopo, &
 524                                                                        ctopo2
 525    real,     dimension( ims:ime )                                            , &
 526              optional                                                        , &
 527              intent(inout)   ::                                        regime
 528 !
 529 ! local vars
 530 !
 531    real,     dimension( its:ite, kts:kte ),                                    &
 532              intent(in)      ::                                         a_u2d, &
 533                                                                         a_v2d, &
 534                                                                         a_t2d, &
 535                                                                         a_q2d, &
 536                                                                         b_u2d, &
 537                                                                         b_v2d, &
 538                                                                         b_t2d, &
 539                                                                         b_q2d, &
 540                                                                         b_e2d, &
 541                                                                         a_e2d, &
 542                                                                         sfk2d, &
 543                                                                         vlk2d, &
 544                                                                         dlu2d, &
 545                                                                         dlg2d
 546
 547    real,     dimension( its:ite ),                                             &
 548              intent(in)      ::                                      frc_urb1d
 549    real ::                                                           bepswitch
 550    real,     dimension( its:ite )            ::                           hol
 551    real,     dimension( its:ite, kts:kte+1 ) ::                            zq
 552 !
 553    real,     dimension( its:ite, kts:kte )   ::                                &
 554                                                                thx,thvx,thlix, &
 555                                                                           del, &
 556                                                                           dza, &
 557                                                                           dzq, &
 558                                                                         xkzom, &
 559                                                                         xkzoh, &
 560                                                                            za
 561 !
 562    real,    dimension( its:ite )             ::                                &
 563                                                                          rhox, &
 564                                                                        govrth, &
 565                                                                   zl1,thermal, &
 566                                                                        wscale, &
 567                                                                   hgamt,hgamq, &
 568                                                                     brdn,brup, &
 569                                                                     phim,phih, &
 570                                                                   dusfc,dvsfc, &
 571                                                                   dtsfc,dqsfc, &
 572                                                                         prpbl, &
 573                                                               wspd1,thermalli
 574 !
 575    real,    dimension( its:ite, kts:kte )    ::                     xkzm,xkzh, &
 576                                                                         f1,f2, &
 577                                                                         r1,r2, &
 578                                                                         ad,au, &
 579                                                                            cu, &
 580                                                                            al, &
 581                                                                          xkzq, &
 582                                                                          zfac, &
 583                                                                         rhox2, &
 584                                                                        hgamt2, &
 585                                                                      ad1, adm
 586 !
 587 !jdf added exch_hx
 588 !
 589    real,    dimension( ims:ime, kms:kme )                                    , &
 590             intent(inout)   ::                                        exch_hx, &
 591                                                                       exch_mx
 592 !
 593    real,    dimension( ims:ime )                                             , &
 594             intent(inout)    ::                                           u10, &
 595                                                                           v10
 596    real,    dimension( ims:ime )                                             , &
 597             intent(in  )    ::                                            uox, &
 598                                                                           vox
 599    real,    dimension( its:ite )    ::                                         &
 600                                                                          brcr, &
 601                                                                         sflux, &
 602                                                                          zol1, &
 603                                                                     brcr_sbro
 604 !
 605    real,    dimension( its:ite, kts:kte, ndiff)  ::                     r3,f3
 606    integer, dimension( its:ite )             ::                  kpbl,kpblold
 607 !
 608    logical, dimension( its:ite )             ::                        pblflg, &
 609                                                                        sfcflg, &
 610                                                                        stable, &
 611                                                                      cloudflg
 612
 613    logical                                   ::                     definebrup
 614 !
 615    integer ::  n,i,k,l,ic,is,kk
 616    integer ::  klpbl, ktrace1, ktrace2, ktrace3
 617 !
 618 !
 619    real    ::  dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0
 620    real    ::  ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri
 621    real    ::  brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz
 622    real    ::  utend,vtend,ttend,qtend
 623    real    ::  dtstep,govrthv
 624    real    ::  cont, conq, conw, conwrc
 625 !
 626
 627    real, dimension( its:ite, kts:kte )     ::                wscalek,wscalek2
 628    real, dimension( ims:ime )              ::                           wstar
 629    real, dimension( ims:ime )              ::                           delta
 630    real, dimension( its:ite, kts:kte )     ::                     xkzml,xkzhl, &
 631                                                                zfacent,entfac
 632    real, dimension( its:ite )              ::                            ust3, &
 633                                                                        wstar3, &
 634                                                                      wstar3_2, &
 635                                                                   hgamu,hgamv, &
 636                                                                       wm2, we, &
 637                                                                        bfxpbl, &
 638                                                                 hfxpbl,qfxpbl, &
 639                                                                 ufxpbl,vfxpbl, &
 640                                                                         dthvx
 641    real    ::  prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx,                           &
 642                dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr,     &
 643                prfac,prfac2,phim8z,radsum,tmp1,templ,rvls,temps,ent_eff,    &
 644                rcldb,bruptmp,radflux,vconvlim,vconvnew,fluxc,vconvc,vconv
 645 !topo-corr
 646    real,    dimension( ims:ime, kms:kme )    ::                          fric, &
 647                                                                        tke_ysu,&
 648                                                                         el_ysu,&
 649                                                                      shear_ysu,&
 650                                                                      buoy_ysu
 651    real,    dimension( ims:ime )             ::                       pblh_ysu,&
 652                                                                       vconvfx
 653 !
 654 !-------------------------------------------------------------------------------
 655 !
 656    klpbl = kte
 657 !
 658    cont=cp/g
 659    conq=xlv/g
 660    conw=1./g
 661    conwrc = conw*sqrt(rcl)
 662    conpr = bfac*karman*sfcfrac
 663 !
 664 !  k-start index for tracer diffusion
 665 !
 666    ktrace1 = 0
 667    ktrace2 = 0 + kte
 668    ktrace3 = 0 + kte*2
 669 !
 670    do k = kts,kte
 671      do i = its,ite
 672        thx(i,k) = tx(i,k)/pi2d(i,k)
 673        thlix(i,k) = (tx(i,k)-xlv*qx(i,ktrace2+k)/cp-2.834E6*qx(i,ktrace3+k)/cp)/pi2d(i,k)
 674      enddo
 675    enddo
 676 !
 677    do k = kts,kte
 678      do i = its,ite
 679        tvcon = (1.+ep1*qx(i,k))
 680        thvx(i,k) = thx(i,k)*tvcon
 681      enddo
 682    enddo
 683 !
 684    do i = its,ite
 685      tvcon = (1.+ep1*qx(i,1))
 686      rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon)
 687      govrth(i) = g/thx(i,1)
 688    enddo
 689 !
 690 !-----compute the height of full- and half-sigma levels above ground
 691 !     level, and the layer thicknesses.
 692 !
 693    do i = its,ite
 694      zq(i,1) = 0.
 695    enddo
 696 !
 697    do k = kts,kte
 698      do i = its,ite
 699        zq(i,k+1) = dz8w2d(i,k)+zq(i,k)
 700        tvcon = (1.+ep1*qx(i,k))
 701        rhox2(i,k) = p2d(i,k)/(rd*tx(i,k)*tvcon)
 702      enddo
 703    enddo
 704 !
 705    do k = kts,kte
 706      do i = its,ite
 707        za(i,k) = 0.5*(zq(i,k)+zq(i,k+1))
 708        dzq(i,k) = zq(i,k+1)-zq(i,k)
 709        del(i,k) = p2di(i,k)-p2di(i,k+1)
 710      enddo
 711    enddo
 712 !
 713    do i = its,ite
 714      dza(i,1) = za(i,1)
 715    enddo
 716 !
 717    do k = kts+1,kte
 718      do i = its,ite
 719        dza(i,k) = za(i,k)-za(i,k-1)
 720      enddo
 721    enddo
 722 !
 723 !
 724 !-----initialize vertical tendencies and
 725 !
 726    utnp(its:ite,:) = 0.
 727    vtnp(its:ite,:) = 0.
 728    ttnp(its:ite,:) = 0.
 729    qtnp(its:ite,:) = 0.
 730 !
 731    do i = its,ite
 732      wspd1(i) = sqrt( (ux(i,1)-uox(i))*(ux(i,1)-uox(i)) + (vx(i,1)-vox(i))*(vx(i,1)-vox(i)) )+1.e-9
 733    enddo
 734 !
 735 !---- compute vertical diffusion
 736 !
 737 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 738 !     compute preliminary variables
 739 !
 740    dtstep = dt
 741    dt2 = 2.*dtstep
 742    rdt = 1./dt2
 743 !
 744    do i = its,ite
 745      bfxpbl(i) = 0.0
 746      hfxpbl(i) = 0.0
 747      qfxpbl(i) = 0.0
 748      ufxpbl(i) = 0.0
 749      vfxpbl(i) = 0.0
 750      hgamu(i)  = 0.0
 751      hgamv(i)  = 0.0
 752      delta(i)  = 0.0
 753      wstar3_2(i) =  0.0
 754    enddo
 755 !
 756    do k = kts,klpbl
 757      do i = its,ite
 758        wscalek(i,k) = 0.0
 759        wscalek2(i,k) = 0.0
 760      enddo
 761    enddo
 762 !
 763    do k = kts,klpbl
 764      do i = its,ite
 765        zfac(i,k) = 0.0
 766      enddo
 767    enddo
 768    do k = kts,klpbl-1
 769      do i = its,ite
 770        xkzom(i,k) = xkzminm
 771        xkzoh(i,k) = xkzminh
 772      enddo
 773    enddo
 774 !
 775    do i = its,ite
 776      dusfc(i) = 0.
 777      dvsfc(i) = 0.
 778      dtsfc(i) = 0.
 779      dqsfc(i) = 0.
 780    enddo
 781 !
 782    do i = its,ite
 783      hgamt(i)  = 0.
 784      hgamq(i)  = 0.
 785      wscale(i) = 0.
 786      kpbl(i)   = 1
 787      hpbl(i)   = zq(i,1)
 788      zl1(i)    = za(i,1)
 789      thermal(i)= thvx(i,1)
 790      thermalli(i) = thlix(i,1)
 791      pblflg(i) = .true.
 792      sfcflg(i) = .true.
 793      sflux(i) = hfx(i)/rhox(i)/cp + qfx(i)/rhox(i)*ep1*thx(i,1)
 794      if(br(i).gt.0.0) sfcflg(i) = .false.
 795    enddo
 796 !
 797 !     compute the first guess of pbl height
 798 !
 799    do i = its,ite
 800      stable(i) = .false.
 801      brup(i) = br(i)
 802      brcr(i) = brcr_ub
 803    enddo
 804 !
 805    do k = 2,klpbl
 806      do i = its,ite
 807        if(.not.stable(i))then
 808          brdn(i) = brup(i)
 809          spdk2   = max(ux(i,k)**2+vx(i,k)**2,1.)
 810          brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
 811          kpbl(i) = k
 812          stable(i) = brup(i).gt.brcr(i)
 813        endif
 814      enddo
 815    enddo
 816 !
 817    do i = its,ite
 818      k = kpbl(i)
 819      if(brdn(i).ge.brcr(i))then
 820        brint = 0.
 821      elseif(brup(i).le.brcr(i))then
 822        brint = 1.
 823      else
 824        brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
 825      endif
 826      hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
 827      if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
 828      if(kpbl(i).le.1) pblflg(i) = .false.
 829    enddo
 830 !
 831    do i = its,ite
 832      fm = psim(i)
 833      fh = psih(i)
 834      zol1(i) = max(br(i)*fm*fm/fh,rimin)
 835      if(sfcflg(i))then
 836        zol1(i) = min(zol1(i),-zfmin)
 837      else
 838        zol1(i) = max(zol1(i),zfmin)
 839      endif
 840      hol1 = zol1(i)*hpbl(i)/zl1(i)*sfcfrac
 841      if(sfcflg(i))then
 842        phim(i) = (1.-aphi16*hol1)**(-1./4.)
 843        phih(i) = (1.-aphi16*hol1)**(-1./2.)
 844        bfx0 = max(sflux(i),0.)
 845        hfx0 = max(hfx(i)/rhox(i)/cp,0.)
 846        qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.)
 847        wstar3(i) = (govrth(i)*bfx0*hpbl(i))
 848        wstar(i) = (wstar3(i))**h1
 849      else
 850        phim(i) = (1.+aphi5*hol1)
 851        phih(i) = phim(i)
 852        wstar(i)  = 0.
 853        wstar3(i) = 0.
 854      endif
 855      ust3(i)   = ust(i)**3.
 856      wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1
 857      wscale(i) = min(wscale(i),ust(i)*aphi16)
 858      wscale(i) = max(wscale(i),ust(i)/aphi5)
 859    enddo
 860 !
 861 !     compute the surface variables for pbl height estimation
 862 !     under unstable conditions
 863 !
 864    do i = its,ite
 865      if(sfcflg(i).and.sflux(i).gt.0.0)then
 866        gamfac   = bfac/rhox(i)/wscale(i)
 867        hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
 868        hgamq(i) = min(gamfac*qfx(i),gamcrq)
 869        vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac
 870        thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
 871        thermalli(i)= thermalli(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
 872        hgamt(i) = max(hgamt(i),0.0)
 873        hgamq(i) = max(hgamq(i),0.0)
 874        brint    = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.)
 875        hgamu(i) = brint*ux(i,1)
 876        hgamv(i) = brint*vx(i,1)
 877      else
 878        pblflg(i) = .false.
 879      endif
 880    enddo
 881 !
 882 !     enhance the pbl height by considering the thermal
 883 !
 884    do i = its,ite
 885      if(pblflg(i))then
 886        kpbl(i) = 1
 887        hpbl(i) = zq(i,1)
 888      endif
 889    enddo
 890 !
 891    do i = its,ite
 892      if(pblflg(i))then
 893        stable(i) = .false.
 894        brup(i) = br(i)
 895        brcr(i) = brcr_ub
 896      endif
 897    enddo
 898 !
 899    do k = 2,klpbl
 900      do i = its,ite
 901        if(.not.stable(i).and.pblflg(i))then
 902          brdn(i) = brup(i)
 903          spdk2   = max(ux(i,k)**2+vx(i,k)**2,1.)
 904          brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
 905          kpbl(i) = k
 906          stable(i) = brup(i).gt.brcr(i)
 907        endif
 908      enddo
 909    enddo
 910 !
 911 !     enhance pbl by theta-li
 912 !
 913    if (ysu_topdown_pblmix.eq.1)then
 914      do i = its,ite
 915         kpblold(i) = kpbl(i)
 916         definebrup=.false.
 917         do k = kpblold(i), kte-1
 918            spdk2   = max(ux(i,k)**2+vx(i,k)**2,1.)
 919            bruptmp = (thlix(i,k)-thermalli(i))*(g*za(i,k)/thlix(i,1))/spdk2
 920            stable(i) = bruptmp.ge.brcr(i)
 921            if (definebrup) then
 922            kpbl(i) = k
 923            brup(i) = bruptmp
 924            definebrup=.false.
 925            endif
 926            if (.not.stable(i)) then !overwrite brup brdn values
 927            brdn(i)=bruptmp
 928            definebrup=.true.
 929            pblflg(i)=.true.
 930            endif
 931         enddo
 932      enddo
 933    endif
 934
 935    do i = its,ite
 936      if(pblflg(i)) then
 937        k = kpbl(i)
 938        if(brdn(i).ge.brcr(i))then
 939          brint = 0.
 940        elseif(brup(i).le.brcr(i))then
 941          brint = 1.
 942        else
 943          brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
 944        endif
 945        hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
 946        if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
 947        if(kpbl(i).le.1) pblflg(i) = .false.
 948      endif
 949    enddo
 950 !
 951 !     stable boundary layer
 952 !
 953    do i = its,ite
 954      if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
 955        brup(i) = br(i)
 956        stable(i) = .false.
 957      else
 958        stable(i) = .true.
 959      endif
 960    enddo
 961 !
 962    do i = its,ite
 963      if((.not.stable(i)).and.((xland(i)-1.5).ge.0))then
 964        wspd10 = u10(i)*u10(i) + v10(i)*v10(i)
 965        wspd10 = sqrt(wspd10)
 966        ross = wspd10 / (cori*znt(i))
 967        brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3)
 968      endif
 969    enddo
 970 !
 971    do i = its,ite
 972      if(.not.stable(i))then
 973        if((xland(i)-1.5).ge.0)then
 974          brcr(i) = brcr_sbro(i)
 975        else
 976          brcr(i) = brcr_sb
 977        endif
 978      endif
 979    enddo
 980 !
 981    do k = 2,klpbl
 982      do i = its,ite
 983        if(.not.stable(i))then
 984          brdn(i) = brup(i)
 985          spdk2   = max(ux(i,k)**2+vx(i,k)**2,1.)
 986          brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
 987          kpbl(i) = k
 988          stable(i) = brup(i).gt.brcr(i)
 989        endif
 990      enddo
 991    enddo
 992 !
 993    do i = its,ite
 994      if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
 995        k = kpbl(i)
 996        if(brdn(i).ge.brcr(i))then
 997          brint = 0.
 998        elseif(brup(i).le.brcr(i))then
 999          brint = 1.
1000        else
1001          brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
1002        endif
1003        hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
1004        if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
1005        if(kpbl(i).le.1) pblflg(i) = .false.
1006      endif
1007    enddo
1008 !
1009 !     estimate the entrainment parameters
1010 !
1011    do i = its,ite
1012      cloudflg(i)=.false.
1013      if(pblflg(i)) then
1014        k = kpbl(i) - 1
1015        wm3       = wstar3(i) + 5. * ust3(i)
1016        wm2(i)    = wm3**h2
1017        bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
1018        dthvx(i)  = max(thvx(i,k+1)-thvx(i,k),tmin)
1019        we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
1020        if((qx(i,ktrace2+k)+qx(i,ktrace3+k)).gt.0.01e-3.and.ysu_topdown_pblmix.eq.1)then
1021            if ( kpbl(i) .ge. 2) then
1022                 cloudflg(i)=.true.
1023                 templ=thlix(i,k)*(p2di(i,k+1)/100000)**rovcp
1024                 !rvls is ws at full level
1025                 rvls=100.*6.112*EXP(17.67*(templ-273.16)/(templ-29.65))*(ep2/p2di(i,k+1))
1026                 temps=templ + ((qx(i,k)+qx(i,ktrace2+k))-rvls)/(cp/xlv  + &
1027                 ep2*xlv*rvls/(rd*templ**2))
1028                 rvls=100.*6.112*EXP(17.67*(temps-273.15)/(temps-29.65))*(ep2/p2di(i,k+1))
1029                 rcldb=max((qx(i,k)+qx(i,ktrace2+k))-rvls,0.)
1030                 !entrainment efficiency
1031                 dthvx(i)  = (thlix(i,k+2)+thx(i,k+2)*ep1*(qx(i,k+2)+qx(i,ktrace2+k+2))) &
1032                           - (thlix(i,k) + thx(i,k)  *ep1*(qx(i,k)  +qx(i,ktrace2+k)))
1033                 dthvx(i)  = max(dthvx(i),0.1)
1034                 tmp1      = xlv/cp * rcldb/(pi2d(i,k)*dthvx(i))
1035                 ent_eff   = 0.2 * 8. * tmp1 +0.2
1036
1037                 radsum=0.
1038                 do kk = 1,kpbl(i)-1
1039                    radflux=rthraten(i,kk)*pi2d(i,kk) !converts theta/s to temp/s
1040                    radflux=radflux*cp/g*(p2diORG(i,kk)-p2diORG(i,kk+1)) ! converts temp/s to W/m^2
1041                    if (radflux < 0.0 ) radsum=abs(radflux)+radsum
1042                 enddo
1043                 radsum=max(radsum,0.0)
1044
1045                 !recompute entrainment from sfc thermals
1046                 bfx0 = max(max(sflux(i),0.0)-radsum/rhox2(i,k)/cp,0.)
1047                 bfx0 = max(sflux(i),0.0)
1048                 wm3 = (govrth(i)*bfx0*hpbl(i))+5. * ust3(i)
1049                 wm2(i)    = wm3**h2
1050                 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
1051                 dthvx(i)  = max(thvx(i,k+1)-thvx(i,k),tmin)
1052                 we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
1053
1054                 !entrainment from PBL top thermals
1055                 bfx0 = max(radsum/rhox2(i,k)/cp-max(sflux(i),0.0),0.)
1056                 bfx0 = max(radsum/rhox2(i,k)/cp,0.)
1057                 wm3       = (g/thvx(i,k)*bfx0*hpbl(i)) ! this is wstar3(i)
1058                 wm2(i)    = wm2(i)+wm3**h2
1059                 bfxpbl(i) = - ent_eff * bfx0
1060                 dthvx(i)  = max(thvx(i,k+1)-thvx(i,k),0.1)
1061                 we(i) = we(i) + max(bfxpbl(i)/dthvx(i),-sqrt(wm3**h2))
1062
1063                 !wstar3_2
1064                 bfx0 = max(radsum/rhox2(i,k)/cp,0.)
1065                 wstar3_2(i) =  (g/thvx(i,k)*bfx0*hpbl(i))
1066                 !recompute hgamt
1067                 wscale(i) = (ust3(i)+phifac*karman*(wstar3(i)+wstar3_2(i))*0.5)**h1
1068                 wscale(i) = min(wscale(i),ust(i)*aphi16)
1069                 wscale(i) = max(wscale(i),ust(i)/aphi5)
1070                 gamfac   = bfac/rhox(i)/wscale(i)
1071                 hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
1072                 hgamq(i) = min(gamfac*qfx(i),gamcrq)
1073                 gamfac   = bfac/rhox2(i,k)/wscale(i)
1074                 hgamt2(i,k) = min(gamfac*radsum/cp,gamcrt)
1075                 hgamt(i) = max(hgamt(i),0.0) + max(hgamt2(i,k),0.0)
1076                 brint    = -15.9*ust(i)*ust(i)/wspd(i)*(wstar3(i)+wstar3_2(i))/(wscale(i)**4.)
1077                 hgamu(i) = brint*ux(i,1)
1078                 hgamv(i) = brint*vx(i,1)
1079            endif
1080        endif
1081        prpbl(i) = 1.0
1082        dthx  = max(thx(i,k+1)-thx(i,k),tmin)
1083        dqx   = min(qx(i,k+1)-qx(i,k),0.0)
1084        hfxpbl(i) = we(i)*dthx
1085        qfxpbl(i) = we(i)*dqx
1086 !
1087        dux = ux(i,k+1)-ux(i,k)
1088        dvx = vx(i,k+1)-vx(i,k)
1089        if(dux.gt.tmin) then
1090          ufxpbl(i) = max(prpbl(i)*we(i)*dux,-ust(i)*ust(i))
1091        elseif(dux.lt.-tmin) then
1092          ufxpbl(i) = min(prpbl(i)*we(i)*dux,ust(i)*ust(i))
1093        else
1094          ufxpbl(i) = 0.0
1095        endif
1096        if(dvx.gt.tmin) then
1097          vfxpbl(i) = max(prpbl(i)*we(i)*dvx,-ust(i)*ust(i))
1098        elseif(dvx.lt.-tmin) then
1099          vfxpbl(i) = min(prpbl(i)*we(i)*dvx,ust(i)*ust(i))
1100        else
1101          vfxpbl(i) = 0.0
1102        endif
1103        delb  = govrth(i)*d3*hpbl(i)
1104        delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.)
1105      endif
1106    enddo
1107 !
1108    do k = kts,klpbl
1109      do i = its,ite
1110        if(pblflg(i).and.k.ge.kpbl(i))then
1111          entfac(i,k) = ((zq(i,k+1)-hpbl(i))/delta(i))**2.
1112        else
1113          entfac(i,k) = 1.e30
1114        endif
1115      enddo
1116    enddo
1117 !
1118 !     compute diffusion coefficients below pbl
1119 !
1120    do k = kts,klpbl
1121      do i = its,ite
1122        if(k.lt.kpbl(i)) then
1123          zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.)
1124          zfacent(i,k) = (1.-zfac(i,k))**3.
1125          wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
1126          wscalek2(i,k) = (phifac*karman*wstar3_2(i)*(zfac(i,k)))**h1
1127          if(sfcflg(i)) then
1128            prfac = conpr
1129            prfac2 = 15.9*(wstar3(i)+wstar3_2(i))/ust3(i)/(1.+4.*karman*(wstar3(i)+wstar3_2(i))/ust3(i))
1130            prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2.
1131          else
1132            prfac = 0.
1133            prfac2 = 0.
1134            prnumfac = 0.
1135            phim8z = 1.+aphi5*zol1(i)*zq(i,k+1)/zl1(i)
1136            wscalek(i,k) = ust(i)/phim8z
1137            wscalek(i,k) = max(wscalek(i,k),0.001)
1138          endif
1139          prnum0 = (phih(i)/phim(i)+prfac)
1140          prnum0 = max(min(prnum0,prmax),prmin)
1141            xkzm(i,k) = wscalek(i,k) *karman*    zq(i,k+1)      *    zfac(i,k)**pfac+ &
1142                        wscalek2(i,k)*karman*(hpbl(i)-zq(i,k+1))*(1-zfac(i,k))**pfac
1143          !Do not include xkzm at kpbl-1 since it changes entrainment
1144          if (k.eq.kpbl(i)-1.and.cloudflg(i).and.we(i).lt.0.0) then
1145            xkzm(i,k) = 0.0
1146          endif
1147          prnum =  1. + (prnum0-1.)*exp(prnumfac)
1148          xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac)
1149          prnum0 = prnum0/(1.+prfac2*karman*sfcfrac)
1150          prnum =  1. + (prnum0-1.)*exp(prnumfac)
1151          xkzh(i,k) = xkzm(i,k)/prnum
1152          xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
1153          xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
1154          xkzq(i,k) = xkzq(i,k)+xkzoh(i,k)
1155          xkzm(i,k) = min(xkzm(i,k),xkzmax)
1156          xkzh(i,k) = min(xkzh(i,k),xkzmax)
1157          xkzq(i,k) = min(xkzq(i,k),xkzmax)
1158        endif
1159      enddo
1160    enddo
1161 !
1162 !     compute diffusion coefficients over pbl (free atmosphere)
1163 !
1164    do k = kts,kte-1
1165      do i = its,ite
1166        if(k.ge.kpbl(i)) then
1167          ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k))                         &
1168               +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k)))                        &
1169               /(dza(i,k+1)*dza(i,k+1))+1.e-9
1170          govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
1171          ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1))
1172          if(imvdif.eq.1.and.ndiff.ge.3)then
1173            if((qx(i,ktrace2+k)+qx(i,ktrace3+k)).gt.0.01e-3.and.(qx(i           &
1174              ,ktrace2+k+1)+qx(i,ktrace3+k+1)).gt.0.01e-3)then
1175 !      in cloud
1176              qmean = 0.5*(qx(i,k)+qx(i,k+1))
1177              tmean = 0.5*(tx(i,k)+tx(i,k+1))
1178              alph  = xlv*qmean/rd/tmean
1179              chi   = xlv*xlv*qmean/cp/rv/tmean/tmean
1180              ri    = (1.+alph)*(ri-g*g/ss/tmean/cp*((chi-alph)/(1.+chi)))
1181            endif
1182          endif
1183          zk = karman*zq(i,k+1)
1184          rlamdz = min(max(0.1*dza(i,k+1),rlam),300.)
1185          rlamdz = min(dza(i,k+1),rlamdz)
1186          rl2 = (zk*rlamdz/(rlamdz+zk))**2
1187          dk = rl2*sqrt(ss)
1188          if(ri.lt.0.)then
1189 ! unstable regime
1190            ri = max(ri, rimin)
1191            sri = sqrt(-ri)
1192            xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
1193            xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
1194          else
1195 ! stable regime
1196            xkzh(i,k) = dk/(1+5.*ri)**2
1197            prnum = 1.0+2.1*ri
1198            prnum = min(prnum,prmax)
1199            xkzm(i,k) = xkzh(i,k)*prnum
1200          endif
1201 !
1202          xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
1203          xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
1204          xkzm(i,k) = min(xkzm(i,k),xkzmax)
1205          xkzh(i,k) = min(xkzh(i,k),xkzmax)
1206          xkzml(i,k) = xkzm(i,k)
1207          xkzhl(i,k) = xkzh(i,k)
1208        endif
1209      enddo
1210    enddo
1211 !
1212 !     compute tridiagonal matrix elements for heat
1213 !
1214    do k = kts,kte
1215      do i = its,ite
1216        au(i,k) = 0.
1217        al(i,k) = 0.
1218        ad(i,k) = 0.
1219        f1(i,k) = 0.
1220      enddo
1221    enddo
1222 !
1223    do i = its,ite
1224      ad(i,1) = 1.
1225      f1(i,1) = thx(i,1)-300.+(1.0-bepswitch)*hfx(i)/cont/del(i,1)*dt2
1226    enddo
1227 !
1228    do k = kts,kte-1
1229      do i = its,ite
1230        dtodsd = sfk2d(i,k)*dt2/del(i,k)
1231        dtodsu = sfk2d(i,k)*dt2/del(i,k+1)
1232        dsig   = p2d(i,k)-p2d(i,k+1)
1233        rdz    = 1./dza(i,k+1)
1234        tem1   = dsig*xkzh(i,k)*rdz
1235        if(pblflg(i).and.k.lt.kpbl(i)) then
1236          dsdzt = tem1*(-hgamt(i)/hpbl(i)-hfxpbl(i)*zfacent(i,k)/xkzh(i,k))
1237          f1(i,k)   = f1(i,k)+dtodsd*dsdzt
1238          f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt
1239        elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1240          xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
1241          xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
1242          xkzh(i,k) = max(xkzh(i,k),xkzoh(i,k))
1243          xkzh(i,k) = min(xkzh(i,k),xkzmax)
1244          f1(i,k+1) = thx(i,k+1)-300.
1245        else
1246          f1(i,k+1) = thx(i,k+1)-300.
1247        endif
1248        tem1   = dsig*xkzh(i,k)*rdz
1249        dsdz2     = tem1*rdz
1250        au(i,k)   = -dtodsd*dsdz2/vlk2d(i,k)
1251        al(i,k)   = -dtodsu*dsdz2/vlk2d(i,k)
1252        ad(i,k)   = ad(i,k)-au(i,k)
1253        ad(i,k+1) = 1.-al(i,k)
1254        exch_hx(i,k+1) = xkzh(i,k)
1255      enddo
1256    enddo
1257 !
1258 ! add bep/bep+bem forcing for heat if flag_bep=.true.
1259 !
1260    do k = kts,kte
1261      do i = its,ite
1262       ad(i,k) = ad(i,k) - a_t2d(i,k)*dt2
1263       f1(i,k) = f1(i,k) + b_t2d(i,k)*dt2
1264      enddo
1265    enddo
1266 !
1267 ! copies here to avoid duplicate input args for tridin
1268 !
1269    do k = kts,kte
1270      do i = its,ite
1271        cu(i,k) = au(i,k)
1272        r1(i,k) = f1(i,k)
1273      enddo
1274    enddo
1275 !
1276    call tridin_ysu(al,ad,cu,r1,au,f1,its,ite,kts,kte,1)
1277 !
1278 !     recover tendencies of heat
1279 !
1280    do k = kte,kts,-1
1281      do i = its,ite
1282        ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k)
1283        ttnp(i,k) = ttnp(i,k)+ttend
1284        dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)/pi2d(i,k)
1285      enddo
1286    enddo
1287 !
1288 !     compute tridiagonal matrix elements for moisture, clouds, and gases
1289 !
1290    do k = kts,kte
1291      do i = its,ite
1292        au(i,k) = 0.
1293        al(i,k) = 0.
1294        ad(i,k) = 0.
1295      enddo
1296    enddo
1297 !
1298    do ic = 1,ndiff
1299      do i = its,ite
1300        do k = kts,kte
1301          f3(i,k,ic) = 0.
1302        enddo
1303      enddo
1304    enddo
1305 !
1306    do i = its,ite
1307      ad(i,1) = 1.
1308      f3(i,1,1) = qx(i,1)+(1.0-bepswitch)*qfx(i)*g/del(i,1)*dt2
1309    enddo
1310 !
1311    if(ndiff.ge.2) then
1312      do ic = 2,ndiff
1313        is = (ic-1) * kte
1314        do i = its,ite
1315          f3(i,1,ic) = qx(i,1+is)
1316        enddo
1317      enddo
1318    endif
1319 !
1320    do k = kts,kte-1
1321      do i = its,ite
1322        if(k.ge.kpbl(i)) then
1323          xkzq(i,k) = xkzh(i,k)
1324        endif
1325      enddo
1326    enddo
1327 !
1328    do k = kts,kte-1
1329      do i = its,ite
1330        dtodsd = sfk2d(i,k)*dt2/del(i,k)
1331        dtodsu = sfk2d(i,k)*dt2/del(i,k+1)
1332        dsig   = p2d(i,k)-p2d(i,k+1)
1333        rdz    = 1./dza(i,k+1)
1334        tem1   = dsig*xkzq(i,k)*rdz
1335        if(pblflg(i).and.k.lt.kpbl(i)) then
1336          dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzq(i,k))
1337          f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq
1338          f3(i,k+1,1) = qx(i,k+1)-dtodsu*dsdzq
1339        elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1340          xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
1341          xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k))
1342          xkzq(i,k) = max(xkzq(i,k),xkzoh(i,k))
1343          xkzq(i,k) = min(xkzq(i,k),xkzmax)
1344          f3(i,k+1,1) = qx(i,k+1)
1345        else
1346          f3(i,k+1,1) = qx(i,k+1)
1347        endif
1348        tem1   = dsig*xkzq(i,k)*rdz
1349        dsdz2     = tem1*rdz
1350        au(i,k)   = -dtodsd*dsdz2/vlk2d(i,k)
1351        al(i,k)   = -dtodsu*dsdz2/vlk2d(i,k)
1352        ad(i,k)   = ad(i,k)-au(i,k)
1353        ad(i,k+1) = 1.-al(i,k)
1354 !      exch_hx(i,k+1) = xkzh(i,k)
1355      enddo
1356    enddo
1357 !
1358    if(ndiff.ge.2) then
1359      do ic = 2,ndiff
1360        is = (ic-1) * kte
1361        do k = kts,kte-1
1362          do i = its,ite
1363            f3(i,k+1,ic) = qx(i,k+1+is)
1364          enddo
1365        enddo
1366      enddo
1367    endif
1368 !
1369 ! add bep/bep+bem forcing for water vapor if flag_bep=.true.
1370 !
1371    do k = kts,kte
1372      do i = its,ite
1373       ad(i,k) = ad(i,k) - a_q2d(i,k)*dt2
1374       f3(i,k,1) = f3(i,k,1) + b_q2d(i,k)*dt2
1375      enddo
1376    enddo
1377 !
1378 ! copies here to avoid duplicate input args for tridin
1379 !
1380    do k = kts,kte
1381      do i = its,ite
1382        cu(i,k) = au(i,k)
1383      enddo
1384    enddo
1385 !
1386    do ic = 1,ndiff
1387      do k = kts,kte
1388        do i = its,ite
1389          r3(i,k,ic) = f3(i,k,ic)
1390        enddo
1391      enddo
1392    enddo
1393 !
1394 !     solve tridiagonal problem for moisture, clouds, and gases
1395 !
1396    call tridin_ysu(al,ad,cu,r3,au,f3,its,ite,kts,kte,ndiff)
1397 !
1398 !     recover tendencies of heat and moisture
1399 !
1400    do k = kte,kts,-1
1401      do i = its,ite
1402        qtend = (f3(i,k,1)-qx(i,k))*rdt
1403        qtnp(i,k) = qtnp(i,k)+qtend
1404        dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k)
1405      enddo
1406    enddo
1407 !
1408    if(ndiff.ge.2) then
1409      do ic = 2,ndiff
1410        is = (ic-1) * kte
1411        do k = kte,kts,-1
1412          do i = its,ite
1413            qtend = (f3(i,k,ic)-qx(i,k+is))*rdt
1414            qtnp(i,k+is) = qtnp(i,k+is)+qtend
1415          enddo
1416        enddo
1417      enddo
1418    endif
1419 !
1420 !     compute tridiagonal matrix elements for momentum
1421 !
1422    do i = its,ite
1423      do k = kts,kte
1424        au(i,k) = 0.
1425        al(i,k) = 0.
1426        ad(i,k) = 0.
1427        f1(i,k) = 0.
1428        f2(i,k) = 0.
1429      enddo
1430    enddo
1431 !
1432 ! paj: ctopo=1 if topo_wind=0 (default)
1433 !raquel---paj tke code (could be replaced with shin-hong tke in future
1434    do i = its,ite
1435       do k= kts, kte-1
1436         shear_ysu(i,k)=xkzm(i,k)*((-hgamu(i)/hpbl(i)+(ux(i,k+1)-ux(i,k))/dza(i,k+1))*(ux(i,k+1)-ux(i,k))/dza(i,k+1) &
1437         + (-hgamv(i)/hpbl(i)+(vx(i,k+1)-vx(i,k))/dza(i,k+1))*(vx(i,k+1)-vx(i,k))/dza(i,k+1))
1438          buoy_ysu(i,k)=xkzh(i,k)*g*(1.0/thx(i,k))*(-hgamt(i)/hpbl(i)+(thx(i,k+1)-thx(i,k))/dza(i,k+1))
1439
1440        zk = karman*zq(i,k+1)
1441  !over pbl
1442        if (k.ge.kpbl(i)) then
1443         rlamdz = min(max(0.1*dza(i,k+1),rlam),300.)
1444         rlamdz = min(dza(i,k+1),rlamdz)
1445        else
1446  !in pbl
1447         rlamdz = 150.0
1448        endif
1449        el_ysu(i,k) = zk*rlamdz/(rlamdz+zk)
1450        tke_ysu(i,k)=16.6*el_ysu(i,k)*(shear_ysu(i,k)-buoy_ysu(i,k)+b_e2d(i,k))
1451  !q2 when q3 positive
1452        if(tke_ysu(i,k).le.0) then
1453         tke_ysu(i,k)=0.0
1454        else
1455         tke_ysu(i,k)=(tke_ysu(i,k))**0.66
1456        endif
1457       enddo
1458  !Hybrid pblh of MYNN
1459  !tke is q2
1460       CALL GET_PBLH(KTS,KTE,pblh_ysu(i),thvx(i,kts:kte),&
1461       &    tke_ysu(i,kts:kte),zq(i,kts:kte+1),dzq(i,kts:kte),xland(i))
1462 !--- end of paj tke
1463 ! compute vconv
1464 !      Use Beljaars over land
1465         if (xland(i).lt.1.5) then
1466         fluxc = max(sflux(i),0.0)
1467         vconvc=1.
1468         VCONV = vconvc*(g/thvx(i,1)*pblh_ysu(i)*fluxc)**.33
1469         else
1470 ! for water there is no topo effect so vconv not needed
1471         VCONV = 0.
1472         endif
1473         vconvfx(i) = vconv
1474 !raquel
1475 !ctopo stability correction
1476       fric(i,1)=ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2         &
1477         *(wspd1(i)/wspd(i))**2
1478       if(present(ctopo)) then
1479         vconvnew=0.9*vconvfx(i)+1.5*(max((pblh_ysu(i)-500)/1000.0,0.0))
1480         vconvlim = min(vconvnew,1.0)
1481        ad(i,1) = 1.+(1.0-bepswitch*frc_urb1d(i))* &
1482                      (fric(i,1)*vconvlim+ctopo(i)*fric(i,1)*(1-vconvlim)) - &
1483                      fric(i,1)*bepswitch*(1-frc_urb1d(i))
1484       else
1485        ad(i,1) = 1.+(1.0-bepswitch)*fric(i,1)
1486       endif
1487      f1(i,1) = ux(i,1)+uox(i)*ust(i)**2*rhox(i)*g/del(i,1)*dt2/wspd1(i)*(wspd1(i)/wspd(i))**2
1488      f2(i,1) = vx(i,1)+vox(i)*ust(i)**2*rhox(i)*g/del(i,1)*dt2/wspd1(i)*(wspd1(i)/wspd(i))**2
1489    enddo
1490 !
1491    do k = kts,kte-1
1492      do i = its,ite
1493        dtodsd = sfk2d(i,k)*dt2/del(i,k)
1494        dtodsu = sfk2d(i,k)*dt2/del(i,k+1)
1495        dsig   = p2d(i,k)-p2d(i,k+1)
1496        rdz    = 1./dza(i,k+1)
1497        tem1   = dsig*xkzm(i,k)*rdz
1498        if(pblflg(i).and.k.lt.kpbl(i))then
1499          dsdzu     = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k))
1500          dsdzv     = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k))
1501          f1(i,k)   = f1(i,k)+dtodsd*dsdzu
1502          f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu
1503          f2(i,k)   = f2(i,k)+dtodsd*dsdzv
1504          f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv
1505        elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1506          xkzm(i,k) = prpbl(i)*xkzh(i,k)
1507          xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
1508          xkzm(i,k) = max(xkzm(i,k),xkzom(i,k))
1509          xkzm(i,k) = min(xkzm(i,k),xkzmax)
1510          f1(i,k+1) = ux(i,k+1)
1511          f2(i,k+1) = vx(i,k+1)
1512        else
1513          f1(i,k+1) = ux(i,k+1)
1514          f2(i,k+1) = vx(i,k+1)
1515        endif
1516        tem1   = dsig*xkzm(i,k)*rdz
1517        dsdz2     = tem1*rdz
1518        au(i,k)   = -dtodsd*dsdz2/vlk2d(i,k)
1519        al(i,k)   = -dtodsu*dsdz2/vlk2d(i,k)
1520        ad(i,k)   = ad(i,k)-au(i,k)
1521        ad(i,k+1) = 1.-al(i,k)
1522        exch_mx(i,k+1) = xkzm(i,k)
1523      enddo
1524    enddo
1525 !
1526 ! add bep/bep+bem forcing for momentum if flag_bep=.true.
1527 !
1528    do k = kts,kte
1529      do i = its,ite
1530        ad1(i,k) = ad(i,k)
1531      end do
1532    end do
1533    do k = kts,kte
1534      do i = its,ite
1535        ad(i,k) = ad(i,k) - a_u2d(i,k)*dt2
1536        ad1(i,k) = ad1(i,k) - a_v2d(i,k)*dt2
1537        f1(i,k) = f1(i,k) + b_u2d(i,k)*dt2
1538        f2(i,k) = f2(i,k) + b_v2d(i,k)*dt2
1539      enddo
1540    enddo
1541 !
1542 ! copies here to avoid duplicate input args for tridin
1543 !
1544    do k = kts,kte
1545      do i = its,ite
1546        cu(i,k) = au(i,k)
1547        r1(i,k) = f1(i,k)
1548        r2(i,k) = f2(i,k)
1549      enddo
1550    enddo
1551 !
1552 !     solve tridiagonal problem for momentum
1553 !
1554     call tridi2n(al,ad,ad1,cu,r1,r2,au,f1,f2,its,ite,kts,kte,1)
1555 !
1556 !     recover tendencies of momentum
1557 !
1558    do k = kte,kts,-1
1559      do i = its,ite
1560        utend = (f1(i,k)-ux(i,k))*rdt
1561        vtend = (f2(i,k)-vx(i,k))*rdt
1562        utnp(i,k) = utnp(i,k)+utend
1563        vtnp(i,k) = vtnp(i,k)+vtend
1564        dusfc(i) = dusfc(i) + utend*conwrc*del(i,k)
1565        dvsfc(i) = dvsfc(i) + vtend*conwrc*del(i,k)
1566      enddo
1567    enddo
1568 !
1569 ! paj: ctopo2=1 if topo_wind=0 (default)
1570 !
1571    do i = its,ite
1572      if(present(ctopo).and.present(ctopo2)) then ! mchen for NMM
1573        u10(i) = ctopo2(i)*u10(i)+(1-ctopo2(i))*ux(i,1)
1574        v10(i) = ctopo2(i)*v10(i)+(1-ctopo2(i))*vx(i,1)
1575      endif !mchen
1576    enddo
1577 !
1578 !---- end of vertical diffusion
1579 !
1580    do i = its,ite
1581      kpbl1d(i) = kpbl(i)
1582    enddo
1583 !
1584    end subroutine ysu2d
1585 !-------------------------------------------------------------------------------
1586 !
1587 !-------------------------------------------------------------------------------
1588    subroutine tridi2n(cl,cm,cm1,cu,r1,r2,au,f1,f2,its,ite,kts,kte,nt)
1589 !-------------------------------------------------------------------------------
1590    implicit none
1591 !-------------------------------------------------------------------------------
1592 !
1593    integer, intent(in )      ::     its,ite, kts,kte, nt
1594 !
1595    real, dimension( its:ite, kts+1:kte+1 )                                   , &
1596          intent(in   )  ::                                                 cl
1597 !
1598    real, dimension( its:ite, kts:kte )                                       , &
1599          intent(in   )  ::                                                 cm, &
1600                                                                           cm1, &
1601                                                                            r1
1602    real, dimension( its:ite, kts:kte,nt )                                    , &
1603          intent(in   )  ::                                                 r2
1604 !
1605    real, dimension( its:ite, kts:kte )                                       , &
1606          intent(inout)  ::                                                 au, &
1607                                                                            cu, &
1608                                                                            f1
1609    real, dimension( its:ite, kts:kte,nt )                                    , &
1610          intent(inout)  ::                                                 f2
1611 !
1612    real    :: fk
1613    integer :: i,k,l,n,it
1614 !
1615 !-------------------------------------------------------------------------------
1616 !
1617    l = ite
1618    n = kte
1619 !
1620    do i = its,l
1621      fk = 1./cm(i,1)
1622      au(i,1) = fk*cu(i,1)
1623      f1(i,1) = fk*r1(i,1)
1624    enddo
1625 !
1626    do it = 1,nt
1627      do i = its,l
1628        fk = 1./cm1(i,1)
1629        f2(i,1,it) = fk*r2(i,1,it)
1630      enddo
1631    enddo
1632
1633    do k = kts+1,n-1
1634      do i = its,l
1635        fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1636        au(i,k) = fk*cu(i,k)
1637        f1(i,k) = fk*(r1(i,k)-cl(i,k)*f1(i,k-1))
1638      enddo
1639    enddo
1640 !
1641    do it = 1,nt
1642      do k = kts+1,n-1
1643        do i = its,l
1644          fk = 1./(cm1(i,k)-cl(i,k)*au(i,k-1))
1645          f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
1646        enddo
1647      enddo
1648    enddo
1649 !
1650    do i = its,l
1651      fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1652      f1(i,n) = fk*(r1(i,n)-cl(i,n)*f1(i,n-1))
1653    enddo
1654 !
1655    do it = 1,nt
1656      do i = its,l
1657        fk = 1./(cm1(i,n)-cl(i,n)*au(i,n-1))
1658        f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
1659      enddo
1660    enddo
1661 !
1662    do k = n-1,kts,-1
1663      do i = its,l
1664        f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1)
1665      enddo
1666    enddo
1667 !
1668    do it = 1,nt
1669      do k = n-1,kts,-1
1670        do i = its,l
1671          f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1672        enddo
1673      enddo
1674    enddo
1675 !
1676    end subroutine tridi2n
1677 !-------------------------------------------------------------------------------
1678 !
1679 !-------------------------------------------------------------------------------
1680    subroutine tridin_ysu(cl,cm,cu,r2,au,f2,its,ite,kts,kte,nt)
1681 !-------------------------------------------------------------------------------
1682    implicit none
1683 !-------------------------------------------------------------------------------
1684 !
1685    integer, intent(in )      ::     its,ite, kts,kte, nt
1686 !
1687    real, dimension( its:ite, kts+1:kte+1 )                                   , &
1688          intent(in   )  ::                                                 cl
1689 !
1690    real, dimension( its:ite, kts:kte )                                       , &
1691          intent(in   )  ::                                                 cm
1692    real, dimension( its:ite, kts:kte,nt )                                    , &
1693          intent(in   )  ::                                                 r2
1694 !
1695    real, dimension( its:ite, kts:kte )                                       , &
1696          intent(inout)  ::                                                 au, &
1697                                                                            cu
1698    real, dimension( its:ite, kts:kte,nt )                                    , &
1699          intent(inout)  ::                                                 f2
1700 !
1701    real    :: fk
1702    integer :: i,k,l,n,it
1703 !
1704 !-------------------------------------------------------------------------------
1705 !
1706    l = ite
1707    n = kte
1708 !
1709    do it = 1,nt
1710      do i = its,l
1711        fk = 1./cm(i,1)
1712        au(i,1) = fk*cu(i,1)
1713        f2(i,1,it) = fk*r2(i,1,it)
1714      enddo
1715    enddo
1716 !
1717    do it = 1,nt
1718      do k = kts+1,n-1
1719        do i = its,l
1720          fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1721          au(i,k) = fk*cu(i,k)
1722          f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
1723        enddo
1724      enddo
1725    enddo
1726 !
1727    do it = 1,nt
1728      do i = its,l
1729        fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1730        f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
1731      enddo
1732    enddo
1733 !
1734    do it = 1,nt
1735      do k = n-1,kts,-1
1736        do i = its,l
1737          f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1738        enddo
1739      enddo
1740    enddo
1741 !
1742    end subroutine tridin_ysu
1743 !-------------------------------------------------------------------------------
1744 !
1745 !-------------------------------------------------------------------------------
1746    subroutine ysuinit(rublten,rvblten,rthblten,rqvblten,                       &
1747                       rqcblten,rqiblten,p_qi,p_first_scalar,                   &
1748                       restart, allowed_to_read,                                &
1749                       ids, ide, jds, jde, kds, kde,                            &
1750                       ims, ime, jms, jme, kms, kme,                            &
1751                       its, ite, jts, jte, kts, kte                 )
1752 !-------------------------------------------------------------------------------
1753    implicit none
1754 !-------------------------------------------------------------------------------
1755 !
1756    logical , intent(in)          :: restart, allowed_to_read
1757    integer , intent(in)          ::  ids, ide, jds, jde, kds, kde,             &
1758                                      ims, ime, jms, jme, kms, kme,             &
1759                                      its, ite, jts, jte, kts, kte
1760    integer , intent(in)          ::  p_qi,p_first_scalar
1761    real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) ::             &
1762                                                                       rublten, &
1763                                                                       rvblten, &
1764                                                                      rthblten, &
1765                                                                      rqvblten, &
1766                                                                      rqcblten, &
1767                                                                      rqiblten
1768    integer :: i, j, k, itf, jtf, ktf
1769 !
1770    jtf = min0(jte,jde-1)
1771    ktf = min0(kte,kde-1)
1772    itf = min0(ite,ide-1)
1773 !
1774    if(.not.restart)then
1775      do j = jts,jtf
1776        do k = kts,ktf
1777          do i = its,itf
1778             rublten(i,k,j) = 0.
1779             rvblten(i,k,j) = 0.
1780             rthblten(i,k,j) = 0.
1781             rqvblten(i,k,j) = 0.
1782             rqcblten(i,k,j) = 0.
1783          enddo
1784        enddo
1785      enddo
1786    endif
1787 !
1788    if (p_qi .ge. p_first_scalar .and. .not.restart) then
1789      do j = jts,jtf
1790        do k = kts,ktf
1791          do i = its,itf
1792            rqiblten(i,k,j) = 0.
1793          enddo
1794        enddo
1795      enddo
1796    endif
1797 !
1798    end subroutine ysuinit
1799 !-------------------------------------------------------------------------------
1800 ! ==================================================================
1801
1802       SUBROUTINE GET_PBLH(KTS,KTE,zi,thetav1D,qke1D,zw1D,dz1D,landsea)
1803 ! Copied from MYNN PBL
1804
1805       !---------------------------------------------------------------
1806       !             NOTES ON THE PBLH FORMULATION
1807       !
1808       !The 1.5-theta-increase method defines PBL heights as the level at
1809       !which the potential temperature first exceeds the minimum potential
1810       !temperature within the boundary layer by 1.5 K. When applied to
1811       !observed temperatures, this method has been shown to produce PBL-
1812       !height estimates that are unbiased relative to profiler-based
1813       !estimates (Nielsen-Gammon et al. 2008). However, their study did not
1814       !include LLJs. Banta and Pichugina (2008) show that a TKE-based
1815       !threshold is a good estimate of the PBL height in LLJs. Therefore,
1816       !a hybrid definition is implemented that uses both methods, weighting
1817       !the TKE-method more during stable conditions (PBLH < 400 m).
1818       !A variable tke threshold (TKEeps) is used since no hard-wired
1819       !value could be found to work best in all conditions.
1820       !---------------------------------------------------------------
1821
1822       INTEGER,INTENT(IN) :: KTS,KTE
1823       REAL, INTENT(OUT) :: zi
1824       REAL, INTENT(IN) :: landsea
1825       REAL, DIMENSION(KTS:KTE), INTENT(IN) :: thetav1D, qke1D, dz1D
1826       REAL, DIMENSION(KTS:KTE+1), INTENT(IN) :: zw1D
1827       !LOCAL VARS
1828       REAL ::  PBLH_TKE,qtke,qtkem1,wt,maxqke,TKEeps,minthv
1829       REAL :: delt_thv   !delta theta-v; dependent on land/sea point
1830       REAL, PARAMETER :: sbl_lim  = 200. !Theta-v PBL lower limit of trust (m).
1831       REAL, PARAMETER :: sbl_damp = 400. !Damping range for averaging with TKE-based PBLH (m).
1832       INTEGER :: I,J,K,kthv,ktke
1833
1834       !FIND MAX TKE AND MIN THETAV IN THE LOWEST 500 M
1835       k = kts+1
1836       kthv = 1
1837       ktke = 1
1838       maxqke = 0.
1839       minthv = 9.E9
1840
1841       DO WHILE (zw1D(k) .LE. 500.)
1842         qtke  =MAX(Qke1D(k),0.)   ! maximum QKE
1843          IF (maxqke < qtke) then
1844             maxqke = qtke
1845             ktke = k
1846          ENDIF
1847          IF (minthv > thetav1D(k)) then
1848              minthv = thetav1D(k)
1849              kthv = k
1850          ENDIF
1851          k = k+1
1852       ENDDO
1853       !TKEeps = maxtke/20. = maxqke/40.
1854       TKEeps = maxqke/40.
1855       TKEeps = MAX(TKEeps,0.025)
1856       TKEeps = MIN(TKEeps,0.25)
1857
1858       !FIND THETAV-BASED PBLH (BEST FOR DAYTIME).
1859       zi=0.
1860       k = kthv+1
1861       IF((landsea-1.5).GE.0)THEN
1862       ! WATER
1863           delt_thv = 0.75
1864       ELSE
1865       ! LAND
1866           delt_thv = 1.5
1867       ENDIF
1868
1869       zi=0.
1870       k = kthv+1
1871       DO WHILE (zi .EQ. 0.)
1872          IF (thetav1D(k) .GE. (minthv + delt_thv))THEN
1873             zi = zw1D(k) - dz1D(k-1)* &
1874               & MIN((thetav1D(k)-(minthv + delt_thv))/MAX(thetav1D(k)-thetav1D(k-1),1E-6),1.0)
1875         ENDIF
1876         k = k+1
1877          IF (k .EQ. kte-1) zi = zw1D(kts+1) !EXIT SAFEGUARD
1878       ENDDO
1879
1880       !print*,"IN GET_PBLH:",thsfc,zi
1881       !FOR STABLE BOUNDARY LAYERS, USE TKE METHOD TO COMPLEMENT THE
1882       !THETAV-BASED DEFINITION (WHEN THE THETA-V BASED PBLH IS BELOW ~0.5 KM).
1883       !THE TANH WEIGHTING FUNCTION WILL MAKE THE TKE-BASED DEFINITION NEGLIGIBLE
1884       !WHEN THE THETA-V-BASED DEFINITION IS ABOVE ~1 KM.
1885       !FIND TKE-BASED PBLH (BEST FOR NOCTURNAL/STABLE CONDITIONS).
1886
1887       PBLH_TKE=0.
1888       k = ktke+1
1889      DO WHILE (PBLH_TKE .EQ. 0.)
1890         !QKE CAN BE NEGATIVE (IF CKmod == 0)... MAKE TKE NON-NEGATIVE.
1891          qtke  =MAX(Qke1D(k)/2.,0.)      ! maximum TKE
1892          qtkem1=MAX(Qke1D(k-1)/2.,0.)
1893          IF (qtke .LE. TKEeps) THEN
1894                PBLH_TKE = zw1D(k) - dz1D(k-1)* &
1895                & MIN((TKEeps-qtke)/MAX(qtkem1-qtke, 1E-6), 1.0)
1896              !IN CASE OF NEAR ZERO TKE, SET PBLH = LOWEST LEVEL.
1897              PBLH_TKE = MAX(PBLH_TKE,zw1D(kts+1))
1898              !print *,"PBLH_TKE:",i,j,PBLH_TKE, Qke1D(k)/2., zw1D(kts+1)
1899          ENDIF
1900          k = k+1
1901          IF (k .EQ. kte-1) PBLH_TKE = zw1D(kts+1) !EXIT SAFEGUARD
1902       ENDDO
1903
1904     !BLEND THE TWO PBLH TYPES HERE:
1905
1906       wt=.5*TANH((zi - sbl_lim)/sbl_damp) + .5
1907       zi=PBLH_TKE*(1.-wt) + zi*wt
1908
1909    END SUBROUTINE GET_PBLH
1910 ! ==================================================================
1911
1912 end module module_bl_ysu
1913 !-------------------------------------------------------------------------------