phys/module_bl_boulac.F

   1 MODULE module_bl_boulac
   2
   3 !USE module_model_constants
   4
   5
   6 !------------------------------------------------------------------------
   7 !    Calculation of the tendency due to momentum, heat
   8 !    and moisture turbulent fluxes follwing the approach
   9 !    of Bougeault and Lacarrere, 1989 (MWR, 117, 1872-1890).
  10 !    The scheme computes a prognostic ecuation for TKE and derives
  11 !    dissipation and turbulent coefficients using length scales.
  12 !
  13 !    Subroutine written by Alberto Martilli, CIEMAT, Spain,
  14 !    e-mail:alberto_martilli@ciemat.es
  15 !    August 2006.
  16 !------------------------------------------------------------------------
  17 ! IN THIS VERSION TKE IS NOT ADVECTED!!!!
  18 ! TO BE CHANGED IN THE FUTURE
  19 !
  20 ! -----------------------------------------------------------------------
  21 !  Constant used in the module
  22 !  ck_b=constant used in the compuation of diffusion coefficients
  23 !  ceps_b=constant used inthe computation of dissipation
  24 !  temin= minimum value allowed for TKE
  25 !  vk=von karman constant
  26 ! -----------------------------------------------------------------------
  27
  28       real ck_b,ceps_b,vk,temin    ! constant for Bougeault and Lacarrere
  29       parameter(ceps_b=1/1.4,ck_b=0.4,temin=0.0001,vk=0.4) ! impose minimum values for tke similar to those of MYJ
  30 ! -----------------------------------------------------------------------
  31
  32
  33    CONTAINS
  34
  35       subroutine boulac(frc_urb2d,idiff,flag_bep,dz8w,dt,u_phy,v_phy   &
  36                       ,th_phy,rho,qv_curr,qc_curr,hfx                                  &
  37                       ,qfx,ustar,cp,g                                          &
  38                       ,rublten,rvblten,rthblten                                &
  39                       ,rqvblten,rqcblten                        &
  40                       ,tke,dlk,wu,wv,wt,wq,exch_h,exch_m,pblh        &
  41                       ,a_u_bep,a_v_bep,a_t_bep,a_q_bep          &
  42                       ,a_e_bep,b_u_bep,b_v_bep                  &
  43                       ,b_t_bep,b_q_bep,b_e_bep,dlg_bep          &
  44                       ,dl_u_bep,sf_bep,vl_bep                   &
  45                       ,ids,ide, jds,jde, kds,kde                &
  46                       ,ims,ime, jms,jme, kms,kme                &
  47                       ,its,ite, jts,jte, kts,kte)
  48
  49       implicit none
  50
  51
  52
  53 !-----------------------------------------------------------------------
  54 !     Input
  55 !------------------------------------------------------------------------
  56    INTEGER::                        ids,ide, jds,jde, kds,kde,  &
  57                                     ims,ime, jms,jme, kms,kme,  &
  58                                     its,ite, jts,jte, kts,kte
  59
  60    integer, INTENT(IN) :: idiff
  61    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   DZ8W      !vertical resolution
  62    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   qv_curr   !moisture
  63    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   qc_curr   !liquid water
  64    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   RHO       !air density
  65    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   TH_PHY    !potential temperature
  66    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   U_PHY     !x-component of wind
  67    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   V_PHY     !y-component of wind
  68    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   ustar              !friction velocity
  69    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   hfx                !sensible heat flux (W/m2) at surface
  70    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   qfx                !moisture flux at surface
  71    real,  INTENT(IN   )    :: g,cp                                              !gravity and Cp
  72    REAL, INTENT(IN )::   DT                                                      ! Time step
  73
  74    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   FRC_URB2D          !fraction cover urban
  75    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)    ::   PBLH          !PBL height
  76 !
  77 ! variable added for urban
  78    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_u_bep        ! Implicit component for the momemtum in X-direction
  79    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_v_bep        ! Implicit component for the momemtum in Y-direction
  80    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_t_bep        ! Implicit component for the Pot. Temp.
  81    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_q_bep        ! Implicit component for Moisture
  82    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_e_bep        ! Implicit component for the TKE
  83    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_u_bep        ! Explicit component for the momemtum in X-direction
  84    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_v_bep        ! Explicit component for the momemtum in Y-direction
  85    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_t_bep        ! Explicit component for the Pot. Temp.
  86    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_q_bep        ! Explicit component for Moisture
  87    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_e_bep        ! Explicit component for the TKE
  88
  89    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT)    ::dlg_bep        ! Height above ground (L_ground in formula (24) of the BLM paper).
  90    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::dl_u_bep        ! Length scale (lb in formula (22) ofthe BLM paper).
  91 ! urban surface and volumes
  92    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::sf_bep           ! surface of the urban grid cells
  93    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::vl_bep             ! volume of the urban grid cells
  94    LOGICAL, INTENT(IN) :: flag_bep                                             !flag for BEP
  95
  96 !
  97 !-----------------------------------------------------------------------
  98 !     Local, carried on from one timestep to the other
  99 !------------------------------------------------------------------------
 100 !      real, save, allocatable, dimension (:,:,:)::TKE  ! Turbulent kinetic energy
 101       real,  dimension (ims:ime, kms:kme, jms:jme)  ::th_0 ! reference state for potential temperature
 102
 103 !------------------------------------------------------------------------
 104 !     Output
 105 !------------------------------------------------------------------------
 106         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  exch_h ! exchange coefficient for heat
 107         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  exch_m ! exchange coefficient for momentum
 108         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT   )    ::  tke  ! Turbulence Kinetic Energy
 109         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wu  ! Turbulent flux of momentum (x)
 110         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wv  ! Turbulent flux of momentum (y)
 111         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wt  ! Turbulent flux of temperature
 112         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wq  ! Turbulent flux of water vapor
 113         real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  dlk  ! Turbulent flux of water vapor
 114 ! only if idiff not equal 1:
 115         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RUBLTEN  !tendency for U_phy
 116         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RVBLTEN  !tendency for V_phy
 117         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RTHBLTEN !tendency for TH_phy
 118         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RQVBLTEN !tendency for QV_curr
 119         REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RQCBLTEN !tendency for QV_curr
 120
 121 !--------------------------------------------------------------
 122 ! Local
 123 !--------------------------------------------------------------
 124 ! 1D array used for the input and output of the routine boulac1D
 125
 126       real z1D(kms:kme)               ! vertical coordinates (faces of the grid)
 127       real dz1D(kms:kme)              ! vertical resolution
 128       real u1D(kms:kme)                 ! wind speed in the x directions
 129       real v1D(kms:kme)                 ! wind speed in the y directions
 130       real th1D(kms:kme)                ! potential temperature
 131       real q1D(kms:kme)                 ! moisture
 132       real qc1D(kms:kme)                 ! liquid water
 133       real rho1D(kms:kme)               ! air density
 134       real rhoz1D(kms:kme)            ! air density at the faces
 135       real tke1D(kms:kme)               ! air pressure
 136       real th01D(kms:kme)               ! reference potential temperature
 137       real dlk1D(kms:kme)               ! dlk
 138       real dls1D(kms:kme)               ! dls
 139       real exch1D(kms:kme)            ! exch
 140       real sf1D(kms:kme)              ! surface of the grid cells
 141       real vl1D(kms:kme)                ! volume of the  grid cells
 142       real a_u1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the X-direction
 143       real a_v1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the Y-direction
 144       real a_t1D(kms:kme)               ! Implicit component of the heat sources or sinks
 145       real a_q1D(kms:kme)               ! Implicit component of the moisture sources or sinks
 146       real a_qc1D(kms:kme)               ! Implicit component of the liquid water sources or sinks
 147       real a_e1D(kms:kme)               ! Implicit component of the TKE sources or sinks
 148       real b_u1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the X-direction
 149       real b_v1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the Y-direction
 150       real b_t1D(kms:kme)               ! Explicit component of the heat sources or sinks
 151       real b_q1D(kms:kme)               ! Explicit component of the moisture sources or sinks
 152       real b_qc1D(kms:kme)               ! Explicit component of the liquid water sources or sinks
 153       real b_e1D(kms:kme)               ! Explicit component of the TKE sources or sinks
 154       real dlg1D(kms:kme)               ! Height above ground (L_ground in formula (24) of the BLM paper).
 155       real dl_u1D(kms:kme)              ! Length scale (lb in formula (22) ofthe BLM paper)
 156       real sh1D(kms:kme)              ! shear
 157       real bu1D(kms:kme)              ! buoyancy
 158       real wu1D(kms:kme)              ! turbulent flux of momentum (x component)
 159       real wv1D(kms:kme)              ! turbulent flux of momentum (y component)
 160       real wt1D(kms:kme)              ! turbulent flux of temperature
 161       real wq1D(kms:kme)              ! turbulent flux of water vapor
 162       real wqc1D(kms:kme)              ! turbulent flux of liquid water
 163       real gamma1D(kms:kme)              ! non local term
 164       real t2_1D(kms:kme)              ! temperature variance
 165       real w2_1D(kms:kme)              ! vertical velocity variance
 166 ! local added only for diagnostic output
 167       real a_e(ims:ime,kms:kme,jms:jme) ! implicit term in TKE
 168       real b_e(ims:ime,kms:kme,jms:jme) ! explicit term in TKE
 169       real bu(ims:ime,kms:kme,jms:jme) ! buoyancy term in TKE
 170       real sh(ims:ime,kms:kme,jms:jme) ! shear term in TKE
 171       real wrk(ims:ime) ! working array
 172       integer ix,iy,iz,id,iz_u,iw_u,ig,ir_u,ix1,iy1,igamma
 173       real ufrac_int                                              ! urban fraction
 174       real vect,time_tke,hour,zzz
 175       real ustarf,wstar,wts,t2,w2,tstar_w,zzi
 176       real summ1,summ2,summ3
 177       save time_tke,hour
 178 !
 179 !
 180
 181 !here I fix the value of the reference state equal to the value of the potnetial temperature
 182 ! the only use of this variable in the code is to compute the paramter BETA = g/th0
 183 ! I fix it to 300K.
 184
 185         do ix=its,ite
 186         do iy=jts,jte
 187         do iz=kts,kte
 188 !         th_0(ix,iz,iy)=th_phy(ix,iz,iy)
 189          th_0(ix,iz,iy)=300.
 190         enddo
 191         enddo
 192         enddo
 193 ! initialization
 194        z1D=0.
 195        dz1D=0.
 196        u1D =0.
 197        v1D =0.
 198        th1D=0.
 199        q1D=0.
 200        rho1D=0.
 201        rhoz1D=0.
 202        tke1D =0.
 203        th01D =0.
 204        dlk1D =0.
 205        dls1D =0.
 206        exch1D=0.
 207        sf1D  =1.
 208        vl1D  =1.
 209        a_u1D =0.
 210        a_v1D =0.
 211        a_t1D =0.
 212        a_q1D =0.
 213        a_qc1D =0.
 214        a_e1D =0.
 215        b_u1D =0.
 216        b_v1D =0.
 217        b_t1D =0.
 218        b_q1D =0.
 219        b_qc1D =0.
 220        b_e1D =0.
 221        dlg1D =0.
 222        dl_u1D=0.
 223        sh1D  =0.
 224        bu1D  =0.
 225        wu1D  =0.
 226        wv1D  =0.
 227        wt1D =0.
 228        wq1D =0.
 229
 230 ! flag to choose the method for the calcaulation of the gamma non local term:
 231 ! igamma=0 - no term
 232 ! igamma=1 Troen and Mahrt
 233 ! igamma=2 Deardroff and Therry-Lacarrere
 234 ! igamma=3 Holstag and Moeng
 235
 236       igamma=1
 237 ! loop over the columns.
 238 ! put variables in 1D temporary arrays
 239 !
 240
 241        do ix=its,ite
 242        do iy=jts,jte
 243          z1d(kts)=0.
 244          do iz= kts,kte
 245           u1D(iz)=u_phy(ix,iz,iy)
 246           v1D(iz)=v_phy(ix,iz,iy)
 247           th1D(iz)=th_phy(ix,iz,iy)
 248           q1D(iz)=qv_curr(ix,iz,iy)
 249           qc1D(iz)=qc_curr(ix,iz,iy)
 250           tke1D(iz)=tke(ix,iz,iy)
 251           rho1D(iz)=rho(ix,iz,iy)
 252           th01D(iz)=th_0(ix,iz,iy)
 253           dz1D(iz)=dz8w(ix,iz,iy)
 254           z1D(iz+1)=z1D(iz)+dz1D(iz)
 255          enddo
 256         rhoz1D(kts)=rho1D(kts)
 257         do iz=kts+1,kte
 258          rhoz1D(iz)=(rho1D(iz)*dz1D(iz-1)+rho1D(iz-1)*dz1D(iz))/(dz1D(iz-1)+dz1D(iz))
 259         enddo
 260         rhoz1D(kte+1)=rho1D(kte)
 261         if(flag_bep)then
 262          do iz=kts,kte
 263           a_e1D(iz)=a_e_bep(ix,iz,iy)
 264           b_e1D(iz)=b_e_bep(ix,iz,iy)
 265           dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.*(1.-frc_urb2d(ix,iy))+dlg_bep(ix,iz,iy)*frc_urb2d(ix,iy)
 266           dl_u1D(iz)=dl_u_bep(ix,iz,iy)
 267           if((1.-frc_urb2d(ix,iy)).lt.1.)dl_u1D(iz)=dl_u1D(iz)/frc_urb2d(ix,iy)
 268           vl1D(iz)=vl_bep(ix,iz,iy)
 269           sf1D(iz)=sf_bep(ix,iz,iy)
 270          enddo
 271          ufrac_int=frc_urb2d(ix,iy)
 272          sf1D(kte+1)=sf_bep(ix,1,iy)
 273         else
 274          do iz=kts,kte
 275           a_e1D(iz)=0.
 276           b_e1D(iz)=0.
 277           dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.
 278           dl_u1D(iz)=0.
 279           vl1D(iz)=1.
 280           sf1D(iz)=1.
 281          enddo
 282          ufrac_int=0.
 283          sf1D(kte+1)=1.
 284         endif
 285
 286 ! compute the pbl_height
 287         call pbl_height(kms,kme,kts,kte,dz1d,z1d,th1D,q1D,pblh(ix,iy))
 288
 289 ! compute the values of wstar
 290         wts=max(0.,hfx(ix,iy)/rho1D(1)/cp)
 291         wstar=(g*wts*pblh(ix,iy)/th01D(1))**(1./3.)
 292         if (wts .ne. 0.0) then
 293           tstar_w=wts/wstar
 294         else
 295           tstar_w=0.0
 296         endif
 297         t2_1D=0.
 298         w2_1D=0.
 299         gamma1D=0.
 300 ! compute the variances
 301          do iz=kts+1,kte
 302            zzi=z1D(iz)/pblh(ix,iy)
 303            t2_1D(iz)=1.8*(zzi**(-2./3.))*(tstar_w**2.)
 304            w2_1D(iz)=1.8*(zzi**(2./3.))*((1.-0.8*zzi)**2.)*(wstar**2.)
 305         enddo
 306
 307 ! compute gamma
 308
 309
 310        if(igamma.eq.1)then
 311 ! (Troen and Mahrt)
 312         do iz=kts+1,kte
 313          if(z1D(iz).le.1.0*pblh(ix,iy).and.wts.gt.0.)then
 314           gamma1D(iz)=10.*wts/wstar/pblh(ix,iy)
 315          else
 316           gamma1D(iz)=0.
 317          endif
 318         enddo
 319
 320        elseif(igamma.eq.2)then
 321 ! Deardorff, and Therry -Lacarrere
 322          do iz=kts+1,kte
 323           if(wts.gt.0)then
 324            if(z1D(iz).le.(1.0*pblh(ix,iy)).and.z1D(iz).gt.(0.1*pblh(ix,iy)))then
 325             gamma1D(iz)=g/th01D(iz)*t2_1D(iz)/w2_1D(iz)
 326            else
 327             gamma1D(iz)=0.
 328            endif
 329           endif
 330          enddo
 331
 332        elseif(igamma.eq.3)then! (Holtslag and Moeng)
 333          do iz=kts+1,kte
 334           if(z1D(iz).le.(1.0*pblh(ix,iy)).and.wts.gt.0)then
 335            gamma1D(iz)=2.*wstar*wts/w2_1D(iz)/pblh(ix,iy)
 336           else
 337            gamma1D(iz)=0.
 338           endif
 339          enddo
 340
 341
 342        endif
 343
 344
 345          call boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz1d,z1D,dt,u1D,v1D,th1D,rho1D,rhoz1D,q1D,th01D,&
 346                    tke1D,ustar(ix,iy),hfx(ix,iy),qfx(ix,iy),cp,g,        &
 347                    a_e1D,b_e1D,                              &
 348                    dlg1D,dl_u1D,sf1D,vl1D,dlk1D,dls1D,exch1D,sh1D,bu1D,gamma1D)
 349
 350
 351
 352 ! store turbulent exchange coefficients, TKE, and other variables
 353
 354          do iz= kts,kte
 355           a_e(ix,iz,iy)=a_e1D(iz)
 356           b_e(ix,iz,iy)=b_e1D(iz)
 357           if(flag_bep)then
 358           dlg_bep(ix,iz,iy)=dlg1D(iz)
 359           endif
 360           tke(ix,iz,iy)=tke1D(iz)
 361           dlk(ix,iz,iy)=dlk1D(iz)
 362           sh(ix,iz,iy)=sh1D(iz)
 363           bu(ix,iz,iy)=bu1D(iz)
 364           exch_h(ix,iz,iy)=exch1D(iz)
 365           exch_m(ix,iz,iy)=exch1D(iz)
 366          enddo
 367
 368          if(idiff.ne.1)then
 369
 370 ! estimate the tendencies
 371
 372         if(flag_bep)then
 373          do iz=kts,kte
 374           a_t1D(iz)=a_t_bep(ix,iz,iy)
 375           b_t1D(iz)=b_t_bep(ix,iz,iy)
 376           a_u1D(iz)=a_u_bep(ix,iz,iy)
 377           b_u1D(iz)=b_u_bep(ix,iz,iy)
 378           a_v1D(iz)=a_v_bep(ix,iz,iy)
 379           b_v1D(iz)=b_v_bep(ix,iz,iy)
 380           a_q1D(iz)=a_q_bep(ix,iz,iy)
 381           b_q1D(iz)=b_q_bep(ix,iz,iy)
 382          enddo
 383         else
 384          do iz=kts,kte
 385           a_t1D(iz)=0.
 386           b_t1D(iz)=0.
 387           a_u1D(iz)=0.
 388           b_u1D(iz)=0.
 389           a_v1D(iz)=0.
 390           b_v1D(iz)=0.
 391           a_q1D(iz)=0.
 392           b_q1D(iz)=0.
 393          enddo
 394           b_t1D(1)=hfx(ix,iy)/dz1D(1)/rho1D(1)/cp
 395           b_q1D(1)=qfx(ix,iy)/dz1D(1)/rho1D(1)
 396           a_u1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5))
 397           a_v1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5))
 398         endif
 399
 400
 401 !
 402
 403 ! compute the value of the extra term that will be added to b_t1D
 404         do iz=kts+1,kte
 405          if(z1D(iz).le.1.0*pblh(ix,iy).and.wts.gt.0.)then
 406           b_t1D(iz)=b_t1D(iz)-(exch1D(iz+1)*gamma1D(iz+1)-exch1D(iz)*gamma1D(iz))/dz1D(iz)
 407          endif
 408         enddo
 409
 410
 411 !
 412
 413
 414
 415 ! solve diffusion equation for momentum x component
 416        call diff(kms,kme,kts,kte,1,1,dt,u1D,rho1D,rhoz1D,exch1D,a_u1D,b_u1D,sf1D,vl1D,dz1D,wu1D)
 417
 418 ! solve diffusion equation for momentum y component
 419        call diff(kms,kme,kts,kte,1,1,dt,v1D,rho1D,rhoz1D,exch1D,a_v1D,b_v1D,sf1D,vl1D,dz1D,wv1D)
 420
 421 ! solve diffusion equation for potential temperature
 422        call diff(kms,kme,kts,kte,1,1,dt,th1D,rho1D,rhoz1D,exch1D,a_t1D,b_t1D,sf1D,vl1D,dz1D,wt1D)
 423
 424 ! solve diffusion equation for water vapor mixing ratio
 425        call diff(kms,kme,kts,kte,1,1,dt,q1D,rho1D,rhoz1D,exch1D,a_q1D,b_q1D,sf1D,vl1D,dz1D,wq1D)
 426
 427 ! solve diffusion equation for liquid water mixing ratio
 428        call diff(kms,kme,kts,kte,1,1,dt,qc1D,rho1D,rhoz1D,exch1D,a_qc1D,b_qc1D,sf1D,vl1D,dz1D,wqc1D)
 429
 430 ! compute the tendencies
 431
 432          do iz= kts,kte
 433           rthblten(ix,iz,iy)=(th1D(iz)-th_phy(ix,iz,iy))/dt
 434           rqvblten(ix,iz,iy)=(q1D(iz)-qv_curr(ix,iz,iy))/dt
 435           rqcblten(ix,iz,iy)=(qc1D(iz)-qc_curr(ix,iz,iy))/dt
 436           rublten(ix,iz,iy)=(u1D(iz)-u_phy(ix,iz,iy))/dt
 437           rvblten(ix,iz,iy)=(v1D(iz)-v_phy(ix,iz,iy))/dt
 438           wu(ix,iz,iy)=wu1D(iz)
 439           wv(ix,iz,iy)=wv1D(iz)
 440           wt(ix,iz,iy)=wt1D(iz)
 441           wq(ix,iz,iy)=wq1D(iz)
 442         enddo
 443       endif
 444
 445       enddo  ! iy
 446       enddo  ! ix
 447
 448
 449       return
 450       end subroutine boulac
 451
 452
 453 ! ===6=8===============================================================72
 454
 455 ! ===6=8===============================================================72
 456
 457       subroutine boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz,z,dt,u,v,th,rho,rhoz,qa,th0,te,    &
 458                    ustar,hfx,qfx,cp,g,                               &
 459                    a_e,b_e,                        &
 460                    dlg,dl_u,sf,vl,dlk,dls,exch,sh,bu,gamma)
 461
 462 ! ----------------------------------------------------------------------
 463 ! 1D resolution of TKE following Bougeault and Lacarrere
 464 ! ----------------------------------------------------------------------
 465
 466       implicit none
 467
 468       integer iz,ix,iy
 469
 470 ! ----------------------------------------------------------------------
 471 ! INPUT:
 472 ! ----------------------------------------------------------------------
 473
 474
 475       integer kms,kme,kts,kte
 476       real z(kms:kme)               ! Altitude above the ground of the cell interfaces.
 477       real dz(kms:kme)                ! vertical resolution
 478       real u(kms:kme)                ! Wind speed in the x direction
 479       real v(kms:kme)                ! Wind speed in the y direction
 480       real th(kms:kme)                ! Potential temperature
 481       real rho(kms:kme)                ! Air density
 482       real g                     ! gravity
 483       real cp                    !
 484       real te(kms:kme)                ! turbulent kinetic energy
 485       real qa(kms:kme)                ! air humidity
 486       real th0(kms:kme)               ! Reference potential temperature
 487       real dt                    ! Time step
 488       real ustar                 ! ustar
 489       real hfx                   ! sensbile heat flux
 490       real qfx                   ! kinematic latent heat flux
 491       real sf(kms:kme)              ! surface of the urban grid cells
 492       real vl(kms:kme)                ! volume of the urban grid cells
 493       real a_e(kms:kme)               ! Implicit component of the TKE sources or sinks
 494       real b_e(kms:kme)               ! Explicit component of the TKE sources or sinks
 495       real dlg(kms:kme)               ! Height above ground (L_ground in formula (24) of the BLM paper).
 496       real dl_u(kms:kme)              ! Length scale (lb in formula (22) ofthe BLM paper)
 497       real ufrac_int             ! urban fraction
 498 ! local variables not needed in principle, but that can be used to estimate the budget and turbulent fluxes
 499
 500       real we(kms:kme),dwe(kms:kme)
 501
 502 ! local variables
 503       real sh(kms:kme)    ! shear term in TKE eqn.
 504       real bu(kms:kme)    ! buoyancy term in TKE eqn.
 505       real gamma(kms:kme)    ! gamma term
 506       real td(kms:kme)    ! dissipation term in TKE eqn.
 507       real exch(kms:kme) ! turbulent diffusion coefficients (defined at the faces)
 508       real dls(kms:kme)   ! dissipation length scale
 509       real dlk(kms:kme)   ! length scale used to estimate exch
 510       real dlu(kms:kme)   ! l_up
 511       real dld(kms:kme)   ! l_down
 512       real rhoz(kms:kme) !air density at the faces of the cell
 513       real tstar     ! derived from hfx and ustar
 514       real beta
 515       real summ1,summ2,summ3,summ4
 516 ! interpolate air density at the faces
 517
 518
 519
 520 ! estimation of tstar
 521
 522         tstar=-hfx/rho(1)/cp/ustar
 523
 524 ! first compute values of dlu and dld (length scales up and down).
 525
 526        call dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
 527
 528 !then average them to obtain dls and dlk (length scales for dissipation and eddy coefficients)
 529
 530        call length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
 531
 532 ! compute the turbulent diffusion coefficients exch
 533
 534        call cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
 535
 536 ! compute source and sink terms in the TKE equation (shear, buoyancy and dissipation)
 537
 538        call tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch,dls,td,sh,bu,gamma,b_e,a_e,sf,ufrac_int)
 539
 540 ! solve for tke
 541
 542        call diff(kms,kme,kts,kte,1,1,dt,te,rho,rhoz,exch,a_e,b_e,sf,vl,dz,we)
 543
 544 ! avoid negative values for tke
 545
 546        do iz=kts,kte
 547         if(te(iz).lt.temin) te(iz)=temin
 548        enddo
 549
 550        return
 551        end subroutine boulac1d
 552 !
 553 ! ===6=8===============================================================72
 554
 555 ! ===6=8===============================================================72
 556          subroutine dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
 557 ! compute the length scales up and down
 558          implicit none
 559          integer kms,kme,kts,kte,iz,izz,ix,iy
 560          real dzt,zup,beta,zup_inf,bbb,tl,zdo,zdo_sup,zzz,g
 561          real te(kms:kme),dlu(kms:kme),dld(kms:kme),dz(kms:kme)
 562          real z(kms:kme),th(kms:kme),th0(kms:kme)
 563
 564          do iz=kts,kte
 565           zup=0.
 566           dlu(iz)=z(kte+1)-z(iz)-dz(iz)/2.
 567           zzz=0.
 568           zup_inf=0.
 569           beta=g/th0(iz)      !Buoyancy coefficient
 570           do izz=iz,kte-1
 571            dzt=(dz(izz+1)+dz(izz))/2.
 572            zup=zup-beta*th(iz)*dzt
 573            zup=zup+beta*(th(izz+1)+th(izz))*dzt/2.
 574            zzz=zzz+dzt
 575            if(te(iz).lt.zup.and.te(iz).ge.zup_inf)then
 576             bbb=(th(izz+1)-th(izz))/dzt
 577             if(bbb.ne.0)then
 578              tl=(-beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zup_inf))))/bbb/beta
 579             else
 580              if(th(izz).ne.th(iz))then
 581               tl=(te(iz)-zup_inf)/(beta*(th(izz)-th(iz)))
 582              else
 583               tl=0.
 584              endif
 585             endif
 586             dlu(iz)=zzz-dzt+tl
 587            endif
 588            zup_inf=zup
 589           enddo
 590
 591           zdo=0.
 592           zdo_sup=0.
 593           dld(iz)=z(iz)+dz(iz)/2.
 594           zzz=0.
 595           do izz=iz,kts+1,-1
 596            dzt=(dz(izz-1)+dz(izz))/2.
 597            zdo=zdo+beta*th(iz)*dzt
 598            zdo=zdo-beta*(th(izz-1)+th(izz))*dzt/2.
 599            zzz=zzz+dzt
 600            if(te(iz).lt.zdo.and.te(iz).ge.zdo_sup)then
 601             bbb=(th(izz)-th(izz-1))/dzt
 602             if(bbb.ne.0.)then
 603              tl=(beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zdo_sup))))/bbb/beta
 604             else
 605              if(th(izz).ne.th(iz))then
 606               tl=(te(iz)-zdo_sup)/(beta*(th(izz)-th(iz)))
 607              else
 608               tl=0.
 609              endif
 610             endif
 611
 612             dld(iz)=zzz-dzt+tl
 613            endif
 614            zdo_sup=zdo
 615           enddo
 616          enddo
 617
 618
 619          end subroutine dissip_bougeault
 620 !
 621 ! ===6=8===============================================================72
 622 ! ===6=8===============================================================72
 623          subroutine length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
 624 ! compute the length scales for dissipation and turbulent coefficients
 625          implicit none
 626          integer kms,kme,kts,kte,iz,ix,iy
 627          real dlu(kms:kme),dld(kms:kme),dl_u(kms:kme)
 628          real dls(kms:kme),dlk(kms:kme),dlg(kms:kme)
 629
 630          do iz=kts,kte
 631           dld(iz)=min(dld(iz),dlg(iz))
 632           dls(iz)=sqrt(dlu(iz)*dld(iz))
 633           dlk(iz)=min(dlu(iz),dld(iz))
 634
 635          if(dl_u(iz).gt.0.)then
 636            dls(iz)=1./(1./dls(iz)+1./dl_u(iz))
 637            dlk(iz)=1./(1./dlk(iz)+1./dl_u(iz))
 638           endif
 639          enddo
 640
 641          return
 642          end subroutine length_bougeault
 643 !
 644
 645 ! ===6=8===============================================================72
 646 ! ===6=8===============================================================72
 647
 648        subroutine cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
 649 ! compute turbulent coefficients
 650        implicit none
 651        integer iz,kms,kme,kts,kte,ix,iy
 652        real te_m,dlk_m
 653        real te(kms:kme),exch(kms:kme)
 654        real dz(kms:kme),z(kms:kme)
 655        real dlk(kms:kme)
 656        real fact
 657
 658        exch(kts)=0.
 659
 660 !       do iz=2,nz-1
 661        do iz=kts+1,kte
 662         te_m=(te(iz-1)*dz(iz)+te(iz)*dz(iz-1))/(dz(iz)+dz(iz-1))
 663         dlk_m=(dlk(iz-1)*dz(iz)+dlk(iz)*dz(iz-1))/(dz(iz)+dz(iz-1))
 664         exch(iz)=ck_b*dlk_m*sqrt(te_m)
 665 !        exch(iz)=max(exch(iz),0.0001)
 666         exch(iz)=max(exch(iz),0.1)
 667        enddo
 668
 669        exch(kte+1)=0.1
 670
 671        return
 672        end subroutine cdtur_bougeault
 673
 674
 675 ! ===6=8===============================================================72
 676 ! ===6=8===============================================================72
 677
 678        subroutine diff(kms,kme,kts,kte,iz1,izf,dt,co,rho,rhoz,cd,aa,bb,sf,vl,dz,fc)
 679
 680
 681 !------------------------------------------------------------------------
 682 !           Calculation of the diffusion in 1D
 683 !------------------------------------------------------------------------
 684 !  - Input:
 685 !       nz    : number of points
 686 !       iz1   : first calculated point
 687 !       co    : concentration of the variable of interest
 688 !       dz    : vertical levels
 689 !       cd    : diffusion coefficients
 690 !       dtext : external time step
 691 !       rho    : density of the air at the center
 692 !       rhoz   : density of the air at the face
 693 !       itest : if itest eq 1 then update co, else store in a flux array
 694 !  - Output:
 695 !       co :concentration of the variable of interest
 696
 697 !  - Internal:
 698 !       cddz  : constant terms in the equations
 699 !       dt    : diffusion time step
 700 !       nt    : number of the diffusion time steps
 701 !       cstab : ratio of the stability condition for the time step
 702 !---------------------------------------------------------------------
 703
 704         implicit none
 705         integer iz,iz1,izf
 706         integer kms,kme,kts,kte
 707         real dt,dzv
 708         real co(kms:kme),cd(kms:kme),dz(kms:kme)
 709         real rho(kms:kme),rhoz(kms:kme)
 710         real cddz(kms:kme+1),fc(kms:kme),df(kms:kme)
 711         real a(kms:kme,3),c(kms:kme)
 712         real sf(kms:kme),vl(kms:kme)
 713         real aa(kms:kme),bb(kms:kme)
 714
 715
 716 ! Compute cddz=2*cd/dz
 717
 718         cddz(kts)=sf(kts)*rhoz(kts)*cd(kts)/dz(kts)
 719         do iz=kts+1,kte
 720          cddz(iz)=2.*sf(iz)*rhoz(iz)*cd(iz)/(dz(iz)+dz(iz-1))
 721         enddo
 722         cddz(kte+1)=sf(kte+1)*rhoz(kte+1)*cd(kte+1)/dz(kte)
 723
 724          do iz=kts,iz1-1
 725           a(iz,1)=0.
 726           a(iz,2)=1.
 727           a(iz,3)=0.
 728           c(iz)=co(iz)
 729          enddo
 730
 731           do iz=iz1,kte-izf
 732            dzv=vl(iz)*dz(iz)
 733            a(iz,1)=-cddz(iz)*dt/dzv/rho(iz)
 734            a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dzv/rho(iz)-aa(iz)*dt
 735            a(iz,3)=-cddz(iz+1)*dt/dzv/rho(iz)
 736            c(iz)=co(iz)+bb(iz)*dt
 737           enddo
 738
 739           do iz=kte-(izf-1),kte
 740            a(iz,1)=0.
 741            a(iz,2)=1
 742            a(iz,3)=0.
 743            c(iz)=co(iz)
 744           enddo
 745
 746           call invert (kms,kme,kts,kte,a,c,co)
 747
 748           do iz=kts,iz1
 749            fc(iz)=0.
 750           enddo
 751
 752           do iz=iz1+1,kte
 753            fc(iz)=-(cddz(iz)*(co(iz)-co(iz-1)))/rho(iz)
 754           enddo
 755
 756 !          do iz=1,iz1
 757 !           df(iz)=0.
 758 !          enddo
 759 !
 760 !          do iz=iz1+1,nz-izf
 761 !           dzv=vl(iz)*dz(iz)
 762 !           df(iz)=+(co(iz-1)*cddz(iz)-co(iz)*(cddz(iz)+cddz(iz+1))+co(iz+1)*cddz(iz+1))/dzv/rho(iz)
 763 !          enddo
 764 !
 765 !          do iz=nz-izf,nz
 766 !           df(iz)=0.
 767 !          enddo
 768
 769        return
 770        end subroutine diff
 771
 772 ! ===6=8===============================================================72
 773 ! ===6=8===============================================================72
 774
 775        subroutine buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,gamma,ustar,tstar,ufrac_int)
 776 ! compute buoyancy term
 777        implicit none
 778        integer kms,kme,kts,kte,iz,ix,iy
 779        real dtdz1,dtdz2,cdm,dtmdz,g
 780        real th(kms:kme),exch(kms:kme),dz(kms:kme),bu(kms:kme),gamma(kms:kme)
 781        real th0(kms:kme),ustar,tstar,ufrac_int,gammam
 782
 783 !       bu(1)=-ustar*tstar*g/th0(1)*(1.-ufrac_int)
 784       bu(kts)=0.
 785
 786
 787        do iz=kts+1,kte-1
 788         dtdz1=2.*(th(iz)-th(iz-1))/(dz(iz-1)+dz(iz))
 789         dtdz2=2.*(th(iz+1)-th(iz))/(dz(iz+1)+dz(iz))
 790         dtmdz=0.5*(dtdz1+dtdz2)
 791         cdm=0.5*(exch(iz+1)+exch(iz))
 792         gammam=0.5*(gamma(iz+1)+gamma(iz))
 793         bu(iz)=-cdm*(dtmdz-gammam)*g/th0(iz)
 794        enddo
 795 !
 796
 797        bu(kte)=0.
 798
 799        return
 800        end subroutine buoy
 801
 802 ! ===6=8===============================================================72
 803 ! ===6=8===============================================================72
 804
 805        subroutine shear(ix,iy,g,kms,kme,kts,kte,u,v,cdua,dz,sh,ustar,tstar,th,ufrac_int)
 806 ! compute shear term
 807        implicit none
 808        integer kms,kme,kts,kte,iz,ix,iy
 809        real dudz1,dudz2,dvdz1,dvdz2,cdm,dumdz,ustar
 810        real tstar,th,al,phim,g
 811        real u(kms:kme),v(kms:kme),cdua(kms:kme),dz(kms:kme),sh(kms:kme)
 812        real u1,u2,v1,v2,ufrac_int
 813
 814 !       al=vk*g*tstar/(th*(ustar**2.))
 815 !       if(al.ge.0.)phim=1.+4.7*dz(1)/2.*al
 816 !       if(al.lt.0.)phim=(1.-15*dz(1)/2.*al)**(-0.25)
 817 !
 818 !       sh(1)=(ustar**3.)/vk/(dz(1)/2.)*(1.-ufrac_int)
 819        sh(kts)=0.
 820        do iz=kts+1,kte-1
 821         u2=(dz(iz+1)*u(iz)+dz(iz)*u(iz+1))/(dz(iz)+dz(iz+1))
 822         u1=(dz(iz)*u(iz-1)+dz(iz-1)*u(iz))/(dz(iz-1)+dz(iz))
 823         v2=(dz(iz+1)*v(iz)+dz(iz)*v(iz+1))/(dz(iz)+dz(iz+1))
 824         v1=(dz(iz)*v(iz-1)+dz(iz-1)*v(iz))/(dz(iz-1)+dz(iz))
 825         cdm=0.5*(cdua(iz)+cdua(iz+1))
 826         dumdz=((u2-u1)/dz(iz))**2.+((v2-v1)/dz(iz))**2.
 827         sh(iz)=cdm*dumdz
 828        enddo
 829
 830 !!!!!!!
 831        sh(kte)=0.
 832
 833        return
 834        end subroutine shear
 835
 836 ! ===6=8===============================================================72
 837 ! ===6=8===============================================================72
 838
 839        subroutine invert(kms,kme,kts,kte,a,c,x)
 840
 841 !ccccccccccccccccccccccccccccccc
 842 ! Aim: Inversion and resolution of a tridiagonal matrix
 843 !          A X = C
 844 ! Input:
 845 !  a(*,1) lower diagonal (Ai,i-1)
 846 !  a(*,2) principal diagonal (Ai,i)
 847 !  a(*,3) upper diagonal (Ai,i+1)
 848 !  c
 849 ! Output
 850 !  x     results
 851 !ccccccccccccccccccccccccccccccc
 852
 853        implicit none
 854        integer in
 855        integer kts,kte,kms,kme
 856        real a(kms:kme,3),c(kms:kme),x(kms:kme)
 857
 858         do in=kte-1,kts,-1
 859          c(in)=c(in)-a(in,3)*c(in+1)/a(in+1,2)
 860          a(in,2)=a(in,2)-a(in,3)*a(in+1,1)/a(in+1,2)
 861         enddo
 862
 863         do in=kts+1,kte
 864          c(in)=c(in)-a(in,1)*c(in-1)/a(in-1,2)
 865         enddo
 866
 867         do in=kts,kte
 868
 869          x(in)=c(in)/a(in,2)
 870
 871         enddo
 872
 873         return
 874         end subroutine invert
 875
 876 ! ===6=8===============================================================72
 877 ! ===6=8===============================================================72
 878
 879        subroutine tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch,   &
 880                          dls,td,sh,bu,gamma,b_e,a_e,sf,ufrac_int)
 881 ! in this routine the shear, buoyancy and part of the dissipation terms
 882 ! of the TKE equation are computed
 883
 884        implicit none
 885        integer kms,kme,kts,kte,iz,ix,iy
 886        real g,ustar,tstar,ufrac_int
 887        real z(kms:kme),dz(kms:kme),u(kms:kme),v(kms:kme),th(kms:kme),th0(kms:kme),te(kms:kme)
 888        real exch(kms:kme),dls(kms:kme),td(kms:kme),sh(kms:kme),bu(kms:kme),gamma(kms:kme)
 889        real a_e(kms:kme),b_e(kms:kme)
 890        real vl(kms:kme),sf(kms:kme)
 891        real te1,dl1
 892
 893        call shear(ix,iy,g,kms,kme,kts,kte,u,v,exch,dz,sh,ustar,tstar,th(kts),ufrac_int)
 894
 895        call buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,gamma,ustar,tstar,ufrac_int)
 896
 897        do iz=kts,kte
 898         te1=max(te(iz),temin)
 899         dl1=max(dls(iz),0.1)
 900         td(iz)=-ceps_b*sqrt(te1)/dl1
 901         sh(iz)=sh(iz)*sf(iz)
 902         bu(iz)=bu(iz)*sf(iz)
 903         a_e(iz)=a_e(iz)+td(iz)
 904         b_e(iz)=b_e(iz)+sh(iz)+bu(iz)
 905        enddo
 906
 907
 908        return
 909        end subroutine tke_bougeault
 910
 911 ! ===6=8===============================================================72
 912       SUBROUTINE BOULACINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN,          &
 913      &                      TKE_PBL,EXCH_H,RESTART,ALLOWED_TO_READ,     &
 914      &                      IDS,IDE,JDS,JDE,KDS,KDE,                    &
 915      &                      IMS,IME,JMS,JME,KMS,KME,                    &
 916      &                      ITS,ITE,JTS,JTE,KTS,KTE                 )
 917 !-----------------------------------------------------------------------
 918       IMPLICIT NONE
 919 !-----------------------------------------------------------------------
 920       LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
 921       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,                    &
 922      &                      IMS,IME,JMS,JME,KMS,KME,                    &
 923      &                      ITS,ITE,JTS,JTE,KTS,KTE
 924
 925       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) ::    EXCH_H, &
 926      &                                                         RUBLTEN, &
 927      &                                                         RVBLTEN, &
 928      &                                                        RTHBLTEN, &
 929      &                                                        RQVBLTEN, &
 930      &                                                        RQCBLTEN, &
 931      &                                                         TKE_PBL
 932       INTEGER :: I,J,K,ITF,JTF,KTF
 933 !-----------------------------------------------------------------------
 934 !-----------------------------------------------------------------------
 935
 936       JTF=MIN0(JTE,JDE-1)
 937       KTF=MIN0(KTE,KDE-1)
 938       ITF=MIN0(ITE,IDE-1)
 939
 940       IF(.NOT.RESTART)THEN
 941         DO J=JTS,JTF
 942         DO K=KTS,KTF
 943         DO I=ITS,ITF
 944           TKE_PBL(I,K,J)=0.0001
 945           RUBLTEN(I,K,J)=0.
 946           RVBLTEN(I,K,J)=0.
 947           RTHBLTEN(I,K,J)=0.
 948           RQVBLTEN(I,K,J)=0.
 949           RQCBLTEN(I,K,J)=0.
 950           EXCH_H(I,K,J)=0.
 951         ENDDO
 952         ENDDO
 953         ENDDO
 954       ENDIF
 955
 956       END SUBROUTINE BOULACINIT
 957 !######################################################################
 958        subroutine pbl_height(kms,kme,kts,kte,dz,z,th,q,pblh)
 959
 960 ! this routine computes the PBL height
 961 ! with an approach similar to MYNN
 962        implicit none
 963        integer kms,kme,kts,kte,iz
 964        real z(kms:kme),dz(kms:kme),th(kms:kme),q(kms:kme)
 965        real pblh
 966 !Local
 967        real thv(kms:kme),zc(kms:kme)
 968        real thsfc
 969
 970 ! compute the height of the center of the grid cells
 971       do iz=kts,kte
 972        zc(iz)=z(iz)+dz(iz)/2.
 973       enddo
 974
 975 ! compute the virtual potential temperature
 976
 977        do iz=kts,kte
 978         thv(iz)=th(iz)*(1.+0.61*q(iz))
 979        enddo
 980 ! now compute the PBL height
 981
 982        pblh=0.
 983        thsfc=thv(kts)+0.5
 984        do iz=kts+1,kte
 985         if(pblh.eq.0.and.thv(iz).gt.thsfc)then
 986          pblh=zc(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(zc(iz)-zc(iz-1))
 987 !         pblh=z(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(z(iz)-z(iz-1))
 988         endif
 989        enddo
 990
 991        return
 992        end subroutine pbl_height
 993
 994
 995 ! ===6=8===============================================================72
 996
 997
 998
 999 END MODULE module_bl_boulac
1000