!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    SUBROUTINE HDIFFS

!     ******************************************************************
!$$$  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
!   13-09-09  MESINGER, SUEIRO, GOMES, MORELLI
!               Code changed to account for slopes (Lines marked by Cfm
!                 and CGCM)
!               Changes due to diffusion limited not to change the sign
!                 of the Laplacian, in order to prevent instability due
!                 to large values of DEF

! 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
!$$$
!***********************************************************************
! m Calculaton of hdiff v at points that have a neighboring "blocked" v
!-- switched on (loops 410 and 420), using velocities at points on
!-- slopes, but only half of the diffusion coefficient (note Fig. 2 of
!-- the "upgraded Eta" paper)
!-- fm and Sandra Morelli, June-July 2013
!-----------------------------------------------------------------------
!     ******************************************************************
    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_PHYS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_VRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_PVRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_CLDWTR.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_INDX.f90"
! SM-----------------------------------------------------------------------
    INCLUDE "COMM_SLOPES.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= .TRUE. 
    LUL=UL(1)
!-----------------------------------------------------------------------
! m Fill the v points at slopes with the wind above --------------------
    DO 205 L=2,LM
        DO 205 J=MYJS_P1,MYJE_P1
            DO 205 I=MYIS_P1,MYIE_P1
                IF (VTMS(I,J,L) == 1) THEN
                    U(I,J,L)=U(I,J,L-1)
                    V(I,J,L)=V(I,J,L-1)
                ENDIF
    205 END DO
!-----------------------------------------------------------------------
!---------------------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)
    ! SM
        CALL ZERO2(CKE)
    ! SM
        CALL ZERO2(TDIF)
        CALL ZERO2(QDIF)
        CALL ZERO2(UDIF)
        CALL ZERO2(VDIF)
        CALL ZERO2(Q2DIF)
        DEFS=-1
    !-----------------------------------------------------------------------
        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
            
                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)
                DEF (I,J)=DEFTK*DEFTK+DEFSK*DEFSK
                DEF (I,J)=SQRT(DEF(I,J)+DEF(I,J))*HBM2(I,J)

                if (DEF(I,J) > DEFS) then
                    DEFS=DEF(I,J)
                    IDefMax=I
                    JDefMax=J
                    LDefMax=L
                    mypeDefMax=mype
                endif

            ! m
            !   Maximazing DEF so as not to be less than DEFC, the line below,
            !   that someone uncommented, led to the crash of a 1 km run at CPTEC.
            !   Commenting the line back, as it was some years ago, removed the
            !   problem.  Apparently, the line led to the diffusion coefficient
            !   greater that it is allowed for stability.
            !      DEF(I,J)=AMAX1(DEF(I,J),DEFC)
            !      DEF(I,J)=AMIN1(DEF(I,J),DEFM)
            ! SM      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)
            ! SM     2          *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)
            ! SM     2          *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
                
                ! m Do not permit sign change due to horizontal diffusion.  This will
                !     prevent instability due to large DEF as well.
                    TFPAV=.25*(T(I+IHE(J),J+1,L)+T(I+IHW(J),J+1,L) &
                    +T(I+IHW(J),J-1,L)+T(I+IHE(J),J-1,L))
                    IF (ABS(TDIF(I,J)) > ABS(TFPAV-T(I,J,L))) THEN
                        TDIF(I,J)=TFPAV-T(I,J,L)
                    ENDIF
                
                    T (I,J,L)=T (I,J,L)+TDIF (I,J)
                
                    QFPAV=.25*(Q(I+IHE(J),J+1,L)+Q(I+IHW(J),J+1,L) &
                    +Q(I+IHW(J),J-1,L)+Q(I+IHE(J),J-1,L))
                    IF (ABS(QDIF(I,J)) > ABS(QFPAV-Q(I,J,L))) THEN
                        QDIF(I,J)=QFPAV-Q(I,J,L)
                    ENDIF
                
                    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
                    ! GSM          Q2(I,J,L)=Q2(I,J,L)+Q2DIF(I,J)*HTM(I,J,L+1)
                        Q2(I,J,L)=AMAX1(Q2(I,J,L),EPSQ2)
                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))*MAX(VTM(I+IVE(J),J+1,L), &
            ! SM     1          *VTM(I,J,L)*VTM(I+IVE(J),J+1,L)
            ! SM     2          *HMASK(I+IVE(J),J)*HMASK(I,J+1)
			    VTMS(I+IVE(J),J+1,L))
                                 
                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))*MAX(VTM(I+IVE(J),J-1,L), &
            ! SM     1          *VTM(I+IVE(J),J-1,L)*VTM(I,J,L)
            ! SM     2          *HMASK(I,J-1)*HMASK(I+IVE(J),J)
                VTMS(I+IVE(J),J-1,L))
                      
                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
                ! SM        U(I,J,L)=U(I,J,L)+UDIF(I,J)
                ! SM        V(I,J,L)=V(I,J,L)+VDIF(I,J)
                
                ! m Do not permit sign change due to horizontal diffusion. This will
                !     prevent instability due to large DEF as well.
                    UFPAV=.25*(U(I+IVE(J),J+1,L)+U(I+IVW(J),J+1,L) &
                    +U(I+IVW(J),J-1,L)+U(I+IVE(J),J-1,L))
                    IF (ABS(UDIF(I,J)) > ABS(UFPAV-U(I,J,L))) THEN
                        UDIF(I,J)=UFPAV-U(I,J,L)
                    ENDIF
                    VFPAV=.25*(V(I+IVE(J),J+1,L)+V(I+IVW(J),J+1,L) &
                    +V(I+IVW(J),J-1,L)+V(I+IVE(J),J-1,L))
                    IF (ABS(VDIF(I,J)) > ABS(VFPAV-V(I,J,L))) THEN
                        VDIF(I,J)=VFPAV-V(I,J,L)
                    ENDIF
                
                    UTK=U(I,J,L)
                    VTK=V(I,J,L)
                    NNTMP=  VTM(I+IVW(J),J+1,L)*VTM(I+IVE(J),J+1,L) &
                    *VTM(I+IVW(J),J-1,L)*VTM(I+IVE(J),J-1,L)
                    DCMD=NNTMP+0.5*MOD(NNTMP+1,2)
                    U(I,J,L)=U(I,J,L)+UDIF(I,J)*VTM(I,J,L)*DCMD
                    V(I,J,L)=V(I,J,L)+VDIF(I,J)*VTM(I,J,L)*DCMD
                    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)
                ! m
                !   Let individual CKE be positive, but not the four point average used
                !     to increase T
                !        IF(CKE(I,J).GT.0)THEN
                !          CKE(I,J)=0.
                ! m      ENDIF
                
            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
        ENDIF
    !-----------------------------------------------------------------------
        IF(HEAT)THEN
            DO 480 J=MYJS2,MYJE2
                DO 480 I=MYIS1,MYIE1
                    CKET=CKE(I+IHE(J),J)+CKE(I,J+1)+CKE(I+IHW(J),J)+CKE(I,J-1)
                ! SM           IF ( ABS(CKET).gt.10) then
                ! SM               write(201,*) I,J,L,CKET,
                ! SM     &         CKE(I+IHE(J),J),CKE(I,J+1),CKE(I+IHW(J),J),CKE(I,J-1)
                ! SM           ENDIF
                    IF (CKET > 0.) CKET=0.
                    T(I,J,L)=-RFCP*CKET*HBM2(I,J) +T(I,J,L)
            480 END DO
        ENDIF
    !-----------------------------------------------------------------------
    ! m Fill the v points at slopes back with zeros ------------------------
        IF (L > 1) THEN
            DO 485 J=MYJS_P1,MYJE_P1
                DO 485 I=MYIS_P1,MYIE_P1
                    IF (VTMS(I,J,L) == 1) THEN
                        U(I,J,L)=0.
                        V(I,J,L)=0.
                    ENDIF
            485 END DO
        ENDIF
    !-----------------------------------------------------------------------
    500 END DO
    write(200+mypeDefMax,*) mypeDefMax,NTSD,IDefMax,JDefMax,LDefMax,DEFS
!-----------------------------------------------------------------------
    RETURN
    END SUBROUTINE HDIFFS
