!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    SUBROUTINE VDIFH(LMHK,KTM,DTQ2,THZ0,QZ0,AKHS,CT,CKLQ &
! Dule     &,                T,Q,AKH,APE,Z)
    ,                T,Q,AKH,APE,Z,EVPR)
!     ******************************************************************
!     *                                                                *
!     *            VERTICAL DIFFUSION OF MASS VARIABLES                *
!     *                                                                *
!     ******************************************************************
!-----------------------------------------------------------------------
    INCLUDE "parmeta.f90"
    INCLUDE "mpp.h"
#include "sp.h"
!-----------------------------------------------------------------------
    PARAMETER &
    (LP1=LM+1,LM1=LM-1,JAM=6+2*(JM-10))
!-----------------------------------------------------------------------
    INCLUDE "COMM_DYNAM.f90"
!-----------------------------------------------------------------------
    DIMENSION &
    T     (LM),Q     (LM)
    DIMENSION &
    AKH   (LM1) &
    ,APE   (LM) &
    ,Z     (LP1)
    DIMENSION &
    CM    (LM1),CR    (LM1),RST   (LM1),RSQ   (LM1) &
    ,DTOZ  (LM1),AKCT  (LM1)
!-----------------------------------------------------------------------
!***********************************************************************
    DTDIF=DTQ2/FLOAT(KTM)
    LMHM=LMHK-1
    LMHP=LMHK+1
!-----------------------------------------------------------------------
    EVPR=0.E0                                                         !dule
!-----------------------------------------------------------------------
    DO 100 L=1,LMHM
        DTOZ(L)=DTDIF/(Z(L)-Z(L+1))
        CR(L)=-DTOZ(L)*AKH(L)
        AKCT(L)=AKH(L)*(Z(L)-Z(L+2))*0.5*CT
    100 END DO

    CM(1)=DTOZ(1)*AKH(1)+1.
!-----------------------------------------------------------------------
    DO 300 KT=1,KTM
    !-----------------------------------------------------------------------
        RST(1)=-AKCT(1)*DTOZ(1)+T(1)*APE(1)
        RSQ(1)=Q(1)
    !-----------------------------------------------------------------------
        DO 110 L=2,LMHM
            DTOZL=DTOZ(L)
            CF=-DTOZL*AKH(L-1)/CM(L-1)
            CM(L)=-CR(L-1)*CF+(AKH(L-1)+AKH(L))*DTOZL+1.
            RST(L)=-RST(L-1)*CF+(AKCT(L-1)-AKCT(L))*DTOZL+T(L)*APE(L)
            RSQ(L)=-RSQ(L-1)*CF+Q(L)
        110 END DO
    !-----------------------------------------------------------------------
        DTOZS=DTDIF/(Z(LMHK)-Z(LMHP))
        AKHH=AKH(LMHM)
    
        CF=-DTOZS*AKHH/CM(LMHM)
        AKQS=AKHS*CKLQ
    
        CMB=CR(LMHM)*CF
        CMTB=-CMB+(AKHH+AKHS)*DTOZS+1.
        CMQB=-CMB+(AKHH+AKQS)*DTOZS+1.
    
        RSTB=-RST(LMHM)*CF+(AKCT(LMHM)-AKHS*CT)*DTOZS &
        +T(LMHK)*APE(LMHK)
        RSQB=-RSQ(LMHM)*CF+Q(LMHK)
    !-----------------------------------------------------------------------
        T(LMHK)=(DTOZS*AKHS*THZ0+RSTB)/(APE(LMHK)*CMTB)
        EVPR=EVPR-Q(LMHK)*DETA(LMHK)                                      !dule
        Q(LMHK)=(DTOZS*AKQS*QZ0 +RSQB)/CMQB
        EVPR=EVPR+Q(LMHK)*DETA(LMHK)                                      !dule
    !-----------------------------------------------------------------------
        DO 120 L=LMHM,1,-1
            RCML=1./CM(L)
            T(L)=(-CR(L)*T(L+1)*APE(L+1)+RST(L))*RCML/APE(L)
            EVPR=EVPR-Q(L)*DETA(L)                                            !dule
            Q(L)=(-CR(L)*Q(L+1)         +RSQ(L))*RCML
            EVPR=EVPR+Q(L)*DETA(L)                                            !dule
        120 END DO
    !-----------------------------------------------------------------------
    300 END DO
!-----------------------------------------------------------------------
    RETURN
    END SUBROUTINE VDIFH