C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
                             SUBROUTINE HDIFFS
C
C     ******************************************************************
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    HDIFF       HORIZONTAL DIFFUSION
C   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 93-11-17
C     
C ABSTRACT:
C     HDIFF CALCULATES THE CONTRIBUTION OF THE HORIZONTAL DIFFUSION
C     TO THE TENDENCIES OF TEMPERATURE, SPECIFIC HUMIDITY, WIND
C     COMPONENTS, AND TURBULENT KINETIC ENERGY AND THEN UPDATES THOSE
C     VARIABLES.  A SECOND-ORDER NONLINEAR SCHEME SIMILAR TO
C     SMAGORINSKYS IS USED WHERE THE DIFFUSION COEFFICIENT IS
C     A FUNCTION OF THE DEFORMATION FIELD AND OF THE TURBULENT
C     KINETIC ENERGY.
C     
C PROGRAM HISTORY LOG:
C   87-06-??  JANJIC     - ORIGINATOR
C   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
C   96-03-28  BLACK      - ADDED EXTERNAL EDGE
C   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
C   13-09-09  MESINGER, SUEIRO, GOMES, MORELLI
C               Code changed to account for slopes (Lines marked by Cfm
C                 and CGCM) 
C               Changes due to diffusion limited not to change the sign
C                 of the Laplacian, in order to prevent instability due
C                 to large values of DEF  
C     
C USAGE: CALL HDIFF FROM MAIN PROGRAM EBU
C
C   INPUT ARGUMENT LIST:
C       NONE     
C  
C   OUTPUT ARGUMENT LIST: 
C     NONE
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C  
C     UNIQUE: NONE
C  
C     LIBRARY: NONE
C  
C   COMMON BLOCKS: CTLBLK
C                  MASKS
C                  PHYS
C                  VRBLS
C                  PVRBLS
C                  INDX
C   
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE : IBM SP
C$$$  
C***********************************************************************
Cfm Calculaton of hdiff v at points that have a neighboring "blocked" v
C-- switched on (loops 410 and 420), using velocities at points on
C-- slopes, but only half of the diffusion coefficient (note Fig. 2 of
C-- the "upgraded Eta" paper) 
C-- fm and Sandra Morelli, June-July 2013
C-----------------------------------------------------------------------
C     ******************************************************************
                             P A R A M E T E R
     & (DEFC=8.0,DEFM=32.0,SCQ2=50.0
     &, EPSQ2=0.2,FCDIF=1.0,RFCP=.25/1004.6)
C----------------------------------------------------------------------
      INCLUDE "parmeta"
      INCLUDE "parm.tbl"
      INCLUDE "mpp.h"
#include "sp.h"
C----------------------------------------------------------------------
                             P A R A M E T E R
     & (IMJM=IM*JM-JM/2,LP1=LM+1,KSMUD=1)
                             P A R A M E T E R
     &(JAM=6+2*(JM-10),JAMD=(JAM*2-10)*3)
C-----------------------------------------------------------------------
                             L O G I C A L
     & RUN,FIRST,RESTRT,SIGMA
     &,SECOND,HEAT,STTDF
C----------------------------------------------------------------------
      INCLUDE "CTLBLK.comm"
C-----------------------------------------------------------------------
      INCLUDE "MASKS.comm"
C-----------------------------------------------------------------------
      INCLUDE "PHYS.comm"
C-----------------------------------------------------------------------
      INCLUDE "VRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "PVRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "CLDWTR.comm"
C-----------------------------------------------------------------------
      INCLUDE "INDX.comm"
CGSM-----------------------------------------------------------------------
      INCLUDE "SLOPES.comm"
C-----------------------------------------------------------------------
                             D I M E N S I O N
     & 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)
C
                             D I M E N S I O N
     & 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)
C-----------------------------------------------------------------------
C***
C***  DIFFUSING Q2 AT GROUND LEVEL DOESNT MATTER, USTAR2 IS RECALCULATED
C***
C-----------------------------------------------------------------------
      SECOND=.TRUE.
      HEAT=.TRUE.
      LUL=UL(1)
C-----------------------------------------------------------------------
C---------------------MAIN VERTICAL INTEGRATION LOOP--------------------
C-----------------------------------------------------------------------
!$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)
C-----------------------------------------------------------------------
C
Cfm Fill the v points at slopes with the wind above --------------------
      IF (L.GT.1) THEN
       DO 205 J=MYJS_P1,MYJE_P1
       DO 205 I=MYIS_P1,MYIE_P1
        IF (VTMS(I,J,L).EQ.1) THEN
         U(I,J,L)=U(I,J,L-1)
         V(I,J,L)=V(I,J,L-1)
        ENDIF
  205 CONTINUE
      END IF
C-----------------------------------------------------------------------
                             DO 500 L=1,LM
C-----------------------------------------------------------------------
      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)
CGSM
      CALL ZERO2(CKE)
CGSM
      CALL ZERO2(TDIF)
      CALL ZERO2(QDIF)
      CALL ZERO2(UDIF)
      CALL ZERO2(VDIF)
      CALL ZERO2(Q2DIF)
C-----------------------------------------------------------------------
      DO 210 J=MYJS_P1,MYJE_P1
      DO 210 I=MYIS_P1,MYIE_P1
      Q2L(I,J)=AMAX1(Q2(I,J,L),EPSQ2)
  210 CONTINUE
C--------------DEFORMATIONS---------------------------------------------
      DO 220 J=MYJS1_P1,MYJE1_P1
      DO 220 I=MYIS_P1,MYIE1_P1
C
      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).gt.DEFS) then
         DEFS=DEF(I,J)
         IDefMax=I
         JDefMax=J
         LDefMax=L
         mypeDefMax=mype
       endif

Cfm
C   Maximazing DEF so as not to be less than DEFC, the line below,
C   that someone uncommented, led to the crash of a 1 km run at CPTEC.
C   Commenting the line back, as it was some years ago, removed the
C   problem.  Apparently, the line led to the diffusion coefficient
C   greater that it is allowed for stability.
C      DEF(I,J)=AMAX1(DEF(I,J),DEFC)
c      DEF(I,J)=AMIN1(DEF(I,J),DEFM)
CGSM      DEF(I,J)=DEF(I,J)*HMSKL(I,J) 
 220  CONTINUE
C--------------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))
     1          *HTM(I,J,L)*HTM(I+IHE(J),J+1,L)
CGSM     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 CONTINUE
C
      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))
     1          *HTM(I+IHE(J),J-1,L)*HTM(I,J,L)
CGSM     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 CONTINUE
C-----------------------------------------------------------------------
      DO 270 J=MYJS1,MYJE1
      DO 270 I=MYIS1,MYIE
      TDIF (I,J)=(TNE (I,J)-TNE (I+IHW(J),J-1)
     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)
     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)
     1           +Q2SE(I,J)-Q2SE(I+IHW(J),J+1))*HDAC(I,J)
  270 CONTINUE
C--------------2-ND ORDER DIFFUSION-------------------------------------
      IF(SECOND)THEN
        DO 280 J=MYJS2,MYJE2
        DO 280 I=MYIS1,MYIE1
C
Cfm Do not permit sign change due to horizontal diffusion.  This will
C     prevent instability due to large DEF as well.
            TFPAV=.25*(T(I+IHE(J),J+1,L)+T(I+IHW(J),J+1,L)
     1       +T(I+IHW(J),J-1,L)+T(I+IHE(J),J-1,L))
        IF (ABS(TDIF(I,J)).GT.ABS(TFPAV-T(I,J,L))) THEN
             TDIF(I,J)=TFPAV-T(I,J,L)
            ENDIF
C
        T (I,J,L)=T (I,J,L)+TDIF (I,J)
C
                   QFPAV=.25*(Q(I+IHE(J),J+1,L)+Q(I+IHW(J),J+1,L)
     1        +Q(I+IHW(J),J-1,L)+Q(I+IHE(J),J-1,L))
             IF (ABS(QDIF(I,J)).GT.ABS(QFPAV-Q(I,J,L))) THEN
             QDIF(I,J)=QFPAV-Q(I,J,L)
            ENDIF
C
        Q (I,J,L)=Q (I,J,L)+QDIF (I,J)
C
  280   CONTINUE
C
C-----------------------------------------------------------------------
        IF(L.NE.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     CONTINUE
        ENDIF
C
        GO TO 360
      ENDIF
C--------------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 CONTINUE
C
      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 CONTINUE
C-----------------------------------------------------------------------
      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)
     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)
     1                  +QSE (I,J)-QSE (I+IHW(J),J+1))*HDAC(I,J)
     2                  *FCDIF
  330 CONTINUE
C
C-----------------------------------------------------------------------
      IF(L.NE.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)
     1                      +Q2SE(I,J)-Q2SE(I+IHW(J),J+1))*HDAC(I,J)
     2                                                  *HTM(I,J,L+1)
  340   CONTINUE
      ENDIF
C--------------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))
CGSM     1          *VTM(I,J,L)*VTM(I+IVE(J),J+1,L)
CGSM     2          *HMASK(I+IVE(J),J)*HMASK(I,J+1)
     1          *MAX(VTM(I+IVE(J),J+1,L),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 CONTINUE
C
      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))
CGSM     1          *VTM(I+IVE(J),J-1,L)*VTM(I,J,L)
CGSM     2          *HMASK(I,J-1)*HMASK(I+IVE(J),J)
     1          *MAX(VTM(I+IVE(J),J-1,L),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 CONTINUE
C-----------------------------------------------------------------------
      DO 430 J=MYJS1,MYJE1
      DO 430 I=MYIS,MYIE1
      UDIF(I,J)=(UNE(I,J)-UNE(I+IVW(J),J-1)
     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)
     1          +VSE(I,J)-VSE(I+IVW(J),J+1))*HDACV(I,J)
  430 CONTINUE
C--------------2-ND ORDER DIFFUSION-------------------------------------
      IF(SECOND)THEN
        DO 440 J=MYJS2,MYJE2
        DO 440 I=MYIS1,MYIE1
CGSM        U(I,J,L)=U(I,J,L)+UDIF(I,J)
CGSM        V(I,J,L)=V(I,J,L)+VDIF(I,J)
C  
Cfm Do not permit sign change due to horizontal diffusion. This will
C     prevent instability due to large DEF as well.
                  UFPAV=.25*(U(I+IVE(J),J+1,L)+U(I+IVW(J),J+1,L)
     1        +U(I+IVW(J),J-1,L)+U(I+IVE(J),J-1,L))
              IF (ABS(UDIF(I,J)).GT.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)
     1                +V(I+IVW(J),J-1,L)+V(I+IVE(J),J-1,L))
            IF (ABS(VDIF(I,J)).GT.ABS(VFPAV-V(I,J,L))) THEN
                VDIF(I,J)=VFPAV-V(I,J,L)
            ENDIF
C
      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)
     1           *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
     1             +V(I,J,L)*V(I,J,L)-VTK*VTK)
Cfm
C   Let individual CKE be positive, but not the four point average used 
C     to increase T
C        IF(CKE(I,J).GT.0)THEN
C          CKE(I,J)=0.
Cfm      ENDIF
C
  440   CONTINUE
      ELSE
c  GO TO 500
c      ENDIF
C--------------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 CONTINUE
C
      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 CONTINUE
C-----------------------------------------------------------------------
      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)
     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)
     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
     1             +V(I,J,L)*V(I,J,L)-VTK*VTK)
  470 CONTINUE
      ENDIF
C-----------------------------------------------------------------------
       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)
           IF ( ABS(CKET).gt.10) then
											write(201,*) I,J,L,CKET,
     &         CKE(I+IHE(J),J),CKE(I,J+1),CKE(I+IHW(J),J),CKE(I,J-1)       
           ENDIF
              IF (CKET.GT.0.) CKET=0.
              T(I,J,L)=-RFCP*CKET*HBM2(I,J) +T(I,J,L)        
  480    CONTINUE
       ENDIF
C-----------------------------------------------------------------------
Cfm Fill the v points at slopes back with zeros ------------------------
      IF (L.GT.1) THEN
      DO 485 J=MYJS_P1,MYJE_P1
       DO 485 I=MYIS_P1,MYIE_P1
        IF (VTMS(I,J,L).EQ.1) THEN
           U(I,J,L)=0.
           V(I,J,L)=0.
        ENDIF
  485   CONTINUE
       ENDIF
C-----------------------------------------------------------------------
  500                        CONTINUE
C       write(200+mypeDefMax,*) mypeDefMax,NTSD,IDefMax,JDefMax,LDefMax,DEFS
C-----------------------------------------------------------------------
                             RETURN
                             END
