wrfv2_fire/phys/module_bl_acm.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3 MODULE module_bl_acm
   4
   5 !  USE module_model_constants
   6
   7   REAL, PARAMETER      :: RIC    = 0.25                ! critical Richardson number
   8   REAL, PARAMETER      :: CRANKP = 0.5                 ! CRANK-NIC PARAMETER
   9
  10 CONTAINS
  11
  12 !-------------------------------------------------------------------
  13    SUBROUTINE ACMPBL(XTIME,    DTPBL,    ZNW, SIGMAH,                &
  14                      U3D,      V3D,      PP3D,  DZ8W, TH3D, T3D,        &
  15                      QV3D,     QC3D,     QI3D, RR3D,                    &       !!  Moisture
  16                      UST,      HFX,  QFX,     TSK,                      &
  17                      PSFC,     EP1,      G,                             &
  18                      ROVCP,    RD,       CPD,                           &
  19                      PBLH,      KPBL2D,   REGIME,                       &
  20                      GZ1OZ0,   WSPD,     PSIM, MUT,                     &
  21                      RUBLTEN,  RVBLTEN,  RTHBLTEN,                      &
  22                      RQVBLTEN, RQCBLTEN, RQIBLTEN,                      &
  23                      ids,ide, jds,jde, kds,kde,                         &
  24                      ims,ime, jms,jme, kms,kme,                         &
  25                      its,ite, jts,jte, kts,kte)
  26 !**********************************************************************
  27
  28 !   THIS MODULE COMPUTES VERTICAL MIXING IN AND ABOVE THE PBL ACCORDING TO
  29 !   THE ASYMMETRICAL CONVECTIVE MODEL, VERSION 2  (ACM2), WHICH IS A COMBINED
  30 !   LOCAL NON-LOCAL CLOSURE SCHEME BASED ON THE ORIGINAL ACM (PLEIM AND CHANG 1992)
  31 !
  32 !   REFERENCES:
  33 !   Pleim (2007) A combined local and non-local closure model for the atmospheric
  34 !                boundary layer.  Part1: Model description and testing.
  35 !                JAMC, 46, 1383-1395
  36 !   Pleim (2007) A combined local and non-local closure model for the atmospheric
  37 !                boundary layer.  Part2: Application and evaluation in a mesoscale
  38 !                meteorology model.  JAMC, 46, 1396-1409
  39 !
  40 !  REVISION HISTORY:
  41 !     AX     3/2005 - developed WRF version based on the MM5 PX LSM
  42 !     RG and JP 7/2006 - Finished WRF adaptation
  43 !
  44 !**********************************************************************
  45 !   ARGUMENT LIST:
  46 !
  47 !... Inputs:
  48 !-- XTIME           Time since simulation start (min)
  49 !-- DTPBL           PBL time step
  50 !-- ZNW             Sigma at full layer
  51 !-- SIGMAH          Sigma at half layer
  52 !-- U3D             3D u-velocity interpolated to theta points (m/s)
  53 !-- V3D             3D v-velocity interpolated to theta points (m/s)
  54 !-- PP3D            Pressure at half levels (Pa)
  55 !-- DZ8W            dz between full levels (m)
  56 !-- TH3D            Potential Temperature (K)
  57 !-- T3D             Temperature (K)
  58 !-- QV3D            3D water vapor mixing ratio (Kg/Kg)
  59 !-- QC3D            3D cloud mixing ratio (Kg/Kg)
  60 !-- QI3D            3D ice mixing ratio (Kg/Kg)
  61 !-- RR3D            3D dry air density (kg/m^3)
  62 !-- UST             Friction Velocity (m/s)
  63 !-- HFX             Upward heat flux at the surface (w/m^2)
  64 !-- QFX             Upward moisture flux at the surface (Kg/m^2/s)
  65 !-- TSK             Surface temperature (K)
  66 !-- PSFC            Pressure at the surface (Pa)
  67 !-- EP1             Constant for virtual temperature (r_v/r_d-1) (dimensionless)
  68 !-- G               Gravity (m/s^2)
  69 !-- ROVCP           r/cp
  70 !-- RD              gas constant for dry air (j/kg/k)
  71 !-- CPD             heat capacity at constant pressure for dry air (j/kg/k)
  72 !-- GZ1OZ0          log(z/z0) where z0 is roughness length
  73 !-- WSPD            wind speed at lowest model level (m/s)
  74 !-- PSIM            similarity stability function for momentum
  75 !-- MUT             Total Mu : Psfc - Ptop
  76
  77 !-- ids             start index for i in domain
  78 !-- ide             end index for i in domain
  79 !-- jds             start index for j in domain
  80 !-- jde             end index for j in domain
  81 !-- kds             start index for k in domain
  82 !-- kde             end index for k in domain
  83 !-- ims             start index for i in memory
  84 !-- ime             end index for i in memory
  85 !-- jms             start index for j in memory
  86 !-- jme             end index for j in memory
  87 !-- kms             start index for k in memory
  88 !-- kme             end index for k in memory
  89 !-- jts             start index for j in tile
  90 !-- jte             end index for j in tile
  91 !-- kts             start index for k in tile
  92 !-- kte             end index for k in tile
  93 !
  94 !... Outputs:
  95 !-- PBLH            PBL height (m)
  96 !-- KPBL2D          K index for PBL layer
  97 !-- REGIME          Flag indicating PBL regime (stable, unstable, etc.)
  98 !-- RUBLTEN         U tendency due to PBL parameterization (m/s^2)
  99 !-- RVBLTEN         V tendency due to PBL parameterization (m/s^2)
 100 !-- RTHBLTEN        Theta tendency due to PBL parameterization (K/s)
 101 !-- RQVBLTEN        Qv tendency due to PBL parameterization (kg/kg/s)
 102 !-- RQCBLTEN        Qc tendency due to PBL parameterization (kg/kg/s)
 103 !-- RQIBLTEN        Qi tendency due to PBL parameterization (kg/kg/s)
 104 !
 105 !-------------------------------------------------------------------
 106      IMPLICIT NONE
 107
 108 !.......Arguments
 109 ! DECLARATIONS - INTEGER
 110     INTEGER,  INTENT(IN   )   ::      ids,ide, jds,jde, kds,kde, &
 111                                       ims,ime, jms,jme, kms,kme, &
 112                                       its,ite, jts,jte, kts,kte, XTIME
 113
 114 ! DECLARATIONS - REAL
 115     REAL,                                INTENT(IN)  ::  DTPBL, EP1,   &
 116                                                         G, ROVCP, RD, CPD
 117     REAL,    DIMENSION( kms:kme ),       INTENT(IN)  :: ZNW, SIGMAH
 118     REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ),                         &
 119              INTENT(IN) ::                              U3D, V3D,            &
 120                                                         PP3D, DZ8W, T3D,     &
 121                                                         QV3D, QC3D, QI3D,    &
 122                                                         RR3D, TH3D
 123     REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::  UST
 124     REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN) ::                    &
 125                                                         HFX, QFX, TSK,       &
 126                                                         PSFC
 127     INTEGER, DIMENSION( ims:ime, jms:jme ),                                  &
 128              INTENT(OUT  )   ::                         KPBL2D
 129     REAL,    DIMENSION( ims:ime, jms:jme ),                                  &
 130              INTENT(INOUT) ::                           PBLH
 131     REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)   ::  REGIME
 132     REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)  ::                   &
 133                                                         psim,                &
 134                                                         gz1oz0,              &
 135                                                         wspd ,   mut
 136     REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ),                         &
 137              INTENT(INOUT)   ::                         RUBLTEN, RVBLTEN,    &
 138                                                         RTHBLTEN, RQVBLTEN,  &
 139                                                         RQCBLTEN, RQIBLTEN
 140
 141 !... Local Variables
 142
 143 !... Integer
 144       INTEGER :: J, K
 145 !... Real
 146       REAL, DIMENSION( kts:kte ) :: DSIGH, DSIGHI, DSIGFI
 147       REAL, DIMENSION( 0:kte )   :: SIGMAF
 148       REAL  RDT
 149       REAL, PARAMETER :: KARMAN = 0.4
 150 !...
 151
 152    RDT = 1.0 / DTPBL
 153
 154    DO K = 1, kte
 155      SIGMAF(K-1) = ZNW(K)
 156    ENDDO
 157    SIGMAF(kte) = 0.0
 158
 159    DO K = 1, kte
 160      DSIGH(K)  = SIGMAF(K) - SIGMAF(K-1)
 161      DSIGHI(K) = 1.0 / DSIGH(K)
 162    ENDDO
 163
 164    DO K = kts,kte-1
 165      DSIGFI(K) = 1.0 / (SIGMAH(K+1) - SIGMAH(K))
 166    ENDDO
 167
 168    DSIGFI(kte) = DSIGFI(kte-1)
 169
 170     do j = jts,jte
 171       CALL ACM2D(j=J,xtime=XTIME, dtpbl=DTPBL, sigmaf=SIGMAF, sigmah=SIGMAH    &
 172               ,dsigfi=DSIGFI,dsighi=DSIGHI,dsigh=DSIGH             &
 173               ,us=u3d(ims,kms,j),vs=v3d(ims,kms,j)                 &
 174               ,theta=th3d(ims,kms,j),tt=t3d(ims,kms,j)             &
 175               ,qvs=qv3d(ims,kms,j),qcs=qc3d(ims,kms,j)             &
 176               ,qis=qi3d(ims,kms,j),dzf=DZ8W(ims,kms,j)             &
 177               ,densx=RR3D(ims,kms,j)                               &
 178               ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j)     &
 179               ,ttnp=rthblten(ims,kms,j),qvtnp=rqvblten(ims,kms,j)  &
 180               ,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) &
 181               ,cpd=cpd,g=g,rovcp=rovcp,rd=rd,rdt=rdt               &
 182               ,psfcpa=psfc(ims,j),ust=ust(ims,j)                   &
 183               ,pbl=pblh(ims,j)                                     &
 184               ,regime=regime(ims,j),psim=psim(ims,j)               &
 185               ,hfx=hfx(ims,j),qfx=qfx(ims,j)                       &
 186               ,tg=tsk(ims,j),gz1oz0=gz1oz0(ims,j)                  &
 187               ,wspd=wspd(ims,j) ,klpbl=kpbl2d(ims,j)               &
 188               ,mut=mut(ims,j)                                      &
 189               ,ep1=ep1,karman=karman                               &
 190               ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde   &
 191               ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme   &
 192               ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte   )
 193     enddo
 194
 195    END SUBROUTINE ACMPBL
 196
 197 !--------------------------------------------------------------------
 198    SUBROUTINE ACM2D(j,XTIME, DTPBL, sigmaf, sigmah          &
 199               ,dsigfi,dsighi,dsigh                          &
 200               ,us,vs,theta,tt,qvs,qcs,qis                   &
 201               ,dzf,densx,utnp,vtnp,ttnp,qvtnp,qctnp,qitnp   &
 202               ,cpd,g,rovcp,rd,rdt,psfcpa,ust                &
 203               ,pbl,regime,psim                              &
 204               ,hfx,qfx,tg,gz1oz0,wspd ,klpbl                &
 205               ,mut                                          &
 206               ,ep1,karman                                   &
 207               ,ids,ide, jds,jde, kds,kde   &
 208               ,ims,ime, jms,jme, kms,kme   &
 209               ,its,ite, jts,jte, kts,kte   )
 210
 211 !            USE module_wrf_error
 212             IMPLICIT NONE
 213
 214 !.......Arguments
 215
 216 !... Real
 217       REAL, DIMENSION( 0:kte ),             INTENT(IN)  :: SIGMAF
 218       REAL, DIMENSION( kms:kme ),           INTENT(IN)  :: SIGMAH
 219       REAL, DIMENSION( kts:kte ),           INTENT(IN)  :: DSIGH, DSIGHI, DSIGFI
 220       REAL ,                                INTENT(IN)  :: DTPBL, G, RD,ep1,karman,CPD,ROVCP,RDT
 221       REAL , DIMENSION( ims:ime ),          INTENT(INOUT)  :: PBL, UST
 222
 223       REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN)  :: US,VS, THETA, TT,   &
 224                                                            QVS, QCS, QIS, DENSX
 225       REAL,  DIMENSION( ims:ime, kms:kme ), intent(in)  :: DZF
 226       REAL,  DIMENSION( ims:ime, kms:kme ), intent(inout)   ::  utnp, &
 227                                                                 vtnp, &
 228                                                                 ttnp, &
 229                                                                 qvtnp, &
 230                                                                 qctnp, &
 231                                                                 qitnp
 232       real,     dimension( ims:ime ), intent(in   )   ::   psfcpa
 233       real,     dimension( ims:ime ), intent(in   )   ::   tg
 234       real,     dimension( ims:ime ), intent(inout)   ::   regime
 235       real,     dimension( ims:ime ), intent(in)      ::   wspd, psim, gz1oz0
 236       real,     dimension( ims:ime ), intent(in)      ::   hfx, qfx
 237       real,     dimension( ims:ime ), intent(in)      ::   mut
 238
 239 !... Integer
 240       INTEGER, DIMENSION( ims:ime ),       INTENT(OUT):: KLPBL
 241       INTEGER,  INTENT(IN)      ::      XTIME
 242       integer,  intent(in   )   ::      ids,ide, jds,jde, kds,kde, &
 243                                         ims,ime, jms,jme, kms,kme, &
 244                                         its,ite, jts,jte, kts,kte, j
 245 !--------------------------------------------------------------------
 246 !--Local
 247       INTEGER I, K
 248       INTEGER :: KPBLHT
 249       INTEGER, DIMENSION( its:ite ) :: KPBLH, NOCONV
 250
 251 !... Real
 252       REAL    ::  TVCON, WSS, TCONV, TH1, TOG, DTMP, WSSQ
 253       REAL    ::  psix, THV1
 254       REAL, DIMENSION( its:ite )          :: FINT, PSTAR, CPAIR
 255       REAL, DIMENSION( its:ite, kts:kte ) :: THETAV, RIB,             &
 256                                              EDDYZ, UX, VX, THETAX,   &
 257                                              QVX, QCX, QIX, ZA
 258       REAL, DIMENSION( its:ite, 0:kte )   :: ZF
 259       REAL,    DIMENSION( its:ite)           :: WST, TST, QST, USTM, TSTV
 260       REAL,    DIMENSION( its:ite )          :: PBLSIG, MOL
 261       REAL    ::  FINTT, ZMIX, UMIX, VMIX
 262       REAL    ::  TMPFX, TMPVTCON, TMPP, TMPTHS, TMPTH1, TMPVCONV, WS1, DTH
 263         REAL    ::  A,TST12,RL,ZFUNC
 264 !    REAL, PARAMETER :: KARMAN = 0.4
 265
 266 !... Integer
 267       INTEGER :: KL, jtf, ktf, itf, KMIX
 268 !...
 269         character*512 :: message
 270 !-----initialize vertical tendencies and
 271 !
 272       do i = its,ite
 273         do k = kts,kte
 274           utnp(i,k) = 0.
 275           vtnp(i,k) = 0.
 276           ttnp(i,k) = 0.
 277         enddo
 278       enddo
 279 !
 280       do k = kts,kte
 281         do i = its,ite
 282           qvtnp(i,k) = 0.
 283         enddo
 284       enddo
 285 !
 286       do k = kts,kte
 287         do i = its,ite
 288           qctnp(i,k) = 0.
 289           qitnp(i,k) = 0.
 290         enddo
 291       enddo
 292 !
 293 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 294 !!!  Compute Micromet Scaling variables, not availiable in WRF for ACM
 295      DO I = its,ite
 296            cpair(i)     = CPD * (1.0 + 0.84 * QVS(I,1))           ! J/(K KG)
 297            TMPFX        = HFX(I)  / (cpair(i) * DENSX(I,1))
 298            TMPVTCON     = 1.0 + EP1 * QVS(I,1)              ! COnversion factor for virtual temperature
 299            WS1          = SQRT(US(I,1)**2 + VS(I,1)**2)     ! Level 1 wind speed
 300            TST(I)     = -TMPFX / UST(I)
 301            QST(I)     = -QFX(I) / (UST(I)*DENSX(I,1))
 302            USTM(I)    = UST(I) * WS1 / wspd(i)
 303            THV1         = TMPVTCON * THETA(I,1)
 304            TSTV(I)      = TST(I)*TMPVTCON + THV1*EP1*QST(I)
 305            IF(ABS(TSTV(I)).LT.1.0E-6) THEN
 306              TSTV(I) = SIGN(1.0E-6,TSTV(I))
 307            ENDIF
 308            MOL(I) = THV1 * UST(i)**2/(KARMAN*G*TSTV(I))
 309            WST(I) = UST(I) * (PBL(I)/(KARMAN*ABS(MOL(I)))) ** 0.333333
 310            PSTAR(I)   =  MUT(I)/1000.                        ! P* in cb
 311      ENDDO
 312 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 313
 314 !... Compute PBL height
 315
 316 !... compute the height of full- and half-sigma level above ground level
 317      DO I = its,ite
 318        ZF(I,0)    = 0.0
 319        KLPBL(I)   = 1
 320      ENDDO
 321
 322      DO K = kts, kte
 323        DO I = its,ite
 324          ZF(I,K) = DZF(I,K) + ZF(I,K-1)
 325          ZA(I,K) = 0.5 * (ZF(I,K) + ZF(I,K-1))
 326        ENDDO
 327      ENDDO
 328
 329      DO K = kts, kte
 330        DO I = its,ite
 331          TVCON       = 1.0 + EP1 * QVS(I,K)
 332          THETAV(I,K) = THETA(I,K) * TVCON
 333        ENDDO
 334      ENDDO
 335
 336
 337 !...  COMPUTE PBL WHERE RICHARDSON NUMBER = RIC (0.25) HOLTSLAG ET AL 1990
 338      DO 100 I = its,ite
 339        if(MOL(I).LT.0.0 .AND. XTIME.GT.1) then
 340          WSS   = (UST(I) ** 3 + 0.6 * WST(I) ** 3) ** 0.33333
 341          TCONV = -8.5 * UST(I) * TSTV(I) / WSS
 342          TH1   = THETAV(I,1) + TCONV
 343        else
 344          TH1   = THETAV(I,1)
 345        endif
 346
 347 99     KMIX = 1
 348        DO K = 1,kte
 349          DTMP   = THETAV(I,K) - TH1
 350          IF (DTMP.LT.0.0) KMIX = K
 351        ENDDO
 352        IF(KMIX.GT.1) THEN
 353          FINTT = (TH1 - THETAV(I,KMIX)) / (THETAV(I,KMIX+1)               &
 354                - THETAV(I,KMIX))
 355          ZMIX = FINTT * (ZA(I,KMIX+1)-ZA(I,KMIX)) + ZA(I,KMIX)
 356          UMIX = FINTT * (US(I,KMIX+1)-US(I,KMIX)) + US(I,KMIX)
 357          VMIX = FINTT * (VS(I,KMIX+1)-VS(I,KMIX)) + VS(I,KMIX)
 358        ELSE
 359          ZMIX = ZA(I,1)
 360          UMIX = US(I,1)
 361          VMIX = VS(I,1)
 362        ENDIF
 363        DO K = KMIX,kte
 364          DTMP   = THETAV(I,K) - TH1
 365          TOG = 0.5 * (THETAV(I,K) + TH1) / G
 366          WSSQ = (US(I,K)-UMIX)**2                                     &
 367               + (VS(I,K)-VMIX)**2
 368          IF (KMIX == 1) WSSQ = WSSQ + 100.*UST(I)*UST(I)
 369          WSSQ = MAX( WSSQ, 0.1 )
 370          RIB(I,K) = ABS(ZA(I,K)-ZMIX) * DTMP / (TOG * WSSQ)
 371          IF (RIB(I,K) .GE. RIC) GO TO 201
 372        ENDDO
 373
 374        write (message, *)' RIBX never exceeds RIC, RIB(i,kte) = ',rib(i,5),        &
 375                ' THETAV(i,1) = ',thetav(i,1),' MOL=',mol(i),            &
 376                ' TCONV = ',TCONV,' WST = ',WST(I),                      &
 377                ' KMIX = ',kmix,' UST = ',UST(I),                       &
 378                ' TST = ',TST(I),' U,V = ',US(I,1),VS(I,1),              &
 379                ' I,J=',I,J
 380        CALL wrf_error_fatal ( message )
 381 201    CONTINUE
 382
 383        KPBLH(I) = K
 384
 385 100  CONTINUE
 386
 387      DO I = its,ite
 388        IF (KPBLH(I) .NE. 1) THEN
 389 !---------INTERPOLATE BETWEEN LEVELS -- jp 7/93
 390          FINT(I) = (RIC - RIB(I,KPBLH(I)-1)) / (RIB(I,KPBLH(I)) -       &
 391                     RIB(I,KPBLH(I)-1))
 392          IF (FINT(I) .GT. 0.5) THEN
 393            KPBLHT  = KPBLH(I)
 394            FINT(I) = FINT(I) - 0.5
 395          ELSE
 396            KPBLHT  = KPBLH(I) - 1
 397            FINT(I) = FINT(I) + 0.5
 398          ENDIF
 399          PBL(I)  = FINT(I) * (ZF(I,KPBLHT) - ZF(I,KPBLHT-1)) +          &
 400                      ZF(I,KPBLHT-1)
 401          KLPBL(I) = KPBLHT
 402          PBLSIG(I)   = FINT(I) * DSIGH(KPBLHT) + SIGMAF(KPBLHT-1)    ! sigma at PBL height
 403        ELSE
 404          KLPBL(I) = 1
 405          PBL(I)    = ZF(I,1)
 406          PBLSIG(I)   = SIGMAF(1)
 407        ENDIF
 408
 409      ENDDO
 410
 411      DO I = its,ite
 412        NOCONV(I) = 0
 413
 414 ! Check for CBL and identify conv. vs. non-conv cells
 415        IF (PBL(I) / MOL(I) .LT. -0.02 .AND. KLPBL(I) .GT. 3        &
 416            .AND. THETAV(I,1) .GT. THETAV(I,2) .AND. XTIME .GT. 1) THEN
 417           NOCONV(I)   = 1
 418           REGIME(I) = 4.0                     ! FREE CONVECTIVE - ACM
 419        ENDIF
 420      ENDDO
 421
 422 !... Calculate Kz
 423      CALL EDDYX(DTPBL, ZF,  ZA,     MOL, PBL,  UST,                &
 424                 US,    VS,  TT,  THETAV, DENSX, PSTAR,              &
 425                 QVS,   QCS, QIS, DSIGFI, G, RD, CPAIR,              &
 426                 EDDYZ, its,ite, kts,kte,ims,ime, kms,kme)
 427
 428      CALL ACM(DTPBL, PSTAR,  NOCONV, SIGMAF, DSIGH, DSIGHI, J,      &
 429                  KLPBL, PBL,   PBLSIG, MOL,  UST,                  &
 430                  TST, QST,  USTM,   EDDYZ,  DENSX,                  &
 431                  US,    VS,     THETA,  QVS,    QCS,    QIS,        &
 432                  UX,    VX,     THETAX, QVX,    QCX,    QIX,        &
 433                  ids,ide, jds,jde, kds,kde,                         &
 434                  ims,ime, jms,jme, kms,kme,                         &
 435                  its,ite, jts,jte, kts,kte)
 436
 437 !... Calculate tendency due to PBL parameterization
 438
 439      DO K = kts, kte
 440        DO I = its, ite
 441          UTNP(I,K)  = UTNP(I,K) + (UX(I,K) - US(I,K)) * RDT
 442          VTNP(I,K)  = VTNP(I,K) + (VX(I,K) - VS(I,K)) * RDT
 443          TTNP(I,K)  = TTNP(I,K) + (THETAX(I,K) - THETA(I,K)) * RDT
 444          QVTNP(I,K) = QVTNP(I,K) + (QVX(I,K) - QVS(I,K)) * RDT
 445          QCTNP(I,K) = QCTNP(I,K) + (QCX(I,K) - QCS(I,K)) * RDT
 446          QITNP(I,K) = QITNP(I,K) + (QIX(I,K) - QIS(I,K)) * RDT
 447        ENDDO
 448      ENDDO
 449 !     IF(J.EQ.36)PRINT *,' PBL,THETA,THETAX,QVS,QVX=',PBL(20),THETA(20,1),THETAX(20,1),QVS(20,1),QVX(20,1)
 450 !     IF(J.EQ.36)PRINT *,' UST,ustm,TG,TST=',UST(20),ustm(20),Tg(20),TST(20)
 451 !     IF(J.EQ.36)PRINT *,' HFX,MOL,WST,WSPD=',HFX(20),MOL(20),WST(20),wspd(20)
 452 !     IF(J.EQ.36)then
 453 !       i=20
 454 !       do k=1,kte
 455 !         PRINT *,' qvten,uten,vten=',QVTNP(I,K),UTNP(I,K),VTNP(I,K)
 456 !         PRINT *,' k,thten,th,thx,edy=',k,TTNP(I,K),THETA(20,k),THETAX(20,k),EDDYZ(20,K)
 457 !       enddo
 458 !     ENDIF
 459    END SUBROUTINE ACM2D
 460
 461 !-------------------------------------------------------------------
 462    SUBROUTINE ACMINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,           &
 463                       RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR,       &
 464                       restart, allowed_to_read ,                   &
 465                       ids, ide, jds, jde, kds, kde,                &
 466                       ims, ime, jms, jme, kms, kme,                &
 467                       its, ite, jts, jte, kts, kte                 )
 468 !**********************************************************************
 469 !
 470 !    This subroutine is for preparing ACM PBL variables.
 471 !    Called from module_physics_init.F
 472 !
 473 !  REVISION HISTORY:
 474 !     AX     3/2005 - Originally developed
 475 !**********************************************************************
 476 !   ARGUMENT LIST:
 477 !
 478 !-------------------------------------------------------------------
 479      IMPLICIT NONE
 480 !
 481    LOGICAL , INTENT(IN)          :: restart , allowed_to_read
 482
 483    INTEGER , INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
 484                                      ims, ime, jms, jme, kms, kme, &
 485                                      its, ite, jts, jte, kts, kte
 486
 487    INTEGER , INTENT(IN)          ::  P_QI,P_FIRST_SCALAR
 488
 489 !   REAL , DIMENSION( kms:kme ), INTENT(IN)  :: SHALF
 490    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::         &
 491                                                          RUBLTEN, &
 492                                                          RVBLTEN, &
 493                                                          RTHBLTEN, &
 494                                                          RQVBLTEN, &
 495                                                          RQCBLTEN, &
 496                                                          RQIBLTEN
 497
 498 !... Local Variables
 499    INTEGER :: i, j, k, itf, jtf, ktf
 500
 501 !
 502    jtf=min0(jte,jde-1)
 503    ktf=min0(kte,kde-1)
 504    itf=min0(ite,ide-1)
 505
 506    IF(.not.restart)THEN
 507      DO j=jts,jtf
 508      DO k=kts,ktf
 509      DO i=its,itf
 510         RUBLTEN(i,k,j)=0.
 511         RVBLTEN(i,k,j)=0.
 512         RTHBLTEN(i,k,j)=0.
 513         RQVBLTEN(i,k,j)=0.
 514         RQCBLTEN(i,k,j)=0.
 515      ENDDO
 516      ENDDO
 517      ENDDO
 518    ENDIF
 519
 520    IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
 521       DO j=jts,jtf
 522       DO k=kts,ktf
 523       DO i=its,itf
 524          RQIBLTEN(i,k,j)=0.
 525       ENDDO
 526       ENDDO
 527       ENDDO
 528    ENDIF
 529
 530
 531    END SUBROUTINE acminit
 532
 533 !-------------------------------------------------------------------
 534    SUBROUTINE EDDYX(DTPBL, ZF,  ZA,     MOL, PBL,  UST,               &
 535                     US,    VS,  TT,  THETAV, DENSX, PSTAR,            &
 536                     QVS,   QCS, QIS, DSIGFI, G, RD, CPAIR,            &
 537                     EDDYZ, its,ite, kts,kte,ims,ime,kms,kme )
 538
 539
 540 !**********************************************************************
 541 !   Two methods for computing Kz:
 542 !   1.  Boundary scaling similar to Holtslag and Boville (1993)
 543 !   2.  Local Kz computed as function of local Richardson # and vertical
 544 !       wind shear, similar to LIU & CARROLL (1996)
 545 !
 546 !**********************************************************************
 547 !
 548 !-- DTPBL           time step of the minor loop for the land-surface/pbl model
 549 !-- ZF              height of full sigma level
 550 !-- ZA              height of half sigma level
 551 !-- MOL             Monin-Obukhov length in 1D form
 552 !-- PBL             PBL height in 1D form
 553 !-- UST             friction velocity U* in 1D form (m/s)
 554 !-- US              U wind
 555 !-- VS              V wind
 556 !-- TT              temperature
 557 !-- THETAV          potential virtual temperature
 558 !-- DENSX           dry air density (kg/m^3)
 559 !-- PSTAR           P*=Psfc-Ptop
 560 !-- QVS             water vapor mixing ratio (Kg/Kg)
 561 !-- QCS             cloud mixing ratio (Kg/Kg)
 562 !-- QIS             ice mixing ratio (Kg/Kg)
 563 !-- DSIGFI          inverse of sigma layer delta
 564 !-- G               gravity
 565 !-- RD              gas constant for dry air (j/kg/k)
 566 !-- CPAIR           specific heat of moist air (M^2 S^-2 K^-1)
 567 !-- EDDYZ           eddy diffusivity KZ
 568 !-----------------------------------------------------------------------
 569
 570       IMPLICIT NONE
 571
 572 !.......Arguments
 573
 574 !... Integer
 575       INTEGER,  INTENT(IN)   ::    its,ite, kts,kte,ims,ime,kms,kme
 576 !... Real
 577       REAL , DIMENSION( ims:ime ),          INTENT(IN)  :: PBL, UST
 578       REAL ,                                INTENT(IN)  :: DTPBL, G, RD
 579       REAL , DIMENSION( kts:kte ),          INTENT(IN)  :: DSIGFI
 580       REAL , DIMENSION( its:ite ),          INTENT(IN)  :: MOL, PSTAR, CPAIR
 581
 582       REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN)  :: US,VS, TT,   &
 583                                                            QVS, QCS, QIS, DENSX
 584       REAL, DIMENSION( its:ite, kts:kte ), INTENT(IN) :: ZA, THETAV
 585       REAL, DIMENSION( its:ite, 0:kte )  , INTENT(IN) :: ZF
 586
 587       REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: EDDYZ
 588
 589 !.......Local variables
 590
 591 !... Integer
 592       INTEGER  :: ILX, KL, KLM, K, I
 593
 594 !... Real
 595       REAL     :: ZOVL, PHIH, WT, ZSOL, ZFUNC, DZF, SS, GOTH, EDYZ
 596       REAL     :: RI, QMEAN, TMEAN, XLV, ALPH, CHI, ZK, SQL, DENSF,kzo
 597
 598 !... Parameters
 599       REAL, PARAMETER :: RV     = 461.5
 600       REAL, PARAMETER :: RC     = 0.25
 601       REAL, PARAMETER :: RLAM   = 80.0
 602       REAL, PARAMETER :: GAMH   = 16.0 !15.0  !  Holtslag and Boville (1993)
 603       REAL, PARAMETER :: BETAH  = 5.0   !  Holtslag and Boville (1993)
 604       REAL, PARAMETER :: KARMAN = 0.4
 605       REAL, PARAMETER :: EDYZ0  = 0.01  ! New Min Kz
 606 !      REAL, PARAMETER :: EDYZ0  = 0.1
 607 !-- IMVDIF      imvdif=1 for moist adiabat vertical diffusion
 608       INTEGER, PARAMETER :: imvdif = 1
 609 !
 610       ILX = ite
 611       KL  = kte
 612       KLM = kte - 1
 613
 614       DO K = kts,KLM
 615         DO I = its,ILX
 616           EDYZ = 0.0
 617           ZOVL = 0.0
 618           DZF  = ZA(I,K+1) - ZA(I,K)
 619 !          kzo = min(0.001 * dzf,EDYZ0)
 620           kzo = 0.001 * dzf
 621 !          kzo = EDYZ0
 622 !-------------------------------------------------
 623           IF (ZF(I,K) .LT. PBL(I)) THEN
 624             ZOVL = ZF(I,K) / MOL(I)
 625             IF (ZOVL .LT. 0.0) THEN
 626               IF (ZF(I,K) .LT. 0.1 * PBL(I)) THEN
 627                 PHIH = 1.0 / SQRT(1.0 - GAMH * ZOVL)
 628                 WT   = UST(I) / PHIH
 629               ELSE
 630                 ZSOL = 0.1 * PBL(I) / MOL(I)
 631                 PHIH = 1.0 / SQRT(1.0 - GAMH * ZSOL)
 632                 WT   = UST(I) / PHIH
 633               ENDIF
 634             ELSE IF (ZOVL .LT. 1.0) THEN
 635               PHIH = 1.0 + BETAH * ZOVL
 636               WT   = UST(I) / PHIH
 637             ELSE
 638               PHIH = BETAH + ZOVL
 639               WT   = UST(I) / PHIH
 640             ENDIF
 641             ZFUNC      = ZF(I,K) * (1.0 - ZF(I,K) / PBL(I)) ** 2
 642             EDYZ = KARMAN * WT * ZFUNC
 643 !            EDYZ = AMAX1(EDYZ,EDYZ0)
 644          ENDIF
 645  !            PRINT *,I,K,TT(I,K)
 646  !--------------------------------------------------------------------------
 647 !
 648 !          RC = 0.257 * DZF ** 0.175
 649             SS   = ((US(I,K+1) - US(I,K)) ** 2 + (VS(I,K+1) - VS(I,K)) ** 2)   &
 650                   / (DZF * DZF) + 1.0E-9
 651             GOTH = 2.0 * G / (THETAV(I,K+1) + THETAV(I,K))
 652             RI   = GOTH * (THETAV(I,K+1) - THETAV(I,K)) / (DZF * SS)
 653 !            PRINT *,I,K,TT(I,K), RI, IMVDIF
 654 !
 655 !-- Adjustment to vert diff in Moist air
 656           if(imvdif.eq.1)then
 657             IF ((QCS(I,K)+QIS(I,K)) .GT. 0.01E-3 .OR. (QCS(I,K+1)+             &
 658                  QIS(I,K+1)) .GT. 0.01E-3) THEN
 659               QMEAN = 0.5 * (QVS(I,K) + QVS(I,K+1))
 660               TMEAN = 0.5 * (TT(I,K) + TT(I,K+1))
 661               XLV   = (2.501 - 0.00237 * (TMEAN - 273.15)) * 1.E6
 662               ALPH  =  XLV * QMEAN / RD / TMEAN
 663               CHI   =  XLV * XLV * QMEAN / CPAIR(I) / RV / TMEAN / TMEAN
 664               RI    = (1.0 + ALPH) * (RI -G * G / SS / TMEAN / CPAIR(I) *       &
 665                       ((CHI - ALPH) / (1.0 + CHI)))
 666             ENDIF
 667           endif
 668
 669             ZK  = 0.4 * ZF(I,K)
 670             SQL = (ZK * RLAM / (RLAM + ZK)) ** 2
 671
 672             IF (RI .GE. RC) THEN
 673               EDDYZ(I,K) = kzo !EDYZ0
 674             ELSE IF (RI .GE. 0.0) THEN
 675               EDDYZ(I,K) = kzo + SQRT(SS) * (1.- RI / RC) ** 2 * SQL
 676             ELSE
 677               EDDYZ(I,K) = kzo + SQRT(SS * (1.0 - 25.0 * RI)) * SQL
 678             ENDIF
 679
 680           IF(ZOVL.LT.0.0.AND.EDYZ.GT.EDDYZ(I,K)) THEN
 681             EDDYZ(I,K) = EDYZ
 682           ENDIF
 683
 684           EDDYZ(I,K) = MIN(1000.0,EDDYZ(I,K))
 685           EDDYZ(I,K) = MAX(kzo,EDDYZ(I,K))
 686
 687           DENSF     = 0.5 * (DENSX(I,K+1) + DENSX(I,K))
 688
 689           EDDYZ(I,K) = EDDYZ(I,K) * (DENSF * G / PSTAR(I)) ** 2 *       &
 690                        DTPBL * DSIGFI(K)*1E-6
 691         ENDDO             ! for I loop
 692       ENDDO               ! for k loop
 693 !
 694       DO I = its,ILX
 695         EDDYZ(I,KL) = 0.0 ! EDDYZ(I,KLM) -- changed jp 3/08
 696       ENDDO
 697    END SUBROUTINE EDDYX
 698
 699 !-------------------------------------------------------------------
 700    SUBROUTINE ACM (DTPBL, PSTAR,  NOCONV, SIGMAF, DSIGH, DSIGHI, JX, &
 701                    KLPBL, PBL,   PBLSIG, MOL,  UST,                  &
 702                    TST, QST,  USTM,   EDDYZ,  DENSX,               &
 703                    US,    VS,     THETA,  QVS,    QCS,    QIS,     &
 704                    UX,    VX,     THETAX, QVX,    QCX,    QIX,     &
 705                    ids,ide, jds,jde, kds,kde,                      &
 706                    ims,ime, jms,jme, kms,kme,                      &
 707                    its,ite, jts,jte, kts,kte)
 708 !**********************************************************************
 709 !   PBL model called the Asymmetric Convective Model, Version 2 (ACM2)
 710 !   -- See top of module for summary and references
 711 !
 712 !---- REVISION HISTORY:
 713 !   AX     3/2005 - developed WRF version based on ACM2 in the MM5 PX LSM
 714 !   JP and RG 8/2006 - updates
 715 !
 716 !**********************************************************************
 717 !  ARGUMENTS:
 718 !-- DTPBL           PBL time step
 719 !-- PSTAR           Psurf - Ptop in cb
 720 !-- NOCONV          If free convection =0, no; =1, yes
 721 !-- SIGMAF          Sigma for full layer
 722 !-- DSIGH           Sigma thickness
 723 !-- DSIGHI          Inverse of sigma thickness
 724 !-- JX              N-S index
 725 !-- KLPBL           PBL level at K index
 726 !-- PBL             PBL height in m
 727 !-- PBLSIG          Sigma level for PBL
 728 !-- MOL             Monin-Obukhov length in 1D form
 729 !-- UST             U* in 1D form
 730 !-- TST             Theta* in 1D form
 731 !-- QST             Q* in 1D form
 732 !-- USTM            U* for computation of momemtum flux
 733 !-- EDDYZ           eddy diffusivity KZ
 734 !-- DENSX           dry air density (kg/m^3)
 735 !-- US              U wind
 736 !-- VS              V wind
 737 !-- THETA           potential temperature
 738 !-- QVS             water vapor mixing ratio (Kg/Kg)
 739 !-- QCS             cloud mixing ratio (Kg/Kg)
 740 !-- QIS             ice mixing ratio (Kg/Kg)
 741 !-- UX              new U wind
 742 !-- VX              new V wind
 743 !-- THETAX          new potential temperature
 744 !-- QVX             new water vapor mixing ratio (Kg/Kg)
 745 !-- QCX             new cloud mixing ratio (Kg/Kg)
 746 !-- QIX             new ice mixing ratio (Kg/Kg)
 747 !-----------------------------------------------------------------------
 748
 749       IMPLICIT NONE
 750
 751 !.......Arguments
 752
 753 !... Integer
 754       INTEGER,  INTENT(IN)      ::      ids,ide, jds,jde, kds,kde, &
 755                                         ims,ime, jms,jme, kms,kme, &
 756                                         its,ite, jts,jte, kts,kte, JX
 757       INTEGER,  DIMENSION( its:ite ), INTENT(IN)  :: NOCONV
 758       INTEGER,  DIMENSION( ims:ime ), INTENT(IN)  :: KLPBL
 759
 760 !... Real
 761       REAL , DIMENSION( ims:ime ),          INTENT(IN)  :: PBL, UST
 762       REAL ,                                INTENT(IN)  :: DTPBL
 763       REAL , DIMENSION( its:ite ),          INTENT(IN)  :: PSTAR, PBLSIG,  &
 764                                                            MOL, TST, &
 765                                                            QST, USTM
 766       REAL , DIMENSION( kts:kte ),          INTENT(IN)  :: DSIGHI, DSIGH
 767       REAL , DIMENSION( 0:kte ),            INTENT(IN)  :: SIGMAF
 768       REAL , DIMENSION( its:ite, kts:kte ), INTENT(INOUT)  :: EDDYZ
 769       REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN)  :: US,VS, THETA,   &
 770                                                            QVS, QCS, QIS, DENSX
 771       REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: UX, VX, THETAX,      &
 772                                                            QVX, QCX, QIX
 773
 774 !.......Local variables
 775
 776 !... Parameters
 777       INTEGER, PARAMETER :: NSP   = 6
 778 !
 779 !......ACM2 Parameters
 780 !     INTEGER, PARAMETER :: IFACM = 0
 781 !
 782       REAL,    PARAMETER :: G1000 = 9.8 * 1.0E-3
 783       REAL,    PARAMETER :: XX    = 0.5          ! FACTOR APPLIED TO CONV MIXING TIME STEP
 784       REAL,    PARAMETER :: KARMAN = 0.4
 785
 786 !... Integer
 787       INTEGER :: ILX, KL, KLM, I, K, NSPX, NLP, NL, JJ, L
 788       INTEGER :: KCBLMX
 789       INTEGER, DIMENSION( its:ite ) :: KCBL
 790
 791 !... Real
 792       REAL                               :: G1000I, MBMAX, HOVL, MEDDY, MBAR
 793       REAL                               :: EKZ, RZ, FM, WSPD, DTS, DTRAT, F1
 794       REAL, DIMENSION( its:ite )         :: PSTARI, FSACM, RAH, DTLIM
 795       REAL, DIMENSION( kts:kte, its:ite) :: MBARKS, MDWN
 796       REAL, DIMENSION( 1:NSP, its:ite )  :: FS, BCBOTN
 797       REAL, DIMENSION( kts:kte )         :: XPLUS, XMINUS
 798       REAL  DELC
 799       REAL, DIMENSION( 1:NSP,its:ite,kts:kte  ) :: VCI
 800
 801       REAL, DIMENSION( kts:kte )               :: AI, BI, CI, EI !, Y
 802       REAL, DIMENSION( 1:NSP,kts:kte )         :: DI, UI
 803 !
 804 !--Start Exicutable ----
 805
 806       ILX = ite
 807       KL  = kte
 808       KLM = kte - 1
 809
 810       G1000I = 1.0 / G1000
 811       KCBLMX = 0
 812       MBMAX  = 0.0
 813
 814 !---COMPUTE ACM MIXING RATE
 815       DO I = its, ILX
 816         DTLIM(I)  = DTPBL
 817         PSTARI(I) = 1.0 / PSTAR(I)
 818         KCBL(I)   = 1
 819         FSACM(I)  = 0.0
 820
 821      IF (NOCONV(I) .EQ. 1) THEN
 822           KCBL(I) = KLPBL(I)
 823 !-------MBARKS IS UPWARD MIXING RATE; MDWN IS DOWNWARD MIXING RATE
 824 !--New couple ACM & EDDY-------------------------------------------------------------
 825           HOVL     = -PBL(I) / MOL(I)
 826           FSACM(I) = 1./(1.+((KARMAN/(HOVL))**0.3333)/(0.72*KARMAN))
 827           MEDDY    = EDDYZ(I,1) / (DTPBL * (PBLSIG(I) - SIGMAF(1)))
 828           MBAR     = MEDDY * FSACM(I)
 829           DO K = kts,KCBL(I)-1
 830             EDDYZ(I,K) = EDDYZ(I,K) * (1.0 - FSACM(I))
 831           ENDDO
 832 !          if(i.eq.100.and.jx.eq.43) PRINT *,' Edy,meddy,mbar=', EDDYZ(I,1),MEDDY,MBAR
 833           MBMAX = AMAX1(MBMAX,MBAR)
 834           DO K = kts+1,KCBL(I)
 835             MBARKS(K,I) = MBAR
 836             MDWN(K,I)   = MBAR * (PBLSIG(I) - SIGMAF(K-1)) * DSIGHI(K)
 837           ENDDO
 838           MBARKS(1,I) = MBAR
 839           MBARKS(KCBL(I),I) = MDWN(KCBL(I),I)
 840           MDWN(KCBL(I)+1,I) = 0.0
 841         ENDIF
 842       ENDDO                              ! end of I loop
 843
 844 !      if(NOCONV(20).EQ.1) print *,' KCBL,PBLSIG=',KCBL(20),PBLSIG(20),SIGMAF(KCBL(20)-1),DSIGHI(KCBL(20))
 845 !      do k = kts,klm
 846 !        if(NOCONV(20).EQ.1.and.MBAR.lt.1.0e-3)print *,' k,MBARKS,MDWN=',k,MBARKS(k,20),MDWN(k,20)
 847 !      enddo
 848
 849        DO K = kts,KLM
 850          DO I = its,ILX
 851           EKZ   = EDDYZ(I,K) / DTPBL * DSIGHI(K)
 852           DTLIM(I) = AMIN1(0.75 / EKZ,DTLIM(I))
 853          ENDDO
 854        ENDDO
 855
 856       DO I = its,ILX
 857         IF (NOCONV(I) .EQ. 1) THEN
 858           KCBLMX = AMAX0(KLPBL(I),KCBLMX)
 859           RZ     = (SIGMAF(KCBL(I)) - SIGMAF(1)) * DSIGHI(1)
 860           DTLIM(I)  = AMIN1(XX / (MBARKS(1,I) * RZ),DTLIM(I))
 861         ENDIF
 862       ENDDO
 863
 864       DO I = its,ILX
 865       ENDDO
 866 !
 867       DO K = kts,KL
 868         DO I = its,ILX
 869           VCI(1,I,K) = THETA(I,K)
 870           VCI(2,I,K) = QVS(I,K)
 871           VCI(3,I,K) = US(I,K)
 872           VCI(4,I,K) = VS(I,K)
 873           ! -- Also mix cloud water and ice if necessary
 874           ! IF (IMOISTX.NE.1.AND.IMOISTX.NE.3) THEN  !!! Check other PBL models
 875           VCI(5,I,K) = QCS(I,K)
 876           VCI(6,I,K) = QIS(I,K)
 877         ENDDO
 878       ENDDO
 879
 880       NSPX=6
 881
 882       DO I = its,ILX
 883 !       RAH(I)  = RA(I,J) + 3.976 / UST(I)
 884         FS(1,I) = -UST(I) * TST(I) * DENSX(I,1) * PSTARI(I)
 885         FS(2,I) = -UST(I) * QST(I) * DENSX(I,1) * PSTARI(I)
 886         FM      = -USTM(I) * USTM(I) * DENSX(I,1) * PSTARI(I)
 887         WSPD    = SQRT(US(I,1) * US(I,1) + VS(I,1) * VS(I,1)) + 1.E-9
 888         FS(3,I) = FM * US(I,1) / WSPD
 889         FS(4,I) = FM * VS(I,1) / WSPD
 890         FS(5,I) = 0.0
 891         FS(6,I) = 0.0                      ! SURFACE FLUXES OF CLOUD WATER AND ICE = 0
 892       ENDDO
 893 !
 894 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 895       DO I = its,ILX
 896 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 897         NLP   = INT(DTPBL / DTLIM(I) + 1.0)
 898         DTS   = (DTPBL / NLP)
 899         DTRAT = DTS / DTPBL
 900         DO NL = 1,NLP           ! LOOP OVER SUB TIME LOOP
 901 !
 902 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 903 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 904
 905 !-- COMPUTE ARRAY ELEMENTS THAT ARE INDEPENDANT OF SPECIES
 906
 907           DO K = kts,KL
 908             AI(K) = 0.0
 909             BI(K) = 0.0
 910             CI(K) = 0.0
 911             EI(K) = 0.0
 912           ENDDO
 913
 914           DO K = 2, KCBL(I)
 915             EI(K-1) = -CRANKP * MDWN(K,I) * DTS * DSIGH(K) * DSIGHI(K-1)
 916             BI(K)   = 1.0 + CRANKP * MDWN(K,I) * DTS
 917             AI(K)   = -CRANKP * MBARKS(K,I) * DTS
 918           ENDDO
 919
 920           EI(1) = EI(1) -EDDYZ(I,1) * CRANKP * DSIGHI(1 )* DTRAT
 921           AI(2) = AI(2) -EDDYZ(I,1) * CRANKP * DSIGHI(2) * DTRAT
 922
 923           DO K =  KCBL(I)+1, KL
 924             BI(K) = 1.0
 925           ENDDO
 926
 927           DO K = 2,KL
 928             XPLUS(K)  = EDDYZ(I,K) * DSIGHI(K) * DTRAT
 929             XMINUS(K) = EDDYZ(I,K-1) * DSIGHI(K) * DTRAT
 930             CI(K)     = - XMINUS(K) * CRANKP
 931             EI(K)     = EI(K) - XPLUS(K) * CRANKP
 932             BI(K)     = BI(K) + XPLUS(K) * CRANKP + XMINUS(K) * CRANKP
 933           ENDDO
 934
 935           IF (NOCONV(I) .EQ. 1) THEN
 936             BI(1) = 1.0 + CRANKP * MBARKS(1,I) * (PBLSIG(I) - SIGMAF(1)) *    &
 937                     DTS * DSIGHI(1) + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
 938       ELSE
 939         BI(1) = 1.0  + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
 940       ENDIF
 941
 942
 943       DO K = 1,KL
 944         DO L = 1,NSPX
 945            DI(L,K) = 0.0
 946         ENDDO
 947      ENDDO
 948 !
 949 !**   COMPUTE TENDENCY OF CBL CONCENTRATIONS - SEMI-IMPLICIT SOLUTION
 950      DO K = 2,KCBL(I)
 951        DO L = 1,NSPX
 952           DELC = DTS * (MBARKS(K,I) * VCI(L,I,1) - MDWN(K,I) *          &
 953                  VCI(L,I,K) + DSIGH(K+1) * DSIGHI(K) *                  &
 954                         MDWN(K+1,I) * VCI(L,I,K+1))
 955           DI(L,K)   = VCI(L,I,K) + (1.0 - CRANKP) * DELC
 956        ENDDO
 957      ENDDO
 958
 959      DO K = KCBL(I)+1, KL
 960        DO L = 1,NSPX
 961          DI(L,K) = VCI(L,I,K)
 962        ENDDO
 963      ENDDO
 964
 965      DO K = 2,KL
 966        IF (K .EQ. KL) THEN
 967          DO L = 1,NSPX
 968            DI(L,K) = DI(L,K)  - (1.0 - CRANKP) * XMINUS(K) *                  &
 969                      (VCI(L,I,K) - VCI(L,I,K-1))
 970          ENDDO
 971        ELSE
 972          DO L = 1,NSPX
 973            DI(L,K) = DI(L,K) + (1.0 - CRANKP) * XPLUS(K) *                   &
 974                       (VCI(L,I,K+1) - VCI(L,I,K))  -                         &
 975                       (1.0 - CRANKP) * XMINUS(K) *                           &
 976                       (VCI(L,I,K) - VCI(L,I,K-1))
 977          ENDDO
 978        ENDIF
 979      ENDDO
 980
 981      IF (NOCONV(I) .EQ. 1) THEN
 982        DO L = 1,NSPX
 983          F1    = -G1000I * (MBARKS(1,I) *                                &
 984                  (PBLSIG(I) - SIGMAF(1)) * VCI(L,I,1) -                  &
 985                   MDWN(2,I) * VCI(L,I,2) * DSIGH(2))
 986
 987          DI(L,1) = VCI(L,I,1) - G1000 * (FS(L,I) - (1.0 - CRANKP)        &
 988                  * F1) * DSIGHI(1) * DTS
 989        ENDDO
 990       ELSE
 991         DO L = 1,NSPX
 992           DI(L,1) = VCI(L,I,1) - G1000 * FS(L,I) * DSIGHI(1) * DTS
 993         ENDDO
 994       ENDIF
 995       DO L = 1,NSPX
 996         DI(L,1) = DI(L,1) + (1.0 - CRANKP) * EDDYZ(I,1) * DSIGHI(1)      &
 997                 * DTRAT * (VCI(L,I,2) - VCI(L,I,1))
 998       ENDDO
 999       IF ( NOCONV(I) .EQ. 1 ) THEN
1000         CALL MATRIX (AI, BI, CI, DI, EI, UI, KL, NSPX)
1001       ELSE
1002         CALL TRI (CI, BI, EI, DI, UI, KL, NSPX)
1003       END IF
1004 !
1005 !-- COMPUTE NEW THETAV AND Q
1006          DO K = 1,KL
1007            DO L = 1,NSPX
1008               VCI(L,I,K) = UI(L,K)
1009            ENDDO
1010          ENDDO
1011 !
1012         ENDDO                   ! END I LOOP
1013       ENDDO                     ! END SUB TIME LOOP
1014 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1015 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1016
1017 !
1018       DO K = kts,KL
1019         DO I = its,ILX
1020           THETAX(I,K) = VCI(1,I,K)
1021           QVX(I,K)    = VCI(2,I,K)
1022           UX(I,K)     = VCI(3,I,K)
1023           VX(I,K)     = VCI(4,I,K)
1024 !            IF( (I.EQ.61 .OR. I.EQ.15) .AND. JX.EQ.22) THEN
1025 !               PRINT *,'TEST -->',I,K,THETAX(I,K)
1026 !            ENDIF
1027         ENDDO
1028       ENDDO
1029
1030       DO K = kts,KL
1031         DO I = its,ILX
1032           QCX(I,K) = VCI(5,I,K)
1033           QIX(I,K) = VCI(6,I,K)
1034         ENDDO
1035       ENDDO
1036 !      ENDIF
1037
1038    END SUBROUTINE ACM
1039 !--------------------------------------------------------
1040
1041    SUBROUTINE MATRIX(A,B,C,D,E,X,KL,NSP)
1042
1043    IMPLICIT NONE
1044 !
1045 !-- Bordered band diagonal matrix solver for ACM2
1046
1047 !-- ACM2 Matrix is in this form:
1048 !   B1 E1
1049 !   A2 B2 E2
1050 !   A3 C3 B3 E3
1051 !   A4    C4 B4 E4
1052 !   A5       C5 B5 E5
1053 !   A6          C6 B6
1054
1055 !--Upper Matrix is
1056 !  U11 U12
1057 !      U22 U23
1058 !          U33 U34
1059 !              U44 U45
1060 !                  U55 U56
1061 !                      U66
1062
1063 !--Lower Matrix is:
1064 !  1
1065 ! L21  1
1066 ! L31 L32  1
1067 ! L41 L42 L43  1
1068 ! L51 L52 L53 L54  1
1069 ! L61 L62 L63 L64 L65 1
1070 !---------------------------------------------------------
1071 !...Arguments
1072       INTEGER, INTENT(IN)   :: KL
1073       INTEGER, INTENT(IN)   :: NSP
1074       REAL A(KL),B(KL),E(KL)
1075       REAL C(KL),D(NSP,KL),X(NSP,KL)
1076
1077 !...Locals
1078       REAL Y(NSP,KL),AIJ,SUM
1079       REAL L(KL,KL),UII(KL),UIIP1(KL),RUII(KL)
1080       INTEGER I,J,V
1081
1082 !-- Define Upper and Lower matrices
1083       L(1,1) = 1.
1084 !      U(1,1) = B(1)
1085       UII(1) = B(1)
1086       RUII(1) = 1./UII(1)
1087       DO I = 2, KL
1088         L(I,I) = 1.
1089         L(I,1) = A(I)/B(1)
1090 !        U(I-1,I) =E(I-1)
1091         UIIP1(I-1)=E(I-1)
1092         IF(I.GE.3) THEN
1093            DO J = 2,I-1
1094              IF(I.EQ.J+1) THEN
1095                 AIJ = C(I)
1096              ELSE
1097                 AIJ = 0.
1098              ENDIF
1099              L(I,J) = (AIJ-L(I,J-1)*E(J-1))/      &
1100                       (B(J)-L(J,J-1)*E(J-1))
1101            ENDDO
1102           ENDIF
1103       ENDDO
1104
1105       DO I = 2,KL
1106 !        U(I,I) = B(I)-L(I,I-1)*E(I-1)
1107         UII(I) = B(I)-L(I,I-1)*E(I-1)
1108         RUII(I) = 1./UII(I)
1109       ENDDO
1110
1111 !-- Forward sub for Ly=d
1112       DO V= 1, NSP
1113         Y(V,1) = D(V,1)
1114         DO I=2,KL
1115           SUM = D(V,I)
1116           DO J=1,I-1
1117             SUM = SUM - L(I,J)*Y(V,J)
1118           ENDDO
1119           Y(V,I) = SUM
1120         ENDDO
1121       ENDDO
1122
1123 !-- Back sub for Ux=y
1124
1125       DO V= 1, NSP
1126 !        X(V,KL) = Y(V,KL)/U(KL,KL)
1127         X(V,KL) = Y(V,KL)*RUII(KL)
1128       ENDDO
1129       DO I = KL-1,1,-1
1130         DO V= 1, NSP
1131 !         X(V,I) = (Y(V,I)-U(I,I+1)*X(V,I+1))/U(I,I)
1132          X(V,I) = (Y(V,I)-UIIP1(I)*X(V,I+1))*RUII(I)
1133         ENDDO
1134       ENDDO
1135
1136    END SUBROUTINE MATRIX
1137 !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1138       SUBROUTINE TRI ( L, D, U, B, X,KL,NSP)
1139 !-----------------------------------------------------------------------
1140
1141 !  FUNCTION:
1142 !    Solves tridiagonal system by Thomas algorithm.
1143 ! The associated tri-diagonal system is stored in 3 arrays
1144 !   D : diagonal
1145 !   L : sub-diagonal
1146 !   U : super-diagonal
1147 !   B : right hand side function
1148 !   X : return solution from tridiagonal solver
1149
1150 !     [ D(1) U(1) 0    0    0 ...       0     ]
1151 !     [ L(2) D(2) U(2) 0    0 ...       .     ]
1152 !     [ 0    L(3) D(3) U(3) 0 ...       .     ]
1153 !     [ .       .     .     .           .     ] X(i) = B(i)
1154 !     [ .             .     .     .     0     ]
1155 !     [ .                   .     .     .     ]
1156 !     [ 0                           L(n) D(n) ]
1157
1158 !-----------------------------------------------------------------------
1159
1160       IMPLICIT NONE
1161
1162 ! Arguments:
1163
1164       INTEGER, INTENT(IN)   :: KL
1165       INTEGER, INTENT(IN)   :: NSP
1166
1167       REAL        L( KL )               ! subdiagonal
1168       REAL        D(KL)   ! diagonal
1169       REAL        U( KL )               ! superdiagonal
1170       REAL        B(NSP,KL )   ! R.H. side
1171       REAL        X( NSP,KL )   ! solution
1172
1173 ! Local Variables:
1174
1175       REAL        GAM( KL )
1176       REAL        BET
1177       INTEGER     V, K
1178
1179 ! Decomposition and forward substitution:
1180       BET = 1.0 / D( 1 )
1181       DO V = 1, NSP
1182          X( V,1 ) = BET * B(V,1 )
1183       ENDDO
1184
1185       DO K = 2, KL
1186          GAM(K ) = BET * U( K-1 )
1187          BET = 1.0 / ( D( K ) - L( K ) * GAM( K ) )
1188          DO V = 1, NSP
1189             X( V, K ) = BET * ( B( V,K ) - L( K ) * X( V,K-1 ) )
1190          ENDDO
1191       END DO
1192
1193 ! Back-substitution:
1194
1195       DO K = KL - 1, 1, -1
1196          DO V = 1, NSP
1197             X( V,K ) = X( V,K ) - GAM( K+1 ) * X( V,K+1 )
1198          END DO
1199       ENDDO
1200
1201   END SUBROUTINE TRI
1202
1203 END MODULE module_bl_acm
1204