!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    SUBROUTINE HDIFF
!     ******************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    HDIFF       HORIZONTAL DIFFUSION
!   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 93-11-17

! ABSTRACT:
!     HDIFF CALCULATES THE CONTRIBUTION OF THE HORIZONTAL DIFFUSION
!     TO THE TENDENCIES OF TEMPERATURE, SPECIFIC HUMIDITY, WIND
!     COMPONENTS, AND TURBULENT KINETIC ENERGY AND THEN UPDATES THOSE
!     VARIABLES.  A SECOND-ORDER NONLINEAR SCHEME SIMILAR TO
!     SMAGORINSKYS IS USED WHERE THE DIFFUSION COEFFICIENT IS
!     A FUNCTION OF THE DEFORMATION FIELD AND OF THE TURBULENT
!     KINETIC ENERGY.

! PROGRAM HISTORY LOG:
!   87-06-??  JANJIC     - ORIGINATOR
!   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
!   96-03-28  BLACK      - ADDED EXTERNAL EDGE
!   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY

! USAGE: CALL HDIFF 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
!                  PHYS
!                  VRBLS
!                  PVRBLS
!                  INDX

! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$
!     ******************************************************************
    PARAMETER &
    (DEFC=8.0,DEFM=32.0,SCQ2=50.0 &
    , EPSQ2=0.2,FCDIF=1.0,RFCP=.25/1004.6)
!----------------------------------------------------------------------
    INCLUDE "parmeta.f90"
    INCLUDE "parm.tbl.f90"
    INCLUDE "mpp.h"
#include "sp.h"
!----------------------------------------------------------------------
    PARAMETER &
    (IMJM=IM*JM-JM/2,LP1=LM+1,KSMUD=1)
    PARAMETER &
    (JAM=6+2*(JM-10),JAMD=(JAM*2-10)*3)
!-----------------------------------------------------------------------
    LOGICAL :: &
    RUN,FIRST,RESTRT,SIGMA &
    ,SECOND,HEAT,STTDF
!----------------------------------------------------------------------
    INCLUDE "COMM_CTLBLK.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_MASKS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_DYNAMD.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_PHYS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_VRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_PVRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_CLDWTR.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_INDX.f90"
!-----------------------------------------------------------------------
    DIMENSION &
    Q2L  (IDIM1:IDIM2,JDIM1:JDIM2),UT   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,HKNE (IDIM1:IDIM2,JDIM1:JDIM2),HKSE (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,VKNE (IDIM1:IDIM2,JDIM1:JDIM2),VKSE (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,HMASK(IDIM1:IDIM2,JDIM1:JDIM2),HMSKL(IDIM1:IDIM2,JDIM1:JDIM2)

    DIMENSION &
    TNE  (IDIM1:IDIM2,JDIM1:JDIM2),TSE  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,QNE  (IDIM1:IDIM2,JDIM1:JDIM2),QSE  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,Q2NE (IDIM1:IDIM2,JDIM1:JDIM2),Q2SE (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,UNE  (IDIM1:IDIM2,JDIM1:JDIM2),USE  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,VNE  (IDIM1:IDIM2,JDIM1:JDIM2),VSE  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,TDIF (IDIM1:IDIM2,JDIM1:JDIM2),QDIF (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,UDIF (IDIM1:IDIM2,JDIM1:JDIM2),VDIF (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,Q2DIF(IDIM1:IDIM2,JDIM1:JDIM2) &
    ,DEF  (IDIM1:IDIM2,JDIM1:JDIM2),CKE  (IDIM1:IDIM2,JDIM1:JDIM2)
!-----------------------------------------------------------------------
!***
!***  DIFFUSING Q2 AT GROUND LEVEL DOESNT MATTER, USTAR2 IS RECALCULATED
!***
!-----------------------------------------------------------------------
    SECOND= .TRUE. 
    HEAT= .FALSE. 
    LUL=UL(1)
!-----------------------------------------------------------------------
    DO J=MYJS_P1,MYJE_P2
    ! 011115      DO I=MYIS_P1,MYIE_P2
        DO I=MYIS_P2,MYIE_P2
            HMASK(I,J)=1.
            HMSKL(I,J)=1.
        ENDDO
    ENDDO
!-----------------------------------------------------------------------
    IF(SIGMA)THEN
        DO 100 J=MYJS2_P1,MYJE2_P1
        
            DO I=MYIS1_P1,MYIE1_P1
                DH1=ABS(FIS(I+IHW(J),J-1)-FIS(I,J))
                DH2=ABS(FIS(I+IHE(J),J-1)-FIS(I,J))
                DH3=ABS(FIS(I+IHW(J),J+1)-FIS(I,J))
                DH4=ABS(FIS(I+IHE(J),J+1)-FIS(I,J))
            
                DHM=AMAX1(DH1,DH2,DH3,DH4)/DY
            
                IF(DHM > 0.100)THEN
                    HMASK(I,J)=0.
                    HMSKL(I,J)=0.
                ENDIF
            ENDDO
        
        100 END DO
    ENDIF
!-----------------------------------------------------------------------
    DO 600 KS=1,KSMUD
    !-----------------------------------------------------------------------
    !---------------------MAIN VERTICAL INTEGRATION LOOP--------------------
    !-----------------------------------------------------------------------
    ! omp parallel do
    ! omp& private(cke,def,defsk,deftk,hkne,hkse,hmskl,q2dif,q2l,q2ne,q2se,
    ! omp&         qdif,qne,qse,tdif,tne,tse,udif,une,use,utk,vdif,vkne,
    ! omp&         vkse,vne,vse,vtk)
    !-----------------------------------------------------------------------
        DO 500 L=1,LM
        !-----------------------------------------------------------------------
            CALL ZERO2(DEF)
            CALL ZERO2(Q2NE)
            CALL ZERO2(Q2SE)
            CALL ZERO2(QNE)
            CALL ZERO2(QSE)
            CALL ZERO2(TNE)
            CALL ZERO2(TSE)
            CALL ZERO2(UNE)
            CALL ZERO2(USE)
            CALL ZERO2(VSE)
            CALL ZERO2(VNE)
            CALL ZERO2(VSE)
            CALL ZERO2(TDIF)
            CALL ZERO2(QDIF)
            CALL ZERO2(UDIF)
            CALL ZERO2(VDIF)
            CALL ZERO2(Q2DIF)
        !-----------------------------------------------------------------------
            DO 210 J=MYJS_P1,MYJE_P1
                DO 210 I=MYIS_P1,MYIE_P1
                    Q2L(I,J)=AMAX1(Q2(I,J,L),EPSQ2)
            210 END DO
        !--------------DEFORMATIONS---------------------------------------------
            DO 220 J=MYJS1_P1,MYJE1_P1
                DO 220 I=MYIS_P1,MYIE1_P1
                
                    IF(L < LUL)THEN
                    ! j     HMSKL(I,J)=1.
                        HMSKL(I,J)=HMASK(I,J)
                    ELSE
                        HMSKL(I,J)=HMASK(I,J)
                    ENDIF
                
                    DEFTK  =U(I+IHE(J),J,L)-U(I+IHW(J),J,L)-V(I,J+1,L)+V(I,J-1,L)
                    DEFSK  =U(I,J+1,L)-U(I,J-1,L)+V(I+IHE(J),J,L)-V(I+IHW(J),J,L)
                ! SM      DEF (I,J)=DEFTK*DEFTK+DEFSK*DEFSK+SCQ2*Q2L(I,J)/(1.+Q2(I,J,L))
                    DEF (I,J)=DEFTK*DEFTK+DEFSK*DEFSK
                    DEF (I,J)=SQRT(DEF(I,J)+DEF(I,J))*HBM2(I,J)
                    DEF(I,J)=AMAX1(DEF(I,J),DEFC)
                !     DEF(I,J)=AMIN1(DEF(I,J),DEFM)
                    DEF(I,J)=DEF(I,J)*HMSKL(I,J)
            220 END DO
        !--------------T,Q, Q2 DIAGONAL CONTRIBUTIONS---------------------------
            DO 250 J=MYJS_P1,MYJE1_P1
                DO 250 I=MYIS_P1,MYIE1_P1
                    HKNE(I,J)=(DEF(I,J)+DEF(I+IHE(J),J+1)) &
                    *HTM(I,J,L)*HTM(I+IHE(J),J+1,L) &
                    *HMSKL(I,J)*HMSKL(I+IHE(J),J+1)
                    TNE (I,J)=(T (I+IHE(J),J+1,L)-T (I,J,L))*HKNE(I,J)
                    QNE (I,J)=(Q (I+IHE(J),J+1,L)-Q (I,J,L))*HKNE(I,J)
                    Q2NE(I,J)=(Q2(I+IHE(J),J+1,L)-Q2(I,J,L))*HKNE(I,J)
            250 END DO
        
            DO 260 J=MYJS1_P1,MYJE_P1
                DO 260 I=MYIS_P1,MYIE1_P1
                    HKSE(I,J)=(DEF(I+IHE(J),J-1)+DEF(I,J)) &
                    *HTM(I+IHE(J),J-1,L)*HTM(I,J,L) &
                    *HMSKL(I+IHE(J),J-1)*HMSKL(I,J)
                    TSE (I,J)=(T (I+IHE(J),J-1,L)-T (I,J,L))*HKSE(I,J)
                    QSE (I,J)=(Q (I+IHE(J),J-1,L)-Q (I,J,L))*HKSE(I,J)
                    Q2SE(I,J)=(Q2(I+IHE(J),J-1,L)-Q2(I,J,L))*HKSE(I,J)
            260 END DO
        !-----------------------------------------------------------------------
            DO 270 J=MYJS1,MYJE1
                DO 270 I=MYIS1,MYIE
                    TDIF (I,J)=(TNE (I,J)-TNE (I+IHW(J),J-1) &
                    +TSE (I,J)-TSE (I+IHW(J),J+1))*HDAC(I,J)
                    QDIF (I,J)=(QNE (I,J)-QNE (I+IHW(J),J-1) &
                    +QSE (I,J)-QSE (I+IHW(J),J+1))*HDAC(I,J)*FCDIF
                    Q2DIF(I,J)=(Q2NE(I,J)-Q2NE(I+IHW(J),J-1) &
                    +Q2SE(I,J)-Q2SE(I+IHW(J),J+1))*HDAC(I,J)
            270 END DO
        !--------------2-ND ORDER DIFFUSION-------------------------------------
            IF(SECOND)THEN
                DO 280 J=MYJS2,MYJE2
                    DO 280 I=MYIS1,MYIE1
                        T (I,J,L)=T (I,J,L)+TDIF (I,J)
                        Q (I,J,L)=Q (I,J,L)+QDIF (I,J)
                280 END DO
            
            !-----------------------------------------------------------------------
                IF(L /= LM)THEN
                    DO 290 J=MYJS2,MYJE2
                        DO 290 I=MYIS1,MYIE1
                            Q2(I,J,L)=Q2(I,J,L)+Q2DIF(I,J)*HTM(I,J,L+1)
                    290 END DO
                ENDIF
            
                GO TO 360
            ENDIF
        !--------------4-TH ORDER DIAGONAL CONTRIBUTIONS------------------------
            DO 310 J=MYJS,MYJE1
                DO 310 I=MYIS,MYIE1
                    TNE (I,J)=(TDIF (I+IHE(J),J+1)-TDIF (I,J))*HKNE(I,J)
                    QNE (I,J)=(QDIF (I+IHE(J),J+1)-QDIF (I,J))*HKNE(I,J)
                    Q2NE(I,J)=(Q2DIF(I+IHE(J),J+1)-Q2DIF(I,J))*HKNE(I,J)
            310 END DO
        
            DO 320 J=MYJS1,MYJE
                DO 320 I=MYIS,MYIE1
                    TSE (I,J)=(TDIF (I+IHE(J),J-1)-TDIF (I,J))*HKSE(I,J)
                    QSE (I,J)=(QDIF (I+IHE(J),J-1)-QDIF (I,J))*HKSE(I,J)
                    Q2SE(I,J)=(Q2DIF(I+IHE(J),J-1)-Q2DIF(I,J))*HKSE(I,J)
            320 END DO
        !-----------------------------------------------------------------------
            DO 330 J=MYJS2,MYJE2
                DO 330 I=MYIS1,MYIE1
                    T(I,J,L)=T(I,J,L)-(TNE (I,J)-TNE (I+IHW(J),J-1) &
                    +TSE (I,J)-TSE (I+IHW(J),J+1))*HDAC(I,J)
                    Q(I,J,L)=Q(I,J,L)-(QNE (I,J)-QNE (I+IHW(J),J-1) &
                    +QSE (I,J)-QSE (I+IHW(J),J+1))*HDAC(I,J) &
                    *FCDIF
            330 END DO
        
        !-----------------------------------------------------------------------
            IF(L /= LM)THEN
                DO 340 J=MYJS2,MYJE2
                    DO 340 I=MYIS1,MYIE1
                        Q2(I,J,L)=Q2(I,J,L)-(Q2NE(I,J)-Q2NE(I+IHW(J),J-1) &
                        +Q2SE(I,J)-Q2SE(I+IHW(J),J+1))*HDAC(I,J) &
                        *HTM(I,J,L+1)
                340 END DO
            ENDIF
        !--------------U,V, DIAGONAL CONTRIBUTIONS------------------------------
            360 DO 410 J=MYJS_P1,MYJE1_P1
                DO 410 I=MYIS_P1,MYIE1_P1
                    VKNE(I,J)=(DEF(I+IVE(J),J)+DEF(I,J+1)) &
                    *VTM(I,J,L)*VTM(I+IVE(J),J+1,L) &
                    *HMASK(I+IVE(J),J)*HMASK(I,J+1)
                    UNE(I,J)=(U(I+IVE(J),J+1,L)-U(I,J,L))*VKNE(I,J)
                    VNE(I,J)=(V(I+IVE(J),J+1,L)-V(I,J,L))*VKNE(I,J)
            410 END DO
        
            DO 420 J=MYJS1_P1,MYJE_P1
                DO 420 I=MYIS_P1,MYIE1_P1
                    VKSE(I,J)=(DEF(I,J-1)+DEF(I+IVE(J),J)) &
                    *VTM(I+IVE(J),J-1,L)*VTM(I,J,L) &
                    *HMASK(I,J-1)*HMASK(I+IVE(J),J)
                    USE(I,J)=(U(I+IVE(J),J-1,L)-U(I,J,L))*VKSE(I,J)
                    VSE(I,J)=(V(I+IVE(J),J-1,L)-V(I,J,L))*VKSE(I,J)
            420 END DO
        !-----------------------------------------------------------------------
            DO 430 J=MYJS1,MYJE1
                DO 430 I=MYIS,MYIE1
                    UDIF(I,J)=(UNE(I,J)-UNE(I+IVW(J),J-1) &
                    +USE(I,J)-USE(I+IVW(J),J+1))*HDACV(I,J)
                    VDIF(I,J)=(VNE(I,J)-VNE(I+IVW(J),J-1) &
                    +VSE(I,J)-VSE(I+IVW(J),J+1))*HDACV(I,J)
            430 END DO
        !--------------2-ND ORDER DIFFUSION-------------------------------------
            IF(SECOND)THEN
                DO 440 J=MYJS2,MYJE2
                    DO 440 I=MYIS1,MYIE1
                        U(I,J,L)=U(I,J,L)+UDIF(I,J)
                        V(I,J,L)=V(I,J,L)+VDIF(I,J)
                440 END DO
            ELSE
            !  GO TO 500
            !      ENDIF
            !--------------4-TH ORDER DIAGONAL CONTRIBUTIONS------------------------
                DO 450 J=MYJS,MYJE1
                    DO 450 I=MYIS,MYIE1
                        UNE(I,J)=(UDIF(I+IVE(J),J+1)-UDIF(I,J))*VKNE(I,J)
                        VNE(I,J)=(VDIF(I+IVE(J),J+1)-VDIF(I,J))*VKNE(I,J)
                450 END DO
            
                DO 460 J=MYJS1,MYJE
                    DO 460 I=MYIS,MYIE1
                        USE(I,J)=(UDIF(I+IVE(J),J-1)-UDIF(I,J))*VKSE(I,J)
                        VSE(I,J)=(VDIF(I+IVE(J),J-1)-VDIF(I,J))*VKSE(I,J)
                460 END DO
            !-----------------------------------------------------------------------
                DO 470 J=MYJS2,MYJE2
                    DO 470 I=MYIS1,MYIE1
                        UTK=U(I,J,L)
                        VTK=V(I,J,L)
                        U(I,J,L)=U(I,J,L)-(UNE(I,J)-UNE(I+IVW(J),J-1) &
                        +USE(I,J)-USE(I+IVW(J),J+1))*HDACV(I,J)
                        V(I,J,L)=V(I,J,L)-(VNE(I,J)-VNE(I+IVW(J),J-1) &
                        +VSE(I,J)-VSE(I+IVW(J),J+1))*HDACV(I,J)
                        CKE(I,J)=0.5*(U(I,J,L)*U(I,J,L)-UTK*UTK &
                        +V(I,J,L)*V(I,J,L)-VTK*VTK)
                470 END DO
            !-----------------------------------------------------------------------
                IF(HEAT)THEN
                    DO 480 J=MYJS2,MYJE2
                        DO 480 I=MYIS1,MYIE1
                            T(I,J,L)=-RFCP*(CKE(I+IHE(J),J)+CKE(I,J+1) &
                            +CKE(I+IHW(J),J)+CKE(I,J-1))*HBM2(I,J) &
                            +T(I,J,L)
                    480 END DO
                ENDIF
            
            ENDIF
        !-----------------------------------------------------------------------
        500 END DO
    600 END DO
!-----------------------------------------------------------------------
    RETURN
    END SUBROUTINE HDIFF

