phys/module_bl_mrf.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3 MODULE module_bl_mrf
   4
   5 CONTAINS
   6
   7 !-------------------------------------------------------------------
   8    SUBROUTINE MRF(U3D,V3D,TH3D,T3D,QV3D,QC3D,P3D,PI3D,             &
   9                   RUBLTEN,RVBLTEN,RTHBLTEN,                        &
  10                   RQVBLTEN,RQCBLTEN,                               &
  11                   CP,G,ROVCP,R,ROVG,                               &
  12                   dz8w,z,XLV,RV,PSFC,                              &
  13                   p1000mb,                                         &
  14                   ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH,                   &
  15                   XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR,                &
  16                   DT,DTMIN,KPBL2D,                                 &
  17                   SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
  18                   flag_qi,                                         &
  19                   ids,ide, jds,jde, kds,kde,                       &
  20                   ims,ime, jms,jme, kms,kme,                       &
  21                   its,ite, jts,jte, kts,kte,                       &
  22                ! Optional
  23                   QI3D,RQIBLTEN,                                   &
  24                   regime                                           )
  25 !-------------------------------------------------------------------
  26       IMPLICIT NONE
  27 !-------------------------------------------------------------------
  28 !-- U3D         3D u-velocity interpolated to theta points (m/s)
  29 !-- V3D         3D v-velocity interpolated to theta points (m/s)
  30 !-- TH3D        3D potential temperature (K)
  31 !-- T3D         temperature (K)
  32 !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
  33 !-- QC3D        3D cloud mixing ratio (Kg/Kg)
  34 !-- QI3D        3D ice mixing ratio (Kg/Kg)
  35 !-- P3D         3D pressure (Pa)
  36 !-- PI3D        3D exner function (dimensionless)
  37 !-- rr3D        3D dry air density (kg/m^3)
  38 !-- RUBLTEN     U tendency due to
  39 !               PBL parameterization (m/s^2)
  40 !-- RVBLTEN     V tendency due to
  41 !               PBL parameterization (m/s^2)
  42 !-- RTHBLTEN    Theta tendency due to
  43 !               PBL parameterization (K/s)
  44 !-- RQVBLTEN    Qv tendency due to
  45 !               PBL parameterization (kg/kg/s)
  46 !-- RQCBLTEN    Qc tendency due to
  47 !               PBL parameterization (kg/kg/s)
  48 !-- RQIBLTEN    Qi tendency due to
  49 !               PBL parameterization (kg/kg/s)
  50 !-- CP          heat capacity at constant pressure for dry air (J/kg/K)
  51 !-- G           acceleration due to gravity (m/s^2)
  52 !-- ROVCP       R/CP
  53 !-- R           gas constant for dry air (J/kg/K)
  54 !-- ROVG        R/G
  55 !-- dz8w        dz between full levels (m)
  56 !-- z           height above sea level (m)
  57 !-- XLV         latent heat of vaporization (J/kg)
  58 !-- RV          gas constant for water vapor (J/kg/K)
  59 !-- PSFC        pressure at the surface (Pa)
  60 !-- ZNT         roughness length (m)
  61 !-- UST         u* in similarity theory (m/s)
  62 !-- ZOL         z/L height over Monin-Obukhov length
  63 !-- HOL         PBL height over Monin-Obukhov length
  64 !-- PBL         PBL height (m)
  65 !-- REGIME      flag indicating PBL regime (stable, unstable, etc.)
  66 !-- PSIM        similarity stability function for momentum
  67 !-- PSIH        similarity stability function for heat
  68 !-- XLAND       land mask (1 for land, 2 for water)
  69 !-- HFX         upward heat flux at the surface (W/m^2)
  70 !-- QFX         upward moisture flux at the surface (kg/m^2/s)
  71 !-- TSK         surface temperature (K)
  72 !-- GZ1OZ0      log(z/z0) where z0 is roughness length
  73 !-- WSPD        wind speed at lowest model level (m/s)
  74 !-- BR          bulk Richardson number in surface layer
  75 !-- DT          time step (s)
  76 !-- DTMIN       time step (minute)
  77 !-- rvovrd      R_v divided by R_d (dimensionless)
  78 !-- SVP1        constant for saturation vapor pressure (kPa)
  79 !-- SVP2        constant for saturation vapor pressure (dimensionless)
  80 !-- SVP3        constant for saturation vapor pressure (K)
  81 !-- SVPT0       constant for saturation vapor pressure (K)
  82 !-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
  83 !-- EP2         constant for specific humidity calculation
  84 !-- KARMAN      Von Karman constant
  85 !-- EOMEG       angular velocity of earth's rotation (rad/s)
  86 !-- STBOLT      Stefan-Boltzmann constant (W/m^2/K^4)
  87 !-- ids         start index for i in domain
  88 !-- ide         end index for i in domain
  89 !-- jds         start index for j in domain
  90 !-- jde         end index for j in domain
  91 !-- kds         start index for k in domain
  92 !-- kde         end index for k in domain
  93 !-- ims         start index for i in memory
  94 !-- ime         end index for i in memory
  95 !-- jms         start index for j in memory
  96 !-- jme         end index for j in memory
  97 !-- kms         start index for k in memory
  98 !-- kme         end index for k in memory
  99 !-- its         start index for i in tile
 100 !-- ite         end index for i in tile
 101 !-- jts         start index for j in tile
 102 !-- jte         end index for j in tile
 103 !-- kts         start index for k in tile
 104 !-- kte         end index for k in tile
 105 !-------------------------------------------------------------------
 106
 107       INTEGER,  INTENT(IN   )   ::      ids,ide, jds,jde, kds,kde, &
 108                                         ims,ime, jms,jme, kms,kme, &
 109                                         its,ite, jts,jte, kts,kte
 110
 111 !
 112       REAL,     INTENT(IN   )   ::      P1000mb
 113       REAL,     INTENT(IN   )   ::      DT,DTMIN,CP,G,ROVCP,       &
 114                                         ROVG,R,XLV,RV
 115
 116       REAL,     INTENT(IN )     ::     SVP1,SVP2,SVP3,SVPT0
 117       REAL,     INTENT(IN )     ::     EP1,EP2,KARMAN,EOMEG,STBOLT
 118
 119       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 120                 INTENT(IN   )   ::                           QV3D, &
 121                                                              QC3D, &
 122                                                               P3D, &
 123                                                              PI3D, &
 124                                                              TH3D, &
 125                                                               T3D, &
 126                                                              dz8w, &
 127                                                                 z
 128 !
 129       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 130                 INTENT(INOUT)   ::                        RUBLTEN, &
 131                                                           RVBLTEN, &
 132                                                          RTHBLTEN, &
 133                                                          RQVBLTEN, &
 134                                                          RQCBLTEN
 135
 136       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 137                 INTENT(IN   )   ::                          XLAND, &
 138                                                               HFX, &
 139                                                               QFX
 140
 141       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 142                 INTENT(INOUT)   ::                            HOL, &
 143                                                               UST, &
 144                                                               PBL, &
 145                                                               ZNT
 146
 147       LOGICAL,  INTENT(IN)      ::                        FLAG_QI
 148 !
 149 !m The following 5 variables are changed to memory size from tile size--
 150 !
 151      REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  ::    &
 152                                                              PSIM, &
 153                                                              PSIH
 154
 155      REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)   ::   &
 156                                                              WSPD
 157
 158      REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(IN   )   ::   &
 159                                                            GZ1OZ0, &
 160                                                                BR
 161
 162       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 163                 INTENT(IN   )               ::               PSFC
 164
 165       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 166                 INTENT(IN   )   ::                            TSK
 167
 168       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 169                 INTENT(INOUT)   ::                            ZOL
 170
 171       INTEGER,  DIMENSION( ims:ime, jms:jme )                    , &
 172                 INTENT(OUT  )   ::                         KPBL2D
 173
 174       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 175                 INTENT(IN   )   ::                            U3D, &
 176                                                               V3D
 177 !
 178 ! Optional
 179 !
 180       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 181                 OPTIONAL                                         , &
 182                 INTENT(INOUT)   ::                         REGIME
 183
 184       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 185                 INTENT(INOUT)   ::                       RQIBLTEN
 186
 187       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 188                 OPTIONAL                                         , &
 189                 INTENT(IN   )   ::                           QI3D
 190
 191 ! LOCAL VARS
 192    REAL,       DIMENSION( its:ite, kts:kte )          ::   dz8w2d, &
 193                                                               z2d
 194
 195
 196    INTEGER ::  I,J,K,NK
 197
 198 !
 199    DO J=jts,jte
 200       DO k=kts,kte
 201       NK=kme-k
 202       DO i=its,ite
 203          dz8w2d(I,K) = dz8w(i,NK,j)
 204          z2d(I,K) = z(i,NK,j)
 205       ENDDO
 206       ENDDO
 207
 208
 209       CALL MRF2D(J,U3D(ims,kms,j),V3D(ims,kms,j),T3D(ims,kms,j),    &
 210                QV3D(ims,kms,j),QC3D(ims,kms,j),                     &
 211                P3D(ims,kms,j),RUBLTEN(ims,kms,j),RVBLTEN(ims,kms,j),&
 212                RTHBLTEN(ims,kms,j),RQVBLTEN(ims,kms,j),             &
 213                RQCBLTEN(ims,kms,j),                                 &
 214                p1000mb,                                             &
 215                CP,G,ROVCP,R,ROVG,                                   &
 216                dz8w2d,z2d,XLV,Rv,                                   &
 217                PSFC(ims,j),ZNT(ims,j),                              &
 218                UST(ims,j),ZOL(ims,j),                               &
 219                HOL(ims,j),PBL(ims,j),PSIM(ims,j),                   &
 220                PSIH(ims,j),XLAND(ims,j),HFX(ims,j),QFX(ims,j),      &
 221                TSK(ims,j),GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),      &
 222                DT,DTMIN,KPBL2D(ims,j),                              &
 223                SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,    &
 224                flag_qi,                                             &
 225                ids,ide, jds,jde, kds,kde,                           &
 226                ims,ime, jms,jme, kms,kme,                           &
 227                its,ite, jts,jte, kts,kte,                           &
 228             !optional
 229                QI2DTEN=RQIBLTEN(ims,kms,j),                         &
 230                REGIME=REGIME(ims,j),QI2D=QI3D(ims,kms,j)            )
 231
 232
 233       DO k=kts,kte
 234       DO i=its,ite
 235          RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/PI3D(I,K,J)
 236       ENDDO
 237       ENDDO
 238
 239     ENDDO
 240
 241    END SUBROUTINE MRF
 242
 243 !-------------------------------------------------------------------
 244    SUBROUTINE MRF2D(J,U2D,V2D,T2D,QV2D,QC2D,     P2D,              &
 245                   U2DTEN,V2DTEN,T2DTEN,                            &
 246                   QV2DTEN,QC2DTEN,                                 &
 247                   p1000mb,                                         &
 248                   CP,G,ROVCP,R,ROVG,                               &
 249                   dz8w2d,z2d,XLV,RV,PSFCPA,                        &
 250                   ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH,                   &
 251                   XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR,                &
 252                   DT,DTMIN,KPBL1D,                                 &
 253                   SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
 254                   flag_qi,                                         &
 255                   ids,ide, jds,jde, kds,kde,                       &
 256                   ims,ime, jms,jme, kms,kme,                       &
 257                   its,ite, jts,jte, kts,kte,                       &
 258                ! optional
 259                   regime, qi2d, QI2DTEN                            )
 260 !-------------------------------------------------------------------
 261       IMPLICIT NONE
 262 !-------------------------------------------------------------------
 263 !     BASED ON THE "COUNTERGRADIENT" TRANSPORT TERM OF TROEN
 264 !     AND MAHRT (1986) FOR THE UNSTABLE PBL.
 265 !     THIS ROUTINE USES AN IMPLICIT APPROACH FOR VERTICAL FLUX
 266 !     DIVERGENCE AND DOES NOT REQUIRE "MITER" TIMESTEPS.
 267 !     IT INCLUDES VERTICAL DIFFUSION IN THE STABLE ATMOSPHERE
 268 !     AND MOIST VERTICAL DIFFUSION IN CLOUDS.
 269 !     SURFACE FLUXES CALCULATED AS IN HIRPBL.
 270 !     5-LAYER SOIL MODEL OPTION REQUIRED IN SLAB DUE TO LONG TIMESTEP
 271 !
 272 !     CODED BY SONG-YOU HONG (NCEP), IMPLEMENTED BY JIMY DUDHIA (NCAR)
 273 !     FALL 1996
 274 !
 275 !     REFERENCES:
 276 !
 277 !        HONG AND PAN (1996), MON. WEA. REV.
 278 !        TROEN AND MAHRT (1986), BOUNDARY LAYER MET.
 279 !
 280 !     CHANGES:
 281 !        INCREASE RLAM FROM 30 TO 150, AND CHANGE FREE ATMOSPHERE
 282 !        STABILITY FUNCTION TO INCREASE VERTICAL DIFFUSION
 283 !        (HONG, JUNE 1997)
 284 !
 285 !        PUT LOWER LIMIT ON PSI FOR STABLE CONDITIONS. THIS WILL
 286 !        PREVENT FLUXES FROM BECOMING TOO SMALL (DUDHIA, OCTOBER 1997)
 287 !
 288 !        CORRECTION TO REGIME CALCULATION. THIS WILL ALLOW POINTS IN
 289 !        REGIME 4 MUCH MORE FREQUENTLY GIVING LARGER SURFACE FLUXES
 290 !        REGIME 3 NO LONGER USES HOL < 1.5 OR THVX LAPSE-RATE CHECK
 291 !        IN MRF SCHEME. THIS WILL MAKE REGIME 3 MUCH LESS FREQUENT.
 292 !
 293 !        ADD SURFACE PRESSURE, PS(I), ARRAY FOR EFFICIENCY
 294 !
 295 !        FIX FOR PROBLEM WITH THIN LAYERS AND HIGH ROUGHNESS
 296 !
 297 !        CHARNOCK CONSTANT NOW COMES FROM NAMELIST (DEFAULT SAME)
 298 !
 299 !-------------------------------------------------------------------
 300
 301       REAL      RLAM,PRMIN,PRMAX,XKZMIN,XKZMAX,RIMIN,BRCR,         &
 302                 CFAC,PFAC,SFCFRAC,CKZ,ZFMIN,APHI5,APHI16,GAMCRT,   &
 303                 GAMCRQ,XKA,PRT
 304
 305       PARAMETER (RLAM=150.,PRMIN=0.5,PRMAX=4.)
 306       PARAMETER (XKZMIN=0.01,XKZMAX=1000.,RIMIN=-100.)
 307       PARAMETER (BRCR=0.5,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1)
 308       PARAMETER (CKZ=0.001,ZFMIN=1.E-8,APHI5=5.,APHI16=16.)
 309       PARAMETER (GAMCRT=3.,GAMCRQ=2.E-3)
 310       PARAMETER (XKA=2.4E-5)
 311       PARAMETER (PRT=1.)
 312 !
 313       INTEGER,  INTENT(IN   )   ::      ids,ide, jds,jde, kds,kde, &
 314                                         ims,ime, jms,jme, kms,kme, &
 315                                         its,ite, jts,jte, kts,kte, &
 316                                         J
 317 !
 318       LOGICAL,  INTENT(IN)      ::                        FLAG_QI
 319 !
 320       REAL,     INTENT(IN   )   ::      P1000mb
 321       REAL,     INTENT(IN   )   ::      DT,DTMIN,CP,G,ROVCP,       &
 322                                         ROVG,R,XLV,RV
 323
 324       REAL,     INTENT(IN )     ::     SVP1,SVP2,SVP3,SVPT0
 325       REAL,     INTENT(IN )     ::     EP1,EP2,KARMAN,EOMEG,STBOLT
 326
 327       REAL,     DIMENSION( ims:ime, kms:kme )                    , &
 328                 INTENT(IN   )   ::                           QV2D, &
 329                                                              QC2D, &
 330                                                               P2D, &
 331                                                               T2D
 332 !
 333       REAL,     DIMENSION( ims:ime, kms:kme )                    , &
 334                 INTENT(INOUT)   ::                         U2DTEN, &
 335                                                            V2DTEN, &
 336                                                            T2DTEN, &
 337                                                           QV2DTEN, &
 338                                                           QC2DTEN
 339
 340       REAL,     DIMENSION( ims:ime )                             , &
 341                 INTENT(INOUT)   ::                            HOL, &
 342                                                               UST, &
 343                                                               PBL, &
 344                                                               ZNT
 345
 346       REAL,     DIMENSION( ims:ime )                             , &
 347                 INTENT(IN   )   ::                          XLAND, &
 348                                                               HFX, &
 349                                                               QFX
 350 !
 351 !m The following 5 are changed to memory size---
 352 !
 353      REAL,     DIMENSION( ims:ime ), INTENT(IN   )   ::      PSIM, &
 354                                                              PSIH
 355
 356      REAL,     DIMENSION( ims:ime ), INTENT(INOUT)   ::      WSPD
 357
 358      REAL,     DIMENSION( ims:ime ), INTENT(IN   )   ::    GZ1OZ0, &
 359                                                                BR
 360
 361       REAL,     DIMENSION( ims:ime )                             , &
 362                 INTENT(IN   )               ::             PSFCPA
 363
 364       REAL,     DIMENSION( ims:ime )                             , &
 365                 INTENT(IN   )   ::                            TSK
 366
 367       REAL,     DIMENSION( ims:ime )                             , &
 368                 INTENT(INOUT)   ::                            ZOL
 369
 370       INTEGER,  DIMENSION( ims:ime )                             , &
 371                 INTENT(OUT  )   ::                         KPBL1D
 372
 373       REAL,     DIMENSION( ims:ime, kms:kme )                    , &
 374                 INTENT(IN   )   ::                            U2D, &
 375                                                               V2D
 376
 377 ! MODULE-LOCAL VARIABLES (DEFINED IN SUBROUTINE MRF)
 378 !
 379       REAL,     DIMENSION( its:ite, kts:kte ) ,                    &
 380                 INTENT(IN)      ::                         dz8w2d, &
 381                                                               z2d
 382 !
 383 !
 384 ! Optional
 385 !
 386       REAL,     DIMENSION( ims:ime )                             , &
 387                 OPTIONAL                                         , &
 388                 INTENT(INOUT)   ::                         REGIME
 389
 390       REAL,     DIMENSION( ims:ime, kms:kme )                    , &
 391                 OPTIONAL                                         , &
 392                 INTENT(IN   )   ::                           QI2D
 393
 394       REAL,     DIMENSION( ims:ime, kms:kme )                    , &
 395                 OPTIONAL                                         , &
 396                 INTENT(INOUT)   ::                        QI2DTEN
 397
 398 ! LOCAL VARS
 399
 400       REAL,     DIMENSION( its:ite, kts:kte+1 ) ::             ZQ
 401
 402       REAL,     DIMENSION( its:ite, kts:kte )   ::                 &
 403                                                          UX,VX,QX, &
 404                                                      QCX,THX,THVX, &
 405                                                           DZQ,DZA, &
 406                                                         TTNP,QTNP, &
 407                                                          QCTNP,ZA, &
 408                                                           UXS,VXS, &
 409                                                          THXS,QXS, &
 410                                                          QCXS,QIX, &
 411                                                        QITNP,QIXS, &
 412                                                         UTNP,VTNP
 413 !
 414       REAL,    DIMENSION( its:ite )             ::     QIXSV,RHOX, &
 415                                                      WSPD1,GOVRTH, &
 416                                                        PBL0,THXSV, &
 417                                                         UXSV,VXSV, &
 418                                                        QXSV,QCXSV, &
 419                                                      QGH,TGDSA,PS
 420
 421       INTEGER                                   ::   ILXM,JLXM,KL, &
 422                                                    KLM,KLP1,KLPBL
 423 !
 424       INTEGER, DIMENSION( its:ite )             ::     KPBL,KPBL0
 425 !
 426       REAL,    DIMENSION( its:ite, kts:kte )    ::      SCR3,SCR4
 427 !
 428       REAL,    DIMENSION( its:ite )             ::           DUM1, &
 429                                                            XKZMKL
 430 !
 431       REAL,    DIMENSION( its:ite )             ::    ZL1,THERMAL, &
 432                                                      WSCALE,HGAMT, &
 433                                                        HGAMQ,BRDN, &
 434                                                         BRUP,PHIM, &
 435                                                          PHIH,CPM, &
 436                                                       DUSFC,DVSFC, &
 437                                                       DTSFC,DQSFC
 438
 439 !
 440       REAL,    DIMENSION( its:ite, kts:kte )    ::      XKZM,XKZH, &
 441                                                             A1,A2, &
 442                                                             AD,AU, &
 443                                                                TX
 444 !
 445       REAL,    DIMENSION( its:ite, kts:kte )  ::               AL
 446 !
 447       LOGICAL, DIMENSION( its:ite )             ::         PBLFLG, &
 448                                                            SFCFLG, &
 449                                                            STABLE
 450 !
 451       REAL, DIMENSION( its:ite )                ::           THGB
 452
 453       INTEGER ::  N,I,K,KK,L,NZOL,IMVDIF
 454
 455       INTEGER ::  JBGN,JEND,IBGN,IEND,NK
 456
 457       REAL    ::  ZOLN,X,Y,CONT,CONQ,CONW,PL,THCON,TVCON,E1,DTSTEP
 458       REAL    ::  ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL
 459       REAL    ::  DTTHX,PSIX,DTG,PSIQ,USTM
 460       REAL    ::  DT4,RDT,SPDK2,FM,FH,HOL1,GAMFAC,VPERT,PRNUM
 461       REAL    ::  ZFAC,XKZO,SS,RI,QMEAN,TMEAN,ALPH,CHI,ZK,RL2,DK,SRI
 462       REAL    ::  BRINT,DTODSD,DSIG,RDZ,DSDZT,DSDZQ,DSDZ2,TTEND,QTEND
 463       REAL    ::  UTEND,VTEND,QCTEND,QITEND,TGC,DTODSU
 464
 465 !----------------------------------------------------------------------
 466
 467       KLPBL=1
 468       KL=kte
 469       ILXM=ite-1
 470       JLXM=jte-1
 471       KLM=kte-1
 472       KLP1=kte+1
 473 !
 474       CONT=1000.*CP/G
 475       CONQ=1000.*XLV/G
 476       CONW=1000./G
 477
 478 !-- IMVDIF      imvdif=1 for moist adiabat vertical diffusion
 479
 480       IMVDIF=1
 481
 482 !     DO i=its,ite
 483 !!PS PSFC cmb
 484 !        PSFC(I)=PSFCPA(I)/1000.
 485 !     ENDDO
 486
 487
 488 !
 489 !----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
 490 !
 491       DO 5 I=its,ite
 492         TGDSA(I)=TSK(I)
 493 ! PS PSFC cmb
 494         PS(I)=PSFCPA(I)/1000.
 495 !        THGB(I)=TSK(I)*(100./PS(I))**ROVCP
 496         THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP
 497     5 CONTINUE
 498 !
 499 !-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
 500 !     T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
 501 !
 502 !     *** NOTE ***
 503 !         THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
 504 !         SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
 505 !         TENDENCIES.
 506 !
 507       DO 24 K=kts,kte
 508         NK=kme-K
 509         DO 24 I=its,ite
 510           UX(I,K)=U2D(I,NK)
 511           VX(I,K)=V2D(I,NK)
 512    24 CONTINUE
 513 !
 514 !.....SCR3(I,K) STORE TEMPERATURE,
 515 !     SCR4(I,K) STORE VIRTUAL TEMPERATURE.
 516 !
 517       DO 30 K=kts,kte
 518         NK=kme-K
 519         DO 30 I=its,ite
 520 ! PL cmb
 521           PL=P2D(I,NK)/1000.
 522           SCR3(I,K)=T2D(I,NK)
 523 !          THCON=(100./PL)**ROVCP
 524           THCON=(P1000mb/(PL*1000.))**ROVCP
 525           THX(I,K)=SCR3(I,K)*THCON
 526           TX(I,K)=SCR3(I,K)
 527           SCR4(I,K)=SCR3(I,K)
 528           THVX(I,K)=THX(I,K)
 529           QX(I,K)=0.
 530    30 CONTINUE
 531 !
 532       DO I=its,ite
 533          QGH(i)=0.
 534          CPM(i)=CP
 535       ENDDO
 536 !
 537 !     IF(IDRY.EQ.1)GOTO 80
 538       DO 50 K=kts,kte
 539         NK=kme-K
 540         DO 50 I=its,ite
 541           QX(I,K)=QV2D(I,NK)
 542           TVCON=(1.+EP1*QX(I,K))
 543           THVX(I,K)=THX(I,K)*TVCON
 544           SCR4(I,K)=SCR3(I,K)*TVCON
 545    50 CONTINUE
 546 !
 547       DO 60 I=its,ite
 548         E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
 549         QGH(I)=EP2*E1/(PS(I)-E1)
 550         CPM(I)=CP*(1.+0.8*QX(I,KL))
 551    60 CONTINUE
 552 !
 553 !     IF(IMOIST.EQ.1)GOTO 80
 554       DO 70 K=kts,kte
 555         NK=kme-K
 556         DO 70 I=its,ite
 557           QCX(I,K)=QC2D(I,NK)
 558    70 CONTINUE
 559
 560       IF (flag_QI .AND. PRESENT( QI2D ) ) THEN
 561          DO K=kts,kte
 562             NK=kme-K
 563             DO I=its,ite
 564                QIX(I,K)=QI2D(I,NK)
 565             ENDDO
 566          ENDDO
 567       ELSE
 568          DO K=kts,kte
 569             NK=kme-K
 570             DO I=its,ite
 571                QIX(I,K)=0.
 572             ENDDO
 573          ENDDO
 574       ENDIF
 575
 576    80 CONTINUE
 577
 578 !
 579 !-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
 580 !     LEVEL, AND THE LAYER THICKNESSES.
 581 !
 582       DO 90 I=its,ite
 583         ZQ(I,KLP1)=0.
 584         RHOX(I)=PS(I)*1000./(R*SCR4(I,KL))
 585    90 CONTINUE
 586 !
 587       DO 110 KK=kts,kte
 588         K=kme-KK
 589         DO 100 I=its,ite
 590           DUM1(I)=ZQ(I,K+1)
 591   100   CONTINUE
 592 !
 593         DO 110 I=its,ite
 594            ZQ(I,K)=dz8w2d(I,K)+DUM1(I)
 595   110   CONTINUE
 596 !
 597       DO 120 K=kts,kte
 598         DO 120 I=its,ite
 599           ZA(I,K)=0.5*(ZQ(I,K)+ZQ(I,K+1))
 600           DZQ(I,K)=ZQ(I,K)-ZQ(I,K+1)
 601   120 CONTINUE
 602 !
 603       DO 130 K=kts,kte-1
 604         DO 130 I=its,ite
 605           DZA(I,K)=ZA(I,K)-ZA(I,K+1)
 606   130 CONTINUE
 607
 608       DTSTEP=DT
 609 !
 610       DO 160 I=its,ite
 611         GOVRTH(I)=G/THX(I,KL)
 612   160 CONTINUE
 613 !
 614 !-----INITIALIZE VERTICAL TENDENCIES AND
 615 !
 616       DO I=its,ite
 617       DO K=kts,kte
 618          UTNP(i,k)=0.
 619          VTNP(i,k)=0.
 620          TTNP(i,k)=0.
 621       ENDDO
 622       ENDDO
 623 !
 624 !     IF(IDRY.EQ.1)GOTO 250
 625       DO 230 K=kts,kte
 626         DO 230 I=its,ite
 627           QTNP(I,K)=0.
 628   230 CONTINUE
 629 !
 630 !     IF(IMOIST.EQ.1)GOTO 250
 631       DO 240 K=kts,kte
 632         DO 240 I=its,ite
 633           QCTNP(I,K)=0.
 634           QITNP(I,K)=0.
 635   240 CONTINUE
 636
 637   250 CONTINUE
 638 !
 639 !-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
 640 !     AKB(1976), EQ(12).
 641
 642 !     DO 260 I=its,ite
 643 !       GZ1OZ0(I)=ALOG(ZA(I,KL)/ZNT(I))
 644 !       IF((XLAND(I)-1.5).GE.0)THEN
 645 !         ZL=ZNT(I)
 646 !       ELSE
 647 !         ZL=0.01
 648 !       ENDIF
 649 !       WSPD(I)=SQRT(UX(I,KL)*UX(I,KL)+VX(I,KL)*VX(I,KL))
 650 !       TSKV=THGB(I)*(1.+EP1*QGH(I)*MAVAIL(I))
 651 !       DTHVDZ=(THVX(I,KL)-TSKV)
 652 !       IF(-DTHVDZ.GE.0)THEN
 653 !         DTHVM=-DTHVDZ
 654 !       ELSE
 655 !         DTHVM=0.
 656 !       ENDIF
 657 !       VCONV=VCONVC*SQRT(DTHVM)
 658 !       WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV)
 659 !       WSPD(I)=AMAX1(WSPD(I),1.)
 660 !       BR(I)=GOVRTH(I)*ZA(I,KL)*DTHVDZ/(WSPD(I)*WSPD(I))
 661 ! 260 CONTINUE
 662
 663 !!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
 664 !!
 665 !!     THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
 666 !!     AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
 667 !!
 668 !!     CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
 669 !!
 670 !!        1. BR .GE. 0.2;
 671 !!               REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
 672 !!
 673 !!        2. BR .LT. 0.2 .AND. BR .GT. 0.0;
 674 !!               REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
 675 !!               (REGIME=2),
 676 !!
 677 !!        3. BR .EQ. 0.0
 678 !!               REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
 679 !!
 680 !!        4. BR .LT. 0.0
 681 !!               REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
 682 !!
 683 !!-----
 684 !
 685 !      DO 320 I=its,ite
 686 !!----
 687 !!--     REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT
 688 !!--          IF(BR(I).LT.0..AND.HOL(I).GT.1.5)GOTO 310
 689 !
 690 !        IF(BR(I).LT.0.)GOTO 310
 691 !!
 692 !!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
 693 !!
 694 !        IF(BR(I).LT.0.2)GOTO 270
 695 !        REGIME(I)=1.
 696 !        PSIM(I)=-10.*GZ1OZ0(I)
 697 !!    LOWER LIMIT ON PSI IN STABLE CONDITIONS
 698 !        PSIM(I)=AMAX1(PSIM(I),-10.)
 699 !        PSIH(I)=PSIM(I)
 700 !        HOL(I)=0.0
 701 !        PBL(I)=0.0
 702 !        GOTO 320
 703 !!
 704 !!-----CLASS 2; DAMPED MECHANICAL TURBULENCE:
 705 !!
 706 !  270   IF(BR(I).EQ.0.0)GOTO 280
 707 !        REGIME(I)=2.
 708 !        PSIM(I)=-5.0*BR(I)*GZ1OZ0(I)/(1.1-5.0*BR(I))
 709 !!    LOWER LIMIT ON PSI IN STABLE CONDITIONS
 710 !        PSIM(I)=AMAX1(PSIM(I),-10.)
 711 !!.....AKB(1976), EQ(16).
 712 !        PSIH(I)=PSIM(I)
 713 !        HOL(I)=0.0
 714 !        PBL(I)=0.0
 715 !        GOTO 320
 716 !!
 717 !!-----CLASS 3; FORCED CONVECTION:
 718 !!
 719 !  280   REGIME(I)=3.
 720 !        PSIM(I)=0.0
 721 !        PSIH(I)=PSIM(I)
 722 !
 723 !! special use kte instead of kme
 724 !
 725 !        DO 290 KK=kts,kte-1
 726 !          K=kte-KK
 727 !          IF(THVX(I,K).GT.THVX(I,KL))GOTO 300
 728 !  290   CONTINUE
 729 !        STOP 290
 730 !  300   PBL(I)=ZQ(I,K+1)
 731 !        IF(UST(I).LT.0.01)THEN
 732 !          ZOL(I)=BR(I)*GZ1OZ0(I)
 733 !        ELSE
 734 !          ZOL(I)=KARMAN*GOVRTH(I)*ZA(I,KL)*MOL(I,J)/(UST(I)*UST(I))
 735 !        ENDIF
 736 !        HOL(I)=-ZOL(I)*PBL(I)/ZA(I,KL)
 737 !        GOTO 320
 738 !
 739 !!-----CLASS 4; FREE CONVECTION:
 740 !
 741 !! 310     IF(THVX(I,KLM).GT.THVX(I,KL))GOTO 280
 742 !
 743 !  310   CONTINUE
 744 !        REGIME(I)=4.
 745 !        IF(UST(I).LT.0.01)THEN
 746 !          ZOL(I)=BR(I)*GZ1OZ0(I)
 747 !        ELSE
 748 !          ZOL(I)=KARMAN*GOVRTH(I)*ZA(I,KL)*MOL(I,J)/(UST(I)*UST(I))
 749 !        ENDIF
 750 !        ZOL(I)=AMIN1(ZOL(I),0.)
 751 !        ZOL(I)=AMAX1(ZOL(I),-9.9999)
 752 !        NZOL=INT(-ZOL(I)*100.)
 753 !        RZOL=-ZOL(I)*100.-NZOL
 754 !        PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))
 755 !        PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))
 756 !!---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
 757 !!---  THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
 758 !        PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
 759 !        PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
 760 !  320 CONTINUE
 761
 762 !-----COMPUTE THE FRICTIONAL VELOCITY:
 763 !     ZA(1982) EQS(2.60),(2.61).
 764
 765       DO 330 I=its,ite
 766         DTG=THX(I,KL)-THGB(I)
 767         PSIX=GZ1OZ0(I)-PSIM(I)
 768         IF((XLAND(I)-1.5).GE.0)THEN
 769           ZL=ZNT(I)
 770         ELSE
 771           ZL=0.01
 772         ENDIF
 773         PSIQ=ALOG(KARMAN*UST(I)*ZA(I,KL)/XKA+ZA(I,KL)/ZL)-PSIH(I)
 774         UST(I)=KARMAN*WSPD(I)/PSIX
 775 !
 776         USTM=AMAX1(UST(I),0.1)
 777         IF((XLAND(I)-1.5).GE.0)THEN
 778           UST(I)=UST(I)
 779         ELSE
 780           UST(I)=USTM
 781         ENDIF
 782 !       MOL(I,J)=KARMAN*DTG/(GZ1OZ0(I)-PSIH(I))/PRT
 783   330 CONTINUE
 784 !
 785       DO 420 I=its,ite
 786         WSPD1(I)=SQRT(UX(I,KL)*UX(I,KL)+VX(I,KL)*VX(I,KL))+1.E-9
 787   420 CONTINUE
 788 !
 789 !---- COMPUTE VERTICAL DIFFUSION
 790 !
 791 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 792 !     COMPUTE PRELIMINARY VARIABLES
 793 !
 794 !
 795       DT4=2.*DTSTEP
 796       RDT=1./DT4
 797 !
 798       DO I=its,ite
 799         HGAMT(I)=0.
 800         HGAMQ(I)=0.
 801         WSCALE(I)=0.
 802         KPBL(I)=KL
 803         PBL(I)=ZQ(I,KL)
 804         KPBL0(I)=KL
 805         PBL0(I)=ZQ(I,KL)
 806         PBLFLG(I)=.TRUE.
 807         SFCFLG(I)=.TRUE.
 808         IF(BR(I).GT.0.0)SFCFLG(I)=.FALSE.
 809         ZL1(I)=ZA(I,KL)
 810         THERMAL(I)=THVX(I,KL)
 811       ENDDO
 812
 813 !     COMPUTE THE FIRST GUESS OF PBL HEIGHT
 814
 815       DO I=its,ite
 816         STABLE(I)=.FALSE.
 817         BRUP(I)=BR(I)
 818       ENDDO
 819       DO K=KLM,KLPBL,-1
 820         DO I=its,ite
 821           IF(.NOT.STABLE(I))THEN
 822             BRDN(I)=BRUP(I)
 823             SPDK2=MAX(UX(I,K)**2+VX(I,K)**2,1.)
 824             BRUP(I)=(THVX(I,K)-THERMAL(I))*(G*ZA(I,K)/THVX(I,KL))/SPDK2
 825             KPBL(I)=K
 826             STABLE(I)=BRUP(I).GT.BRCR
 827           ENDIF
 828         ENDDO
 829       ENDDO
 830 !
 831       DO I=its,ite
 832         K=KPBL(I)
 833         IF(BRDN(I).GE.BRCR)THEN
 834           BRINT=0.
 835         ELSEIF(BRUP(I).LE.BRCR)THEN
 836           BRINT=1.
 837         ELSE
 838           BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))
 839         ENDIF
 840         PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
 841         IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1
 842       ENDDO
 843 !
 844       DO I=its,ite
 845         FM=GZ1OZ0(I)-PSIM(I)
 846         FH=GZ1OZ0(I)-PSIH(I)
 847         HOL(I)=MAX(BR(I)*FM*FM/FH,RIMIN)
 848         IF(SFCFLG(I))THEN
 849           HOL(I)=MIN(HOL(I),-ZFMIN)
 850         ELSE
 851           HOL(I)=MAX(HOL(I),ZFMIN)
 852         ENDIF
 853 !
 854         HOL1=HOL(I)*PBL(I)/ZL1(I)*SFCFRAC
 855         HOL(I)=-HOL(I)*PBL(I)/ZL1(I)
 856         IF(SFCFLG(I))THEN
 857           PHIM(I)=(1.-APHI16*HOL1)**(-1./4.)
 858           PHIH(I)=(1.-APHI16*HOL1)**(-1./2.)
 859         ELSE
 860           PHIM(I)=(1.+APHI5*HOL1)
 861           PHIH(I)=PHIM(I)
 862         ENDIF
 863         WSCALE(I)=UST(I)/PHIM(I)
 864         WSCALE(I)=MIN(WSCALE(I),UST(I)*APHI16)
 865         WSCALE(I)=MAX(WSCALE(I),UST(I)/APHI5)
 866       ENDDO
 867
 868 !     COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION
 869 !     UNDER UNSTABLE CONDITIONS
 870
 871       DO I=its,ite
 872         IF(SFCFLG(I))THEN
 873           GAMFAC=CFAC/RHOX(I)/WSCALE(I)
 874           HGAMT(I)=MIN(GAMFAC*HFX(I)/CPM(I),GAMCRT)
 875           HGAMQ(I)=MIN(GAMFAC*QFX(I),GAMCRQ)
 876           IF((XLAND(I)-1.5).GE.0)HGAMQ(I)=0.
 877           VPERT=HGAMT(I)+EP1*THX(I,KL)*HGAMQ(I)
 878           VPERT=MIN(VPERT,GAMCRT)
 879           THERMAL(I)=THERMAL(I)+MAX(VPERT,0.)
 880           HGAMT(I)=MAX(HGAMT(I),0.0)
 881           HGAMQ(I)=MAX(HGAMQ(I),0.0)
 882         ELSE
 883           PBLFLG(I)=.FALSE.
 884         ENDIF
 885       ENDDO
 886 !
 887       DO I=its,ite
 888         IF(PBLFLG(I))THEN
 889           KPBL(I)=KL
 890           PBL(I)=ZQ(I,KL)
 891         ENDIF
 892       ENDDO
 893 !
 894 !     ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL
 895 !
 896       DO I=its,ite
 897         IF(PBLFLG(I))THEN
 898           STABLE(I)=.FALSE.
 899           BRUP(I)=BR(I)
 900         ENDIF
 901       ENDDO
 902       DO K=KLM,KLPBL,-1
 903         DO I=its,ite
 904           IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN
 905             BRDN(I)=BRUP(I)
 906             SPDK2=MAX((UX(I,K)**2+VX(I,K)**2),1.)
 907             BRUP(I)=(THVX(I,K)-THERMAL(I))*(G*ZA(I,K)/THVX(I,KL))/SPDK2
 908             KPBL(I)=K
 909             STABLE(I)=BRUP(I).GT.BRCR
 910           ENDIF
 911         ENDDO
 912       ENDDO
 913 !
 914       DO I=its,ite
 915         IF(PBLFLG(I))THEN
 916           K=KPBL(I)
 917           IF(BRDN(I).GE.BRCR)THEN
 918             BRINT=0.
 919           ELSEIF(BRUP(I).LE.BRCR)THEN
 920             BRINT=1.
 921           ELSE
 922             BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))
 923           ENDIF
 924           PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
 925           IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1
 926           IF(KPBL(I).LE.1)PBLFLG(I)=.FALSE.
 927         ENDIF
 928       ENDDO
 929 !
 930 !     DIAGNOSTIC PBL HEIGHT WITH BRCR EFFECTIVELY ZERO (PBL0)
 931 !
 932       DO I=its,ite
 933         IF(PBLFLG(I))THEN
 934           STABLE(I)=.FALSE.
 935           BRUP(I)=BR(I)
 936         ENDIF
 937       ENDDO
 938       DO K=KLM,KLPBL,-1
 939         DO I=its,ite
 940           IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN
 941             BRDN(I)=BRUP(I)
 942             SPDK2=MAX((UX(I,K)**2+VX(I,K)**2),1.)
 943             BRUP(I)=(THVX(I,K)-THERMAL(I))*(G*ZA(I,K)/THVX(I,KL))/SPDK2
 944             KPBL0(I)=K
 945             STABLE(I)=BRUP(I).GT.0.0
 946           ENDIF
 947
 948         ENDDO
 949       ENDDO
 950 !
 951       DO I=its,ite
 952         IF(PBLFLG(I))THEN
 953           K=KPBL0(I)
 954           IF(BRDN(I).GE.0.0)THEN
 955             BRINT=0.
 956           ELSEIF(BRUP(I).LE.0.0)THEN
 957             BRINT=1.
 958           ELSE
 959             BRINT=(0.0-BRDN(I))/(BRUP(I)-BRDN(I))
 960           ENDIF
 961           PBL0(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
 962           IF(PBL0(I).LT.ZQ(I,KPBL0(I)+1))KPBL0(I)=KPBL0(I)+1
 963           IF(KPBL0(I).LE.1)PBLFLG(I)=.FALSE.
 964         ENDIF
 965       ENDDO
 966
 967 !
 968 !     COMPUTE DIFFUSION COEFFICIENTS BELOW PBL
 969 !
 970       DO K=kte,KLPBL,-1
 971         DO I=its,ite
 972           IF(KPBL(I).LT.K)THEN
 973             PRNUM=(PHIH(I)/PHIM(I)+CFAC*KARMAN*SFCFRAC)
 974             PRNUM=MIN(PRNUM,PRMAX)
 975             PRNUM=MAX(PRNUM,PRMIN)
 976             ZFAC=MAX((1.-(ZQ(I,K)-ZL1(I))/(PBL(I)-ZL1(I))),ZFMIN)
 977             XKZO=CKZ*DZA(I,K-1)
 978             XKZM(I,K)=XKZO+WSCALE(I)*KARMAN*ZQ(I,K)*ZFAC**PFAC
 979             XKZH(I,K)=XKZM(I,K)/PRNUM
 980             XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)
 981             XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)
 982             XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)
 983             XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)
 984           ENDIF
 985         ENDDO
 986       ENDDO
 987 !
 988 !     COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE)
 989 !
 990       DO K=kts+1,kte
 991         DO I=its,ite
 992           XKZO=CKZ*DZA(I,K-1)
 993           IF(K.LE.KPBL(I))THEN
 994             SS=((UX(I,K-1)-UX(I,K))*(UX(I,K-1)-UX(I,K))+(VX(I,K-1)-   &
 995                VX(I,K))*(VX(I,K-1)-VX(I,K)))/(DZA(I,K-1)*DZA(I,K-1))+ &
 996                1.E-9
 997             RI=GOVRTH(I)*(THVX(I,K-1)-THVX(I,K))/(SS*DZA(I,K-1))
 998             IF(IMVDIF.EQ.1)THEN
 999               IF((QCX(I,K)+QIX(I,K)).GT.0.01E-3.AND.(QCX(I,K-1)+      &
1000                 QIX(I,K-1)).GT.0.01E-3)THEN
1001 !      IN CLOUD
1002                 QMEAN=0.5*(QX(I,K)+QX(I,K-1))
1003                 TMEAN=0.5*(SCR3(I,K)+SCR3(I,K-1))
1004                 ALPH=XLV*QMEAN/R/TMEAN
1005                 CHI=XLV*XLV*QMEAN/CP/RV/TMEAN/TMEAN
1006                 RI=(1.+ALPH)*(RI-G*G/SS/TMEAN/CP*((CHI-ALPH)/(1.+CHI)))
1007               ENDIF
1008             ENDIF
1009             ZK=KARMAN*ZQ(I,K)
1010             RL2=(ZK*RLAM/(RLAM+ZK))**2
1011             DK=RL2*SQRT(SS)
1012             IF(RI.LT.0.)THEN
1013 ! UNSTABLE REGIME
1014               SRI=SQRT(-RI)
1015               XKZM(I,K)=XKZO+DK*(1+8.*(-RI)/(1+1.746*SRI))
1016               XKZH(I,K)=XKZO+DK*(1+8.*(-RI)/(1+1.286*SRI))
1017             ELSE
1018 ! STABLE REGIME
1019               XKZH(I,K)=XKZO+DK/(1+5.*RI)**2
1020               PRNUM=1.0+2.1*RI
1021               PRNUM=MIN(PRNUM,PRMAX)
1022               XKZM(I,K)=(XKZH(I,K)-XKZO)*PRNUM+XKZO
1023             ENDIF
1024 !
1025             XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)
1026             XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)
1027             XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)
1028             XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)
1029           ENDIF
1030 !
1031         ENDDO
1032       ENDDO
1033
1034 !     COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE
1035
1036       DO I=its,ite
1037       DO K=kts,kte
1038          AU(i,k)=0.
1039          AL(i,k)=0.
1040          AD(i,k)=0.
1041          A1(i,k)=0.
1042          A2(i,k)=0.
1043       ENDDO
1044       ENDDO
1045
1046       DO I=its,ite
1047         AD(I,1)=1.
1048         A1(I,1)=SCR3(I,KL)+HFX(I)/(RHOX(I)*CPM(I))/ZQ(I,KL)*DT4
1049         A2(I,1)=QX(I,KL)+QFX(I)/(RHOX(I))/ZQ(I,KL)*DT4
1050       ENDDO
1051 !
1052       DO K=kte,kts+1,-1
1053         KK=kme-K
1054         DO I=its,ite
1055           DTODSD=DT4/dz8w2d(I,K)
1056           DTODSU=DT4/dz8w2d(I,K-1)
1057           DSIG=z2d(I,K)-z2d(I,K-1)
1058           DSIG=-DSIG
1059           RDZ=1./DZA(I,K-1)
1060           IF(PBLFLG(I).AND.KPBL(I).LT.K)THEN
1061             DSDZT=DSIG*XKZH(I,K)*RDZ*(G/CP-HGAMT(I)/PBL(I))
1062             DSDZQ=DSIG*XKZH(I,K)*RDZ*(-HGAMQ(I)/PBL(I))
1063             A2(I,KK)=A2(I,KK)+DTODSD*DSDZQ
1064             A2(I,KK+1)=QX(I,K-1)-DTODSU*DSDZQ
1065           ELSE
1066             DSDZT=DSIG*XKZH(I,K)*RDZ*(G/CP)
1067             A2(I,KK+1)=QX(I,K-1)
1068           ENDIF
1069           DSDZ2=DSIG*XKZH(I,K)*RDZ*RDZ
1070           AU(I,KK)=-DTODSD*DSDZ2
1071           AL(I,KK)=-DTODSU*DSDZ2
1072           AD(I,KK)=AD(I,KK)-AU(I,KK)
1073           AD(I,KK+1)=1.-AL(I,KK)
1074           A1(I,KK)=A1(I,KK)+DTODSD*DSDZT
1075           A1(I,KK+1)=SCR3(I,K-1)-DTODSU*DSDZT
1076         ENDDO
1077       ENDDO
1078
1079 !     SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE
1080
1081       CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2,                 &
1082                   its,ite,kts,kte                          )
1083
1084 !     RECOVER TENDENCIES OF HEAT AND MOISTURE
1085
1086       DO K=kte,kts,-1
1087         KK=kme-K
1088         DO I=its,ite
1089           TTEND=(A1(I,KK)-SCR3(I,K))*RDT
1090           QTEND=(A2(I,KK)-QX(I,K))*RDT
1091           TTNP(I,K)=TTNP(I,K)+TTEND
1092           QTNP(I,K)=QTNP(I,K)+QTEND
1093         ENDDO
1094       ENDDO
1095
1096 !     COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM
1097
1098       DO I=its,ite
1099       DO K=kts,kte
1100          AU(i,k)=0.
1101          AL(i,k)=0.
1102          AD(i,k)=0.
1103          A1(i,k)=0.
1104          A2(i,k)=0.
1105       ENDDO
1106       ENDDO
1107
1108       DO I=its,ite
1109         AD(I,1)=1.
1110         A1(I,1)=UX(I,KL)-UX(I,KL)/WSPD1(I)*UST(I)*UST(I)/ZQ(I,KL)  &
1111                 *DT4*(WSPD1(I)/WSPD(I))**2
1112         A2(I,1)=VX(I,KL)-VX(I,KL)/WSPD1(I)*UST(I)*UST(I)/ZQ(I,KL)  &
1113                 *DT4*(WSPD1(I)/WSPD(I))**2
1114       ENDDO
1115 !
1116       DO K=kte,kts+1,-1
1117         KK=kme-K
1118         DO I=its,ite
1119           DTODSD=DT4/dz8w2d(I,K)
1120           DTODSU=DT4/dz8w2d(I,K-1)
1121           DSIG=z2d(I,K)-z2d(I,K-1)
1122           DSIG=-DSIG
1123           RDZ=1./DZA(I,K-1)
1124           DSDZ2=DSIG*XKZM(I,K)*RDZ*RDZ
1125           AU(I,KK)=-DTODSD*DSDZ2
1126           AL(I,KK)=-DTODSU*DSDZ2
1127           AD(I,KK)=AD(I,KK)-AU(I,KK)
1128           AD(I,KK+1)=1.-AL(I,KK)
1129           A1(I,KK+1)=UX(I,K-1)
1130           A2(I,KK+1)=VX(I,K-1)
1131         ENDDO
1132       ENDDO
1133
1134 !     SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM
1135
1136       CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2,                 &
1137                   its,ite,kts,kte                          )
1138
1139 !     RECOVER TENDENCIES OF MOMENTUM
1140
1141       DO K=kte,kts,-1
1142         KK=kme-K
1143         DO I=its,ite
1144           UTEND=(A1(I,KK)-UX(I,K))*RDT
1145           VTEND=(A2(I,KK)-VX(I,K))*RDT
1146           UTNP(I,K)=UTNP(I,K)+UTEND
1147           VTNP(I,K)=VTNP(I,K)+VTEND
1148         ENDDO
1149       ENDDO
1150
1151 !     COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR CLOUD
1152
1153       DO I=its,ite
1154       DO K=kts,kte
1155          AU(i,k)=0.
1156          AL(i,k)=0.
1157          AD(i,k)=0.
1158          A1(i,k)=0.
1159          A2(i,k)=0.
1160       ENDDO
1161       ENDDO
1162
1163 !     IF(IMOIST.EQ.1)GOTO 690
1164       DO I=its,ite
1165         AD(I,1)=1.
1166         A1(I,1)=QCX(I,KL)
1167         A2(I,1)=QIX(I,KL)
1168       ENDDO
1169 !
1170       DO K=kte,kts+1,-1
1171         KK=kme-K
1172         DO I=its,ite
1173           DTODSD=DT4/dz8w2d(I,K)
1174           DTODSU=DT4/dz8w2d(I,K-1)
1175           DSIG=z2d(I,K)-z2d(I,K-1)
1176           DSIG=-DSIG
1177           RDZ=1./DZA(I,K-1)
1178           A1(I,KK+1)=QCX(I,K-1)
1179           A2(I,KK+1)=QIX(I,K-1)
1180           DSDZ2=DSIG*XKZH(I,K)*RDZ*RDZ
1181           AU(I,KK)=-DTODSD*DSDZ2
1182           AL(I,KK)=-DTODSU*DSDZ2
1183           AD(I,KK)=AD(I,KK)-AU(I,KK)
1184           AD(I,KK+1)=1.-AL(I,KK)
1185         ENDDO
1186       ENDDO
1187
1188 !     SOLVE TRIDIAGONAL PROBLEM FOR CLOUD
1189
1190       CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2,                 &
1191                   its,ite,kts,kte                          )
1192 !
1193       DO K=kte,kts,-1
1194         KK=kme-K
1195         DO I=its,ite
1196           QCTEND=(A1(I,KK)-QCX(I,K))*RDT
1197           QITEND=(A2(I,KK)-QIX(I,K))*RDT
1198           QCTNP(I,K)=QCTNP(I,K)+QCTEND
1199           QITNP(I,K)=QITNP(I,K)+QITEND
1200         ENDDO
1201       ENDDO
1202 !
1203 !---- END OF VERTICAL DIFFUSION
1204 !
1205   690 CONTINUE
1206 !
1207 !-----CALCULATION OF NEW VALUES DUE TO VERTICAL EXCHANGE PROCESSES IS
1208 !     COMPLETED. THE FINAL STEP IS TO ADD THE TENDENCIES CALCULATED
1209 !     IN HIRPBL TO THOSE OF MM4.
1210
1211       DO 820 K=kts,kte
1212         NK=kme-K
1213         DO 820 I=its,ite
1214           U2DTEN(I,NK)=UTNP(I,K)
1215           V2DTEN(I,NK)=VTNP(I,K)
1216   820   CONTINUE
1217 !
1218 !     IF(J.EQ.1.AND.IN.GT.1)GOTO 860
1219 !SUE  JBGN=3
1220 !SUE  JEND=JLXM-1
1221
1222 ! change when nest
1223 !SUE  JBGN=2
1224 !SUE  JEND=JLXM
1225
1226       JBGN=jts
1227       JEND=jte
1228       IBGN=its
1229       IEND=ite
1230
1231 !     IF(J.LT.JBGN.OR.J.GT.JEND)GOTO 860
1232 !SUE  IBGN=3
1233 !SUE  IEND=ILXM-1
1234
1235 ! change when nest
1236 !SUE  IBGN=2
1237 !SUE  IEND=ILXM
1238
1239       DO 830 K=kts,kte
1240         NK=kme-K
1241         DO 830 I=IBGN,IEND
1242           T2DTEN(I,NK)=TTNP(I,K)
1243   830   CONTINUE
1244 !
1245 !     IF(IDRY.EQ.1)GOTO 860
1246       DO 840 K=kts,kte
1247         NK=kme-K
1248         DO 840 I=IBGN,IEND
1249           QV2DTEN(I,NK)=QTNP(I,K)
1250   840   CONTINUE
1251
1252 !     IF(IMOIST.EQ.1)GOTO 860
1253       DO 850 K=kts,kte
1254         NK=kme-K
1255         DO 850 I=IBGN,IEND
1256            QC2DTEN(I,NK)=QCTNP(I,K)
1257   850   CONTINUE
1258
1259       IF(flag_QI .AND. PRESENT( QI2DTEN ) ) THEN
1260         DO K=kts,kte
1261         NK=kme-K
1262           DO I=IBGN,IEND
1263              QI2DTEN(I,NK)=QITNP(I,K)
1264           ENDDO
1265         ENDDO
1266       ENDIF
1267
1268   860 CONTINUE
1269 !
1270 !-----APPLY ASSELIN FILTER TO TGD FOR LARGE TIME STEP:
1271 !
1272 !     DO 885 I=its,ite
1273 !       TSK(I)=TSK(I)*(PS(I)/100.)**ROVCP
1274 ! 885 CONTINUE
1275 !
1276   940 CONTINUE
1277 !
1278 ! KPBL IS NEEDED FOR THE FDDA, AND SINCE THERE IS NO LONGER JUST ONE
1279 ! LARGE "J LOOP" IT MUST BE STORED AS (I,J)...
1280 !
1281 ! USE NEW DIAGNOSED PBL DEPTH (CALCULATED WITH brcr=0.0)
1282 ! PBL IS USED FOR OUTPUT AND NEXT-TIME-STEP BELJAARS CALC IN SFCLAY
1283       DO 950 I=its,ite
1284         KPBL1D(I)=KPBL0(I)
1285         PBL(I)=PBL0(I)
1286   950 CONTINUE
1287
1288    END SUBROUTINE MRF2D
1289
1290 !================================================================
1291    SUBROUTINE TRIDI2(CL,CM,CU,R1,R2,AU,A1,A2,                   &
1292                      its,ite,kts,kte                            )
1293 !----------------------------------------------------------------
1294    IMPLICIT NONE
1295 !----------------------------------------------------------------
1296
1297    INTEGER, INTENT(IN )      ::     its,ite, kts,kte
1298
1299    REAL, DIMENSION( its:ite, kts+1:kte+1 )                    , &
1300          INTENT(IN   )  ::                                  CL
1301
1302    REAL, DIMENSION( its:ite, kts:kte )                        , &
1303          INTENT(IN   )  ::                                  CM, &
1304                                                             R1, &
1305                                                             R2
1306    REAL, DIMENSION( its:ite, kts:kte )                        , &
1307          INTENT(INOUT)  ::                                  AU, &
1308                                                             CU, &
1309                                                             A1, &
1310                                                             A2
1311
1312    REAL    :: FK
1313    INTEGER :: I,K,L,N
1314
1315 !----------------------------------------------------------------
1316
1317    L=ite
1318    N=kte
1319
1320    DO I=its,L
1321      FK=1./CM(I,1)
1322      AU(I,1)=FK*CU(I,1)
1323      A1(I,1)=FK*R1(I,1)
1324      A2(I,1)=FK*R2(I,1)
1325    ENDDO
1326    DO K=2,N-1
1327      DO I=its,L
1328        FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
1329        AU(I,K)=FK*CU(I,K)
1330        A1(I,K)=FK*(R1(I,K)-CL(I,K)*A1(I,K-1))
1331        A2(I,K)=FK*(R2(I,K)-CL(I,K)*A2(I,K-1))
1332      ENDDO
1333    ENDDO
1334    DO I=its,L
1335      FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
1336      A1(I,N)=FK*(R1(I,N)-CL(I,N)*A1(I,N-1))
1337      A2(I,N)=FK*(R2(I,N)-CL(I,N)*A2(I,N-1))
1338
1339    ENDDO
1340    DO K=N-1,kts,-1
1341      DO I=its,L
1342        A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
1343        A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
1344      ENDDO
1345    ENDDO
1346
1347    END SUBROUTINE TRIDI2
1348
1349 !===================================================================
1350    SUBROUTINE mrfinit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,           &
1351                       RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR,       &
1352                       restart, allowed_to_read ,                   &
1353                       ids, ide, jds, jde, kds, kde,                &
1354                       ims, ime, jms, jme, kms, kme,                &
1355                       its, ite, jts, jte, kts, kte                 )
1356 !-------------------------------------------------------------------
1357    IMPLICIT NONE
1358 !-------------------------------------------------------------------
1359    LOGICAL , INTENT(IN)          :: restart , allowed_to_read
1360    INTEGER , INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
1361                                      ims, ime, jms, jme, kms, kme, &
1362                                      its, ite, jts, jte, kts, kte
1363    INTEGER , INTENT(IN)          ::  P_QI,P_FIRST_SCALAR
1364
1365    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::         &
1366                                                          RUBLTEN, &
1367                                                          RVBLTEN, &
1368                                                          RTHBLTEN, &
1369                                                          RQVBLTEN, &
1370                                                          RQCBLTEN, &
1371                                                          RQIBLTEN
1372    INTEGER :: i, j, k, itf, jtf, ktf
1373
1374    jtf=min0(jte,jde-1)
1375    ktf=min0(kte,kde-1)
1376    itf=min0(ite,ide-1)
1377
1378    IF(.not.restart)THEN
1379      DO j=jts,jtf
1380      DO k=kts,ktf
1381      DO i=its,itf
1382         RUBLTEN(i,k,j)=0.
1383         RVBLTEN(i,k,j)=0.
1384         RTHBLTEN(i,k,j)=0.
1385         RQVBLTEN(i,k,j)=0.
1386         RQCBLTEN(i,k,j)=0.
1387      ENDDO
1388      ENDDO
1389      ENDDO
1390    ENDIF
1391
1392    IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
1393       DO j=jts,jtf
1394       DO k=kts,ktf
1395       DO i=its,itf
1396          RQIBLTEN(i,k,j)=0.
1397       ENDDO
1398       ENDDO
1399       ENDDO
1400    ENDIF
1401
1402    END SUBROUTINE mrfinit
1403
1404 !-------------------------------------------------------------------
1405
1406 END MODULE module_bl_mrf
1407