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