SUBROUTINE pfdht
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  ! subprogram:    pfdht       pressure gradient/coriolis calc
  !   prgrmmr: janjic          org: w/np2      date: 93-10-28       
  !     
  ! abstract:
  !     pfdht calculates the pressure gradient force, updates the
  !     velocity components due to the effect of the pressure gradient
  !     and corioils forces, computes the divergence including the
  !     modification preventing gravity wave grid separation, and
  !     calculates the horizontal part of the omega-alpha term
  !     (the part proportional to the advection of mass along
  !     eta surfaces).
  !     
  ! program history log:
  !   87-06-??  janjic     - originator
  !   95-03-25  black      - conversion from 1-d to 2-d in horizontal
  !   96-03-29  black      - added external edge
  !     
  ! usage: call pfdht 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
  !                  dynam
  !                  vrbls
  !                  indx
  !   
  IMPLICIT NONE

  !*** problem size

  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER :: imjm=im*jm-jm/2
  INTEGER, PARAMETER :: jam=6+2*(jm-10)                                 
  INTEGER, PARAMETER :: lp1=lm+1 

  !*** global variables

  INCLUDE "ctlblk.f90" ! first(in); ntsd(in); ntstm(in); sigma(in)
  INCLUDE "masks.f90"  ! htm(in); vbm2(in); vbm3(in); vtm(in); hbm2(in)
  INCLUDE "dynam.f90"  ! aeta(in); cpgfu(in); cpgfv(in); curv(in); deta(in)
  ! dfl(in); dx(in); dy(in); fcp(in); fdiv(in); pt(in); 
  ! r(in); wpdar(in); f(in)
  INCLUDE "vrbls.f90"  ! fis(in); pd(in); q(in); res(in); t(inout); 
  ! u(inout); v(inout)
  INCLUDE "contin.f90" ! div(out); omgalf(out); pdsl(out); rtop(out)
  INCLUDE "indx.f90"   ! ihe(in); ihw(in); ive(in); ivw(in)

  !*** local variables

  REAL, PARAMETER    :: d50=.5e0
  REAL, PARAMETER    :: h1=1.e0
  REAL, PARAMETER    :: h2=2.e0
  REAL, PARAMETER    :: d00=0.e0
  REAL, PARAMETER    :: d608=.608e0                   
  REAL, PARAMETER    :: cp=1004.6e0
  INTEGER            :: i, j, l
  REAL               :: fiupk
  REAL               :: dpnek
  REAL               :: dpsek
  REAL               :: dcnek
  REAL               :: dcsek
  REAL               :: dpfnek
  REAL               :: dpfsek
  REAL               :: f0k
  REAL               :: upk
  REAL               :: vpk
  REAL               :: utk
  REAL               :: vtk
  REAL               :: pvnek
  REAL               :: pvsek
  REAL               :: ttk
  REAL, DIMENSION(im,jm) :: filo
  REAL, DIMENSION(im,jm) :: rdpd
  REAL, DIMENSION(im,jm) :: adpdx
  REAL, DIMENSION(im,jm) :: rdpdx
  REAL, DIMENSION(im,jm) :: adpdy
  REAL, DIMENSION(im,jm) :: rdpdy
  REAL, DIMENSION(im,jm) :: adpdne
  REAL, DIMENSION(im,jm) :: adpdse
  REAL, DIMENSION(im,jm) :: pew
  REAL, DIMENSION(im,jm) :: pns
  REAL, DIMENSION(im,jm) :: pcew
  REAL, DIMENSION(im,jm) :: pcns
  REAL, DIMENSION(im,jm) :: dpfew
  REAL, DIMENSION(im,jm) :: dpfns
  REAL, DIMENSION(im,jm) :: fns
  REAL, DIMENSION(im,jm) :: tns
  REAL, DIMENSION(im,jm) :: vm
  REAL, DIMENSION(im,jm) :: hm
  REAL, DIMENSION(0:im+1,0:jm+1) :: dpde
  REAL, DIMENSION(0:im+1,0:jm+1) :: fim
  REAL, DIMENSION(0:im+1,0:jm+1) :: apel
  REAL, DIMENSION(0:im+1,0:jm+1) :: pcxc
  REAL, DIMENSION(0:im+1,0:jm+1) :: udy
  REAL, DIMENSION(0:im+1,0:jm+1) :: vdx
  REAL, DIMENSION(0:im+1,0:jm+1) :: tew
  REAL, DIMENSION(0:im+1,0:jm+1) :: few
  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) :: fne
  REAL, DIMENSION(0:im+1,0:jm+1) :: fse
  REAL, DIMENSION(0:im+1,0:jm+1) :: pne
  REAL, DIMENSION(0:im+1,0:jm+1) :: pse
  REAL, DIMENSION(0:im+1,0:jm+1) :: cne
  REAL, DIMENSION(0:im+1,0:jm+1) :: cse
  REAL, DIMENSION(0:im+1,0:jm+1) :: ppne
  REAL, DIMENSION(0:im+1,0:jm+1) :: ppse
  REAL, DIMENSION(0:im+1,0:jm+1) :: pcne
  REAL, DIMENSION(0:im+1,0:jm+1) :: pcse


  CALL zero2(dpde)
  CALL zero2(fim)
  CALL zero2(apel)
  CALL zero2(pcxc)
  CALL zero2(udy)
  CALL zero2(vdx)
  CALL zero2(tew)
  CALL zero2(few)
  CALL zero2(tne)
  CALL zero2(tse)
  CALL zero2(fne)
  CALL zero2(fse)
  CALL zero2(pne)
  CALL zero2(cne)
  CALL zero2(pse)
  CALL zero2(cse)
  CALL zero2(ppne)
  CALL zero2(pcne)
  CALL zero2(ppse)
  CALL zero2(pcse)

  !*** preparatory calculations

  IF(sigma)THEN
     DO j=1,jm
        DO i=1,im
           filo(i,j)=fis(i,j)
           pdsl(i,j)=pd(i,j)
        END DO
     END DO
  ELSE
     DO j=1,jm
        DO i=1,im
           filo(i,j)=d00
           pdsl(i,j)=res(i,j)*pd(i,j)
        END DO
     END DO
  ENDIF
  DO l=1,lm
     DO j=1,jm
        DO i=1,im
           div(i,j,l)=d00
           omgalf(i,j,l)=d00
        END DO
     END DO
  END DO
  DO j=1,jm
     DO i=1,im
        adpdx(i,j)=d00
        adpdy(i,j)=d00
     ENDDO
  ENDDO

  !*** main vertical integration loop

  DO l=lm,1,-1
     DO j=1,jm
        DO i=1,im
           dpde(i,j)=deta(l)*pdsl(i,j)
           rdpd(i,j)=h1/dpde(i,j)
        END DO
     END DO
     DO j=2,jm-1
        DO i=1,im-1
           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
     DO j=1,jm
        DO i=1,im
           apel(i,j)=pt+aeta(l)*pdsl(i,j)
           rtop(i,j,l)=r*t(i,j,l)*(h1+d608*q(i,j,l))/apel(i,j)
           fiupk=filo(i,j)+rtop(i,j,l)*dpde(i,j)
           fim(i,j)=filo(i,j)+fiupk
           filo(i,j)=dfl(l)+htm(i,j,l)*(fiupk-dfl(l))
        END DO
     END DO

     !*** diagonal contributions to pressure gradient force

     DO j=1,jm-1
        DO i=1,im-1
           adpdne(i,j)=dpde(i+ihe(j),j+1)+dpde(i,j)
           pne(i,j)=h2*(fim(i+ihe(j),j+1)-fim(i,j))
           ppne(i,j)=pne(i,j)*adpdne(i,j)
           cne(i,j)=h2*(rtop(i+ihe(j),j+1,l)+rtop(i,j,l)) &
                *(apel(i+ihe(j),j+1)-apel(i,j))
           pcne(i,j)=cne(i,j)*adpdne(i,j)
        END DO
     END DO
     DO j=2,jm
        DO i=1,im-1
           adpdse(i,j)=dpde(i+ihe(j),j-1)+dpde(i,j)
           pse(i,j)=h2*(fim(i+ihe(j),j-1)-fim(i,j))
           ppse(i,j)=pse(i,j)*adpdse(i,j)
           cse(i,j)=h2*(rtop(i+ihe(j),j-1,l)+rtop(i,j,l)) &
                *(apel(i+ihe(j),j-1)-apel(i,j))
           pcse(i,j)=cse(i,j)*adpdse(i,j)
        END DO
     END DO

     !*** continuity equation modification

     DO j=2,jm-1
        DO i=1,im
           pcxc(i,j)=vbm3(i,j)*vtm(i,j,l)*(pne(i+ivw(j),j) &
                +cne(i+ivw(j),j)+pse(i+ivw(j),j)+cse(i+ivw(j),j) &
                -pne(i,j-1)-cne(i,j-1)-pse(i,j+1)-cse(i,j+1))
        END DO
     END DO
     DO j=3,jm-2
        DO i=2,im-1
           div(i,j,l)=deta(l)*wpdar(i,j) &
                *(pcxc(i+ihe(j),j)-pcxc(i,j+1) &
                +pcxc(i+ihw(j),j)-pcxc(i,j-1))
        END DO
     END DO

     !*** lat & long pressure force components

     DO j=2,jm-1
        DO i=1,im
           dpnek=pne(i+ivw(j),j)+pne(i,j-1)
           dpsek=pse(i+ivw(j),j)+pse(i,j+1)
           pew(i,j)=dpnek+dpsek
           pns(i,j)=dpnek-dpsek
           dcnek=cne(i+ivw(j),j)+cne(i,j-1)
           dcsek=cse(i+ivw(j),j)+cse(i,j+1)
           pcew(i,j)=(dcnek+dcsek)*adpdx(i,j)
           pcns(i,j)=(dcnek-dcsek)*adpdy(i,j)
        END DO
     END DO

     !*** if not first pass, update u and v (coriolis & pgf)

     IF(.NOT. first) THEN
        DO j=3,jm-2
           DO i=1,im-1
              dpfnek=((ppne(i+ivw(j),j)+ppne(i,j-1)) &
                   +(pcne(i+ivw(j),j)+pcne(i,j-1)))*h2
              dpfsek=((ppse(i+ivw(j),j)+ppse(i,j+1)) &
                   +(pcse(i+ivw(j),j)+pcse(i,j+1)))*h2
              dpfew(i,j)=dpfnek+dpfsek
              dpfns(i,j)=dpfnek-dpfsek
           END DO
        END DO
        DO j=3,jm-2
           DO i=1,im-1
              f0k=u(i,j,l)*curv(i,j)+f(i,j)
              vm(i,j)=vtm(i,j,l)*vbm2(i,j)
              upk=((dpfew(i,j)+pcew(i,j))*rdpdx(i,j) &
                   +pew(i,j))*cpgfu(i,j)+f0k*v(i,j,l)+u(i,j,l)
              vpk=((dpfns(i,j)+pcns(i,j))*rdpdy(i,j) &
                   +pns(i,j))*cpgfv-f0k*u(i,j,l)+v(i,j,l)
              utk=u(i,j,l)
              vtk=v(i,j,l)
              u(i,j,l)=((f0k*vpk+upk)/(f0k*f0k+h1)-u(i,j,l)) &
                   *vm(i,j)+u(i,j,l)
              v(i,j,l)=(vpk-f0k*u(i,j,l)-v(i,j,l)) &
                   *vm(i,j)+v(i,j,l)
           END DO
        END DO
     END IF

     !*** if not last pass 

     IF(ntsd.NE.ntstm) THEN

        !*** lat & lon fluxes & omega-alpha components

        DO j=2,jm-1
           DO i=1,im
              udy(i,j)=dy*u(i,j,l)
              few(i,j)=udy(i,j)*adpdx(i,j)
              tew(i,j)=udy(i,j)*pcew(i,j)
              vdx(i,j)=dx(i,j)*v(i,j,l)
              fns(i,j)=vdx(i,j)*adpdy(i,j)
              tns(i,j)=vdx(i,j)*pcns(i,j)
           END DO
        END DO

        !*** diagonal fluxes and diagonally averaged wind

        DO j=2,jm-2
           DO i=1,im-1
              pvnek=(udy(i+ihe(j),j)+vdx(i+ihe(j),j))+(udy(i,j+1)+vdx(i,j+1))
              fne(i,j)=pvnek*adpdne(i,j)
              tne(i,j)=pvnek*pcne(i,j)*h2
           END DO
        END DO
        !
        DO j=3,jm-1
           DO i=1,im-1
              pvsek=(udy(i+ihe(j),j)-vdx(i+ihe(j),j))+(udy(i,j-1)-vdx(i,j-1))
              fse(i,j)=pvsek*adpdse(i,j)
              tse(i,j)=pvsek*pcse(i,j)*h2
           END DO
        END DO

        !*** horizontal part of omega-alpha & divergence

        DO j=3,jm-2
           DO i=2,im-1
              hm(i,j)=htm(i,j,l)*hbm2(i,j)
              ttk=t(i,j,l)
              omgalf(i,j,l)=(tew(i+ihe(j),j)+tew(i+ihw(j),j)+tns(i,j+1) &
                   +tns(i,j-1)+tne(i,j)+tne(i+ihw(j),j-1)+tse(i,j)&
                   +tse(i+ihw(j),j+1))*rdpd(i,j)*fcp(i,j)*hm(i,j)
              t(i,j,l)=omgalf(i,j,l)+t(i,j,l)
              div(i,j,l)=(((few(i+ihe(j),j)+fns(i,j+1)+fne(i,j)+fse(i,j)) &
                   -(few(i+ihw(j),j)+fns(i,j-1)+fne(i+ihw(j),j-1) &
                   +fse(i+ihw(j),j+1)))*fdiv(i,j)+div(i,j,l))*hm(i,j)
           END DO
        END DO
     END IF
  END DO
END SUBROUTINE pfdht
