!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    SUBROUTINE VDIFQ(LMHK,KTM,DTQ2,Q2,EL,Z)
!     ******************************************************************
!     *                                                                *
!     *  VERTICAL DIFFUSION                                            *
!     *                                                                *
!     ******************************************************************
!-----------------------------------------------------------------------
    INCLUDE "parmeta.f90"
#include "sp.h"
!-----------------------------------------------------------------------
    PARAMETER &
    (LP1=LM+1,LM1=LM-1,LM2=LM-2)
    PARAMETER &
    (ESQ=0.20,ELZ0=0.)
!-----------------------------------------------------------------------
    DIMENSION &
    Q2    (LM)
    DIMENSION &
    EL    (LM1) &
    ,Z     (LP1)
    DIMENSION &
    CM    (LM2),CR    (LM2),RSQ2  (LM2),AKQ   (LM2),DTOZ  (LM2)
!-----------------------------------------------------------------------
!***********************************************************************
    DTDIF=DTQ2/FLOAT(KTM)
    LMHM=LMHK-1
    LMH2=LMHK-2
    LMHP=LMHK+1
!-----------------------------------------------------------------------
    DO 300 KT=1,KTM
    !-----------------------------------------------------------------------
        DO 100 L=1,LMH2
            DTOZ(L)=(DTDIF+DTDIF)/(Z(L)-Z(L+2))
            AKQ(L)=SQRT((Q2(L)+Q2(L+1))*0.5)*(EL(L)+EL(L+1))*0.5*ESQ &
            /(Z(L+1)-Z(L+2))
            CR(L)=-DTOZ(L)*AKQ(L)
        100 END DO
    
        CM(1)=DTOZ(1)*AKQ(1)+1.
        RSQ2(1)=Q2(1)
    !-----------------------------------------------------------------------
        DO 110 L=2,LMH2
            CF=-DTOZ(L)*AKQ(L-1)/CM(L-1)
            CM(L)=-CR(L-1)*CF+(AKQ(L-1)+AKQ(L))*DTOZ(L)+1.
            RSQ2(L)=-RSQ2(L-1)*CF+Q2(L)
        110 END DO
    !-----------------------------------------------------------------------
        DTOZS=(DTDIF+DTDIF)/(Z(LMHM)-Z(LMHP))
        AKQS=SQRT((Q2(LMHM)+Q2(LMHK))*0.5)*(EL(LMHM)+ELZ0)*0.5*ESQ &
        /(Z(LMHK)-Z(LMHP))
    
        CF=-DTOZS*AKQ(LMH2)/CM(LMH2)
    !-----------------------------------------------------------------------
        Q2(LMHM)=(DTOZS*AKQS*Q2(LMHK)-RSQ2(LMH2)*CF+Q2(LMHM)) &
        /((AKQ(LMH2)+AKQS)*DTOZS-CR(LMH2)*CF+1.)
    
        DO 120 IVI=1,LMH2
            L=LMHM-IVI
            Q2(L)=(-CR(L)*Q2(L+1)+RSQ2(L))/CM(L)
        120 END DO
    !-----------------------------------------------------------------------
    300 END DO
!-----------------------------------------------------------------------
    RETURN
    END SUBROUTINE VDIFQ