SUBROUTINE hdiff
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !                .      .    .     
  ! subprogram:    hdiff       horizontal diffusion
  !   prgrmmr: janjic          org: w/np2      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
  !     smagorinsky's 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
  !     
  ! 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
  !                  scrch
  !                  vrbls
  !                  pvrbls
  !                  indx
  !   
  IMPLICIT NONE

  !*** problem size

  INCLUDE "parmeta.f90"
  INCLUDE "parm.tbl.f90"
  INTEGER, PARAMETER :: lp1=lm+1

  !*** global variables

  INCLUDE "masks.f90"   ! hbm2(in); htm(in); vtm(in)
  INCLUDE "phys.f90"    ! hdac(in); hdacv(in)
  INCLUDE "vrbls.f90"   ! t(inout); q(inout); u(inout); v(inout)
  INCLUDE "pvrbls.f90"  ! q2(inout)
  INCLUDE "indx.f90"    ! ihe(in); ihw(in); ive(in); ivw(in)

  !*** local variables

  INTEGER, PARAMETER :: ksmud=01
  REAL   , PARAMETER :: scq2=050.e0
  REAL   , PARAMETER :: epsq2=0.2
  REAL   , PARAMETER :: fcdif=1.0e0
  REAL   , PARAMETER :: d50=.5e0
  REAL   , PARAMETER :: rfcp=.25e0/1004.6e0
  LOGICAL, PARAMETER :: second=.TRUE.
  LOGICAL, PARAMETER :: heat=.FALSE.
  REAL   , DIMENSION(im,jm) :: q2l
  REAL   , DIMENSION(im,jm) :: hkne
  REAL   , DIMENSION(im,jm) :: hkse
  REAL   , DIMENSION(im,jm) :: vkne
  REAL   , DIMENSION(im,jm) :: vkse
  REAL   , DIMENSION(0:im+1,0:jm+1) :: tne
  REAL   , DIMENSION(0:im+1,0:jm+1) :: tse
  REAL   , DIMENSION(0:im+1,0:jm+1) :: qne
  REAL   , DIMENSION(0:im+1,0:jm+1) :: qse
  REAL   , DIMENSION(0:im+1,0:jm+1) :: q2ne
  REAL   , DIMENSION(0:im+1,0:jm+1) :: q2se
  REAL   , DIMENSION(0:im+1,0:jm+1) :: une
  REAL   , DIMENSION(0:im+1,0:jm+1) :: use
  REAL   , DIMENSION(0:im+1,0:jm+1) :: vne
  REAL   , DIMENSION(0:im+1,0:jm+1) :: vse
  REAL   , DIMENSION(0:im+1,0:jm+1) :: tdif
  REAL   , DIMENSION(0:im+1,0:jm+1) :: qdif
  REAL   , DIMENSION(0:im+1,0:jm+1) :: udif
  REAL   , DIMENSION(0:im+1,0:jm+1) :: vdif
  REAL   , DIMENSION(0:im+1,0:jm+1) :: q2dif
  REAL   , DIMENSION(0:im+1,0:jm+1) :: def
  REAL   , DIMENSION(0:im+1,0:jm+1) :: cke
  INTEGER :: i, j, ks, l
  REAL    :: deftk
  REAL    :: defsk
  REAL    :: vtk
  REAL    :: utk

  DO ks=1,ksmud

     !*** main vertical integration loop

     DO l=1,lm
        !*** diffusing q2 at ground level doesn't matter;
        !*** ustar2 is recalculated
        CALL zero2(tne)
        CALL zero2(tse)
        CALL zero2(qne)
        CALL zero2(qse)
        CALL zero2(q2ne)
        CALL zero2(q2se)
        CALL zero2(une)
        CALL zero2(use)
        CALL zero2(vne)
        CALL zero2(vse)
        CALL zero2(tdif)
        CALL zero2(qdif)
        CALL zero2(udif)
        CALL zero2(vdif)
        CALL zero2(q2dif)
        CALL zero2(def)
        CALL zero2(cke)
        DO j=1,jm
           DO i=1,im
              q2l(i,j)=amax1(q2(i,j,l),epsq2)
           END DO
        END DO

        !*** deformations

        DO j=2,jm-1
           DO i=1,im-1
              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  +scq2*q2l(i,j)
              def (i,j)=SQRT(def(i,j)+def(i,j))*hbm2(i,j)
           END DO
        END DO

        !*** t, q, q2 diagonal contributions

        DO j=1,jm-1
           DO i=1,im-1
              hkne(i,j)=(def(i,j)+def(i+ihe(j),j+1)) &
                   *htm(i,j,l)*htm(i+ihe(j),j+1,l)
              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)
           END DO
        END DO
        DO j=2,jm
           DO i=1,im-1
              hkse(i,j)=(def(i+ihe(j),j-1)+def(i,j)) &
                   *htm(i+ihe(j),j-1,l)*htm(i,j,l)
              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)
           END DO
        END DO
        DO j=2,jm-1
           DO i=2,im
              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)
           END DO
        END DO

        !*** t, q, q2 diffusion

        IF(second)THEN

           !*** 2-nd order

           DO j=3,jm-2
              DO i=2,im-1
                 t (i,j,l)=t (i,j,l)+tdif (i,j)
                 q (i,j,l)=q (i,j,l)+qdif (i,j)
              END DO
           END DO
           DO j=3,jm-2
              DO i=2,im-1
                 q2(i,j,l)=q2(i,j,l)+q2dif(i,j)
              END DO
           END DO
        ELSE

           !*** 4-th order

           DO j=1,jm-1
              DO i=1,im-1
                 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)
              END DO
           END DO

           DO j=2,jm
              DO i=1,im-1
                 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)
              END DO
           END DO

           DO j=3,jm-2
              DO i=2,im-1
                 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
              END DO
           END DO
           DO j=3,jm-2
              DO i=2,im-1
                 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)
              END DO
           END DO
        ENDIF

        !*** u,v, diagonal contributions

        DO j=1,jm-1
           DO i=1,im-1
              vkne(i,j)=(def(i+ive(j),j)+def(i,j+1)) &
                   *vtm(i,j,l)*vtm(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)
           END DO
        END DO
        DO j=2,jm
           DO i=1,im-1
              vkse(i,j)=(def(i,j-1)+def(i+ive(j),j)) &
                   *vtm(i+ive(j),j-1,l)*vtm(i,j,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)
           END DO
        END DO
        DO j=2,jm-1
           DO i=1,im-1
              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)
           END DO
        END DO

        !*** u, v, diffusion

        IF(second)THEN

           !*** 2-nd order

           DO j=3,jm-2
              DO i=2,im-1
                 u(i,j,l)=u(i,j,l)+udif(i,j)
                 v(i,j,l)=v(i,j,l)+vdif(i,j)
              END DO
           END DO
        ELSE

           !*** 4-th order

           DO j=1,jm-1
              DO i=1,im-1
                 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)
              END DO
           END DO
           DO j=2,jm
              DO i=1,im-1
                 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)
              END DO
           END DO
           DO j=3,jm-2
              DO i=2,im-1
                 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)=d50*(u(i,j,l)*u(i,j,l)-utk*utk &
                      +v(i,j,l)*v(i,j,l)-vtk*vtk)
              END DO
           END DO
           IF(heat)THEN
              DO j=3,jm-2
                 DO i=2,im-1
                    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)
                 END DO
              END DO
           ENDIF
        ENDIF
     END DO
  END DO
END SUBROUTINE hdiff
