      subroutine vdiffacmx(ns,nddm,nz,deltat,temp0,press0,z0,delz,rkv,
     &                     qv,temp,press,wind,kpbl,fconv,con)
c
c----CAMx v7Beta6 190902
c
c     VDIFFACMX performs vertical diffusion of concentrations using
c     the Asymmetric Convective Model v2 (ACM2/ACM1; Pleim, 2007;2014). 
c
c     ACM2 portion of code is based on: 
c         Models-3/CMAQ vdiffacm2.F,v 1.13 2012/01/19 14:37:47 
c
c     Copyright 1996 - 2018
c     Ramboll
c
c     Modifications:
c        5/27/16   Updated to parallel ACM2 speed improvements in CMAQ v5.1,
c                  remove K-theory solution and dry deposition (both done using
c                  CAMx K-theory solver), and extend to H/DDM
c
c     Input arguments:
c        ns                number of species
c        nddm              number of parameters for which sensitivities are
c                          calculated - should be zero if DDM is not enabled
c        nz                number of layers
c        deltat            timestep (s)
c        temp0             Surface temperature (K)
c        press0            Surface pressure (mb)
c        z0                Surface roughness (m)
c        delz              Layer depth (m)
c        rkv               Kv profile (m2/s)
c        qv                Humidity profile (mixing ratio) 
c        temp              Temperature profile (K)
c        press             Pressure profile (mb)
c        wind              Wind speed profile
c        con               Species concentrations (umol/m3, ug/m3) followed by
c                          sensitivities (umol/m3/parameter unit).
c
c     Output arguments:
c        kpbl              Layer index at the top of the PBL
c        fconv             Convective fraction to scale Kv
c        con               Species concentrations (umol/m3, ug/m3) followed by
c                          sensitivities (umol/m3/parameter unit).
c
c     Routines Called:
c        MICROMET
c        MATRIX
c
c     Called by:
c        DIFFUS
c
      implicit none
      include "camx.prm"
c
      integer ns,nz,nddm,kpbl
      real deltat,temp0,press0,z0,fconv
      real delz(nz),rkv(nz),qv(nz),temp(nz),press(nz),
     &     wind(nz)
      real con(nz+nz*nddm,ns)
c
      real theta,thbar
      parameter (theta = 0.5)         !Crank-Nicholson = 0.5 (semi-implicit)
      parameter (thbar = 1 - theta)
      logical lconv,lstable
      integer k,kk,nstep,n,l,m
      real gamma,vk
      real tv,pbl,zf,critk,ustar,el,psih,wstar,hoverl,dt,mbar,meddy,
     &     delc,lfac1,lfac2,rhopbl,rhosum,roplus,rz,dtacm,
     &     dfacp,dfacq
      real zz(0:MXLAYER)
      real dz(MXLAYER),dzinv(MXLAYER),zm(MXLAYER),thetav(MXLAYER),
     &     rho(MXLAYER),dzm(MXLAYER),dzminv(MXLAYER),seddy(MXLAYER)
      real mbarks(MXLAYER),mdwn(MXLAYER),mfac(MXLAYER),dfsp(MXLAYER),
     &     dfsq(MXLAYER)
      real aa(MXLAYER),bb(MXLAYER),ee(MXLAYER),dd(MXLAYER),uu(MXLAYER)
      real cc(MXLAYER,1+MXTRSP)
c
      data gamma /0.286/, vk /0.4/
c
c-----Entry point
c
c-----Determine vertical grid structure variables and thermodynamic profiles.
c
      zz(0) = 0.
      do k = 1,nz
        dz(k) = delz(k)
        zz(k) = zz(k-1) + dz(k)
        zm(k)  = (zz(k) + zz(k-1))/2.
        tv = temp(k)*(1. + 0.608*qv(k))
        thetav(k) = tv*(1000./press(k))**gamma
        rho(k) = press(k)/tv
        seddy(k) = rkv(k)
      enddo
      do k = 1,nz-1
        dzm(k) = zm(k+1) - zm(k)
      enddo
      dzm(nz) = 0.
c
c-----Get the PBL parameters
c
      kpbl = 1
      pbl = dz(1)
      do k = 2,nz-1
        if (kpbl .ne. k-1) goto 100
        zf = dz(k)/dz(k-1)
        critk = 0.03*dz(k-1)*dzm(k)*(1. + zf)/200.
        if (rkv(k-1) .gt. critk) then
          kpbl = k
          pbl = pbl + dz(k)
        endif
      enddo
 100  continue

      call micromet(temp(1),temp0,press(1),press0,dz(1)/2.,wind(1),z0,
     &              pbl,ustar,el,psih,wstar,lstable)

      lconv = .false.
      hoverl = pbl/el
      if (((thetav(1) - thetav(2)) .gt. 1.e-8 ) .and.
     &    (hoverl .lt. -0.1) .and.
     &    (kpbl .gt. 3)) lconv = .true.
      if (.not.lconv) return
c
c-----Initialize ACM2 arrays
c
      do k = 1,nz - 1
        mbarks(k) = 0.
        mdwn(k) = 0.
      enddo
      mdwn(nz) = 0.
c
c-----Calculate ACM2 non-local mixing rates
c
      do k = 1,nz
        dz(k) = dz(k)*rho(k)
        dzinv(k) = 1./dz(k)
      enddo
      do k = 1,nz-1
        roplus = (rho(k)+rho(k+1))/2.
        dzm(k) = dzm(k)*roplus
        dzminv(k) = 1./dzm(k)
        seddy(k) = seddy(k)*roplus*roplus*dzminv(k)
      enddo
      dzminv(nz) = 0.

      mbar = 0.
      rhopbl = 0.
      do k = 2,kpbl
        rhopbl = rhopbl + dz(k)
      enddo
      meddy = seddy(1)/rhopbl
      fconv = 1./(1. + ((vk/(-hoverl))**0.3333)/(0.72*vk))
      mbar = meddy*fconv
      do k = 1,kpbl-1
        mbarks(k) = mbar
        rhosum = 0.
        do kk = k,kpbl
          rhosum = rhosum + dz(kk)
        enddo
        mdwn(k) = mbar*rhosum*dzinv(k)
      enddo
      mbarks(kpbl) = mbar
      mdwn(kpbl) = mbarks(kpbl)
c
c-----Modify timestep for ACM2
c
      rz = rhopbl*dzinv(1)
      dtacm = 1./(mbar*rz)
      dt = min(0.75*dtacm,deltat)
c
c-----Determine number of steps to take
c
      nstep = int(deltat/dt + 0.99)
      dt   = deltat/real(nstep)
      dfacp = theta*dt
      dfacq = thbar*dt

      do k = 1,nz
        dfsp(k) = dfacp*dzinv(k)
        dfsq(k) = dfacq*dzinv(k)
      enddo
c
c-----Initialize matrix elements
c
      do k = 1,nz
        aa(k) = 0.
        bb(k) = 0.
        ee(k) = 0.
      enddo
      bb(1) = 1. + rhopbl*dfsp(1)*mbarks(1)
      lfac1 = dfsq(1)*rhopbl*mbarks(1)
      lfac2 = dfsq(1)*mdwn(2)*dz(2)
      do k = 2,kpbl
        aa (k) = -dfacp*mbarks(k)
        bb(k) = 1. + dfacp*mdwn(k)
        ee(k) = -dfsp(k-1)*dz(k)*mdwn(k)
        mfac(k) = dz(k+1)*dzinv(k)*mdwn(k+1)
      enddo
c
c-----Loop over species
c
      do 301 l = 1,ns
        do m = 1,1+nddm
          do k = 1,nz
            cc(k,m) = con((m-1)*nz+k,l)/rho(k)
          enddo
        enddo
c
c-----Loop over sub-steps
c
        do 300 n = 1,nstep
          do m = 1,1+nddm
            do k = 1,nz
              dd(k) = 0.
              uu(k) = 0.
            enddo
c
c-----Call MATRIX1 solver for ACM convection
c
            dd(1) = cc(1,m) - lfac1*cc(1,m) + lfac2*cc(2,m)
            do k = 2,kpbl
              delc = mbarks(k)*cc(1,m)
     &               - mdwn(k)*cc(k,m)
     &               + mfac(k)*cc(k+1,m)
              dd(k) = cc(k,m) + dfacq*delc
            enddo
            call matrix1(nz,kpbl,aa,bb,ee,dd,uu)
            do k = 1,kpbl
              cc(k,m) = uu(k)
            enddo
          enddo
300     continue

        do m = 1,1+nddm
          do k = 1,nz
            con((m-1)*nz+k,l) = cc(k,m)*rho(k)
          enddo
        enddo

301   continue

      return
      end
c
c-------------------------------------------------------------------------------
c
      SUBROUTINE MATRIX1 ( nz, KL, A, B, E, D, X )
c
c This routine taken from:
c Models-3/CMAQ matrix.F,v 1.5 2011/10/21 16:11:45
c
C-------------------------------------------------------------------------------
C Rather than solving the ACM2 banded tridiagonal matrix using LU decomposition,
C it is much faster to split the solution into the ACM1 convective solver
C followed by the tridiagonal solver.
C MATRIX1 is the ACM1 solver. When the PBL is convective, this solver is called
C followed by TRI. If not convective, only TRI is called.
C
C-- ACM1 Matrix is in this form (there is no subdiagonal:
C   B1 E2                     <- note E2 (flux from layer above), not E1
C   A2 B2 E3
C   A3    B3 E4
C   A4       B4 E5
C   A5          B5 E6
C   A6             B6
C-------------------------------------------------------------------------------

      IMPLICIT NONE

C Arguments:

      integer nz
      INTEGER KL     ! CBL layer height
      REAL A( nz )   ! matrix column one
      REAL B( nz )   ! diagnonal
      REAL E( nz )   ! superdiagonal
      REAL D( nz )   ! R.H.S.
      REAL X( nz )   ! returned solution

C Locals:

      INTEGER NLAYS
      INTEGER L
      REAL BETA
      REAL ALPHA, GAMA

      NLAYS = nz

C-- ACM1 matrix solver

      BETA = D( 1 )
      GAMA = B( 1 )
      ALPHA = 1.0

      DO L = 2, KL
         ALPHA = -ALPHA * E( L ) / B( L )
         BETA  = ALPHA * D( L ) + BETA
         GAMA = GAMA + ALPHA * A( L )
      END DO

      X( 1 )  = BETA / GAMA
      X( KL ) = ( D( KL ) - A( KL ) * X( 1 ) ) / B( KL )

C-- Back sub for Ux=y

      DO L = KL-1, 2, -1
         X( L ) = ( D( L ) - A( L ) * X( 1 ) - E( L+1 ) * X( L+1 ) )
     &            / B( L )
      END DO

      RETURN
      END
