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