SUBROUTINE vdifh(lmh, ktmq2, dtq2, thz0, qz0, akhs, ct, cklq, &
     t, q, akhwork, ape, zint)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !
  ! vertical diffusion of mass variables                *
  !
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER                                     :: lp1=lm+1
  INTEGER, PARAMETER                                     :: lm1=lm-1
  INTEGER, DIMENSION(  im  ,   jm       ), INTENT(IN   ) :: lmh
  INTEGER,                                 INTENT(IN   ) :: ktmq2
  REAL   ,                                 INTENT(IN   ) :: dtq2
  REAL   , DIMENSION(  im  ,   jm       ), INTENT(IN   ) :: thz0
  REAL   , DIMENSION(  im  ,   jm       ), INTENT(IN   ) :: qz0
  REAL   , DIMENSION(  im  ,   jm       ), INTENT(IN   ) :: akhs
  REAL   , DIMENSION(  im  ,   jm       ), INTENT(IN   ) :: ct
  REAL   , DIMENSION(  im  ,   jm       ), INTENT(IN   ) :: cklq
  REAL   , DIMENSION(0:im+1, 0:jm+1, lm ), INTENT(INOUT) :: t, q
  REAL   , DIMENSION(  im  ,   jm  , lm1), INTENT(IN   ) :: akhwork
  REAL   , DIMENSION(  im,     jm  , lm ), INTENT(IN   ) :: ape
  REAL   , DIMENSION(0:im+1, 0:jm+1, lp1), INTENT(IN   ) :: zint
  !
  !
  !
  REAL   , DIMENSION(  im  ,   jm  , lm1) :: dtoz
  REAL   , DIMENSION(  im  ,   jm  , lm1) :: cm
  REAL   , DIMENSION(  im  ,   jm  , lm1) :: cr
  REAL   , DIMENSION(  im  ,   jm  , lm1) :: rst
  REAL   , DIMENSION(  im  ,   jm  , lm1) :: rsq
  REAL   , DIMENSION(  im  ,   jm  , lm1) :: akct
  INTEGER                 :: i, j, l, kt, lmax
  REAL                    :: dtdif, dtozs, cf 
  REAL                    :: akhh, akqs, cmb, cmtb, cmqb
  REAL                    :: rstb, rsqb, rclm

  dtdif=dtq2/float(ktmq2)

  lmax = MAXVAL (lmh(1:im, 3:jm-2))

  DO l=1,lmax-1
     DO j = 3, jm-2
        DO i = 1, im
           IF (l <= lmh(i,j)-1) THEN
              dtoz(i,j,l)=dtdif/(zint(i,j,l)-zint(i,j,l+1))
              cr(i,j,l)=-dtoz(i,j,l)*akhwork(i,j,l)
              akct(i,j,l)=akhwork(i,j,l)*(zint(i,j,l)-zint(i,j,l+2)) &
                   *0.5*ct(i,j)
           END IF
        END DO
     END DO
  END DO
  DO j = 3, jm-2
     DO i = 1, im
        cm(i,j,1)=dtoz(i,j,1)*akhwork(i,j,1)+1.
     END DO
  END DO

  DO kt=1,ktmq2

     DO j = 3, jm-2
        DO i = 1, im
           rst(i,j,1)=-akct(i,j,1)*dtoz(i,j,1)+t(i,j,1)*ape(i,j,1)
           rsq(i,j,1)=q(i,j,1)
        END DO
     END DO

     DO l = 2, lmax-1
        DO j = 3, jm-2
           DO i = 1, im
              IF (l <= lmh(i,j)-1) THEN
                 cf=-dtoz(i,j,l)*akhwork(i,j,l-1)/cm(i,j,l-1)
                 cm(i,j,l)=-cr(i,j,l-1)*cf+ &
                      (akhwork(i,j,l-1)+akhwork(i,j,l))*dtoz(i,j,l)+1.
                 rst(i,j,l)=-rst(i,j,l-1)*cf+ &
                      (akct(i,j,l-1)-akct(i,j,l))*dtoz(i,j,l)+ &
                      t(i,j,l)*ape(i,j,l)
                 rsq(i,j,l)=-rsq(i,j,l-1)*cf+q(i,j,l)
              END IF
           END DO
        END DO
     END DO

     DO j = 3, jm-2
        DO i = 1, im
           dtozs=dtdif/(zint(i,j,lmh(i,j))-zint(i,j,lmh(i,j)+1))
           akhh=akhwork(i,j,lmh(i,j)-1)
           cf=-dtozs*akhh/cm(i,j,lmh(i,j)-1)
           akqs=akhs(i,j)*cklq(i,j)
           cmb=cr(i,j,lmh(i,j)-1)*cf
           cmtb=-cmb+(akhh+akhs(i,j))*dtozs+1.
           cmqb=-cmb+(akhh+akqs)*dtozs+1.
           rstb=-rst(i,j,lmh(i,j)-1)*cf &
                +(akct(i,j,lmh(i,j)-1)-akhs(i,j)*ct(i,j))*dtozs &
                +t(i,j,lmh(i,j))*ape(i,j,lmh(i,j))
           rsqb=-rsq(i,j,lmh(i,j)-1)*cf+q(i,j,lmh(i,j))
           t(i,j,lmh(i,j))=(dtozs*akhs(i,j)*thz0(i,j)+rstb) &
                /(ape(i,j,lmh(i,j))*cmtb)
           q(i,j,lmh(i,j))=(dtozs*akqs*qz0(i,j) +rsqb)/cmqb
        END DO
     END DO

     DO l = lmax-1, 1, -1
        DO j = 3, jm-2
           DO i = 1, im
              IF (l <= lmh(i,j)-1) THEN
                 rclm=1./cm(i,j,l)
                 t(i,j,l)=(-cr(i,j,l)*t(i,j,l+1)*ape(i,j,l+1)+rst(i,j,l))&
                      *rclm/ape(i,j,l)
                 q(i,j,l)=(-cr(i,j,l)*q(i,j,l+1)             +rsq(i,j,l))&
                      *rclm
              END IF
           END DO
        END DO
     END DO
  END DO
END SUBROUTINE vdifh
