SUBROUTINE vdifq (lmh, ktm, dtq2, q2, el, z)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !
  !  vertical diffusion   
  !
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER :: lp1=lm+1, lm1=lm-1, lm2=lm-2
  REAL   , PARAMETER :: esq=0.20, elz0=0.
  INTEGER, DIMENSION(im,jm   ),  INTENT(IN)    :: lmh
  INTEGER,                       INTENT(IN)    :: ktm
  REAL   ,                       INTENT(IN)    :: dtq2
  REAL   , DIMENSION(im,jm,lm ), INTENT(INOUT) :: q2
  REAL   , DIMENSION(im,jm,lm1), INTENT(IN   ) :: el
  REAL   , DIMENSION(0:im+1,0:jm+1,lp1), INTENT(IN   ) :: z


  INTEGER :: kt, l, i, j
  REAL    :: dtdif, cf, dtozs, akqs
  REAL   , DIMENSION(im,jm,lm2)                :: cm
  REAL   , DIMENSION(im,jm,lm2)                :: cr
  REAL   , DIMENSION(im,jm,lm2)                :: rsq2
  REAL   , DIMENSION(im,jm,lm2)                :: akq
  REAL   , DIMENSION(im,jm,lm2)                :: dtoz

  dtdif=dtq2/float(ktm)

  DO kt=1,ktm
     DO l=1,lm
        DO j = 3, jm-2
           DO i = 1, im
              IF (l <= lmh(i,j)-2) THEN
                 dtoz(i,j,l)=(dtdif+dtdif)/(z(i,j,l)-z(i,j,l+2))
                 akq(i,j,l)=SQRT((q2(i,j,l)+q2(i,j,l+1))*0.5)*&
                      (el(i,j,l)+el(i,j,l+1))*0.5*esq &
                      /(z(i,j,l+1)-z(i,j,l+2))
                 cr(i,j,l)=-dtoz(i,j,l)*akq(i,j,l)
              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)*akq(i,j,1)+1.
           rsq2(i,j,1)=q2(i,j,1)
        END DO
     END DO
     DO l=2,lm
        DO j = 3, jm-2
           DO i = 1, im
              IF (l <= lmh(i,j)-2 ) THEN
                 cf=-dtoz(i,j,l)*akq(i,j,l-1)/cm(i,j,l-1)
                 cm(i,j,l)=-cr(i,j,l-1)*cf+&
                      (akq(i,j,l-1)+akq(i,j,l))*dtoz(i,j,l)+1.
                 rsq2(i,j,l)=-rsq2(i,j,l-1)*cf+q2(i,j,l)
              END IF
           END DO
        END DO
     END DO
     DO j = 3, jm-2
        DO i = 1, im
           dtozs=(dtdif+dtdif)/(z(i,j,lmh(i,j)-1)-z(i,j,lmh(i,j)+1))
           akqs=SQRT((q2(i,j,lmh(i,j)-1)+q2(i,j,lmh(i,j)))*0.5)*&
                (el(i,j,lmh(i,j)-1)+elz0)*0.5*esq &
                /(z(i,j,lmh(i,j))-z(i,j,lmh(i,j)+1))
           cf=-dtozs*akq(i,j,lmh(i,j)-2)/cm(i,j,lmh(i,j)-2)
           q2(i,j,lmh(i,j)-1)=(dtozs*akqs*q2(i,j,lmh(i,j))-&
                rsq2(i,j,lmh(i,j)-2)*cf+q2(i,j,lmh(i,j)-1)) &
                /((akq(i,j,lmh(i,j)-2)+akqs)*dtozs-&
                cr(i,j,lmh(i,j)-2)*cf+1.)
        END DO
     END DO
     DO l=lm, 1,-1
        DO j = 3, jm-2
           DO i = 1, im
              IF (l <= lmh(i,j)-2) THEN
                 q2(i,j,l)=(-cr(i,j,l)*q2(i,j,l+1)+rsq2(i,j,l))/cm(i,j,l)
              END IF
           END DO
        END DO
     END DO
  END DO
END SUBROUTINE vdifq
