SUBROUTINE vtadv
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !
  !  subprogram DOcumentation block
  !                .      .    .     
  ! subprogram:    vtadv       vertical advection
  !   prgrmmr: janjic          org: w/np2      date: 93-11-17
  !     
  ! abstract:
  !     vtadv calculates the contribution of the vertical advection
  !     to the tendencies of temperature, specIFic humidity, wind
  !     components, and turbulent kinetic energy and THEN updates those
  !     variables.  for all variables except specIFic humidity a
  !     simple centered dIFference scheme in space is used in
  !     conjunction with the pure euler-backward time scheme.  
  !     a piecewise linear scheme is used to calculate the vertical
  !     advection of specIFic humidity so that no false maxima or
  !     minima are produced.
  !     
  ! program history log:
  !   87-06-??  janji!     - originator
  !   90-??-??  mesinger   - inserted piecewise linear scheme for
  !                          specIFic humidity
  !   95-03-25  black      - conversion from 1-d to 2-d in horizontal
  !   95-11-20  abeles     - parallel optimization
  !   96-03-29  black      - added external edge; removed scrch common
  !     
  ! usage: call vtadv from main program ebu
  !   input argument list:
  !       none     
  !  
  !   output argument list: 
  !     none
  !     
  !   output files:
  !     none
  !     
  !   subprograms called:
  !  
  !     unique: none
  !  
  !     library: none
  !  
  !   common blocks: ctlblk
  !                  masks
  !                  dynam
  !                  vrbls
  !                  pvrbls
  !                  indx
  !  
  IMPLICIT NONE
  REAL, PARAMETER	:: d00=.0e0,d125=.125e0,edqmx=2.e-5,edqmn=-2.e-5
  REAL, PARAMETER	:: epsq=1.e-12,h1=1.e0,h2=2.e0
  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER	:: lda=lm+9,la=13,ksmud=0
  ! warning   IF lm.gt.16 redefine l1 and l3 for the equivalences
  INTEGER, PARAMETER	:: imjm=im*jm-jm/2,jam=6+2*(jm-10)
  INTEGER, PARAMETER	:: lm1=lm-1,lm2=lm-2,lp1=lm+1
  INTEGER, PARAMETER	:: lscrch=4*lm+1+la +1,l1=lda+1,l3=lda+1+lm
  LOGICAL  :: nosla
  !
  INCLUDE "ctlblk.f90"
  INCLUDE "masks.f90"
  INCLUDE "dynam.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "contin.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "indx.f90"
  !
  INTEGER	:: l,nmsap,j,i,nsmud,ks
  REAL	:: ttak,extrem,msa,dqtik,astik,asbik,qdedtk,qdeutk
  REAL	:: sedbk,dqdek,edbfk,edtdk,tq2ak,vmk,tuak,tvak
  REAL, DIMENSION(lm1)	:: wfa ,wfb
  REAL, DIMENSION(im,jm)	:: etadtl,ttb,tq2b,tub,tvb
  REAL, DIMENSION(im,jm)	:: dqti,dqbi,rpdx,rpdy
  REAL, DIMENSION(im,jm)	:: qdedb,qdeub,edbd
  REAL, DIMENSION(im,jm)	:: dqdeb
  REAL, DIMENSION(0:im+1,0:jm+1)	:: fne,fse
  REAL, DIMENSION(im,jm,lm)	:: tstl,ustl,vstl,sam,qbi,q2st
  REAL, DIMENSION(im,jm,lm1)	:: array1,array2
  ! 
  ! define added upstream advection constants
  ! 
  DO l=1,lm1
     wfa(l)=deta(l  )/(deta(l)+deta(l+1))
     wfb(l)=deta(l+1)/(deta(l)+deta(l+1))
  END DO
  ! no moisture slope adjustment IF not wanted
  nosla=.false.
  !  IF false, number of moisture slope adjustment passes
  nmsap=3
  ! smoothing vertical velocity at h points
  IF(ksmud.gt.0)THEN
     DO l=1,lm1
        DO j=1,jm
           DO i=1,im
              etadt(i,j,l)=etadt(i,j,l)*hbm2(i,j)
           END DO
        END DO
        nsmud=ksmud
        ! 
        !  the fne, fse, etadtl, and etadt arrays
        !  are on or associated with h points
        !
        DO ks=1,nsmud
           DO j=1,jm-1
              DO i=1,im
                 fne(i,j)=(etadt(i+ihe(j),j+1,l)-etadt(i,j,l)) &
                      *htm(i,j,l+1)*htm(i+ihe(j),j+1,l+1)
              END DO
           END DO
           DO j=2,jm
              DO i=1,im
                 fse(i,j)=(etadt(i+ihe(j),j-1,l)-etadt(i,j,l)) &
                      *htm(i+ihe(j),j-1,l+1)*htm(i,j,l+1)
              END DO
           END DO
           DO j=3,jm-2
              DO i=1,im
                 etadtl(i,j)=(fne(i,j)-fne(i+ihw(j),j-1) &
                      +fse(i,j)-fse(i+ihw(j),j+1))*hbm2(i,j)
              END DO
           END DO
           DO j=3,jm-2
              DO i=1,im
                 etadt(i,j,l)=etadtl(i,j)*d125+etadt(i,j,l)
              END DO
           END DO
        END DO
     END DO
  ENDIF
  ! vertical (matsuno) advection of t
  DO j=1,jm
     DO i=1,im
        ttb(i,j)=d00
     END DO
  END DO
  DO l=1,lm1
     DO j=3,jm-2
        DO i=1,im
           ttak  =(t(i,j,l+1)-t(i,j,l))*etadt(i,j,l)*f4d
           tstl(i,j,l)=(ttak  +ttb(i,j))*rdeta(l)+t(i,j,l)
           ttb(i,j)=ttak
        END DO
     END DO
  END DO
  DO j=3,jm-2
     DO i=1,im
        tstl(i,j,lm)=t(i,j,lm)+ttb(i,j)*rdeta(lm)
     END DO
  END DO
  ! second (backward) matsuno step
  DO j=1,jm
     DO i=1,im
        ttb(i,j)=d00
     END DO
  END DO
  DO l=1,lm1
     DO j=3,jm-2
        DO i=1,im
           ttak  =(tstl(i,j,l+1)-tstl(i,j,l))*etadt(i,j,l)*f4d
           t(i,j,l)=(ttak  +ttb(i,j))*rdeta(l)+t(i,j,l)
           ttb(i,j)=ttak
        END DO
     END DO
  END DO
  DO j=3,jm-2
     DO i=1,im
        t(i,j,lm)=t(i,j,lm)+ttb(i,j)*rdeta(lm)
     END DO
  END DO
  ! piecewise linear upstream vertical advection of q
  !       intialize q at the bottom interface and the slope adjustment
  !       mask (sam=1 for sa permitted, 0 for not permitted)
  ! 
  DO l=1,lm
     DO j=3,jm-2
        DO i=1,im
           qbi(i,j,l)=q(i,j,l)
           sam(i,j,l)=h1
        END DO
     END DO
  END DO
  !    IF(nosla) go to 290
  IF(.NOT.nosla) THEN
     !
     ! the slope adjustment code
     ! no slope permitted at the top and at the bottom layer
     ! 
     DO j=3,jm-2
        DO i=1,im
           sam(i,j, 1)=d00
           sam(i,j,lm)=d00  
        END DO
     END DO
     DO l=1,lm1
        DO j=3,jm-2
           DO i=1,im
              sam(i,j,l)=sam(i,j,l)*htm(i,j,l+1)
           END DO
        END DO
     END DO
     ! 
     !   now, search for the maxima and minima of q (at the first
     !   pass) and for layers which had overadjusted (at subsequent
     !   passes) due to round-off errors
     ! 
     DO l=2,lm1
        DO j=3,jm-2
           DO i=1,im
              dqti(i,j)=q(i,j,l)-q(i,j,l-1)
              dqbi(i,j)=q(i,j,l+1)-q(i,j,l)
              extrem=dqti(i,j)*dqbi(i,j)
              IF(extrem.lt.d00)sam(i,j,l)=d00
           END DO
        END DO
     END DO
     DO l=2,lm1
        DO j=3,jm-2
           DO i=1,im
              array1(i,j,l)=wfa(l-1)*(h1-sam(i,j,l-1))+wfb(l-1)
              array2(i,j,l)=wfa(l)+wfb(l)*(h1-sam(i,j,l+1))
           END DO
        END DO
     END DO
     DO msa=1,nmsap
        ! 
        !  calculate dq at interfaces and adjust the slopes where
        !  and to the extent permitted observing the monotonicity
        !  condition (e.g. van leer, j. comp. phys. 1977, 276-299)
        ! 
        DO j=3,jm-2
           DO i=1,im
              dqbi(i,j)=h2*q(i,j,2)-qbi(i,j,2) -qbi(i,j,1)
           END DO
        END DO
        DO l=2,lm1
           DO j=3,jm-2
              DO i=1,im
                 dqtik  =dqbi(i,j)
                 astik  =array1(i,j,l)*dqtik
                 dqbi(i,j)=h2*q(i,j,l+1)-qbi(i,j,l+1)-qbi(i,j,l)
                 asbik  =array2(i,j,l)*dqbi(i,j)
                 qbi(i,j,l)=qbi(i,j,l) &
                      +(astik-sign(h1,astik) &
                      *dim(abs(astik),abs(asbik)))*sam(i,j,l)
              END DO
           END DO
        END DO
     END DO
     ! 
     !  slope adjustment of the layers above that next to the surface
     !  is DOne; now adjust the lowermost layer
     ! 
     DO l=9,lm1
        DO j=3,jm-2
           DO i=1,im
              IF(htm(i,j,l+1).eq.d00)qbi(i,j,l)=h2*q(i,j,l)-qbi(i,j,l-1)
           END DO
        END DO
     END DO
     DO j=3,jm-2
        DO i=1,im
           qbi(i,j,lm)=h2*q(i,j,lm)-qbi(i,j,lm1)
        END DO
     END DO
     ! 
     ! end of the slope adjustment code
     ! 
  ENDIF
  DO j=3,jm-2
     DO i=1,im
        qdedb(i,j)=d00
        qdeub(i,j)=d00
        dqdeb(i,j)=h2*(qbi(i,j,1)-q(i,j,1))*rdeta(1)
        edbd (i,j)=d00
     END DO
  END DO
  DO l=1,lm1
     DO j=3,jm-2
        DO i=1,im
           qdedtk  =qdedb(i,j)
           qdeutk  =qdeub(i,j)
           sedbk   =sign(h1,etadt(i,j,l))
           dqdek   =dqdeb(i,j)
           dqdeb(i,j)=h2*(qbi(i,j,l+1)-q(i,j,l+1))*rdeta(l+1)
           edbfk   =etadt(i,j,l)*f4d
           qdedb(i,j)=(h1+sedbk)*(qbi(i,j,l)+dqdek*edbfk)*(-edbfk)
           qdeub(i,j)=(h1-sedbk)*(q(i,j,l+1)+q(i,j,l+1)-qbi(i,j,l+1) &
                +dqdeb(i,j)*edbfk)*edbfk
           edtdk   =edbd(i,j)
           edbd (i,j)=etadt(i,j,l)*(-f4q)
           q(i,j,l)=q(i,j,l)+(qdedtk-qdeutk-qdedb(i,j)+qdeub(i,j) &
                +q(i,j,l)*(edbd(i,j)-edtdk))*rdeta(l)
        END DO
     END DO
  END DO
  DO j=3,jm-2
     DO i=1,im
        q(i,j,lm)=q(i,j,lm)+(qdedb(i,j)-qdeub(i,j) &
             +q(i,j,lm)*(-edbd(i,j)))*rdeta(lm)
     END DO
  END DO
  !
  ! negative moisture may occur due to violation of the cfl
  !
  DO l=1,lm1
     DO j=3,jm-2
        DO i=1,im
           IF(q(i,j,l).lt.epsq)THEN
              dqbi(i,j)=q(i,j,l)
              q(i,j,l)=epsq
              q(i,j,l+1)=q(i,j,l+1)+deta(l)*rdeta(l+1)*dqbi(i,j)
           ENDIF
        END DO
     END DO
  END DO
  DO j=3,jm-2
     DO i=1,im
        IF(q(i,j,lm).lt.epsq)q(i,j,lm)=epsq
     END DO
  END DO
  !
  ! vertical (matsuno) advection of q2
  !
  DO j=3,jm-2
     DO i=1,im
        tq2b(i,j)=q2(i,j,1)*etadt(i,j,1)*f4q2(1)
     END DO
  END DO
  DO l=1,lm2
     DO j=3,jm-2
        DO i=1,im
           tq2ak=(q2(i,j,l+1)-q2(i,j,l))*(etadt(i,j,l)+etadt(i,j,l+1)) &
                *f4q2(l+1)
           q2st(i,j,l)=tq2ak+tq2b(i,j)+q2(i,j,l)
           tq2b(i,j)=tq2ak
        END DO
     END DO
  END DO
  DO j=3,jm-2
     DO i=1,im
        tq2ak=(q2(i,j,lm)-q2(i,j,lm1))*etadt(i,j,lm1)*f4q2(lm)
        q2st(i,j,lm1)=tq2ak+tq2b(i,j)+q2(i,j,lm1)
        q2st(i,j,lm )=q2(i,j,lm)
     END DO
  END DO
  !
  ! second (backward) matsuno step
  !
  DO j=3,jm-2
     DO i=1,im
        tq2b(i,j)=q2st(i,j,1)*etadt(i,j,1)*f4q2(1)
     END DO
  END DO
  DO l=1,lm2
     DO j=3,jm-2
        DO i=1,im
           tq2ak =(q2st(i,j,l+1)-q2st(i,j,l)) &
                *(etadt(i,j,l)+etadt(i,j,l+1))*f4q2(l+1)
           q2(i,j,l)=tq2ak+tq2b(i,j)+q2(i,j,l)
           tq2b(i,j)=tq2ak
        END DO
     END DO
  END DO
  DO j=3,jm-2
     DO i=1,im
        tq2ak  =(q2st(i,j,lm)-q2st(i,j,lm1))*etadt(i,j,lm1)*f4q2(lm)
        q2(i,j,lm1)=tq2ak+tq2b(i,j)+q2(i,j,lm1)
     END DO
  END DO
  !
  ! definition of variables needed at v points
  !
  DO l=1,lm1
     DO j=1,jm
        DO i=1,im
           etadt(i,j,l)=etadt(i,j,l)*pdsl(i,j)*hbm2(i,j)
        END DO
     END DO
  END DO
  DO j=3,jm-2
     DO i=1,im
        rpdx(i,j)=h1/(pdsl(i+ivw(j),j)+pdsl(i+ive(j),j))
        rpdy(i,j)=h1/(pdsl(i,j-1)+pdsl(i,j+1))
     END DO
  END DO
  !
  ! vertical (matsuno) advection of u & v
  !
  DO j=1,jm
     DO i=1,im
        tub(i,j)=d00
        tvb(i,j)=d00
     END DO
  END DO
  DO l=1,lm1
     DO j=3,jm-2
        DO i=1,im
           vmk   =vtm(i,j,l+1)*vbm2(i,j)
           tuak  =(etadt(i+ivw(j),j,l)+etadt(i+ive(j),j,l)) &
                *(u(i,j,l+1)-u(i,j,l))*rpdx(i,j)*(vmk*f4d)
           ustl(i,j,l)=(tuak+tub(i,j))*rdeta(l)+u(i,j,l)
           tub(i,j)=tuak
           tvak  =(etadt(i,j-1,l)+etadt(i,j+1,l))*(v(i,j,l+1)-v(i,j,l)) &
                *rpdy(i,j)*(vmk*f4d)
           vstl(i,j,l)=(tvak+tvb(i,j))*rdeta(l)+v(i,j,l)
           tvb(i,j)=tvak
        END DO
     END DO
  END DO
  DO j=3,jm-2
     DO i=1,im
        ustl(i,j,lm)=u(i,j,lm)+tub(i,j)*rdeta(lm)
        vstl(i,j,lm)=v(i,j,lm)+tvb(i,j)*rdeta(lm)
     END DO
  END DO
  ! 
  ! second (backward) matsuno step
  !
  DO j=1,jm
     DO i=1,im
        tub(i,j)=d00
        tvb(i,j)=d00
     END DO
  END DO
  DO l=1,lm1
     DO j=3,jm-2
        DO i=1,im
           vmk   =vtm(i,j,l+1)*vbm2(i,j)
           tuak  =(etadt(i+ivw(j),j,l)+etadt(i+ive(j),j,l)) &
                *(ustl(i,j,l+1)-ustl(i,j,l))*rpdx(i,j)*(vmk*f4d)
           u(i,j,l)=(tuak+tub(i,j))*rdeta(l)+u(i,j,l)
           tub(i,j)=tuak
           tvak  =(etadt(i,j-1,l)+etadt(i,j+1,l)) &
                *(vstl(i,j,l+1)-vstl(i,j,l))*rpdy(i,j)*(vmk*f4d)
           v(i,j,l)=(tvak+tvb(i,j))*rdeta(l)+v(i,j,l)
           tvb(i,j)=tvak
        END DO
     END DO
  END DO
  DO j=3,jm-2
     DO i=1,im
        u(i,j,lm)=u(i,j,lm)+tub(i,j)*rdeta(lm)
        v(i,j,lm)=v(i,j,lm)+tvb(i,j)*rdeta(lm)
     END DO
  END DO
END SUBROUTINE vtadv
