phys/module_bl_keps.F

   1 MODULE module_bl_keps
   2
   3       real      vk,temin,dmin,tpemin,c1,c2,c3,c_mu,Pr    ! constant from Zonato et al.,2022
   4       parameter(temin=1E-4,dmin=1E-7,tpemin=1E-7,c_theta=0.09,  &
   5                 vk=0.4,c_mu=0.09,c1=1.44,c2=1.92,c3=1.44,c4=0.27,c5=0.08,sigma_eps=1.3) ! impose minimum values for tke similar to those of MYJ
   6       parameter(temax=50.,dmax=50.,tpemax=5.)
   7       integer,private,parameter:: bl_keps_adv = 0
   8 ! -----------------------------------------------------------------------
   9
  10
  11    CONTAINS
  12
  13       subroutine keps(mol,tsk,xtime,frc_urb2d,flag_bep,dz8w,dt,u_phy,v_phy,th_phy &
  14                       ,PI_PHY,RTHRATEN,P8W                          &
  15                       ,rho,qv_curr,qc_curr,hfx                  &
  16                       ,qfx,ustar,cp,g                           &
  17                       ,rublten,rvblten,rthblten                 &
  18                       ,rqvblten,rqcblten                        &
  19                       ,tke_pbl,diss_pbl,tpe_pbl                 &
  20                       ,tke_adv,diss_adv,tpe_adv                 &
  21                       ,pr_pbl                   &
  22                       ,wu,wv,wt,wq,exch_h,exch_m,pblh           &
  23                       ,a_u_bep,a_v_bep,a_t_bep,a_q_bep          &
  24                       ,b_u_bep,b_v_bep                          &
  25                       ,b_t_bep,b_q_bep                          &
  26                       ,b_e_bep                                  &
  27                       ,sf_bep,vl_bep                            &
  28                       ,br,znt,psim,psih                         &
  29                       ,ids,ide, jds,jde, kds,kde                &
  30                       ,ims,ime, jms,jme, kms,kme                &
  31                       ,its,ite, jts,jte, kts,kte)
  32
  33       implicit none
  34
  35
  36
  37 !-----------------------------------------------------------------------
  38 !     Input
  39 !------------------------------------------------------------------------
  40    INTEGER::                        ids,ide, jds,jde, kds,kde,  &
  41                                     ims,ime, jms,jme, kms,kme,  &
  42                                     its,ite, jts,jte, kts,kte
  43
  44
  45    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   DZ8W      !vertical resolution
  46    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   qv_curr   !moisture
  47    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   qc_curr   !liquid water
  48    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   RHO       !air density
  49    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   TH_PHY    !potential temperature
  50    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   PI_PHY
  51    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   RTHRATEN
  52    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   P8W
  53   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   U_PHY     !x-component of wind
  54    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   V_PHY     !y-component of wind
  55    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   ustar              !friction velocity
  56    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   hfx                !sensible heat flux (W/m2) at surface
  57    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   qfx                !moisture flux at surface
  58    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   tsk,mol                !moisture flux at surface
  59    real,  INTENT(IN   )    :: g,cp                                              !gravity and Cp
  60    REAL, INTENT(IN )::   DT                                                      ! Time step
  61    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    :: br
  62    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    :: znt
  63    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   FRC_URB2D          !fraction cover urban
  64    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)    ::   PBLH          !PBL height
  65    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   psim,psih
  66
  67 ! variable added for urban
  68    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_u_bep        ! Implicit component for the momemtum in X-direction
  69    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_v_bep        ! Implicit component for the momemtum in Y-direction
  70    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_t_bep        ! Implicit component for the Pot. Temp.
  71    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_q_bep        ! Implicit component for Moisture
  72    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_u_bep        ! Explicit component for the momemtum in X-direction
  73    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_v_bep        ! Explicit component for the momemtum in Y-direction
  74    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_t_bep        ! Explicit component for the Pot. Temp.
  75    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_q_bep        ! Explicit component for Moisture
  76    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_e_bep        ! Explicit component for Moisture
  77
  78  ! urban surface and volumes
  79    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::sf_bep           ! surface of the urban grid cells
  80    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::vl_bep             ! volume of the urban grid cells
  81    LOGICAL, INTENT(IN) :: flag_bep                                             !flag for BEP
  82
  83 !
  84 !-----------------------------------------------------------------------
  85 !     Local, carried on from one timestep to the other
  86 !------------------------------------------------------------------------
  87 !      real, save, allocatable, dimension (:,:,:)::TKE  ! Turbulent kinetic energy
  88       real,  dimension (ims:ime, kms:kme, jms:jme)  ::th_0 ! reference state for potential temperature
  89
  90 !------------------------------------------------------------------------
  91 !     Output
  92 !------------------------------------------------------------------------
  93         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  exch_h ! exchange coefficient for heat
  94         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  exch_m ! exchange coefficient for momentum
  95         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT   )  ::  tke_pbl  ! Turbulence Kinetic Energy
  96         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT   )  ::  diss_pbl  ! Dissipation Rate
  97         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT   )  ::  tpe_pbl  ! Dissipation Rate
  98         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT   )  ::  pr_pbl  ! Dissipation Rate
  99         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT   )  ::  tke_adv  ! Dissipation Rate
 100         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT   )  ::  diss_adv  ! Dissipation Rate
 101         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT   )  ::  tpe_adv  ! Dissipation Rate
 102         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wu  ! Turbulent flux of momentum (x)
 103         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wv  ! Turbulent flux of momentum (y)
 104         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wt  ! Turbulent flux of temperature
 105         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wq  ! Turbulent flux of water vapor
 106
 107 ! only if idiff not equal 1:
 108         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RUBLTEN  !tendency for U_phy
 109         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RVBLTEN  !tendency for V_phy
 110         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RTHBLTEN !tendency for TH_phy
 111         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RQVBLTEN !tendency for QV_curr
 112         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RQCBLTEN !tendency for QV_curr
 113
 114 !--------------------------------------------------------------
 115 ! Local
 116 !--------------------------------------------------------------
 117 ! 1D array used for the input and output of the routine boulac1D
 118
 119       real z1D(kms:kme)               ! vertical coordinates (faces of the grid)
 120       real dz1D(kms:kme)              ! vertical resolution
 121       real u1D(kms:kme)                 ! wind speed in the x directions
 122       real v1D(kms:kme)                 ! wind speed in the y directions
 123       real th1D(kms:kme)                ! potential temperature
 124       real q1D(kms:kme)                 ! moisture
 125       real qc1D(kms:kme)                 ! liquid water
 126       real rho1D(kms:kme)               ! air density
 127       real rhoz1D(kms:kme)            ! air density at the faces
 128       real tke1D(kms:kme)               ! Turbulent Kinetic Energy
 129       real diss1D(kms:kme)              ! Dissipation Rate
 130       real tpe1D(kms:kme)
 131       real th01D(kms:kme)               ! reference potential temperature
 132       real exch1D(kms:kme)            ! exch
 133       real sf1D(kms:kme)              ! surface of the grid cells
 134       real vl1D(kms:kme)                ! volume of the  grid cells
 135       real a_u1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the X-direction
 136       real a_v1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the Y-direction
 137       real a_t1D(kms:kme)               ! Implicit component of the heat sources or sinks
 138       real a_q1D(kms:kme)               ! Implicit component of the moisture sources or sinks
 139       real a_qc1D(kms:kme)               ! Implicit component of the liquid water sources or sinks
 140       real b_u1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the X-direction
 141       real b_v1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the Y-direction
 142       real b_t1D(kms:kme)               ! Explicit component of the heat sources or sinks
 143       real b_q1D(kms:kme)               ! Explicit component of the moisture sources or sinks
 144       real b_qc1D(kms:kme)               ! Explicit component of the liquid water sources or sinks
 145       real b_e1D(kms:kme)               ! Explicit component of the moisture sources or sinks
 146       real a_e1D(kms:kme)               ! Explicit component of the moisture sources or sinks
 147       real b_d1D(kms:kme)               ! Explicit component of the moisture sources or sinks
 148       real a_d1D(kms:kme)               ! Explicit component of the moisture sources or sinks
 149       real b_tpe1D(kms:kme)               ! Explicit component of the moisture sources or sinks
 150       real a_tpe1D(kms:kme)               ! Explicit component of the moisture sources or sinks
 151       real sh1D(kms:kme)              ! shear
 152       real bu1D(kms:kme)              ! buoyancy
 153       real wu1D(kms:kme)              ! turbulent flux of momentum (x component)
 154       real wv1D(kms:kme)              ! turbulent flux of momentum (y component)
 155       real wt1D(kms:kme)              ! turbulent flux of temperature
 156       real wq1D(kms:kme)              ! turbulent flux of water vapor
 157       real wqc1D(kms:kme)              ! turbulent flux of liquid water
 158       integer iz_pbl(ims:ime,jms:jme)
 159       real pr1D(kms:kme)
 160       real raten1D(kms:kme)
 161       real p8w1D(kms:kme)
 162       real pi_1D(kms:kme)
 163 ! local added only for diagnostic output
 164       real bu(ims:ime,kms:kme,jms:jme) ! buoyancy term in TKE
 165       real sh(ims:ime,kms:kme,jms:jme) ! shear term in TKE
 166       real wrk(ims:ime) ! working array
 167       integer ix,iy,iz,id
 168       real ufrac_int                                              ! urban fraction
 169       real psim1D,psih1D,znt1D,br1D,sflux
 170       integer xtime
 171
 172 !
 173 !
 174
 175
 176
 177
 178 !here I fix the value of the reference state equal to the value of the potnetial temperature
 179 ! the only use of this variable in the code is to compute the paramter BETA = g/th0
 180 ! I fix it to 300K.
 181          rublten=0.
 182          rvblten=0.
 183          rthblten=0.
 184          rqvblten=0.
 185          rqcblten=0.
 186
 187         do ix=its,ite
 188         do iy=jts,jte
 189         do iz=kts,kte
 190          th_0(ix,iz,iy)=300.
 191         enddo
 192         enddo
 193         enddo
 194
 195
 196        z1D=0.
 197        dz1D=0.
 198        u1D =0.
 199        v1D =0.
 200        th1D=0.
 201        q1D=0.
 202        rho1D=0.
 203        rhoz1D=0.
 204        tke1D =0.
 205        diss1D=0.
 206        tpe1D=0.
 207        raten1D=0.
 208        pi_1D=0.
 209        p8w1D=0.
 210        th01D =0.
 211        exch1D=0.
 212        sf1D  =1.
 213        vl1D  =1.
 214        a_u1D =0.
 215        a_v1D =0.
 216        a_t1D =0.
 217        a_q1D =0.
 218        a_qc1D =0.
 219        b_u1D =0.
 220        b_v1D =0.
 221        b_t1D =0.
 222        b_q1D =0.
 223        b_qc1D =0.
 224        b_e1D=0.
 225        a_e1D=0.
 226        b_d1D=0.
 227        a_d1D=0.
 228        b_tpe1D=0.
 229        a_tpe1D=0.
 230        sh1D  =0.
 231        bu1D  =0.
 232        wu1D  =0.
 233        wv1D  =0.
 234        wt1D =0.
 235        wq1D =0.
 236        iz_pbl=0
 237        psim1D=0.
 238        psih1D=0.
 239
 240
 241        do ix=its,ite
 242        do iy=jts,jte
 243          z1d(kts)=0.
 244          br1D=br(ix,iy)
 245          sflux=hfx(ix,iy)/rho(ix,1,iy)/cp+qfx(ix,iy)/rho(ix,1,iy)*(461./287.-1)*th_phy(ix,1,iy)
 246          znt1D=znt(ix,iy)
 247          psim1D=psim(ix,iy)
 248          psih1D=psih(ix,iy)
 249          do iz= kts,kte
 250          if(xtime.le.2)then
 251           tke_pbl(ix,iz,iy)=temin
 252           diss_pbl(ix,iz,iy)=dmin
 253           tpe_pbl(ix,iz,iy)=tpemin
 254           tke_adv(ix,iz,iy)=temin
 255           diss_adv(ix,iz,iy)=dmin
 256           tpe_adv(ix,iz,iy)=tpemin
 257           pr_pbl(ix,iz,iy)=1.
 258          endif
 259
 260           u1D(iz)=u_phy(ix,iz,iy)
 261           v1D(iz)=v_phy(ix,iz,iy)
 262           th1D(iz)=th_phy(ix,iz,iy)
 263           q1D(iz)=qv_curr(ix,iz,iy)
 264           qc1D(iz)=qc_curr(ix,iz,iy)
 265           rho1D(iz)=rho(ix,iz,iy)
 266           th01D(iz)=th_0(ix,iz,iy)
 267           dz1D(iz)=dz8w(ix,iz,iy)
 268           pi_1D(iz)=pi_phy(ix,iz,iy)
 269           raten1D(iz)=rthraten(ix,iz,iy)
 270           p8w1D(iz)=p8w(ix,iz,iy)
 271           pr1D(iz)=pr_pbl(ix,iz,iy)
 272           z1D(iz+1)=z1D(iz)+dz1D(iz)
 273          if(bl_keps_adv==1)then
 274           tke_pbl(ix,iz,iy)=tke_adv(ix,iz,iy)
 275           diss_pbl(ix,iz,iy)=diss_adv(ix,iz,iy)
 276           tpe_pbl(ix,iz,iy)=tpe_adv(ix,iz,iy)
 277          endif
 278           tke1D(iz)=min(temax,max(temin,tke_pbl(ix,iz,iy)))
 279           diss1D(iz)=min(dmax,max(dmin,diss_pbl(ix,iz,iy)))
 280           tpe1D(iz)=min(tpemax,max(tpemin,tpe_pbl(ix,iz,iy)))
 281          enddo
 282           rhoz1D(kts)=rho1D(kts)
 283          do iz=kts+1,kte
 284           rhoz1D(iz)=(rho1D(iz)*dz1D(iz-1)+rho1D(iz-1)*dz1D(iz))/(dz1D(iz-1)+dz1D(iz))
 285          enddo
 286           rhoz1D(kte+1)=rho1D(kte)
 287          do iz=kts,kte
 288           vl1D(iz)=1.
 289           sf1D(iz)=1.
 290          enddo
 291          ufrac_int=0.
 292          sf1D(kte+1)=1.
 293
 294 ! compute the pbl_height
 295           call GET_PBLH(kts,kte,pblh(ix,iy),th1D,tke1D,z1D,dz1D,q1D,iz_pbl(ix,iy))
 296           !call pbl_height(kms,kme,kts,kte,dz1D,z1D,th1D,q1D,pblh(ix,iy),iz_pbl(ix,iy))
 297
 298
 299 ! estimate the tendencies
 300          do iz=kts,kte
 301           a_t1D(iz)=0.
 302           b_t1D(iz)=0.
 303           a_u1D(iz)=0.
 304           b_u1D(iz)=0.
 305           a_v1D(iz)=0.
 306           b_v1D(iz)=0.
 307           a_q1D(iz)=0.
 308           b_q1D(iz)=0.
 309           b_e1D(iz)=0.
 310           a_e1D(iz)=0.
 311           b_d1D(iz)=0.
 312           a_d1D(iz)=0.
 313           b_tpe1D(iz)=0.
 314           a_tpe1D(iz)=0.
 315          enddo
 316           b_t1D(1)=hfx(ix,iy)/dz1D(1)/rho1D(1)/cp
 317           b_q1D(1)=qfx(ix,iy)/dz1D(1)/rho1D(1)
 318           a_u1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((max(1.E-10,abs(u1D(1)))**2.+max(1.E-10,abs(v1D(1)))**2.)**.5))
 319           a_v1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((max(1.E-10,abs(u1D(1)))**2.+max(1.E-10,abs(v1D(1)))**2.)**.5))
 320
 321         if(xtime.gt.2)then
 322
 323         call keps1D(mol(ix,iy),tsk(ix,iy),xtime,ix,iy,ids,ide,jds,jde,kms,kme,kts,kte,dz1D,z1D,dt,                          &
 324                            u1D,v1D,th1D,pi_1D,raten1D,p8w1D,                                          &
 325                            rho1D,rhoz1D,q1D,qc1D,th01D,tke1D,diss1D,tpe1D,pr1D,                       &
 326                            a_u1D,b_u1D,a_v1D,b_v1D,a_t1D,b_t1D,a_q1D,b_q1D,a_qc1D,b_qc1D,             &
 327                            b_e1D,a_e1D,b_d1D,a_d1D,b_tpe1D,a_tpe1D,psim1D,psih1D,                     &
 328                            ustar(ix,iy),hfx(ix,iy),qfx(ix,iy),sflux,cp,g,                             &
 329                            sf1D,vl1D,exch1D,sh1D,bu1D,br1D,znt1D,                                     &
 330                            pblh(ix,iy),iz_pbl(ix,iy),ufrac_int)
 331         endif
 332
 333
 334          do iz= kts,kte
 335           tke_pbl(ix,iz,iy)=tke1D(iz)
 336           diss_pbl(ix,iz,iy)=diss1D(iz)
 337           tpe_pbl(ix,iz,iy)=tpe1D(iz)
 338        if(bl_keps_adv==1)then
 339         tke_adv(ix,iz,iy)=(tke_pbl(ix,iz,iy))
 340         diss_adv(ix,iz,iy)=diss_pbl(ix,iz,iy)
 341         tpe_adv(ix,iz,iy)=tpe_pbl(ix,iz,iy)
 342        endif
 343           pr_pbl(ix,iz,iy)=pr1D(iz)
 344           bu(ix,iz,iy)=bu1D(iz)
 345           exch_h(ix,iz,iy)=exch1D(iz)/pr1D(iz)
 346           exch_m(ix,iz,iy)=exch1D(iz)
 347         enddo
 348
 349
 350 ! compute the tendencies
 351
 352          do iz= kts,kte
 353           rthblten(ix,iz,iy)=rthblten(ix,iz,iy)+(th1D(iz)-th_phy(ix,iz,iy))/dt
 354           rqvblten(ix,iz,iy)=rqvblten(ix,iz,iy)+(q1D(iz)-qv_curr(ix,iz,iy))/dt
 355           rqcblten(ix,iz,iy)=rqcblten(ix,iz,iy)+(qc1D(iz)-qc_curr(ix,iz,iy))/dt
 356           rublten(ix,iz,iy)=rublten(ix,iz,iy)+(u1D(iz)-u_phy(ix,iz,iy))/dt
 357           rvblten(ix,iz,iy)=rvblten(ix,iz,iy)+(v1D(iz)-v_phy(ix,iz,iy))/dt
 358         enddo
 359      ! endif
 360
 361          wt(ix,kts,iy)=hfx(ix,iy)/rho1D(kts)/cp-rho1D(kts)*(th1D(kts)-th_phy(ix,kts,iy))/dt*dz1D(kts)
 362          wu(ix,kts,iy)=-ustar(ix,iy)*ustar(ix,iy)-rho1D(kts)*(u1D(kts)-u_phy(ix,kts,iy))/dt*dz1D(kts)
 363          wv(ix,kts,iy)=-ustar(ix,iy)*ustar(ix,iy)-rho1D(kts)*(u1D(kts)-u_phy(ix,kts,iy))/dt*dz1D(kts)
 364          do iz=kts+1,kte
 365           wt(ix,iz,iy)=wt(ix,iz-1,iy)-rho1D(iz)*(th1D(iz)-th_phy(ix,iz,iy))/dt*dz1D(iz)
 366           wu(ix,iz,iy)=wu(ix,iz-1,iy)-rho1D(iz)*(u1D(iz)-u_phy(ix,iz,iy))/dt*dz1D(iz)
 367           wv(ix,iz,iy)=wv(ix,iz-1,iy)-rho1D(iz)*(v1D(iz)-v_phy(ix,iz,iy))/dt*dz1D(iz)
 368         enddo
 369
 370
 371       enddo  ! iy
 372
 373
 374       enddo  ! ix
 375
 376
 377       return
 378       end subroutine keps
 379
 380
 381 ! ===6=8===============================================================72
 382
 383 ! ===6=8===============================================================72
 384
 385
 386
 387      subroutine  keps1D(mol,tsk,xtime,ix,iy,ids,ide,jds,jde,kms,kme,kts,kte,dz,z,dt,          &
 388                         u,v,th,pi1d,raten,p8w1D,rho,rhoz,qa,qc,th0,te,d,tpe,Pra,      &
 389                         a_u,b_u,a_v,b_v,a_t,b_t,a_q,b_q,a_qc,b_qc,b_e,a_e,            &
 390                         b_d,a_d,a_tpe,b_tpe,                                          &
 391                         psim,psih,ustar,hfx,qfx,sflux,cp,g,                          &
 392                         sf,vl,exch,sh,bu,b_ric,z0,pblh,iz_pbl,ufrac_int)
 393
 394 ! ----------------------------------------------------------------------
 395 ! 1D resolution of TKE and dissipation rate following k epsilon turbulent closure
 396 ! ----------------------------------------------------------------------
 397
 398       implicit none
 399       integer iz,ix,iy,xtime
 400       integer kms,kme,kts,kte
 401       integer ids,ide,jds,jde
 402       real z(kms:kme)               ! Altitude above the ground of the cell interfaces.
 403       real dz(kms:kme)                ! vertical resolution
 404       real dt,time                         ! Time step
 405       real pblh
 406       real tsk,mol
 407       real u(kms:kme)                ! Wind speed in the x direction
 408       real v(kms:kme)                ! Wind speed in the y direction
 409       real th(kms:kme)                ! Potential temperature
 410       real rho(kms:kme)                ! Air density
 411       real rhoz(kms:kme) !air density at the faces of the cell
 412       real qa(kms:kme)                ! air humidity
 413       real qc(kms:kme)                ! air humidity
 414       real th0(kms:kme)               ! Reference potential temperature
 415       real te(kms:kme)                ! turbulent kinetic energy
 416       real d(kms:kme)                 ! dissipation rate
 417       real tpe(kms:kme)
 418       real pi1d(kms:kme)
 419       real raten(kms:kme)
 420       real p8w1D(kms:kme)
 421       real a_u(kms:kme)
 422       real b_u(kms:kme)
 423       real a_v(kms:kme)
 424       real b_v(kms:kme)
 425       real a_t(kms:kme)
 426       real b_t(kms:kme)
 427       real a_q(kms:kme)
 428       real b_q(kms:kme)
 429       real a_qc(kms:kme)
 430       real b_qc(kms:kme)
 431       real a_e(kms:kme)
 432       real b_e(kms:kme)
 433       real a_d(kms:kme)
 434       real b_d(kms:kme)
 435       real a_tpe(kms:kme)
 436       real b_tpe(kms:kme)
 437       real ufrac_int
 438       real ustar                 ! ustar
 439       real hfx                   ! sensbile heat flux
 440       real qfx                   ! kinematic latent heat flux
 441       real g                     ! gravity
 442       real cp                    !
 443       real sf(kms:kme)              ! surface of the urban grid cells
 444       real vl(kms:kme)                ! volume of the urban grid cells
 445       real exch(kms:kme) ! turbulent diffusion coefficients (defined at the faces)
 446       real sh(kms:kme)    ! shear term in TKE eqn.
 447       real bu(kms:kme)    ! buoyancy term in TKE eqn.
 448       real b_ric,z0,sflux
 449       integer iz_pbl
 450       real we(kms:kme)
 451       real wd(kms:kme)
 452       real wtpe(kms:kme)
 453       real wu(kms:kme)
 454       real wv(kms:kme)
 455       real wt(kms:kme)
 456       real wq(kms:kme)
 457       real wqc(kms:kme)
 458       real Ri(kms:kme)
 459       real Pra(kms:kme),coeff(kms:kme)
 460       real Pra_st(kms:kme)
 461       real exch_h(kms:kme)
 462       real nsquare(kms:kme)
 463       real dls(kms:kme)
 464       real psim,psih
 465       real phim,phih,phieps
 466       real dtmdz(kms:kme)
 467
 468        call shear(kms,kme,kts,kte,u,v,dz,sh,sf)
 469
 470        call c_nsquare(Pra,ustar,sflux,z0,tsk,mol,pi1d,g,kms,kme,kts,kte,th,th0,dz,nsquare,dtmdz)
 471
 472        call surface_bl_pra_ri(kms,kme,kts,kte,dz,z,rho,g,cp,z0,sflux,raten,pi1d,p8w1D,   &
 473                               th,qa,sh,nsquare,Ri,b_ric,psim,psih,pblh,ustar, &
 474                               iz_pbl,Pra,Pra_st,phim,phih,phieps,dtmdz(kts))
 475
 476        call  buoy(g,kms,kme,kts,kte,th,th0,dz,bu,te,tpe,Pra)
 477
 478        call  cdtur_ke(kms,kme,kts,kte,te,d,z,dz,exch,exch_h,Pra_st)
 479
 480        call  length_bougeault(g,kms,kme,kts,kte,z,dz,te,th,th0,dls)
 481 !solve analytically the equation for tke and dissipation,following Zonato et
 482 !al.2022
 483         call thetaphi_alphabeta(kms,kme,kts,kte,z0,dz,sh,bu,te,d,dt,Pra,dls,ustar,phieps,phim,Ri,nsquare,a_d)
 484 ! solve diffusion equation for tke
 485         call diff(kms,kme,kts,kte,2,1,dt,te,rho,rhoz,exch,a_e,b_e,sf,vl,dz,we)
 486 ! solve diffusion equation for dissipation rate
 487         call diff(kms,kme,kts,kte,2,1,dt,d,rho,rhoz,exch/sigma_eps,a_d,b_d,sf,vl,dz,wd)
 488
 489 !calculation of the terms for turbulent potential energy equation and addition
 490 !of the countergradient term to tke
 491        call terms_tpe(dtmdz,g,kms,kme,kts,kte,th0,te,d,tpe,a_tpe,b_tpe,Pra,coeff)
 492        call calc_countergradient(dz,g,kms,kme,kts,kte,th0,te,d,tpe,Pra,b_t)
 493 !solve diffusion equation for turbulent potential energy
 494        call diff(kms,kme,kts,kte,1,1,dt,tpe,rho,rhoz,exch,a_tpe,b_tpe,sf,vl,dz,wtpe)
 495
 496 ! solve diffusion equation for momentum x component
 497         call diff(kms,kme,kts,kte,1,1,dt,u,rho,rhoz,exch,a_u,b_u,sf,vl,dz,wu)
 498
 499 ! solve diffusion equation for momentum y component
 500         call diff(kms,kme,kts,kte,1,1,dt,v,rho,rhoz,exch,a_v,b_v,sf,vl,dz,wv)
 501 ! solve diffusion equation for potential temperature
 502         call diff(kms,kme,kts,kte,1,1,dt,th,rho,rhoz,exch_h,a_t,b_t,sf,vl,dz,wt)
 503
 504 ! solve diffusion equation for water vapor mixing ratio
 505         call diff(kms,kme,kts,kte,1,1,dt,qa,rho,rhoz,exch_h,a_q,b_q,sf,vl,dz,wq)
 506
 507 ! solve diffusion equation for liquid water mixing ratio
 508         call diff(kms,kme,kts,kte,1,1,dt,qc,rho,rhoz,exch_h,a_qc,b_qc,sf,vl,dz,wqc)
 509
 510 do iz=kts,kte
 511 te(iz)=min(max(temin,te(iz)),temax)
 512 d(iz)=min(max(dmin,d(iz)),dmax)
 513 tpe(iz)=min(max(tpemin,tpe(iz)),tpemax)
 514 enddo
 515
 516 return
 517 end subroutine keps1D
 518
 519 ! ===6=8===============================================================72
 520 ! ===6=8===============================================================72
 521
 522 subroutine cdtur_ke(kms,kme,kts,kte,te,d,z,dz,exch,exch_h,Pra_st)
 523
 524 ! compute turbulent coefficients
 525 implicit none
 526 integer iz,kms,kme,kts,kte
 527 real te(kms:kme),exch(kms:kme),exch_h(kms:kme),d(kms:kme)
 528 real dz(kms:kme),z(kms:kme),Pra_st(kms:kme)
 529 real fact
 530
 531 exch(kts)=0.
 532 exch_h(kts)=0.
 533 do iz=kts+1,kte
 534 exch(iz)=(c_mu*te(iz-1)**2/d(iz-1)*dz(iz-1)+c_mu*te(iz)**2/d(iz)*dz(iz))/(dz(iz)+dz(iz-1))
 535 exch(iz)=max(exch(iz),0.09)
 536 exch_h(iz)=exch(iz)/Pra_st(iz)
 537 enddo
 538 exch(kte+1)=exch(kte)
 539 exch_h(kte+1)=exch_h(kte)
 540 return
 541 end subroutine cdtur_ke
 542
 543
 544 ! ===6=8===============================================================72
 545 ! ===6=8===============================================================72
 546
 547 subroutine diff(kms,kme,kts,kte,iz1,izf,dt,co,rho,rhoz,cd,aa,bb,sf,vl,dz,fc)
 548
 549
 550 !------------------------------------------------------------------------
 551 !           Calculation of the diffusion in 1D
 552 !------------------------------------------------------------------------
 553 !  - Input:
 554 !       nz    : number of points
 555 !       iz1   : first calculated point
 556 !       co    : concentration of the variable of interest
 557 !       dz    : vertical levels
 558 !       cd    : diffusion coefficients
 559 !       dtext : external time step
 560 !       rho    : density of the air at the center
 561 !       rhoz   : density of the air at the face
 562 !       itest : if itest eq 1 then update co, else store in a flux array
 563 !  - Output:
 564 !       co :concentration of the variable of interest
 565
 566 !  - Internal:
 567 !       cddz  : constant terms in the equations
 568 !       dt    : diffusion time step
 569 !       nt    : number of the diffusion time steps
 570 !       cstab : ratio of the stability condition for the time step
 571 !--------------------------------------------------------------
 572
 573 implicit none
 574 integer iz,iz1,izf
 575 integer kms,kme,kts,kte
 576 real dt,dzv
 577 real co(kms:kme),cd(kms:kme),dz(kms:kme)
 578 real rho(kms:kme),rhoz(kms:kme)
 579 real cddz(kms:kme+1),fc(kms:kme),df(kms:kme)
 580 real a(kms:kme,3),c(kms:kme)
 581 real sf(kms:kme),vl(kms:kme)
 582 real aa(kms:kme),bb(kms:kme)
 583
 584 ! Compute cddz=2*cd/dz
 585  cddz(kts)=sf(kts)*rhoz(kts)*cd(kts)/dz(kts)
 586 do iz=kts+1,kte
 587  cddz(iz)=2.*sf(iz)*rhoz(iz)*cd(iz)/(dz(iz)+dz(iz-1))
 588 enddo
 589  cddz(kte+1)=sf(kte+1)*rhoz(kte+1)*cd(kte+1)/dz(kte)
 590
 591  do iz=kts,iz1-1
 592   a(iz,1)=0.
 593   a(iz,2)=1.
 594   a(iz,3)=0.
 595   c(iz)=co(iz)
 596  enddo
 597
 598   do iz=iz1,kte-izf
 599    dzv=vl(iz)*dz(iz)
 600    a(iz,1)=-cddz(iz)*dt/dzv/rho(iz)
 601    a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dzv/rho(iz)-aa(iz)*dt
 602    a(iz,3)=-cddz(iz+1)*dt/dzv/rho(iz)
 603    c(iz)=co(iz)+bb(iz)*dt
 604   enddo
 605
 606   do iz=kte-(izf-1),kte
 607    a(iz,1)=0.
 608    a(iz,2)=1
 609    a(iz,3)=0.
 610    c(iz)=co(iz)
 611   enddo
 612
 613   call invert (kms,kme,kts,kte,a,c,co)
 614
 615   do iz=kts,iz1
 616    fc(iz)=0.
 617   enddo
 618
 619   do iz=iz1+1,kte
 620    fc(iz)=-(cddz(iz)*(co(iz)-co(iz-1)))/rho(iz)
 621   enddo
 622
 623 return
 624 end subroutine diff
 625
 626 ! ===6=8===============================================================72
 627 ! ===6=8===============================================================72
 628
 629 subroutine buoy(g,kms,kme,kts,kte,th,th0,dz,bu,te,tpe,Pra)
 630        implicit none
 631        integer kms,kme,kts,kte,iz
 632        real dtdz1,dtdz2,g,dtmdz
 633        real th(kms:kme),dz(kms:kme),bu(kms:kme),sf(kms:kme)
 634        real th0(kms:kme),phi_bu(kms:kme)
 635        real te(kms:kme),tpe(kms:kme),Pra(kms:kme)
 636
 637         bu(kts)=0.
 638        do iz=kts+1,kte-1
 639         phi_bu(iz)=g/th0(iz)*tpe(iz)/te(iz)*Pra(iz)*c_theta/c_mu
 640         dtdz1=2.*(th(iz)-th(iz-1))/(dz(iz-1)+dz(iz))
 641         dtdz2=2.*(th(iz+1)-th(iz))/(dz(iz+1)+dz(iz))
 642         dtmdz=0.5*(dtdz1+dtdz2)
 643         bu(iz)=-dtmdz*g/th0(iz)+(g/th0(iz))*phi_bu(iz)
 644        enddo
 645         bu(kte)=0.
 646
 647 return
 648 end subroutine buoy
 649
 650 ! ===6=8===============================================================72
 651 ! ===6=8===============================================================72
 652
 653 subroutine c_nsquare(Pra,ustar,sflux,z0,tsk,mol,pi,g,kms,kme,kts,kte,th,th0,dz,nsquare,dtmdz)
 654        implicit none
 655        integer kms,kme,kts,kte,iz
 656        real thetaz0,sflux,ustar,z0,dtdz1,dtdz2,g,tsk,mol,th0(kms:kme)
 657        real th(kms:kme),dz(kms:kme),nsquare(kms:kme),sf(kms:kme)
 658        real te(kms:kme),d(kms:kme),dtmdz(kms:kme),Pra(kms:kme)
 659        real pi(kms:kme)
 660          dtmdz(kts)=(-sflux/ustar)/vk/(dz(kts)/2.-z0)
 661          nsquare(kts)=-g/th0(kts)*dtmdz(kts)
 662         do iz=kts+1,kte-1
 663          dtdz1=2.*(th(iz)-th(iz-1))/(dz(iz-1)+dz(iz))
 664          dtdz2=2.*(th(iz+1)-th(iz))/(dz(iz+1)+dz(iz))
 665          dtmdz(iz)=0.5*(dtdz1+dtdz2)
 666          nsquare(iz)=-dtmdz(iz)*g/th0(iz)
 667         enddo
 668          nsquare(kte)=0.
 669
 670 return
 671 end subroutine c_nsquare
 672
 673 ! ===6=8===============================================================72
 674 ! ===6=8===============================================================72
 675
 676 subroutine shear(kms,kme,kts,kte,u,v,dz,sh,sf)
 677 ! compute shear term
 678        implicit none
 679        integer kms,kme,kts,kte,iz,ix,iy
 680        real dudz1,dudz2,dvdz1,dvdz2,dumdz
 681        real u(kms:kme),v(kms:kme),dz(kms:kme),sh(kms:kme),sf(kms:kme)
 682        real u1,u2,v1,v2,ufrac_int
 683
 684         sh(kts)=0.
 685        do iz=kts+1,kte-1
 686         u2=(dz(iz+1)*u(iz)+dz(iz)*u(iz+1))/(dz(iz)+dz(iz+1))
 687         u1=(dz(iz)*u(iz-1)+dz(iz-1)*u(iz))/(dz(iz-1)+dz(iz))
 688         v2=(dz(iz+1)*v(iz)+dz(iz)*v(iz+1))/(dz(iz)+dz(iz+1))
 689         v1=(dz(iz)*v(iz-1)+dz(iz-1)*v(iz))/(dz(iz-1)+dz(iz))
 690         dumdz=((u2-u1)/dz(iz))**2.+((v2-v1)/dz(iz))**2.
 691         sh(iz)=dumdz
 692        enddo
 693         sh(kte)=0.
 694
 695 return
 696 end subroutine shear
 697
 698 ! ===6=8===============================================================72
 699 ! ===6=8===============================================================72
 700
 701 subroutine invert(kms,kme,kts,kte,a,c,x)
 702        implicit none
 703        integer in
 704        integer kts,kte,kms,kme
 705        real a(kms:kme,3),c(kms:kme),x(kms:kme)
 706
 707         do in=kte-1,kts,-1
 708          c(in)=c(in)-a(in,3)*c(in+1)/a(in+1,2)
 709          a(in,2)=a(in,2)-a(in,3)*a(in+1,1)/a(in+1,2)
 710         enddo
 711
 712         do in=kts+1,kte
 713          c(in)=c(in)-a(in,1)*c(in-1)/a(in-1,2)
 714         enddo
 715
 716         do in=kts,kte
 717
 718          x(in)=c(in)/a(in,2)
 719
 720         enddo
 721
 722         return
 723         end subroutine invert
 724
 725
 726
 727 !######################################################################
 728        subroutine pbl_height(kms,kme,kts,kte,dz,z,th,q,pblh,iz_pbl)
 729
 730 ! this routine computes the PBL height
 731 ! with an approach similar to MYNN
 732        implicit none
 733        integer kms,kme,kts,kte,iz
 734        real z(kms:kme),dz(kms:kme),th(kms:kme),q(kms:kme)
 735        real pblh
 736 !Local
 737        real thv(kms:kme),zc(kms:kme)
 738        real thsfc
 739        integer iz_pbl
 740 ! compute the height of the center of the grid cells
 741       do iz=kts,kte
 742        zc(iz)=z(iz)+dz(iz)/2.
 743       enddo
 744
 745 ! compute the virtual potential temperature
 746
 747        do iz=kts,kte
 748         thv(iz)=th(iz)*(1.+0.61*q(iz))
 749        enddo
 750 ! now compute the PBL height
 751
 752        pblh=0.
 753        thsfc=thv(kts)+0.5
 754        iz_pbl=1
 755        do iz=kts+1,kte
 756         if(pblh.eq.0.and.thv(iz).gt.thsfc)then
 757          pblh=zc(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(zc(iz)-zc(iz-1))
 758
 759         endif
 760        enddo
 761
 762        do iz=kts,kte
 763        if(zc(iz).le.pblh)then
 764        iz_pbl=iz_pbl+1
 765        endif
 766        enddo
 767        return
 768
 769        end subroutine pbl_height
 770
 771
 772 ! ===6=8===============================================================72
 773 ! ===6=8===============================================================72
 774
 775       SUBROUTINE GET_PBLH(KTS,KTE,zi,theta1D,tke1D,zw1D,dz1D,q1D,iz_pbl)
 776 ! Copied from MYNN PBL
 777
 778       !---------------------------------------------------------------
 779       !             NOTES ON THE PBLH FORMULATION
 780       !
 781       !The 1.5-theta-increase method defines PBL heights as the level at
 782       !which the potential temperature first exceeds the minimum potential
 783       !temperature within the boundary layer by 1.5 K. When applied to
 784       !observed temperatures, this method has been shown to produce PBL-
 785       !height estimates that are unbiased relative to profiler-based
 786       !estimates (Nielsen-Gammon et al. 2008). However, their study did not
 787       !include LLJs. Banta and Pichugina (2008) show that a TKE-based
 788       !threshold is a good estimate of the PBL height in LLJs. Therefore,
 789       !a hybrid definition is implemented that uses both methods, weighting
 790       !the TKE-method more during stable conditions (PBLH < 400 m).
 791       !A variable tke threshold (TKEeps) is used since no hard-wired
 792       !value could be found to work best in all conditions.
 793       !---------------------------------------------------------------
 794
 795       INTEGER,INTENT(IN) :: KTS,KTE
 796       REAL, INTENT(OUT) :: zi
 797       REAL, DIMENSION(KTS:KTE), INTENT(IN) ::  tke1D, dz1D,theta1D,q1D
 798       REAL, DIMENSION(KTS:KTE) :: thetav1D
 799       REAL, DIMENSION(KTS:KTE+1), INTENT(IN) :: zw1D
 800       !LOCAL VARS
 801       REAL ::  PBLH_TKE,tke,tkem1,wt,maxtke,TKEeps,minthv
 802       REAL :: delt_thv   !delta theta-v; dependent on land/sea point
 803       REAL, PARAMETER :: sbl_lim  = 200. !Theta-v PBL lower limit of trust (m).
 804       REAL, PARAMETER :: sbl_damp = 400. !Damping range for averaging with TKE-based PBLH (m).
 805       INTEGER :: I,J,K,kthv,ktke,iz_pbl
 806
 807         do iz=kts,kte
 808         thetav1D(iz)=theta1D(iz)*(1.+0.61*q1D(iz))
 809        enddo
 810
 811
 812
 813       !FIND MAX TKE AND MIN THETAV IN THE LOWEST 500 M
 814       k = kts+1
 815       kthv = 1
 816       ktke = 1
 817       maxtke = 0.
 818       minthv = 9.E9
 819
 820       DO WHILE (zw1D(k) .LE. 500.)
 821         tke  =MAX(tke1D(k),0.)   ! maximum QKE
 822          IF (maxtke < ttke) then
 823             maxtke = ttke
 824             ktke = k
 825          ENDIF
 826          IF (minthv > thetav1D(k)) then
 827              minthv = thetav1D(k)
 828              kthv = k
 829          ENDIF
 830          k = k+1
 831       ENDDO
 832       TKEeps = maxtke/20.
 833       TKEeps = MAX(TKEeps,0.025/2.)
 834       TKEeps = MIN(TKEeps,0.25/2.)
 835
 836       !FIND THETAV-BASED PBLH (BEST FOR DAYTIME).
 837       zi=0.
 838       k = kthv+1
 839           delt_thv = 1.5
 840       DO WHILE (zi .EQ. 0.)
 841          IF (thetav1D(k) .GE. (minthv + delt_thv))THEN
 842             zi = zw1D(k) - dz1D(k-1)* &
 843               & MIN((thetav1D(k)-(minthv+delt_thv))/MAX(thetav1D(k)-thetav1D(k-1),1E-6),1.0)
 844         ENDIF
 845         k = k+1
 846          IF (k .EQ. kte-1) zi = zw1D(kts+1) !EXIT SAFEGUARD
 847       ENDDO
 848
 849       PBLH_TKE=0.
 850       k = ktke+1
 851      DO WHILE (PBLH_TKE .EQ. 0.)
 852          tke  =MAX(tke1D(k)/2.,0.)      ! maximum TKE
 853          tkem1=MAX(tke1D(k-1)/2.,0.)
 854          IF (tke .LE. TKEeps) THEN
 855                PBLH_TKE = zw1D(k) - dz1D(k-1)* &
 856                & MIN((TKEeps-tke)/MAX(tkem1-tke, 1E-6), 1.0)
 857              !IN CASE OF NEAR ZERO TKE, SET PBLH = LOWEST LEVEL.
 858              PBLH_TKE = MAX(PBLH_TKE,zw1D(kts+1))
 859              !print *,"PBLH_TKE:",i,j,PBLH_TKE, Qke1D(k)/2., zw1D(kts+1)
 860          ENDIF
 861          k = k+1
 862          IF (k .EQ. kte-1) PBLH_TKE = zw1D(kts+1) !EXIT SAFEGUARD
 863       ENDDO
 864
 865     !BLEND THE TWO PBLH TYPES HERE:
 866
 867       wt=.5*TANH((zi - sbl_lim)/sbl_damp) + .5
 868       zi=PBLH_TKE*(1.-wt) + zi*wt
 869
 870
 871        iz_pbl=1
 872        do k=kts,kte
 873        if(zw1D(k).le.zi)then
 874        iz_pbl=iz_pbl+1
 875        endif
 876        enddo
 877
 878
 879    END SUBROUTINE GET_PBLH
 880
 881 ! ===6=8===============================================================72
 882 ! ===6=8===============================================================72
 883 subroutine thetaphi_alphabeta(kms,kme,kts,kte,z0,dz,sh,bu,te,d,dt,Pra,dls,ustar,phieps,phim,Ri,nsquare,a_d)
 884
 885       implicit none
 886       integer iz,kms,kme,kts,kte
 887       real te(kms:kme)                ! turbulent kinetic energy
 888       real d(kms:kme)                 ! dissipation rate
 889       real a_d(kms:kme)
 890       real dt,z0                   ! Time step
 891       real sh(kms:kme)    ! shear term in TKE eqn.
 892       real bu(kms:kme)    ! buoyancy term in TKE eqn.
 893       real A(kms:kme)
 894       real B(kms:kme)
 895       real theta(kms:kme)
 896       real phi(kms:kme)
 897       real theta_new(kms:kme)
 898       real phi_new(kms:kme)
 899       real C(kms:kme)
 900       real alpha
 901       real beta
 902       real F
 903       real dz(kms:kme)
 904       real Pra(kms:kme)
 905       real dls(kms:kme)
 906       real Ri(kms:kme)
 907       real nsquare(kms:kme)
 908       real ustar,phieps,phim
 909        alpha=1.
 910        beta=-1./c2
 911        F=c2-1.
 912
 913         do iz=kts+1,kte
 914         if(abs(te(iz))/abs(d(iz)).gt.1E4.and.d(iz).le.1E-6)then
 915         d(iz)=1/1.4*te(iz)**(3./2.)/dls(iz)
 916         endif
 917         theta(iz)=max(temin,te(iz))/max(dmin,d(iz))
 918         phi(iz)=max(dmin,d(iz))**beta*max(temin,te(iz))**alpha
 919         C(iz)=(c_mu*((c1-1.)*sh(iz)+((c3-1.)*bu(iz))/Pra(iz)))
 920         A(iz)=c_mu*(sh(iz)+bu(iz)/Pra(iz))
 921         B(iz)=c_mu*(c1*sh(iz)+c3*bu(iz)/Pra(iz))
 922         if(C(iz).lt.0.)then
 923         C(iz)=-C(iz)
 924         theta_new(iz)=(c2-1.)/(C(iz)*theta(iz))+((c2-1.)**(1./2.)*(C(iz)*theta(iz)**2-c2+1))/(C(iz)*theta(iz)*((c2-1.)**(1./2.)+C(iz)**(1./2.)*theta(iz)*tanh(C(iz)**(1./2.)*dt*(c2-1)**(1./2.))))
 925         elseif(C(iz).eq.0.)then
 926         theta_new(iz)=theta(iz)+F*dt;
 927         else
 928         theta_new(iz)=((F/C(iz))**(0.5)*(tanh(dt*(C(iz)*F)**(0.5)) +theta(iz)*(C(iz)/F)**(0.5)))/(theta(iz)*tanh(dt*(C(iz)*F)**(0.5))*(C(iz)/F)**(0.5) + 1.)
 929         endif
 930         theta_new(iz)=max(1.,theta_new(iz))
 931         phi_new(iz)=(phi(iz)*exp(theta_new(iz)*dt*(alpha*A(iz)+beta*B(iz))))
 932         te(iz)=((phi_new(iz)*theta_new(iz)**(beta))**(1./(alpha+beta)))
 933         d(iz)=(te(iz)/theta_new(iz))
 934         enddo
 935         do iz=kts,kte
 936        if(Ri(iz).gt.0.)then
 937        a_d(iz)=a_d(iz)+c4*min(1.,sqrt(Ri(iz)/c5))*sqrt(max(0.,-nsquare(iz)))
 938        endif
 939        enddo
 940        d(kts)=(ustar**3./vk/(dz(1)/2.-z0))*phieps*(1.+a_d(kts)*d(kts)*dt)
 941        te(kts)=(ustar**2.)/sqrt(c_mu)*sqrt(phieps/phim)
 942
 943 return
 944 end subroutine thetaphi_alphabeta
 945
 946 ! ===6=8===============================================================72
 947 ! ===6=8===============================================================72
 948
 949 subroutine length_bougeault(g,kms,kme,kts,kte,z,dz,te,th,th0,dls)
 950 ! compute the length scales up and down
 951          implicit none
 952          integer kms,kme,kts,kte,iz,izz,ix,iy
 953          real dzt,zup,beta,zup_inf,bbb,tl,zdo,zdo_sup,zzz,g
 954          real te(kms:kme),dlu(kms:kme),dld(kms:kme),dz(kms:kme)
 955          real z(kms:kme),th(kms:kme),th0(kms:kme),dls(kms:kme)
 956
 957          do iz=kts,kte
 958           zup=0.
 959           dlu(iz)=z(kte+1)-z(iz)-dz(iz)/2.
 960           zzz=0.
 961           zup_inf=0.
 962           beta=g/th0(iz)      !Buoyancy coefficient
 963           do izz=iz,kte-1
 964            dzt=(dz(izz+1)+dz(izz))/2.
 965            zup=zup-beta*th(iz)*dzt
 966            zup=zup+beta*(th(izz+1)+th(izz))*dzt/2.
 967            zzz=zzz+dzt
 968            if(te(iz).lt.zup.and.te(iz).ge.zup_inf)then
 969             bbb=(th(izz+1)-th(izz))/dzt
 970             if(bbb.ne.0)then
 971              tl=(-beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zup_inf))))/bbb/beta
 972             else
 973              if(th(izz).ne.th(iz))then
 974               tl=(te(iz)-zup_inf)/(beta*(th(izz)-th(iz)))
 975              else
 976               tl=0.
 977              endif
 978             endif
 979             dlu(iz)=zzz-dzt+tl
 980            endif
 981            zup_inf=zup
 982           enddo
 983
 984           zdo=0.
 985           zdo_sup=0.
 986           dld(iz)=z(iz)+dz(iz)/2.
 987           zzz=0.
 988           do izz=iz,kts+1,-1
 989            dzt=(dz(izz-1)+dz(izz))/2.
 990            zdo=zdo+beta*th(iz)*dzt
 991            zdo=zdo-beta*(th(izz-1)+th(izz))*dzt/2.
 992            zzz=zzz+dzt
 993            if(te(iz).lt.zdo.and.te(iz).ge.zdo_sup)then
 994             bbb=(th(izz)-th(izz-1))/dzt
 995             if(bbb.ne.0.)then
 996              tl=(beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zdo_sup))))/bbb/beta
 997             else
 998              if(th(izz).ne.th(iz))then
 999               tl=(te(iz)-zdo_sup)/(beta*(th(izz)-th(iz)))
1000              else
1001               tl=0.
1002              endif
1003             endif
1004
1005             dld(iz)=zzz-dzt+tl
1006            endif
1007            zdo_sup=zdo
1008           enddo
1009           enddo
1010          do iz=kts,kte
1011           dld(iz)=min(dld(iz),(z(iz)+dz(iz)/2.))
1012           dls(iz)=sqrt(dlu(iz)*dld(iz))
1013          enddo
1014
1015
1016 return
1017 end subroutine length_bougeault
1018
1019 ! ===6=8===============================================================72
1020 ! ===6=8===============================================================72
1021
1022  subroutine surface_bl_pra_ri(kms,kme,kts,kte,dz,z,rho,g,cp,z0,sflux,raten,pi1d,p8w1D,  &
1023                            th,qa,sh,nsquare,Ri,b_ric,psim,psih,pblh,ustar,      &
1024                            iz_pbl,Pra,Pra_st,phim,phih,phieps,dtdzs)
1025
1026     implicit none
1027     integer kms,kme,kts,kte,iz,iz_pbl
1028     real cp,g,b_ric,psim,psih,pblh,ustar,sflux
1029     real phim,phih,phieps
1030     real radflux,radsum
1031     real z0,zz,dtdzs
1032     real dz(kms:kme),prnumfac(kms:kme),Pra(kms:kme),Pra_st(kms:kme),z(kms:kme+1)
1033     real th(kms:kme),qa(kms:kme),rho(kms:kme)
1034     real raten(kms:kme),pi1d(kms:kme),p8w1D(kms:kme)
1035     real sh(kms:kme),nsquare(kms:kme),Ri(kms:kme)
1036     real,parameter :: prmin=0.25,prmax=4.,sfcfrac=0.1,bfac=6.8
1037     real,parameter ::rimin=-100.,zfmin=1.e-8,aphi5 = 5.,aphi16 = 16.
1038     real prnum,prnum0,conpr,prfac,prfac2,wstar3,wstar3_2
1039     real zol1,zl1,hol1,phim_sl,phih_sl
1040     real bfx0,govrth,govrthv
1041     logical sfcflg
1042      Ri=1.
1043      Pra=4.
1044
1045     do iz=kts,kte
1046     if(iz.eq.kts)then
1047       Ri(iz)=b_ric
1048       else
1049       Ri(iz)=min(-nsquare(iz)/sh(iz),100.)
1050       if(sh(iz).lt.1e-15)then
1051       if(-nsquare(iz).gt.0.)then
1052       Ri(iz)=10.
1053       else
1054       Ri(iz)=-1.
1055       endif
1056       endif
1057     endif
1058     enddo
1059       Ri(kte)=Ri(kte-1)
1060
1061      if(iz_pbl.ge.2)then
1062       do iz=1,iz_pbl-1
1063          radflux=raten(iz)*pi1d(iz)
1064          radflux=radflux*cp/g*(p8w1D(iz)-p8w1D(iz+1))
1065       if (radflux < 0.0 ) radsum=abs(radflux)+radsum
1066       enddo
1067       endif
1068       radsum=max(radsum,0.0)/rho(iz_pbl-1)/cp
1069
1070      sfcflg=.true.
1071     if(b_ric.ge.0)sfcflg=.false.
1072      zl1=dz(1)/2.
1073      zz=(zl1+z0)/zl1
1074      zol1 = max(b_ric*psim*psim/psih,rimin)
1075     if(sfcflg)then
1076      zol1 = min(zol1,-zfmin)
1077     else
1078      zol1 = max(zol1,zfmin)
1079     endif
1080      hol1 = zol1*pblh/zl1*sfcfrac
1081
1082     if(sfcflg)then
1083        phim = (1.-aphi16*zol1*zz)**(-1./4.)
1084        phih = (1.-aphi16*zol1*zz)**(-1./2.)
1085        phieps=1.-zol1*zz
1086        phim_sl = (1.-aphi16*hol1)**(-1./4.)
1087        phih_sl = (1.-aphi16*hol1)**(-1./2.)
1088        bfx0 = max(sflux,0.)
1089        govrth=g/th(1)
1090        govrthv=g/(th(iz_pbl-1)*(1+0.606271777*qa(iz_pbl-1)))
1091        wstar3 =(govrth*bfx0*pblh)
1092        wstar3_2 =(govrthv*radsum*pblh)
1093      else
1094        phim_sl = (1.+aphi5*hol1)
1095        phih_sl= phim_sl
1096        phim=(1.+aphi5*zol1*zz)
1097        phieps=(1+2.5*(zol1*zz)**0.6)**(3./2.)
1098        phih= phim_sl
1099        wstar3=0.
1100        wstar3_2=0.
1101      endif
1102        conpr=bfac*vk*sfcfrac
1103      if(iz_pbl.gt.1)then
1104      do iz=kts,kte+1
1105      if(sfcflg)then
1106        prnumfac(iz) = -3.*(max(z(iz+1)-sfcfrac*pblh,0.))**2./pblh**2.
1107        prfac=conpr
1108        prfac2 = 15.9*(wstar3+wstar3_2)/ustar**3/(1.+4.*vk*(wstar3+wstar3_2)/ustar**3)
1109        prnum0 =(phih_sl/phim_sl+prfac)
1110        prnum0= prnum0/(1.+prfac2*vk*sfcfrac)
1111        prnum0 = max(min(prnum0,prmax),prmin)
1112        Pra(iz)= 1. + (prnum0-1.)*exp(prnumfac(iz))
1113      else
1114        Pra(iz)=max(1.,min(4.,1.+3*Ri(iz)))
1115      endif
1116      enddo
1117      endif
1118       dtdzs=-sflux/ustar/vk*phih/(dz(kts)/2-z0)
1119       Pra_st(kts)=0.
1120      do iz=kts+1,kte
1121       Pra_st(iz)=(Pra(iz-1)*dz(iz)+Pra(iz)*dz(iz-1))/(dz(iz)+dz(iz-1))
1122      enddo
1123       Pra_st(kte+1)=Pra_st(kte)
1124
1125 return
1126 end subroutine surface_bl_pra_ri
1127
1128 !==============================================================================
1129 !==============================================================================
1130
1131 subroutine terms_tpe(dtmdz,g,kms,kme,kts,kte,th0,te,d,tpe,a_tpe,b_tpe,Pra,coeff)
1132
1133      implicit none
1134      integer iz,kms,kme,kts,kte
1135      real g,th0(kms:kme),coeff(kms:kme)
1136      real te(kms:kme),d(kms:kme),tpe(kms:kme)
1137      real dtmdz(kms:kme),A_wt(kms:kme),a_tpe(kms:kme),b_tpe(kms:kme),Pra(kms:kme)
1138      do iz=kts,kte
1139       A_wt(iz)=1./te(iz)/tpe(iz)/2.*(c_mu/Pra(iz)*te(iz)**2./d(iz)*dtmdz(iz))**2.
1140       coeff(iz)= max(0.55,3./2.*(1.+A_wt(iz)))
1141       a_tpe(iz)=-g/th0(iz)*c_theta*te(iz)/d(iz)*dtmdz(iz)-coeff(iz)*d(iz)/te(iz)
1142       b_tpe(iz)=c_mu*te(iz)**2./d(iz)*dtmdz(iz)**2./Pra(iz)
1143      enddo
1144   !     tpe(kts)=(1./coeff(kts)*c_mu*dtmdz(kts)**2.*te(kts)**3./Pra(kts))/d(kts)**2.
1145
1146 return
1147 end subroutine terms_tpe
1148
1149 !==============================================================================
1150 !==============================================================================
1151  subroutine  calc_countergradient(dz,g,kms,kme,kts,kte,th0,te,d,tpe,Pra,b_t)
1152  implicit none
1153 integer iz,kms,kme,kts,kte
1154 real dz(kms:kme),g,te(kms:kme),d(kms:kme),Pra(kms:kme)
1155 real phi(kms:kme),dphidz1,dphidz2,phiz1,dphidz(kms:kme)
1156 real b_t(kms:kme),th0(kms:kme),tpe(kms:kme)
1157        phi=0.
1158        do iz=kts,kte
1159         phi(iz)=g/th0(iz)*te(iz)*tpe(iz)*c_theta/d(iz)
1160       enddo
1161        phiz1=(phi(kts)*dz(kts)+phi(kts+1)*dz(kts+1))/(dz(kts)+dz(kts+1))
1162        dphidz(kts)=phiz1/dz(kts)
1163       do iz=kts+1,kte-1
1164        dphidz1=2.*(phi(iz)-phi(iz-1))/(dz(iz-1)+dz(iz))
1165        dphidz2=2.*(phi(iz+1)-phi(iz))/(dz(iz+1)+dz(iz))
1166        dphidz(iz)=0.5*(dphidz1+dphidz2)
1167      enddo
1168       dphidz(kte)=0.
1169       dphidz(kts)=0.
1170      do iz=kts,kte
1171       b_t(iz)=b_t(iz)-dphidz(iz)
1172      enddo
1173
1174
1175 return
1176 end subroutine calc_countergradient
1177 END MODULE module_bl_keps