SUBROUTINE hzadv
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  ! subprogram:    hzadv       horizontal advection
  !   prgrmmr: janjic          org: w/np2      date: 93-10-28
  ! abstract:
  !     hzadv calculates the contribution of the horizontal advection
  !     to the tendencies of temperature, wind components, and
  !     turbulent kinetic energy and then updates those variables.
  !     the janjic advection scheme for the arakawa e grid is used
  !     for all variables inside the fifth row.  an upstream scheme
  !     is used on all variables in the third, fourth, and fifth
  !     outermost rows.  a modified euler-backward time scheme (heun)
  !     is used.
  ! 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

  IMPLICIT NONE
  INCLUDE "parmeta.f90"

  REAL,PARAMETER :: d00=.0e0,h1=1.e0,d99=.99e0,d10=.1e0,d50=.50
  REAL,PARAMETER :: tlc=1.414213562e0,h4=4.e0,d0625=.0625e0,epsq=-2.e-12
  REAL,PARAMETER :: epsq2=0.2
  INTEGER,PARAMETER :: lda=lm+9,la=13,im1=im-1,jam=6+2*(jm-10),imjm=im*jm-jm/2
  INTEGER,PARAMETER :: lp1=lm+1,jamd=(jam*2-10)*3

  INCLUDE "ctlblk.f90"
  INCLUDE "loops.f90"
  INCLUDE "masks.f90"
  INCLUDE "dynam.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "contin.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "indx.f90"

  REAL,DIMENSION(im,jm) :: hm,vm,rdpd,adpdx,adpdy ,rdpdx ,rdpdy,adu,adv,adt
  REAL,DIMENSION(im,jm) :: adq2m,adq2l,q2mns ,q2lns,udy,vdx   
  REAL,DIMENSION(0:im+1,0:jm+1) :: dpde,temp,tst,ust,vst,q2m,q2l,tew,tns
  REAL,DIMENSION(0:im+1,0:jm+1) :: q2mew,q2lew,tne,tse,q2mne,q2mse,q2lne,q2lse
  REAL,DIMENSION(0:im+1,0:jm+1) :: uew,uns,vew,vns,une,USE,vne,vse,few,fns 
  REAL,DIMENSION(0:im+1,0:jm+1) :: fne,fse  
  REAL,DIMENSION(im,jm,lm) :: adq2hl
  REAL,DIMENSION(im,jm,lm+1) :: q2ml
  REAL,DIMENSION(jamd) :: array0,array1,array2,array3

  INTEGER,DIMENSION(jamd) :: khhas,ihlas,jhlas,kvhas,ivlas,jvlas,ispa,isqa

  INTEGER :: i,j,l,jakone,ja,ihl,ihh,jaktwo,ivl,ivh,isp,isq,ifp,ifq,jak,ipq
  REAL :: f0,f1,f2,f3,tta,ttb,pp,qp
  LOGICAL :: iter2

  !note: for the eta model with tke initialization
  !the following section should be removed
  IF(ntsd == 1)THEN
     DO  l=1,lm
        DO  j=1,jm
           DO  i=1,im
              q2(i,j,l)=amax1(q2(i,j,l)*hbm2(i,j),epsq2)
           END DO
        END DO
     END DO
  END IF
  jakone=0   ! (this need only be done once.)
  DO  ja=1,jam
     ihl=ihla(ja)
     ihh=ihha(ja)
     j=jra(ja)
     DO  i=ihl,ihh
        jakone=jakone+1
        khhas(jakone)=ja
        ihlas(jakone)=i
        jhlas(jakone)=j
     END DO
  END DO
  jaktwo=0
  DO  ja=1,jam
     ivl=ivla(ja)
     ivh=ivha(ja)
     j=jra(ja)
     DO  i=ivl,ivh
        jaktwo=jaktwo+1
        kvhas(jaktwo)=ja
        ivlas(jaktwo)=i
        jvlas(jaktwo)=j
     END DO
  END DO
  DO  j=1,jm
     DO  i=1,im
        q2ml(i,j,1)=d00
     END DO
  END DO
  DO  l=2,lm+1
     DO  j=1,jm
        DO  i=1,im
           q2ml(i,j,l)=q2(i,j,l-1)
        END DO
     END DO
  END DO
  DO  l=1,lm
     CALL zero2(temp)
     CALL zero2(dpde)
     CALL zero2(tst)
     CALL zero2(ust)
     CALL zero2(vst)
     CALL zero2(q2m)
     CALL zero2(q2l)
     CALL zero2(tew)
     CALL zero2(tns)
     CALL zero2(q2mew)
     CALL zero2(q2lew)
     CALL zero2(tne)
     CALL zero2(tse)
     CALL zero2(q2mne)
     CALL zero2(q2mse)
     CALL zero2(q2lne)
     CALL zero2(q2lse)
     CALL zero2(uew)
     CALL zero2(uns)
     CALL zero2(vew)
     CALL zero2(vns)
     CALL zero2(une)
     CALL zero2(USE)
     CALL zero2(vne)
     CALL zero2(vse)
     CALL zero2(few)
     CALL zero2(fns)
     CALL zero2(fne)
     CALL zero2(fse)
     iter2=.FALSE.
     DO  j=1,jm
        DO  i=1,im
           q2m(i,j)=d00
        END DO
     END DO
     DO  j=1,jm
        DO  i=1,im
           q2m(i,j)=q2ml(i,j,l)
        END DO
     END DO
     DO  j=1,jm
        DO  i=1,im
           hm(i,j)=htm(i,j,l)*hbm2(i,j)
           dpde(i,j)=pdsl(i,j)*deta(l)
           rdpd(i,j)=h1/dpde(i,j)
           ust(i,j)=u(i,j,l)
           vst(i,j)=v(i,j,l)
           tst(i,j)=t(i,j,l)
           q2l(i,j)=q2ml(i,j,l+1)
        END DO
     END DO
     DO  j=2,jm-1
        DO  i=1,im
           vm(i,j)=vtm(i,j,l)*vbm2(i,j)
           adpdx(i,j)=dpde(i+ivw(j),j)+dpde(i+ive(j),j)
           adpdy(i,j)=dpde(i,j-1)+dpde(i,j+1)
           rdpdx(i,j)=h1/adpdx(i,j)
           rdpdy(i,j)=h1/adpdy(i,j)
        END DO
     END DO
     !mass fluxes and mass points advection components---------
     !the ns and ew fluxes in the 130 and 140 loops are on v points
     DO 
        DO  j=2,jm-1
           DO  i=1,im
              udy(i,j)=ust(i,j)*dy
              few(i,j)=udy(i,j)*adpdx(i,j)
              tew(i,j)=few(i,j)*(tst(i+ive(j),j)-tst(i+ivw(j),j))
              q2mew(i,j)=few(i,j)*(q2m(i+ive(j),j)-q2m(i+ivw(j),j))
              q2lew(i,j)=few(i,j)*(q2l(i+ive(j),j)-q2l(i+ivw(j),j))
           END DO
        END DO
        DO  j=2,jm-1
           DO  i=1,im
              vdx(i,j)=vst(i,j)*dx(i,j)
              fns(i,j)=vdx(i,j)*adpdy(i,j)
              tns(i,j)=fns(i,j)*(tst(i,j+1)-tst(i,j-1))
              q2mns(i,j)=fns(i,j)*(q2m(i,j+1)-q2m(i,j-1))
              q2lns(i,j)=fns(i,j)*(q2l(i,j+1)-q2l(i,j-1))
           END DO
        END DO
        !diagonal fluxes and diagonally averaged wind-------------
        !the ne and se fluxes in the 150 and 160 loops are on h points
        !(actually just to the ne and se of each h point)
        DO  j=3,jm-2
           DO  i=1,im
              temp(i,j)=udy(i,j)+vdx(i,j)
           END DO
        END DO
        DO  j=3,jm-3
           DO  i=1,im
              fne(i,j)=(temp(i+ihe(j),j)+temp(i,j+1)) *(dpde(i,j)+dpde(i+ihe(j),j+1))
              tne(i,j)=fne(i,j)*(tst(i+ihe(j),j+1)-tst(i,j))
              q2mne(i,j)=fne(i,j)*(q2m(i+ihe(j),j+1)-q2m(i,j))
              q2lne(i,j)=fne(i,j)*(q2l(i+ihe(j),j+1)-q2l(i,j))
           END DO
        END DO
        DO  j=3,jm-2
           DO  i=1,im
              temp(i,j)=udy(i,j)-vdx(i,j)
           END DO
        END DO
        DO  j=4,jm-2
           DO  i=1,im
              fse(i,j)=(temp(i+ihe(j),j)+temp(i,j-1)) *(dpde(i,j)+dpde(i+ihe(j),j-1))
              tse(i,j)=fse(i,j)*(tst(i+ihe(j),j-1)-tst(i,j))
              q2mse(i,j)=fse(i,j)*(q2m(i+ihe(j),j-1)-q2m(i,j))
              q2lse(i,j)=fse(i,j)*(q2l(i+ihe(j),j-1)-q2l(i,j))
           END DO
        END DO
        !thermodynamic equation & moisture
        !the ad arrays in the 170 loop are on h points
        DO  j=6,jm-5
           DO  i=1,im
              adt(i,j)=(tew(i+ihw(j),j)+tew(i+ihe(j),j)+tns(i,j-1)+tns(i,j+1)  &
                   +tne(i+ihw(j),j-1)+tne(i,j)+tse(i,j)+tse(i+ihw(j),j+1))  &
                   *rdpd(i,j)*fad(i,j)
              adq2m(i,j)=(q2mew(i+ihw(j),j)+q2mew(i+ihe(j),j)  &
                   +q2mns(i,j-1)+q2mns(i,j+1) +q2mne(i+ihw(j),j-1)+q2mne(i,j)  &
                   +q2mse(i,j)+q2mse(i+ihw(j),j+1)) *rdpd(i,j)*fad(i,j)
              adq2l(i,j)=(q2lew(i+ihw(j),j)+q2lew(i+ihe(j),j)  &
                   +q2lns(i,j-1)+q2lns(i,j+1) +q2lne(i+ihw(j),j-1)+q2lne(i,j)  &
                   +q2lse(i,j)+q2lse(i+ihw(j),j+1)) *rdpd(i,j)*fad(i,j)
           END DO
        END DO
        !upstream advection of t, q and q2------------------------
        DO  jak=1,jakone
           ja =khhas(jak)
           i  =ihlas(jak)
           j  =jhlas(jak)
           tta=emt(ja)*(ust(i,j-1)+ust(i+ihw(j),j) +ust(i+ihe(j),j)+ust(i,j+1))
           ttb=ent    *(vst(i,j-1)+vst(i+ihw(j),j) +vst(i+ihe(j),j)+vst(i,j+1))
           pp =-tta-ttb
           qp =tta-ttb
           ispa(jak)=INT(SIGN(h1,pp))
           isqa(jak)=INT(SIGN(h1,qp))
           pp       =ABS(pp)
           qp       =ABS(qp)
           array3(jak)=pp*qp
           array0(jak)=array3(jak)-pp-qp
           array1(jak)=pp-array3(jak)
           array2(jak)=qp-array3(jak)
        END DO
        jak=0
        DO  ja=1,jam
           ihl=ihla(ja)
           ihh=ihha(ja)
           j=jra(ja)
           DO  i=ihl,ihh
              jak=jak+1
              isp=ispa(jak)
              isq=isqa(jak)
              ifp=(isp-1)/2
              ifq=(-isq-1)/2
              ipq=(isp-isq)/2
              IF(htm(i+ihe(j)+ifp,j+isp,l)*htm(i+ihe(j)+ifq,j+isq,l)  &
                   *htm(i+ipq,j+isp+isq,l) <= d10) THEN
                 IF(htm(i+ihe(j)+ifp,j+isp,l)+htm(i+ihe(j)+ifq,j+isq,l)  &
                      +htm(i+ipq,j+isp+isq,l) < d10)THEN
                    tst(i+ihe(j)+ifp,j+isp)=tst(i,j)
                    tst(i+ihe(j)+ifq,j+isq)=tst(i,j)
                    tst(i+ipq,j+isp+isq)   =tst(i,j)
                 ELSE IF (htm(i+ihe(j)+ifp,j+isp,l)+htm(i+ipq,j+isp+isq,l)<d99) THEN
                    tst(i+ihe(j)+ifp,j+isp)=tst(i,j)
                    tst(i+ipq,j+isp+isq)   =tst(i+ihe(j)+ifq,j+isq)
                 ELSE IF (htm(i+ihe(j)+ifq,j+isq,l)+htm(i+ipq,j+isp+isq,l)<d99)THEN
                    tst(i+ihe(j)+ifq,j+isq)=tst(i,j)
                    tst(i+ipq,j+isp+isq)   =tst(i+ihe(j)+ifp,j+isp)
                 ELSE IF (htm(i+ihe(j)+ifp,j+isp,l)+htm(i+ihe(j)+ifq,j+isq,l)<d99) THEN
                    tst(i+ihe(j)+ifp,j+isp)= d50*(tst(i,j)+tst(i+ipq,j+isp+isq))
                    tst(i+ihe(j)+ifq,j+isq)=tst(i+ihe(j)+ifp,j+isp)
                 ELSE IF (htm(i+ihe(j)+ifp,j+isp,l) < d99)THEN
                    tst(i+ihe(j)+ifp,j+isp)=  &
                         tst(i,j)+tst(i+ipq,j+isp+isq)-tst(i+ihe(j)+ifq,j+isq)
                 ELSE IF (htm(i+ihe(j)+ifq,j+isq,l) < d99)THEN
                    tst(i+ihe(j)+ifq,j+isq)=  &
                         tst(i,j)+tst(i+ipq,j+isp+isq)-tst(i+ihe(j)+ifp,j+isp)
                 ELSE
                    tst(i+ipq,j+isp+isq)=tst(i+ihe(j)+ifp,j+isp)  &
                         +tst(i+ihe(j)+ifq,j+isq)-tst(i,j)
                 END IF
              END IF
              f0=array0(jak)
              f1=array1(jak)
              f2=array2(jak)
              f3=array3(jak)
              adt(i,j)=f0*tst(i,j)+f1*tst(i+ihe(j)+ifp,j+isp)  &
                   +f2*tst(i+ihe(j)+ifq,j+isq)+f3*tst(i+ipq,j+isp+isq)
           END DO
        END DO
        DO  jak=1,jakone
           i=ihlas(jak)
           j=jhlas(jak)
           isp=ispa(jak)
           isq=isqa(jak)
           ifp=(isp-1)/2
           ifq=(-isq-1)/2
           ipq=(isp-isq)/2
           f0=array0(jak)
           f1=array1(jak)
           f2=array2(jak)
           f3=array3(jak)
           adq2m(i,j)=f0*q2m(i,j)+f1*q2m(i+ihe(j)+ifp,j+isp)  &
                +f2*q2m(i+ihe(j)+ifq,j+isq)+f3*q2m(i+ipq,j+isp+isq)
           adq2l(i,j)=f0*q2l(i,j)+f1*q2l(i+ihe(j)+ifp,j+isp)  &
                +f2*q2l(i+ihe(j)+ifq,j+isq)+f3*q2l(i+ipq,j+isp+isq)
        END DO
        !calculation of momentum advection components-------------
        !the ew and ns arrays in the 180 loop are on h points
        DO  j=5,jm-4
           DO  i=1,im
              uew(i,j)=(few(i+ihw(j),j)+few(i+ihe(j),j))  &
                   *(ust(i+ihe(j),j)-ust(i+ihw(j),j))
              uns(i,j)=(fns(i+ihw(j),j)+fns(i+ihe(j),j)) *(ust(i,j+1)-ust(i,j-1))
              vew(i,j)=(few(i,j-1)+few(i,j+1)) *(vst(i+ihe(j),j)-vst(i+ihw(j),j))
              vns(i,j)=(fns(i,j-1)+fns(i,j+1))*(vst(i,j+1)-vst(i,j-1))
           END DO
        END DO
        !the ne and se arrays in the 190 loop are tied to v points
        DO  j=5,jm-4
           DO  i=1,im
              une(i,j)=(fne(i+ivw(j),j)+fne(i+ive(j),j)) *(ust(i+ive(j),j+1)-ust(i,j))
              USE(i,j)=(fse(i+ivw(j),j)+fse(i+ive(j),j)) *(ust(i+ive(j),j-1)-ust(i,j))
              vne(i,j)=(fne(i,j-1)+fne(i,j+1))*(vst(i+ive(j),j+1)-vst(i,j))
              vse(i,j)=(fse(i,j-1)+fse(i,j+1))*(vst(i+ive(j),j-1)-vst(i,j))
           END DO
        END DO
        !equation of motion---------------------------------------
        !adu and adv are on v points
        DO  j=6,jm-5
           DO  i=1,im
              adu(i,j)=(uew(i+ivw(j),j)+uew(i+ive(j),j)+uns(i,j-1)+uns(i,j+1)  &
                   +une(i+ivw(j),j-1)+une(i,j)+USE(i,j)+USE(i+ivw(j),j+1))  &
                   *rdpdx(i,j)*fad(i+ivw(j),j)
              adv(i,j)=(vew(i+ivw(j),j)+vew(i+ive(j),j)+vns(i,j-1)+vns(i,j+1)  &
                   +vne(i+ivw(j),j-1)+vne(i,j)+vse(i,j)+vse(i+ivw(j),j+1))  &
                   *rdpdy(i,j)*fad(i+ivw(j),j)
           END DO
        END DO
        !upstream advection of velocity components----------------
        DO  jak=1,jaktwo
           ja=kvhas(jak)
           i=ivlas(jak)
           j=jvlas(jak)
           tta=em(ja)*ust(i,j)
           ttb=en    *vst(i,j)
           pp=-tta-ttb
           qp=tta-ttb
           isp=INT(SIGN(h1,pp))
           isq=INT(SIGN(h1,qp))
           ifp=(isp-1)/2
           ifq=(-isq-1)/2
           ipq=(isp-isq)/2
           pp=ABS(pp)
           qp=ABS(qp)
           f3=pp*qp
           f0=f3-pp-qp
           f1=pp-f3
           f2=qp-f3
           adu(i,j)=f0*ust(i,j)+f1*ust(i+ive(j)+ifp,j+isp)  &
                +f2*ust(i+ive(j)+ifq,j+isq)+f3*ust(i+ipq,j+isp+isq)
           adv(i,j)=f0*vst(i,j)+f1*vst(i+ive(j)+ifp,j+isp)  &
                +f2*vst(i+ive(j)+ifq,j+isq)+f3*vst(i+ipq,j+isp+isq)
        END DO
        IF( .NOT. iter2) THEN
           DO  j=3,jm-2
              DO  i=2,im-1
                 tst(i,j)=adt  (i,j)*(hm(i,j)*tlc)+tst(i,j)
                 q2m(i,j)=adq2m(i,j)*(hm(i,j)*tlc)+q2m(i,j)
                 q2l(i,j)=adq2l(i,j)*(hm(i,j)*tlc)+q2l(i,j)
              END DO
           END DO
           DO  j=3,jm-2
              DO  i=2,im-1
                 ust(i,j)=adu(i,j)*vm(i,j)+ust(i,j)
                 vst(i,j)=adv(i,j)*vm(i,j)+vst(i,j)
              END DO
           END DO
           iter2=.TRUE.
           CYCLE
        END IF
        EXIT
     END DO
     DO  j=3,jm-2
        DO  i=2,im-1
           t(i,j,l)=adt(i,j)*(2.0*hm(i,j))+t(i,j,l)
        END DO
     END DO
     DO  j=3,jm-2
        DO  i=2,im-1
           u(i,j,l)=adu(i,j)*(2.0*vm(i,j))+u(i,j,l)
           v(i,j,l)=adv(i,j)*(2.0*vm(i,j))+v(i,j,l)
        END DO
     END DO
     IF(l == 1)THEN
        DO  j=3,jm-2
           DO  i=2,im-1
              adq2hl(i,j,1)=adq2l(i,j)
           END DO
        END DO
     ELSE
        DO  j=3,jm-2
           DO  i=2,im-1
              adq2hl(i,j,l)=adq2l(i,j)
              q2(i,j,l-1)=adq2m(i,j)*hm(i,j)+q2(i,j,l-1)
           END DO
        END DO
     END IF
  END DO
  DO  l=2,lm
     DO  j=3,jm-2
        DO  i=2,im-1
           hm(i,j)=htm(i,j,l)*hbm2(i,j)
           q2(i,j,l-1)=adq2hl(i,j,l-1)*hm(i,j)+q2(i,j,l-1)
        END DO
     END DO
  END DO
END SUBROUTINE hzadv
