SUBROUTINE hzadv2
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  ! subprogram:    hzadv2      horizontal advection of vapor and cloud
  !   prgrmmr: janjic          org: w/np2      date: 96-07-19
  ! abstract:
  !     hzadv2 calculates the contribution of the horizontal advection
  !     to the tendencies of specific humidity and cloud water and
  !     then updates those variables.  an anti-filtering technique
  !     is used.
  ! program history log:
  !   96-07-19  janjic  - originator

  IMPLICIT NONE
  INCLUDE "parmeta.f90"

  REAL, PARAMETER :: epsq=2.e-12,ff1=0.52500,ff2=-0.64813,ff3=0.24520
  REAL, PARAMETER :: ff4=-0.12189
  INTEGER,PARAMETER :: lda=lm+9,la=13,im1=im-1,jam=6+2*(jm-10),imjm=im*jm-jm/2 
  INTEGER,PARAMETER :: jamd=(jam*2-10)*3,lp1=lm+1

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

  INTEGER,DIMENSION(im,jm) :: ifpa ,ifqa ,ifpf ,ifqf,jfpa ,jfqa ,jfpf ,jfqf 
  REAL, DIMENSION(im,jm) :: afp,afq,q1,dqst,w1,dwst,dare,dvol,emh,hbm3

  INTEGER :: i,j,ihl,ihh,l,htmijl,jfp,jfq
  REAL :: enh,sumpq,sumnq,sumpw,sumnw,dvolp,tta,ttb,pp,qp,sign,dqstij,dwstij
  REAL :: rfacq, rfacw,rfqij,rfwij,q1ij,w1ij,d2pqq,d2pqw,qstij,wstij,q00,qp0
  REAL :: q0q,w00,wp0,w0q

  enh=float(idtad)*dt/(08.*dy)
  DO  j=1,jm
     DO  i=1,im
        hbm3(i,j)=0.
     END DO
  END DO
  DO  j=4,jm-3
     ihl=2-ihw(j)
     ihh=im-2
     DO  i=ihl,ihh
        hbm3(i,j)=1.
     END DO
  END DO
  DO  j=1,jm
     DO  i=1,im
        emh (i,j)=float(idtad)*dt/(08.*dx(i,j))
        dare(i,j)=hbm2(i,j)*dx(i,j)*dy
     END DO
  END DO
  DO  l=1,lm
     DO  j=1,jm
        DO  i=1,im
           dvol(i,j)=dare(i,j)*pdsl(i,j)*deta(l)
           htmijl=htm(i,j,l)
           q  (i,j,l)=amax1(q  (i,j,l),epsq)*htmijl
           cwm(i,j,l)=amax1(cwm(i,j,l),epsq)*htmijl
           q1  (i,j)=q  (i,j,l)
           w1  (i,j)=cwm(i,j,l)
        END DO
     END DO
     sumpq=0.
     sumnq=0.
     sumpw=0.
     sumnw=0.
     DO  j=3,jm-2
        DO  i=2,im-1
           dvolp=dvol(i,j)*hbm3(i,j)
           tta=(u(i,j-1,l)+u(i+ihw(j),j,l)+u(i+ihe(j),j,l)+u(i,j+1,l))  &
                *hbm2(i,j)*emh(i,j)
           ttb=(v(i,j-1,l)+v(i+ihw(j),j,l)+v(i+ihe(j),j,l)+v(i,j+1,l))  &
                *hbm2(i,j)*enh
           pp=-tta-ttb
           qp= tta-ttb
           jfp=INT(SIGN(1.,pp))
           jfq=INT(SIGN(1.,qp))
           ifpa(i,j)=ihe(j)+i+( jfp-1  )/2
           ifqa(i,j)=ihe(j)+i+(-jfq-1  )/2
           jfpa(i,j)=       j+jfp
           jfqa(i,j)=       j+jfq
           ifpf(i,j)=ihe(j)+i+(-jfp-1  )/2
           ifqf(i,j)=ihe(j)+i+( jfq-1  )/2
           jfpf(i,j)=       j-jfp
           jfqf(i,j)=       j-jfq
           pp=ABS(pp)*htm(i,j,l)*htm(ifpa(i,j),jfpa(i,j),l)
           qp=ABS(qp)*htm(i,j,l)*htm(ifqa(i,j),jfqa(i,j),l)
           afp (i,j)=(((ff4*pp+ff3)*pp+ff2)*pp+ff1)*pp
           afq (i,j)=(((ff4*qp+ff3)*qp+ff2)*qp+ff1)*qp
           dqstij=(q  (ifpa(i,j),jfpa(i,j),l)-q  (i,j,l))*pp  &
                +(q  (ifqa(i,j),jfqa(i,j),l)-q  (i,j,l))*qp
           dwstij=(cwm(ifpa(i,j),jfpa(i,j),l)-cwm(i,j,l))*pp  &
                +(cwm(ifqa(i,j),jfqa(i,j),l)-cwm(i,j,l))*qp
           dqst(i,j)=dqstij
           dwst(i,j)=dwstij
           dqstij=dqstij*dvolp
           dwstij=dwstij*dvolp
           IF(dqstij > 0.)    THEN
              sumpq =sumpq+dqstij
           ELSE
              sumnq =sumnq+dqstij
           END IF
           IF(dwstij > 0.)    THEN
              sumpw =sumpw+dwstij
           ELSE
              sumnw =sumnw+dwstij
           END IF
        END DO
     END DO
     !first moment conserving factor
     IF(sumpq > 1.)    THEN
        rfacq=-sumnq/sumpq
     ELSE
        rfacq=1.
     END IF
     IF(sumpw > 1.)    THEN
        rfacw=-sumnw/sumpw
     ELSE
        rfacw=1.
     END IF
     IF(rfacq < 0.9.OR.rfacq > 1.1)rfacq=1.
     IF(rfacw < 0.9.OR.rfacw > 1.1)rfacw=1.
     !impose conservation on advection-------------------------
     IF(rfacq < 1.)THEN
        DO j=3,jm-2
           DO i=2,im-1
              dqstij=dqst(i,j)
              rfqij=hbm3(i,j)*(rfacq-1.)+1.
              IF(dqstij < 0.)   dqstij=dqstij/rfqij
              q1(i,j)=q  (i,j,l)+dqstij
           END DO
        END DO
     ELSE
        DO j=3,jm-2
           DO i=2,im-1
              dqstij=dqst(i,j)
              rfqij=hbm3(i,j)*(rfacq-1.)+1.
              IF(dqstij >= 0.)   dqstij=dqstij*rfqij
              q1(i,j)=q  (i,j,l)+dqstij
           END DO
        END DO
     END IF
     IF(rfacw < 1.)THEN
        DO j=3,jm-2
           DO i=2,im-1
              dwstij=dwst(i,j)
              rfwij=hbm3(i,j)*(rfacw-1.)+1.
              IF(dwstij < 0.)   dwstij=dwstij/rfwij
              w1(i,j)=cwm(i,j,l)+dwstij
           END DO
        END DO
     ELSE
        DO j=3,jm-2
           DO i=2,im-1
              dwstij=dwst(i,j)
              rfwij=hbm3(i,j)*(rfacw-1.)+1.
              IF(dwstij >= 0.)   dwstij=dwstij*rfwij
              w1(i,j)=cwm(i,j,l)+dwstij
           END DO
        END DO
     END IF
     !anti-filtering step--------------------------------------
     sumpq=0.
     sumnq=0.
     sumpw=0.
     sumnw=0.
     !anti-filtering limiters----------------------------------
     DO  j=3,jm-2
        DO  i=2,im-1
           dvolp=dvol(i,j)
           q1ij =q1(i,j)
           w1ij =w1(i,j)
           d2pqq=((q1(ifpa(i,j),jfpa(i,j))-q1ij                   )  &
                -(q1ij                   -q1(ifpf(i,j),jfpf(i,j)))  &
                *htm(ifpf(i,j),jfpf(i,j),l))*afp(i,j)  &
                +((q1(ifqa(i,j),jfqa(i,j))-q1ij                   )  &
                -(q1ij                   -q1(ifqf(i,j),jfqf(i,j)))  &
                *htm(ifqf(i,j),jfqf(i,j),l))*afq(i,j)
           d2pqw=((w1(ifpa(i,j),jfpa(i,j))-w1ij                   )  &
                -(w1ij                   -w1(ifpf(i,j),jfpf(i,j)))  &
                *htm(ifpf(i,j),jfpf(i,j),l))*afp(i,j)  &
                +((w1(ifqa(i,j),jfqa(i,j))-w1ij                   )  &
                -(w1ij                   -w1(ifqf(i,j),jfqf(i,j)))  &
                *htm(ifqf(i,j),jfqf(i,j),l))*afq(i,j)
           qstij=q1ij-d2pqq
           wstij=w1ij-d2pqw
           q00=q  (i,j,l)
           qp0=q  (ifpa(i,j),jfpa(i,j),l)
           q0q=q  (ifqa(i,j),jfqa(i,j),l)
           w00=cwm(i,j,l)
           wp0=cwm(ifpa(i,j),jfpa(i,j),l)
           w0q=cwm(ifqa(i,j),jfqa(i,j),l)
           qstij=amax1(qstij,amin1(q00,qp0,q0q))
           qstij=amin1(qstij,amax1(q00,qp0,q0q))
           wstij=amax1(wstij,amin1(w00,wp0,w0q))
           wstij=amin1(wstij,amax1(w00,wp0,w0q))
           dqstij=qstij-q1ij
           dwstij=wstij-w1ij
           dqst(i,j)=dqstij
           dwst(i,j)=dwstij
           dqstij=dqstij*dvolp
           dwstij=dwstij*dvolp
           IF(dqstij > 0.)    THEN
              sumpq =sumpq+dqstij
           ELSE
              sumnq =sumnq+dqstij
           END IF
           IF(dwstij > 0.)    THEN
              sumpw =sumpw+dwstij
           ELSE
              sumnw =sumnw+dwstij
           END IF
        END DO
     END DO
     !first moment conserving factor
     IF(sumpq > 1.)    THEN
        rfacq=-sumnq/sumpq
     ELSE
        rfacq=1.
     END IF
     IF(sumpw > 1.)    THEN
        rfacw=-sumnw/sumpw
     ELSE
        rfacw=1.
     END IF
     IF(rfacq < 0.9.OR.rfacq > 1.1)rfacq=1.
     IF(rfacw < 0.9.OR.rfacw > 1.1)rfacw=1.
     !impose conservation on anti-filtering--------------------
     IF(rfacq < 1.)THEN
        DO j=3,jm-2
           DO i=2,im-1
              dqstij=dqst(i,j)
              rfqij=hbm2(i,j)*(rfacq-1.)+1.
              IF(dqstij >= 0.)   dqstij=dqstij*rfqij
              q  (i,j,l)=q1(i,j)+dqstij
           END DO
        END DO
     ELSE
        DO j=3,jm-2
           DO i=2,im-1
              dqstij=dqst(i,j)
              rfqij=hbm2(i,j)*(rfacq-1.)+1.
              IF(dqstij < 0.)   dqstij=dqstij/rfqij
              q  (i,j,l)=q1(i,j)+dqstij
           END DO
        END DO
     END IF
     IF(rfacw < 1.)THEN
        DO j=3,jm-2
           DO i=2,im-1
              dwstij=dwst(i,j)
              rfwij=hbm2(i,j)*(rfacw-1.)+1.
              IF(dwstij >= 0.)   dwstij=dwstij*rfwij
              cwm(i,j,l)=w1(i,j)+dwstij
           END DO
        END DO
     ELSE
        DO j=3,jm-2
           DO i=2,im-1
              dwstij=dwst(i,j)
              rfwij=hbm2(i,j)*(rfacw-1.)+1.
              IF(dwstij < 0.)   dwstij=dwstij/rfwij
              cwm(i,j,l)=w1(i,j)+dwstij
           END DO
        END DO
     END IF
     DO j=1,jm
        DO i=1,im
           htmijl=htm(i,j,l)
           q  (i,j,l)=amax1(q  (i,j,l),epsq)*htmijl
           cwm(i,j,l)=amax1(cwm(i,j,l),epsq)*htmijl
        END DO
     END DO
  END DO
END SUBROUTINE hzadv2
