!
!NCEP_MESO:MODEL_LAYER: INERTIAL GRAVITY WAVE ADJUSTMENT
!
!#define FROM_PROD
#ifdef FROM_PROD
!
! This is the version that comes from the standalone NMM PROD
! Code. It retains the original loop nesting and only the 
! array indices have been permuted to IKJ. Also K runs from
! bottom to top with ascending index values in vertical.
!
!-----------------------------------------------------------------------
!

      MODULE MODULE_IGWAVE_ADJUST 2
!
!-----------------------------------------------------------------------
      USE MODULE_MODEL_CONSTANTS
      USE MODULE_MPP
      USE MODULE_TIMERS  ! this one creates a name conflict at compile time
!-----------------------------------------------------------------------
!***
!***  SPECIFY THE NUMBER OF TIMES TO SMOOTH THE VERTICAL VELOCITY
!***  AND THE NUMBER OF ROWS FROM THE NORTHERN AND SOUTHERN EDGES
!***  OF THE GLOBAL DOMAIN BEYOND WHICH THE SMOOTHING DOES NOT GO
!***  FOR SUBROUTINE PDTEDT
!
      INTEGER :: KSMUD=0,LNSDT=7
!
!-----------------------------------------------------------------------
!
      CONTAINS
!
!***********************************************************************

      SUBROUTINE PFDHT(NTSD,LAST_TIME,PT,DETA1,DETA2,PDTOP,RES,FIS      & 1,1
     &                ,HYDRO,SIGMA,FIRST,DX,DY                          &
     &                ,HTM,HBM2,VTM,VBM2,VBM3                           &
     &                ,FDIV,FCP,WPDAR,DFL,CPGFU,CPGFV                   &
     &                ,PD,PDSL,T,Q,U,V,CWM,OMGALF,PINT,DWDT             &
     &                ,RTOP,DIV,FEW,FNS,FNE,FSE                         &
     &                ,IHE,IHW,IVE,IVW,INDX3_WRK                        &
     &                ,IDS,IDE,JDS,JDE,KDS,KDE                          &
     &                ,IMS,IME,JMS,JME,KMS,KME                          &
     &                ,ITS,ITE,JTS,JTE,KTS,KTE)
!***********************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    PFDHT       DIVERGENCE/HORIZONTAL OMEGA-ALPHA
!   PRGRMMR: JANJIC          ORG: W/NP22     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
!      COORDINATE 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
!   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
!   02-02-01  BLACK      - REWRITTEN FOR WRF CODING STANDARDS
!
! USAGE: CALL PFDHT FROM MAIN PROGRAM SOLVE_RUNSTREAM
!   INPUT ARGUMENT LIST:
!
!   OUTPUT ARGUMENT LIST:
!
!   OUTPUT FILES:
!     NONE
!
!   SUBPROGRAMS CALLED:
!
!     UNIQUE: NONE
!
!     LIBRARY: NONE
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$
!----------------------------------------------------------------------
!**********************************************************************
!----------------------------------------------------------------------
!----------------------------------------------------------------------
      LOGICAL,INTENT(IN) :: FIRST,HYDRO
      INTEGER,INTENT(IN) :: SIGMA
!
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
                           ,IMS,IME,JMS,JME,KMS,KME                    &
                           ,ITS,ITE,JTS,JTE,KTS,KTE
!
      INTEGER, DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
!
!***  NMM_MAX_DIM is set in configure.wrf and must agree with
!***  the value of dimspec q in the Registry/Registry
!
      INTEGER,DIMENSION(-3:3,NMM_MAX_DIM,0:6),INTENT(IN) :: indx3_wrk
!
      INTEGER,INTENT(IN) :: NTSD
      LOGICAL,INTENT(IN) :: LAST_TIME
!
      REAL,INTENT(IN) :: CPGFV,DY,PDTOP,PT
!
      REAL,DIMENSION(KMS:KME),INTENT(IN) :: DETA1,DETA2
!
      REAL,DIMENSION(KMS:KME),INTENT(IN) :: DFL
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CPGFU,DX,FCP,FDIV  &
                                                   ,PD,FIS,RES,WPDAR   &
                                                   ,HBM2,VBM2,VBM3
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: CWM,DWDT   &
                                                           ,Q,HTM,VTM
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: PINT
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: DIV     &
                                                              ,OMGALF  &
                                                              ,RTOP,T  &
                                                              ,U,V
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: FEW,FNS   &
                                                            ,FNE,FSE
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: PDSL


                             REAL  &
     & FILO  (IMS:IME,JMS:JME),fim   (IMS:IME,JMS:JME)  &
     &,RDPD  (IMS:IME,JMS:JME)  &
     &,ADPDX (IMS:IME,JMS:JME),RDPDX (IMS:IME,JMS:JME)  &
     &,ADPDY (IMS:IME,JMS:JME),RDPDY (IMS:IME,JMS:JME)  &
     &,ADPDNE(IMS:IME,JMS:JME),ADPDSE(IMS:IME,JMS:JME)  &
     &,UDY   (IMS:IME,JMS:JME),VDX   (IMS:IME,JMS:JME)  &
     &,PEW   (IMS:IME,JMS:JME),PNS   (IMS:IME,JMS:JME)  &
     &,PCEW  (IMS:IME,JMS:JME),PCNS  (IMS:IME,JMS:JME)  &
     &,DPFEW (IMS:IME,JMS:JME),DPFNS (IMS:IME,JMS:JME)  &
     &,TNS   (IMS:IME,JMS:JME)  &
     &,HM    (IMS:IME,JMS:JME)
!
                             REAL &
     & DPDE  (IMS:IME,JMS:JME) &
     &,APEL  (IMS:IME,JMS:JME),PCXC  (IMS:IME,JMS:JME) &
     &,DFDZ  (IMS:IME,JMS:JME) &
     &,TEW   (IMS:IME,JMS:JME) &
     &,TNE   (IMS:IME,JMS:JME),TSE   (IMS:IME,JMS:JME) &
     &,PNE   (IMS:IME,JMS:JME),PSE   (IMS:IME,JMS:JME) &
     &,CNE   (IMS:IME,JMS:JME),CSE   (IMS:IME,JMS:JME) &
     &,PPNE  (IMS:IME,JMS:JME),PPSE  (IMS:IME,JMS:JME) &
     &,PCNE  (IMS:IME,JMS:JME),PCSE  (IMS:IME,JMS:JME)

      REAL :: DCNEK,DCSEK          &
             ,DPFNEK,DPFSEK,DPNEK,DPSEK                    &
             ,PVNEK,PVSEK,VM
      REAL :: FIUPP, DFI, RTOPP, APELP
      INTEGER :: L, I, J, K
!-----------------------------------------------------------------------
!--------------PREPARATORY CALCULATIONS---------------------------------
!-----------------------------------------------------------------------

#ifdef DEC_ALPHA
FILO    =   0.  ;   FIM     =   0.  ;   RDPD    =   0.  ;   
ADPDX   =   0.  ;   RDPDX   =   0.  ;   ADPDY   =   0.  ;   
RDPDY   =   0.  ;   ADPDNE  =   0.  ;   ADPDSE  =   0.  ;   
UDY     =   0.  ;   VDX     =   0.  ;   PEW     =   0.  ;   
PNS     =   0.  ;   PCEW    =   0.  ;   PCNS    =   0.  ;   
DPFEW   =   0.  ;   DPFNS   =   0.  ;   TNS     =   0.  ;   
HM      =   0.  ;   DPDE    =   0.  ;   APEL    =   0.  ;   
PCXC    =   0.  ;   DFDZ    =   0.  ;   TEW     =   0.  ;   
TNE     =   0.  ;   TSE     =   0.  ;   PNE     =   0.  ;   
PSE     =   0.  ;   CNE     =   0.  ;   CSE     =   0.  ;   
PPNE    =   0.  ;   PPSE    =   0.  ;   PCNE    =   0.  ;   
PCSE    =   0.  ;   
#endif


      IF(SIGMA.EQ.1)THEN
  !IBM* ASSERT(ITERCNT(100))
        DO 50 J=MYJS_P5,MYJE_P5
  !IBM* ASSERT(ITERCNT(100))
        DO 50 I=MYIS_P5,MYIE_P5
        FILO(I,J)=FIS(I,J)
        PDSL(I,J)=PD(I,J)
   50   CONTINUE
      ELSE
  !IBM* ASSERT(ITERCNT(100))
        DO 100 J=MYJS_P5,MYJE_P5
  !IBM* ASSERT(ITERCNT(100))
        DO 100 I=MYIS_P5,MYIE_P5
        FILO(I,J)=0.0
        PDSL(I,J)=RES(I,J)*PD(I,J)
  100   CONTINUE
      ENDIF
!
  !IBM* ASSERT(ITERCNT(60))
      DO L=KTS,KTE
  !IBM* ASSERT(ITERCNT(100))
        DO J=MYJS_P4,MYJE_P4
  !IBM* ASSERT(ITERCNT(100))
        DO I=MYIS_P4,MYIE_P4
          OMGALF(I,L,J)=0.
          DIV(I,L,J)=0.
        ENDDO
        ENDDO
      ENDDO
!-------------------- MAIN VERTICAL INTEGRATION LOOP -------------------
  !IBM* ASSERT(ITERCNT(60))
      DO 400 L = KTS,KTE
!-----------------------------------------------------------------------
  !IBM* ASSERT(ITERCNT(100))
      DO 205 J=MYJS_P4,MYJE_P4
  !IBM* ASSERT(ITERCNT(100))
      DO 205 I=MYIS_P4,MYIE_P4
!
      HM(I,J)=HTM(I,L,J)*HBM2(I,J)
!
      apelp=(pint(i,l,j)+pint(i,l+1,j))*0.5
      rtopp=(q(i,l,j)*0.608-cwm(i,l,j)+1.)*t(i,l,j)*r_d/apelp
      dfi=rtopp*(deta1(l)*pdtop+deta2(l)*pdsl(i,j))
!
      apel(i,j)=apelp
      rtop(i,l,j)=rtopp
      dfdz(i,j)=rtopp
!
      fiupp=filo(i,j)+dfi
!
      fim(i,j)=filo(i,j)+fiupp
!
      filo(i,j)=(fiupp-dfl(l+1))*htm(i,l,j)+dfl(l+1)
  205 CONTINUE
!
  !IBM* ASSERT(ITERCNT(100))
      DO 210 J=MYJS_P4,MYJE_P4
  !IBM* ASSERT(ITERCNT(100))
      DO 210 I=MYIS_P4,MYIE_P4
      dpde(i,j)=deta1(l)*pdtop+deta2(l)*pdsl(i,j)
  210 CONTINUE
!
  !IBM* ASSERT(ITERCNT(100))
      DO 215 J=MYJS,MYJE
  !IBM* ASSERT(ITERCNT(100))
      DO 215 I=MYIS,MYIE
      RDPD(I,J)=1./DPDE(I,J)
  215 CONTINUE
!
  !IBM* ASSERT(ITERCNT(100))
      DO 220 J=MYJS1_P3,MYJE1_P3
  !IBM* ASSERT(ITERCNT(100))
      DO 220 I=MYIS_P3,MYIE_P3
      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)=1./ADPDX(I,J)
      RDPDY(I,J)=1./ADPDY(I,J)
  220 CONTINUE
!
!--------------DIAGONAL CONTRIBUTIONS TO PRESSURE GRADIENT FORCE--------
!
  !IBM* ASSERT(ITERCNT(100))
      DO 240 J=MYJS_P3,MYJE1_P3
  !IBM* ASSERT(ITERCNT(100))
      DO 240 I=MYIS_P3,MYIE_P3
      ADPDNE(I,J)=DPDE(I+IHE(J),J+1)+DPDE(I,J)
      PNE(I,J)=(FIM (I+IHE(J),J+1)-FIM (I,J))  &
     &        *(DWDT(I+IHE(J),L,J+1)+DWDT(I,L,J))
      PPNE(I,J)=PNE(I,J)*ADPDNE(I,J)
      CNE(I,J)=(DFDZ(I+IHE(J),J+1)+DFDZ(I,J))*2.  &
     &        *(APEL(I+IHE(J),J+1)-APEL(I,J))
      PCNE(I,J)=CNE(I,J)*ADPDNE(I,J)
  240 CONTINUE
!
  !IBM* ASSERT(ITERCNT(100))
      DO 250 J=MYJS1_P3,MYJE_P3
  !IBM* ASSERT(ITERCNT(100))
      DO 250 I=MYIS_P3,MYIE_P3
      ADPDSE(I,J)=DPDE(I+IHE(J),J-1)+DPDE(I,J)
      PSE(I,J)=(FIM (I+IHE(J),J-1)-FIM (I,J))  &
     &        *(DWDT(I+IHE(J),L,J-1)+DWDT(I,L,J))
      PPSE(I,J)=PSE(I,J)*ADPDSE(I,J)
!write(0,*)'i,j,ihe(j)',i,j,ihe(j),i+ihe(j)
      CSE(I,J)=(DFDZ(I+IHE(J),J-1)+DFDZ(I,J))*2.  &
     &        *(APEL(I+IHE(J),J-1)-APEL(I,J))
      PCSE(I,J)=CSE(I,J)*ADPDSE(I,J)
  250 CONTINUE
!
!--------------CONTINUITY EQUATION MODIFICATION-------------------------
!
  !IBM* ASSERT(ITERCNT(100))
      DO 260 J=MYJS1_P1,MYJE1_P1
  !IBM* ASSERT(ITERCNT(100))
      DO 260 I=MYIS_P1,MYIE_P1
      PCXC(I,J)=VBM3(I,J)*VTM(I,L,J)*(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))
  260 CONTINUE
!-----------------------------------------------------------------------
  !IBM* ASSERT(ITERCNT(100))
      DO 270 J=MYJS2,MYJE2
  !IBM* ASSERT(ITERCNT(100))
      DO 270 I=MYIS1,MYIE1
      div(i,l,j)=(deta1(l)*pdtop/(101300.-pt)  &
     &           +deta2(l)*(1.-pdtop/(101300.-pt)))  &
     &         *wpdar(i,j)  &
     &         *(PCXC(I+IHE(J),J)-PCXC(I,J+1)  &
     &          +PCXC(I+IHW(J),J)-PCXC(I,J-1))
  270 CONTINUE
!
!--------------LAT & LONG PRESSURE FORCE COMPONENTS---------------------
!
  !IBM* ASSERT(ITERCNT(100))
      DO 280 J=MYJS1_P2,MYJE1_P2
  !IBM* ASSERT(ITERCNT(100))
      DO 280 I=MYIS_P2,MYIE_P3
      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)
  280 CONTINUE
!
!--------------FIRST PASS SWITCH----------------------------------------
      IF(.NOT.FIRST)THEN
!
!--------------UPDATE U AND V PGF---------------------------------------
!
  !IBM* ASSERT(ITERCNT(100))
      DO 290 J=MYJS2_P2,MYJE2_P2
  !IBM* ASSERT(ITERCNT(100))
      DO 290 I=MYIS_P2,MYIE1_P2
      DPFNEK=((PPNE(I+IVW(J),J)+PPNE(I,J-1)) &
     &       +(PCNE(I+IVW(J),J)+PCNE(I,J-1)))
      dpfnek=dpfnek+dpfnek
      DPFSEK=((PPSE(I+IVW(J),J)+PPSE(I,J+1)) &
     &       +(PCSE(I+IVW(J),J)+PCSE(I,J+1)))
      dpfsek=dpfsek+dpfsek
      DPFEW(I,J)=DPFNEK+DPFSEK
      DPFNS(I,J)=DPFNEK-DPFSEK
  290 CONTINUE
!-----------------------------------------------------------------------
  !IBM* ASSERT(ITERCNT(100))
      DO 300 J=MYJS2_P3,MYJE2_P3
  !IBM* ASSERT(ITERCNT(100))
      DO 300 I=MYIS_P2,MYIE1_P2
        VM=VTM(I,L,J)*VBM2(I,J)
        U(I,L,J)=(((DPFEW(I,J)+PCEW(I,J))*RDPDX(I,J) &
     &             +PEW(I,J))*CPGFU(I,J))*VM+U(I,L,J)
        V(I,L,J)=(((DPFNS(I,J)+PCNS(I,J))*RDPDY(I,J) &
     &             +PNS(I,J))*CPGFV)*VM+V(I,L,J)
  300 CONTINUE
!
!--------------LAST PASS SWITCH-----------------------------------------
!      IF(NTSD.EQ.NTSTM)GO TO 400
      ENDIF
      IF(.NOT. LAST_TIME) THEN
!-----------------------------------------------------------------------
 305  CONTINUE
!
!--------------LAT & LON FLUXES & OMEGA-ALPHA COMPONENTS----------------
!
  !IBM* ASSERT(ITERCNT(100))
      DO 310 J=MYJS1_P2,MYJE1_P2
  !IBM* ASSERT(ITERCNT(100))
      DO 310 I=MYIS_P2,MYIE_P3
      UDY(I,J)=DY*U(I,L,J)
      few(i,l,j)=udy(i,j)*adpdx(i,j)
      TEW(I,J)=UDY(I,J)*PCEW(I,J)
      VDX(I,J)=DX(I,J)*V(I,L,J)
      fns(i,l,j)=vdx(i,j)*adpdy(i,j)
      TNS(I,J)=VDX(I,J)*PCNS(I,J)
  310 CONTINUE
!
!--------------DIAGONAL FLUXES AND DIAGONALLY AVERAGED WIND-------------
!
  !IBM* ASSERT(ITERCNT(100))
      DO 320 J=MYJS1_P1,MYJE2_P1
  !IBM* ASSERT(ITERCNT(100))
      DO 320 I=MYIS_P1,MYIE1_P1
      PVNEK=(UDY(I+IHE(J),J)+VDX(I+IHE(J),J))+(UDY(I,J+1)+VDX(I,J+1))
      fne(i,l,j)=pvnek*adpdne(i,j)
      TNE(I,J)=PVNEK*PCNE(I,J)*2.
  320 CONTINUE
!
  !IBM* ASSERT(ITERCNT(100))
      DO 330 J=MYJS2_P1,MYJE1_P1
  !IBM* ASSERT(ITERCNT(100))
      DO 330 I=MYIS_P1,MYIE1_P1
      PVSEK=(UDY(I+IHE(J),J)-VDX(I+IHE(J),J))+(UDY(I,J-1)-VDX(I,J-1))
      fse(i,l,j)=pvsek*adpdse(i,j)
      TSE(I,J)=PVSEK*PCSE(I,J)*2.
  330 CONTINUE
!
!--------------HORIZONTAL PART OF OMEGA-ALPHA & DIVERGENCE--------------
!
  !IBM* ASSERT(ITERCNT(100))
      DO 340 J=MYJS2,MYJE2
  !IBM* ASSERT(ITERCNT(100))
      DO 340 I=MYIS1,MYIE1
      OMGALF(I,L,J)=(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,L,J)=OMGALF(I,L,J)+T(I,L,J)
      div(i,l,j)=(((few(i+ihe(j),l,j)+fns(i,l,j+1)  &
     &             +fne(i,l,j)+fse(i,l,j))  &
     &            -(few(i+ihw(j),l,j)+fns(i,l,j-1)  &
     &             +fne(i+ihw(j),l,j-1)+fse(i+ihw(j),l,j+1)))*fdiv(i,j)  &
     &           +div(i,l,j))*hm(i,j)

  340 CONTINUE
     ENDIF
!-----------------------------------------------------------------------
  400                        CONTINUE
!------------------------------------------------------------------------
                             RETURN
!
      END SUBROUTINE PFDHT
!
!----------------------------------------------------------------------
!**********************************************************************

      SUBROUTINE PDTEDT(                                               & 1,4
#ifdef DM_PARALLEL
                        grid,                                          &
#endif
                        NTSD,DT,PT,ETA2,RES,HYDRO                      &
                       ,HTM,HBM2,EF4T                                  &
                       ,PD,PDSL,PDSLO,T,DWDT,RTOP,OMGALF               &
                       ,PETDT,PINT,DIV,PSDT                            &
                       ,IHE,IHW,IVE,IVW,INDX3_WRK                      &
                       ,IDS,IDE,JDS,JDE,KDS,KDE                        &
                       ,IMS,IME,JMS,JME,KMS,KME                        &
                       ,ITS,ITE,JTS,JTE,KTS,KTE)
!**********************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    PDTEDT      SURFACE PRESSURE TENDENCY CALC
!   PRGRMMR: JANJIC          ORG: W/NP2      DATE: 96-07-??
!
! ABSTRACT:
!     PDTEDT VERTICALLY INTEGRATES THE MASS FLUX DIVERGENCE TO
!     OBTAIN THE SURFACE PRESSURE TENDENCY AND VERTICAL VELOCITY ON
!     THE LAYER INTERFACES.  THEN IT UPDATES THE HYDROSTATIC SURFACE
!     PRESSURE, THE NONHYDROSTATIC PRESSURE, AND ADDS THE LOCAL TIME
!     DERIVATIVE AND VERTICAL ADVECTION OF NONHYDROSTATIC PRESSURE
!     CONTRIBUTION TO THE OMEGA-ALPHA TERM OF THE THERMODYNAMIC
!     EQUATION.  ALSO, THE OMEGA-ALPHA TERM IS COMPUTED FOR DIAGNOSTICS.
!
! PROGRAM HISTORY LOG:
!   87-06-??  JANJIC     - ORIGINATOR
!   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
!   96-05-??  JANJIC     - ADDED NONHYDROSTATIC EFFECTS & MERGED THE
!                          PREVIOUS SUBROUTINES PDTE & PDNEW
!   00-01-03  BLACK      - DISTRIBUTED MEMORY AND THREADS
!   01-02-23  BLACK      - CONVERTED TO WRF FORMAT
!   01-04-11  BLACK      - REWRITTEN FOR WRF CODING STANDARDS
!
! USAGE: CALL PDTEDT FROM SUBROUTINE SOLVE_RUNSTREAM
!   INPUT ARGUMENT LIST:
!
!   OUTPUT ARGUMENT LIST:
!
!   OUTPUT FILES:
!     NONE
!
!   SUBPROGRAMS CALLED:
!
!     UNIQUE: NONE
!
!     LIBRARY: NONE
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$
!***********************************************************************
#ifdef DM_PARALLEL
      USE module_domain
      USE module_dm
#endif
!-----------------------------------------------------------------------
      IMPLICIT NONE
!-----------------------------------------------------------------------
#ifdef DM_PARALLEL
      INCLUDE "mpif.h"
      TYPE (DOMAIN) :: grid
#endif
!-----------------------------------------------------------------------
      LOGICAL,INTENT(IN) :: HYDRO
!
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                     &
                           ,IMS,IME,JMS,JME,KMS,KME                     &
                           ,ITS,ITE,JTS,JTE,KTS,KTE
!
      INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
!
!***  NMM_MAX_DIM is set in configure.wrf and must agree with
!***  the value of dimspec q in the Registry/Registry
!
      INTEGER,DIMENSION(-3:3,NMM_MAX_DIM,0:6),INTENT(IN) :: INDX3_WRK
!
      INTEGER,INTENT(IN) :: NTSD
!
      REAL,INTENT(IN) :: DT,EF4T,PT
!
      REAL,DIMENSION(KMS:KME),INTENT(IN) :: ETA2
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: RES,HBM2
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DWDT,HTM
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: DIV      &
                                                              ,OMGALF   &
                                                              ,RTOP,T
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: PD,PDSL
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: PINT
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: PETDT
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: PDSLO,PSDT
!
!-----------------------------------------------------------------------
!
!***  LOCAL VARIABLES
!
      INTEGER :: I,IHH,IHL,IX,J,JHH,JHL,JJ,JX,K,KNT,KS,NSMUD
      INTEGER :: J1_00,J1_M1,J2_00,J2_P1
      INTEGER :: MY_IS_GLB,MY_IE_GLB,MY_JS_GLB,MY_JE_GLB
!
      REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: APDT,HBMS,PPDT,PRET,TPM
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,-1:0) :: PNE
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE, 0:1) :: PSE
!
      REAL :: DWDTP,PETDTL,RHS,TPMP
!
!jm      REAL(8) :: timef
!-----------------------------------------------------------------------
!**********************************************************************
!----------------------------------------------------------------------
!
!jw   CALL ZERO(PDSLO,IMS,IME,JMS,JME)
      DO J=JMS,JME
      DO I=IMS,IME
        PDSLO(I,J)=0.
      ENDDO
      ENDDO
!
      MY_IS_GLB=ITS
      MY_IE_GLB=ITE
      MY_JS_GLB=JTS
      MY_JE_GLB=JTE
      MYIS=ITS
      MYIE=ITE
      MYJS=JTS
      MYJE=JTE
!
      MYIS1  =MAX(IDS+1,ITS)
      MYIE1  =MIN(IDE-1,ITE)
      MYJS2  =MAX(JDS+2,JTS)
      MYJE2  =MIN(JDE-2,JTE)
!
      MYIS_P1=MAX(IDS,ITS-1)
      MYIE_P1=MIN(IDE,ITE+1)
      MYIS_P2=MAX(IDS,ITS-2)
      MYIE_P2=MIN(IDE,ITE+2)
      MYJS_P2=MAX(JDS,JTS-2)
      MYJE_P2=MIN(JDE,JTE+2)
!----------------------------------------------------------------------
!***  COMPUTATION OF PRESSURE TENDENCY & PREPARATIONS
!----------------------------------------------------------------------
!
      DO J=MYJS_P2,MYJE_P2
        DO K=KTE-1,KTS,-1
        DO I=MYIS_P2,MYIE_P2
          DIV(I,K,J)=DIV(I,K+1,J)+DIV(I,K,J)
        ENDDO
        ENDDO
      ENDDO
!----------------------------------------------------------------------
      DO J=MYJS_P2,MYJE_P2
      DO I=MYIS_P2,MYIE_P2
        PSDT(I,J)=-DIV(I,KTS,J)
        APDT(I,J)=PSDT(I,J)
        PPDT(I,J)=PSDT(I,J)
        PDSLO(I,J)=PDSL(I,J)
      ENDDO
      ENDDO
!----------------------------------------------------------------------
      DO J=MYJS_P2,MYJE_P2
      DO I=MYIS_P2,MYIE_P2
        PD(I,J)=PSDT(I,J)*DT+PD(I,J)
        PRET(I,J)=PSDT(I,J)*RES(I,J)
        PDSL(I,J)=PD(I,J)*RES(I,J)
!
        PINT(I,KTE+1,J)=PT
!
        TPM(I,J)=PT+PINT(I,KTE,J)
      ENDDO
      ENDDO
!----------------------------------------------------------------------
!***  COMPUTATION OF PETDT
!----------------------------------------------------------------------
      DO J=MYJS_P2,MYJE_P2
        DO K=KTE-1,KTS,-1
        DO I=MYIS_P2,MYIE_P2
          PETDT(I,K,J)=-(PRET(I,J)*ETA2(K+1)+DIV(I,K+1,J))             &
                        *HTM(I,K,J)*HBM2(I,J)
        ENDDO
        ENDDO
      ENDDO
!----------------------------------------------------------------------
!***  KINETIC ENERGY GENERATION TERMS IN T EQUATION
!----------------------------------------------------------------------
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        DWDTP=DWDT(I,KTE,J)
        TPMP=PINT(I,KTE,J)+PINT(I,KTE-1,J)
!
        RHS=-DIV(I,KTE,J)*RTOP(I,KTE,J)*HTM(I,KTE,J)*DWDTP*EF4T
        OMGALF(I,KTE,J)=OMGALF(I,KTE,J)+RHS
        T(I,KTE,J)=RHS*HBM2(I,J)+T(I,KTE,J)
        PINT(I,KTE,J)=PRET(I,J)*(ETA2(KTE+1)+ETA2(KTE))*DWDTP*DT       &
                   +TPM(I,J)-PINT(I,KTE+1,J)
!
        TPM(I,J)=TPMP
      ENDDO
      ENDDO
!----------------------------------------------------------------------
      DO J=MYJS,MYJE
        DO K=KTE-1,KTS+1,-1
        DO I=MYIS,MYIE
          DWDTP=DWDT(I,K,J)
          TPMP=PINT(I,K,J)+PINT(I,K-1,J)
!
          RHS=-(DIV(I,K+1,J)+DIV(I,K,J))*RTOP(I,K,J)*HTM(I,K,J)*DWDTP  &
               *EF4T
          OMGALF(I,K,J)=OMGALF(I,K,J)+RHS
          T(I,K,J)=RHS*HBM2(I,J)+T(I,K,J)
          PINT(I,K,J)=PRET(I,J)*(ETA2(K+1)+ETA2(K))*DWDTP*DT           &
                     +TPM(I,J)-PINT(I,K+1,J)
!
          TPM(I,J)=TPMP
        ENDDO
        ENDDO
      ENDDO
!----------------------------------------------------------------------
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
!
        DWDTP=DWDT(I,KTS,J)
!
        RHS=-(DIV(I,KTS+1,J)+DIV(I,KTS,J))*RTOP(I,KTS,J)*HTM(I,KTS,J)  &
             *DWDTP*EF4T
        OMGALF(I,KTS,J)=OMGALF(I,KTS,J)+RHS
        T(I,KTS,J)=RHS*HBM2(I,J)+T(I,KTS,J)
        PINT(I,KTS,J)=PRET(I,J)*(ETA2(KTS+1)+ETA2(KTS))*DWDTP*DT       &
                       +TPM(I,J)-PINT(I,KTS+1,J)
      ENDDO
      ENDDO
!----------------------------------------------------------------------
!***  SMOOTHING VERTICAL VELOCITY ALONG BOUNDARIES
!----------------------------------------------------------------------
      nonhydrostatic_smoothing: IF(.NOT.HYDRO.AND.KSMUD.GT.0)THEN
!
        NSMUD=KSMUD
!
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          HBMS(I,J)=HBM2(I,J)
        ENDDO
        ENDDO
!
        JHL=LNSDT
        JHH=JDE-JHL+1
!
        DO J=JHL,JHH
          IF(J.GE.MY_JS_GLB.AND.J.LE.MY_JE_GLB)THEN
            IHL=JHL/2+1
            IHH=IDE-IHL+MOD(J,2)
!
            DO I=IHL,IHH
              IF(I.GE.MY_IS_GLB.AND.I.LE.MY_IE_GLB)THEN
                IX=I    ! -MY_IS_GLB+1
                JX=J    ! -MY_JS_GLB+1
                HBMS(IX,JX)=0.
              ENDIF
            ENDDO
!
          ENDIF
        ENDDO
!
!----------------------------------------------------------------------
!***
!***  SMOOTH THE VERTICAL VELOCITY
!***
!----------------------------------------------------------------------
!
        DO KS=1,NSMUD
!
!----------------------------------------------------------------------
!
!***  FILL SOUTHERNMOST SLABS OF THE PNE AND PSE WORKING ARRAYS
!
          JJ=MYJS2-1
          DO K=KTS,KTE-1
          DO I=MYIS_P1,MYIE1_P1
            PNE(I,K,-1)=(PETDT(I+IHE(JJ),K,JJ+1)-PETDT(I,K,JJ))        &
                        *HTM(I,K,JJ)*HTM(I+IHE(JJ),K,JJ+1)
          ENDDO
          ENDDO
!
          DO K=KTS,KTE-1
          DO I=MYIS_P1,MYIE1_P1
            PSE(I,K,0)=(PETDT(I+IHE(JJ+1),K,JJ)-PETDT(I,K,JJ+1))       &
                       *HTM(I+IHE(JJ+1),K,JJ)*HTM(I,K,JJ+1)
          ENDDO
          ENDDO
!
          KNT=0
!
!----------------------------------------------------------------------
!
!***  PROCEED NORTHWARD WITH THE SMOOTHING.
!***  PNE AT H(I,J) LIES BETWEEN (I,J) AND THE H POINT TO THE NE.
!***  PSE AT H(I,J) LIES BETWEEN (I,J) AND THE H POINT TO THE SE.
!
          DO J=MYJS2,MYJE2
!
            KNT=KNT+1
            J1_00=-MOD(KNT+1,2)
            J1_M1=-MOD(KNT,2)
            J2_P1=MOD(KNT,2)
            J2_00=MOD(KNT+1,2)
!
            DO K=KTS,KTE-1
            DO I=MYIS_P1,MYIE1_P1
              PNE(I,K,J1_00)=(PETDT(I+IHE(J),K,J+1)-PETDT(I,K,J))      &
                             *HTM(I,K+1,J)*HTM(I+IHE(J),K+1,J+1)
            ENDDO
            ENDDO
!
            DO K=KTS,KTE-1
            DO I=MYIS_P1,MYIE1_P1
              PSE(I,K,J2_P1)=(PETDT(I+IHE(J+1),K,J)-PETDT(I,K,J+1))    &
                             *HTM(I+IHE(J+1),K+1,J)*HTM(I,K+1,J+1)
            ENDDO
            ENDDO
!
            DO K=KTS,KTE-1
            DO I=MYIS1,MYIE1
              PETDTL=(PNE(I,K,J1_00)-PNE(I+IHW(J),K,J1_M1)             &
                     +PSE(I,K,J2_00)-PSE(I+IHW(J),K,J2_P1))*HBM2(I,J)
              PETDT(I,K,J)=PETDTL*HBMS(I,J)*0.125+PETDT(I,K,J)
            ENDDO
            ENDDO
!
          ENDDO
!
#ifdef DM_PARALLEL
!jm          btim=timef()
!jw need to add the following arrays to HALO:
!jw     3D - PETDT
!jw stencil needs to go at least 2 point in the x-dir and 2 in y
!jw       CALL EXCH(PETDT,KTE-1,2,2,IMS,IME,JMS,JME,KMS,KME)
!jw#         include "HALO_NMM_E.inc" not now -- need to move this
!jw                                   to the mediation layer
!jm          exch_tim=exch_tim+timef()-btim

# include <HALO_NMM_E.inc>


#endif


!----------------------------------------------------------------------
!
        ENDDO  ! End of smoothing loop
!
!----------------------------------------------------------------------
      ENDIF nonhydrostatic_smoothing
!----------------------------------------------------------------------
      END SUBROUTINE PDTEDT
!----------------------------------------------------------------------
!**********************************************************************

      SUBROUTINE DDAMP(NTSD,DT,DETA1,DETA2,PDSL,PDTOP,DIV,HBM2,VTM     & 1,3
                      ,T,U,V,DDMPU,DDMPV                               &
                      ,IHE,IHW,IVE,IVW,INDX3_WRK                       &
                      ,IDS,IDE,JDS,JDE,KDS,KDE                         &
                      ,IMS,IME,JMS,JME,KMS,KME                         &
                      ,ITS,ITE,JTS,JTE,KTS,KTE)
!**********************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    DDAMP       DIVERGENCE DAMPING
!   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 94-03-08
!
! ABSTRACT:
!     DDAMP MODIFIES THE WIND COMPONENTS SO AS TO REDUCE THE
!     HORIZONTAL DIVERGENCE.
!
! PROGRAM HISTORY LOG:
!   87-08-??  JANJIC     - ORIGINATOR
!   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
!   95-03-28  BLACK      - ADDED EXTERNAL EDGE
!   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
!   01-03-12  BLACK      - CONVERTED TO WRF STRUCTURE
!
! USAGE: CALL DDAMP FROM SUBROUTINE SOLVE_RUNSTREAM
!
!   INPUT ARGUMENT LIST:
!
!   OUTPUT ARGUMENT LIST:
!
!   OUTPUT FILES:
!     NONE
!
!   SUBPROGRAMS CALLED:
!
!     UNIQUE: NONE
!
!     LIBRARY: NONE
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$
!**********************************************************************
!----------------------------------------------------------------------
      IMPLICIT NONE
!----------------------------------------------------------------------
!
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
                           ,IMS,IME,JMS,JME,KMS,KME                    &
                           ,ITS,ITE,JTS,JTE,KTS,KTE
!
      INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
!
!***  NMM_MAX_DIM is set in configure.wrf and must agree with
!***  the value of dimspec q in the Registry/Registry
!
      INTEGER,DIMENSION(-3:3,NMM_MAX_DIM,0:6),INTENT(IN) :: INDX3_WRK
!
      INTEGER,INTENT(IN) :: NTSD
!
      REAL,INTENT(IN) :: DT,PDTOP
!
      REAL,DIMENSION(KMS:KME-1),INTENT(IN) :: DETA1,DETA2
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: DDMPU,DDMPV        &
                                                   ,HBM2,PDSL
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: VTM
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: DIV,T   &
                                                              ,U,V
!----------------------------------------------------------------------
!
!***  LOCAL VARIABLES
!
      INTEGER :: I,IER,J,J4_00,J4_M1,J4_P1,JJ,JKNT,JSTART,K,STAT
!
      REAL :: RDPDX,RDPDY
!
!***  TYPE 4 WORKING ARRAY   ! See PFDHT
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,-1:1) :: CKE,DPDE
!
!----------------------------------------------------------------------
!**********************************************************************
!----------------------------------------------------------------------
      MYJS2  =MAX(JDS+2,JTS)
      MYJE2  =MIN(JDE-2,JTE)
!
      MYIS_P1=MAX(IDS,ITS-1)
      MYIE_P1=MIN(IDE,ITE+1)
      MYIS_P2=MAX(IDS,ITS-2)
      MYIE_P2=MIN(IDE,ITE+2)
!
!----------------------------------------------------------------------
!
!***  MARCH NORTHWARD THROUGH THE SOUTHERNMOST SLABS TO BEGIN
!***  FILLING THE WORKING ARRAY NEEDED FOR AVERAGING AND
!***  DIFFERENCING IN J
!
!----------------------------------------------------------------------
      JSTART=MYJS2
!
      DO J=-1,0
        JJ=JSTART+J
!
        DO K=KTS,KTE
        DO I=MYIS_P2,MYIE_P2
          DPDE(I,K,J)=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,JJ)
          DIV(I,K,JJ)=DIV(I,K,JJ)*HBM2(I,JJ)
        ENDDO
        ENDDO
!
      ENDDO
!
      JKNT=0
!----------------------------------------------------------------------
!
      main_integration : DO J=MYJS2,MYJE2
!
!----------------------------------------------------------------------
!***
!***  SET THE 3RD INDEX OF THE WORKING ARRAYS (SEE SUBROUTINE INIT
!***                                           AND PFDHT DIAGRAMS)
!***
!***  J[TYPE]_NN WHERE "TYPE" IS THE WORKING ARRAY TYPE SEEN IN THE
!***  LOCAL DECLARATION ABOVE (DEPENDENT UPON THE J EXTENT) AND
!***  NN IS THE NUMBER OF ROWS NORTH OF THE CENTRAL ROW WHOSE J IS
!***  THE CURRENT VALUE OF THE main_integration LOOP.
!***  (P2 denotes +2, etc.)
!***
      JKNT=JKNT+1
!
      J4_P1=INDX3_WRK(1,JKNT,4)
      J4_00=INDX3_WRK(0,JKNT,4)
      J4_M1=INDX3_WRK(-1,JKNT,4)
!
!----------------------------------------------------------------------
      DO K=KTS,KTE
      DO I=MYIS_P1,MYIE_P1
        DPDE(I,K,J4_P1)=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J+1)
        DIV(I,K,J+1)=DIV(I,K,J+1)*HBM2(I,J+1)
      ENDDO
      ENDDO
!
      DO K=KTS,KTE
      DO I=MYIS1_P1,MYIE1_P1
        RDPDX=VTM(I,K,J)/(DPDE(I+IVW(J),K,J4_00)                       &
                         +DPDE(I+IVE(J),K,J4_00))
        U(I,K,J)=U(I,K,J)+(DIV(I+IVE(J),K,J)-DIV(I+IVW(J),K,J))        &
                          *RDPDX*DDMPU(I,J)
!
        RDPDY=VTM(I,K,J)/(DPDE(I,K,J4_M1)+DPDE(I,K,J4_P1))
        V(I,K,J)=V(I,K,J)+(DIV(I,K,J+1)-DIV(I,K,J-1))                  &
                          *RDPDY*DDMPV(I,J)
      ENDDO
      ENDDO
!
!----------------------------------------------------------------------
!
      ENDDO main_integration
!
!----------------------------------------------------------------------
      END SUBROUTINE DDAMP
!----------------------------------------------------------------------
      END MODULE MODULE_IGWAVE_ADJUST

#else

! This is the version that was developed with the first version of 
! WRF_NMM. It is restructured so that the J-loop is outer.
!
!-----------------------------------------------------------------------
!

      MODULE MODULE_IGWAVE_ADJUST 2
!
!-----------------------------------------------------------------------
      USE MODULE_MODEL_CONSTANTS
      USE MODULE_MPP
      USE MODULE_TIMERS  ! this one creates a name conflict at compile time
!-----------------------------------------------------------------------
!***
!***  SPECIFY THE NUMBER OF TIMES TO SMOOTH THE VERTICAL VELOCITY
!***  AND THE NUMBER OF ROWS FROM THE NORTHERN AND SOUTHERN EDGES
!***  OF THE GLOBAL DOMAIN BEYOND WHICH THE SMOOTHING DOES NOT GO
!***  FOR SUBROUTINE PDTEDT
!
      INTEGER :: KSMUD=0,LNSDT=7
!
!-----------------------------------------------------------------------
!
      CONTAINS
!
!***********************************************************************

      SUBROUTINE PFDHT(NTSD,LAST_TIME,PT,DETA1,DETA2,PDTOP,RES,FIS      & 1,1
     &                ,HYDRO,SIGMA,FIRST,DX,DY                          &
     &                ,HTM,HBM2,VTM,VBM2,VBM3                           &
     &                ,FDIV,FCP,WPDAR,DFL,CPGFU,CPGFV                   &
     &                ,PD,PDSL,T,Q,U,V,CWM,OMGALF,PINT,DWDT             &
     &                ,RTOP,DIV,FEW,FNS,FNE,FSE                         &
     &                ,IHE,IHW,IVE,IVW,INDX3_WRK                        &
     &                ,IDS,IDE,JDS,JDE,KDS,KDE                          &
     &                ,IMS,IME,JMS,JME,KMS,KME                          &
     &                ,ITS,ITE,JTS,JTE,KTS,KTE)
!***********************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    PFDHT       DIVERGENCE/HORIZONTAL OMEGA-ALPHA
!   PRGRMMR: JANJIC          ORG: W/NP22     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
!      COORDINATE 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
!   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
!   02-02-01  BLACK      - REWRITTEN FOR WRF CODING STANDARDS
!
! USAGE: CALL PFDHT FROM MAIN PROGRAM SOLVE_RUNSTREAM
!   INPUT ARGUMENT LIST:
!
!   OUTPUT ARGUMENT LIST:
!
!   OUTPUT FILES:
!     NONE
!
!   SUBPROGRAMS CALLED:
!
!     UNIQUE: NONE
!
!     LIBRARY: NONE
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$  
!----------------------------------------------------------------------
!**********************************************************************
!----------------------------------------------------------------------
      IMPLICIT NONE
!----------------------------------------------------------------------
#ifdef DM_PARALLEL
      INCLUDE "mpif.h"
#endif
!----------------------------------------------------------------------
      LOGICAL,INTENT(IN) :: FIRST,HYDRO
      INTEGER,INTENT(IN) :: SIGMA
!
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
                           ,IMS,IME,JMS,JME,KMS,KME                    &
                           ,ITS,ITE,JTS,JTE,KTS,KTE
!
      INTEGER, DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
!
!***  NMM_MAX_DIM is set in configure.wrf and must agree with
!***  the value of dimspec q in the Registry/Registry
!
      INTEGER,DIMENSION(-3:3,NMM_MAX_DIM,0:6),INTENT(IN) :: indx3_wrk
!
      INTEGER,INTENT(IN) :: NTSD
      LOGICAL,INTENT(IN) :: LAST_TIME
!
      REAL,INTENT(IN) :: CPGFV,DY,PDTOP,PT
!
      REAL,DIMENSION(KMS:KME-1),INTENT(IN) :: DETA1,DETA2
!
      REAL,DIMENSION(KMS:KME),INTENT(IN) :: DFL
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CPGFU,DX,FCP,FDIV  &
                                                   ,PD,FIS,RES,WPDAR   &
                                                   ,HBM2,VBM2,VBM3
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: CWM,DWDT   &
                                                           ,Q,HTM,VTM
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: PINT
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: DIV     &
                                                              ,OMGALF  &
                                                              ,RTOP,T  &
                                                              ,U,V
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: FEW,FNS   &
                                                            ,FNE,FSE
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: PDSL

!----------------------------------------------------------------------
!
!***  LOCAL VARIABLES
!
      INTEGER :: I,J,JJ,JKNT,JSTART,K
      INTEGER :: J1_00,J1_M1,J1_P1,J1_P2
      INTEGER :: J2_00,J2_M1,J2_P1
      INTEGER :: J3_00,J3_P1,J3_P2
      INTEGER :: J4_00,J4_M1,J4_P1
      INTEGER :: J5_00,J5_M1
      INTEGER :: J6_00,J6_P1
!
      REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: ALP1,FILO
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE+1,JTS-5:JTE+5) :: PINTLG
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,JTS-5:JTE+5) :: FIM
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE) :: DIVL,TEW
!
      REAL :: ADPDNE,ADPDSE,ADPDX,ADPDY,APELP,DFI,DCNEK,DCSEK          &
             ,DPFEW,DPFNS,DPFNEK,DPFSEK,DPNEK,DPSEK,EDIV,FIUP          &
             ,HM,PCEW,PCNS,PEW,PNS,PVNEK,PVSEK,RTOPP,VM
!
      REAL :: SLP_STD=101300.0
!
!***  TYPE 1 WORKING ARRAY
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,-2:2) :: APEL,DFDZ,DPDE
!
!***  TYPE 2 WORKING ARRAY
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,-2:1) :: CNE,PCNE,PNE,PPNE
!
!***  TYPE 3 WORKING ARRAY
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,-1:2) :: CSE,PCSE,PPSE,PSE
!
!***  TYPE 4 WORKING ARRAY
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,-1:1) :: PCXC,TNS,UDY,VDX
!
!***  TYPE 5 WORKING ARRAY
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,-1:0) :: TNE
!
!***  TYPE 6 WORKING ARRAY
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE, 0:1) :: TSE
          integer myproc
!----------------------------------------------------------------------
!**********************************************************************
!
!                                       
!                CSE                          CSE            -------  1
!                 *                            *  
!                 *                            *    
!       *******   *                  *******   *   
!      *       *  *                 *       *  *  
!   CNE         * *              CNE         * *       
!               TEW----------OMGALF----------TEW             -------  0
!   CSE         * *              CSE         * *         
!      *       *  *                 *       *  *       
!       *******   *                  *******   *     
!                 *                            *   
!                 *                            * 
!                CNE                          CNE            ------- -1
!                                        
!
!
! 
!**********************************************************************
! 
!                              CSE                           -------  2
!                               *
!                               *
!                               *
!                               *
!                      CNE*****TNS                           -------  1
!                      CSE     | *
!                              | *
!                              | *
!                              | *
!                              | CNE
!                            OMGALF                          -------  0
!                              | CSE
!                              | *
!                              | *
!                              | *
!                      CNE     | *
!                      CSE*****TNS                           ------- -1
!                               *
!                               *
!                               *
!                               *
!                              CNE                           ------- -2
! 
!**********************************************************************
!----------------------------------------------------------------------
!***  PREPARATORY CALCULATIONS
!----------------------------------------------------------------------
! zero out temporaries.
      ALP1=0.;FILO=0.
      PINTLG=0.
      FIM=0.
      DIVL=0.;TEW=0.
      APEL=0.;DFDZ=0.;DPDE=0.
      CNE=0.;PCNE=0.;PNE=0.;PPNE=0.
      CSE=0.;PCSE=0.;PPSE=0.;PSE=0.
      PCXC=0.;TNS=0.;UDY=0.;VDX=0.
      TNE=0.
      TSE=0.
!

        call wrf_get_myproc(myproc)
      PDSL = 0.
      OMGALF = 0.
      DIV = 0.
      IF(SIGMA.EQ.1)THEN
        DO J=MYJS_P5,MYJE_P5
        DO I=MYIS_P5,MYIE_P5
          FILO(I,J)=FIS(I,J)
          PDSL(I,J)=PD(I,J)
        ENDDO
        ENDDO
      ELSE
        DO J=MYJS_P5,MYJE_P5
        DO I=MYIS_P5,MYIE_P5
          FILO(I,J)=0.0
          PDSL(I,J)=RES(I,J)*PD(I,J)
        ENDDO
        ENDDO
      ENDIF
!
      DO J=MYJS_P4,MYJE_P4
        DO K=KTS,KTE
        DO I=MYIS_P4,MYIE_P4
          OMGALF(I,K,J)=0.
          DIV(I,K,J)=0.
        ENDDO
        ENDDO
      ENDDO
!
!----------------------------------------------------------------------
!***
!***  INTEGRATE THE GEOPOTENTIAL
!***
!----------------------------------------------------------------------
!
      DO J=MYJS_P5,MYJE_P5
!
        DO K=KTS,KTE
        DO I=MYIS_P5,MYIE_P5
!
          APELP=(PINT(I,K+1,J)+PINT(I,K,J))*0.5
          RTOPP=(Q(I,K,J)*P608-CWM(I,K,J)+1.)*T(I,K,J)*R_D/APELP

          DFI=RTOPP*(DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J))
!
          RTOP(I,K,J)=RTOPP
          FIUP=FILO(I,J)+DFI
          FIM(I,K,J)=FILO(I,J)+FIUP
          FILO(I,J)=(FIUP-DFL(K+1))*HTM(I,K,J)+DFL(K+1)
        ENDDO
        ENDDO
!
      ENDDO
!
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!***  MARCH NORTHWARD THROUGH THE SOUTHERNMOST SLABS TO BEGIN
!***  FILLING THE MAIN WORKING ARRAYS WHICH ARE MULTI-DIMENSIONED
!***  IN J BECAUSE THEY ARE DIFFERENCED OR AVERAGED IN J
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!
      JSTART=MYJS2_P2
!
      DO J=-2,1
        JJ=JSTART+J
!
        DO K=KTS,KTE
        DO I=MYIS_P5,MYIE_P5
          APELP=0.5*(PINT(I,K+1,JJ)+PINT(I,K,JJ))
          RTOPP=(Q(I,K,JJ)*P608-CWM(I,K,JJ)+1.)*T(I,K,JJ)*R_D/APELP
          APEL(I,K,J)=APELP
          DFDZ(I,K,J)=RTOPP
          DPDE(I,K,J)=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,JJ)
        ENDDO
        ENDDO
!
      ENDDO
!
      DO J=-2,0
        JJ=JSTART+J
!
        DO K=KTS,KTE
        DO I=MYIS_P4,MYIE_P4
             ! 20030508, bit for bit agreement on different P; JM
          CNE(I,K,J)=(DFDZ(I+IHE(JJ),K,J+1)+DFDZ(I,K,J))*2.            &
                    *(APEL(I+IHE(JJ),K,J+1)-APEL(I,K,J))
          CSE(I,K,J+1)=(DFDZ(I+IHE(JJ+1),K,J)+DFDZ(I,K,J+1))*2.        &
                      *(APEL(I+IHE(JJ+1),K,J)-APEL(I,K,J+1))
          PNE(I,K,J)=(FIM(I+IHE(JJ),K,JJ+1)-FIM(I,K,JJ))               &
                    *(DWDT(I+IHE(JJ),K,JJ+1)+DWDT(I,K,JJ))
          PSE(I,K,J+1)=(FIM(I+IHE(JJ+1),K,JJ)-FIM(I,K,JJ+1))           &                
                      *(DWDT(I+IHE(JJ+1),K,JJ)+DWDT(I,K,JJ+1))
          ADPDNE= (DPDE(I+IHE(JJ),K,J+1)+DPDE(I,K,J))
          ADPDSE= (DPDE(I+IHE(JJ+1),K,J)+DPDE(I,K,J+1))
          PCNE(I,K,J)   = CNE(I,K,J)*ADPDNE
          PPNE(I,K,J)   = PNE(I,K,J)*ADPDNE
          PCSE(I,K,J+1) = CSE(I,K,J+1)*ADPDSE
          PPSE(I,K,J+1) = PSE(I,K,J+1)*ADPDSE
        ENDDO
        ENDDO
      ENDDO
!
      IF(.NOT.FIRST)THEN   ! Skip at timestep 0
        J=0
        JJ=JSTART+J
!
        DO K=KTS,KTE
        DO I=MYIS_P2,MYIE1_P2
          DPFNEK=((PPNE(I+IVW(JJ),K,J)+PPNE(I,K,J-1))                  &
                 +(PCNE(I+IVW(JJ),K,J)+PCNE(I,K,J-1)))*2.
          DPFSEK=((PPSE(I+IVW(JJ),K,J)+PPSE(I,K,J+1))                  &
                 +(PCSE(I+IVW(JJ),K,J)+PCSE(I,K,J+1)))*2.
          DPFEW=DPFNEK+DPFSEK
          DPFNS=DPFNEK-DPFSEK
          ADPDX=DPDE(I+IVW(JJ),K,J)+DPDE(I+IVE(JJ),K,J)
          ADPDY=DPDE(I,K,J-1)+DPDE(I,K,J+1)
          DPNEK=PNE(I+IVW(JJ),K,J)+PNE(I,K,J-1)
          DPSEK=PSE(I+IVW(JJ),K,J)+PSE(I,K,J+1)
          PEW=DPNEK+DPSEK
          PNS=DPNEK-DPSEK
          DCNEK=CNE(I+IVW(JJ),K,J)+CNE(I,K,J-1)
          DCSEK=CSE(I+IVW(JJ),K,J)+CSE(I,K,J+1)
          PCEW=(DCNEK+DCSEK)*ADPDX
          PCNS=(DCNEK-DCSEK)*ADPDY
          VM=VTM(I,K,JJ)*VBM2(I,JJ)
          U(I,K,JJ)=(((DPFEW+PCEW)/ADPDX+PEW)*CPGFU(I,JJ))*VM+U(I,K,JJ)
          V(I,K,JJ)=(((DPFNS+PCNS)/ADPDY+PNS)*CPGFV      )*VM+V(I,K,JJ)
        ENDDO
        ENDDO
      ENDIF
!
      DO J=-1,0
        JJ=JSTART+J
!
        DO K=KTS,KTE
        DO I=MYIS_P3,MYIE_P3
          UDY(I,K,J)=DY*U(I,K,JJ)
          VDX(I,K,J)=DX(I,JJ)*V(I,K,JJ)
          DCNEK=CNE(I+IVW(JJ),K,J)+CNE(I,K,J-1)
          DCSEK=CSE(I+IVW(JJ),K,J)+CSE(I,K,J+1)
          ADPDY=DPDE(I,K,J-1)+DPDE(I,K,J+1)
          TNS(I,K,J)=VDX(I,K,J)*((DCNEK-DCSEK)*ADPDY)
          FNS(I,K,JJ)=VDX(I,K,J)*ADPDY
          PCXC(I,K,J)=(PNE(I+IVW(JJ),K,J)-PNE(I,K,J-1)                 &
                      +CNE(I+IVW(JJ),K,J)-CNE(I,K,J-1)                 &
                      +PSE(I+IVW(JJ),K,J)-PSE(I,K,J+1)                 &
                      +CSE(I+IVW(JJ),K,J)-CSE(I,K,J+1))                &
                      *VBM3(I,JJ)*VTM(I,K,JJ)
        ENDDO
        ENDDO
!
      ENDDO
!
      JJ=JSTART
      DO K=KTS,KTE
      DO I=MYIS_P2,MYIE1_P2
        ADPDNE=DPDE(I+IHE(JJ-1),K,0)+DPDE(I,K,-1)
        PVNEK=(UDY(I+IHE(JJ-1),K,-1)+VDX(I+IHE(JJ-1),K,-1))            &
             +(UDY(I,K,0)          +VDX(I,K,0))
        PCNE(I,K,-1)=CNE(I,K,-1)*ADPDNE
        PPNE(I,K,-1)=PNE(I,K,-1)*ADPDNE



        TNE(I,K,-1)=PVNEK*PCNE(I,K,-1)*2.
        FNE(I,K,JJ-1)=PVNEK*ADPDNE
!
        ADPDSE=DPDE(I+IHE(JJ),K,-1)+DPDE(I,K,0)
        PVSEK=(UDY(I+IHE(JJ),K,0)-VDX(I+IHE(JJ),K,0))                  &
             +(UDY(I,K,-1)      -VDX(I,K,-1))
        PCSE(I,K,0)=CSE(I,K,0)*ADPDSE
        PPSE(I,K,0)=PSE(I,K,0)*ADPDSE
        TSE(I,K,0)=PVSEK*PCSE(I,K,0)*2.
        FSE(I,K,JJ)=PVSEK*ADPDSE
      ENDDO
      ENDDO
!
      JKNT=0
!
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!***  MAIN INTEGRATION LOOP
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!
      main_integration : DO J=MYJS2_P2,MYJE2_P2
!
!----------------------------------------------------------------------
!***
!***  SET THE 3RD INDEX IN THE WORKING ARRAYS (SEE SUBROUTINE INIT
!***                                           AND ABOVE DIAGRAMS)
!***
!***  J[TYPE]_NN WHERE "TYPE" IS THE WORKING ARRAY TYPE SEEN IN THE
!***  LOCAL DECLARATION ABOVE (DEPENDENT UPON THE J EXTENT) AND
!***  NN IS THE NUMBER OF ROWS NORTH OF THE CENTRAL ROW WHOSE J IS
!***  THE CURRENT VALUE OF THE main_integration LOOP.
!***  (P2 denotes +2, etc.)
!***
      JKNT=JKNT+1
!
      J1_P2=INDX3_WRK(2,JKNT,1)
      J1_P1=INDX3_WRK(1,JKNT,1)
      J1_00=INDX3_WRK(0,JKNT,1)
      J1_M1=INDX3_WRK(-1,JKNT,1)
!
      J2_P1=INDX3_WRK(1,JKNT,2)
      J2_00=INDX3_WRK(0,JKNT,2)
      J2_M1=INDX3_WRK(-1,JKNT,2)
!
      J3_P2=INDX3_WRK(2,JKNT,3)
      J3_P1=INDX3_WRK(1,JKNT,3)
      J3_00=INDX3_WRK(0,JKNT,3)
!
      J4_P1=INDX3_WRK(1,JKNT,4)
      J4_00=INDX3_WRK(0,JKNT,4)
      J4_M1=INDX3_WRK(-1,JKNT,4)
!
      J5_00=INDX3_WRK(0,JKNT,5)
      J5_M1=INDX3_WRK(-1,JKNT,5)
!
      J6_P1=INDX3_WRK(1,JKNT,6)
      J6_00=INDX3_WRK(0,JKNT,6)
!
!----------------------------------------------------------------------
      DO K=KTS,KTE
      DO I=MYIS_P5,MYIE_P5
        APELP=0.5*(PINT(I,K+1,J+2)+PINT(I,K,J+2))
        RTOPP=(Q(I,K,J+2)*P608-CWM(I,K,J+2)+1.)*T(I,K,J+2)*R_D/APELP
        APEL(I,K,J1_P2)=APELP
        DFDZ(I,K,J1_P2)=RTOPP
        DPDE(I,K,J1_P2)=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J+2)
      ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!***  DIAGONAL CONTRIBUTIONS TO PRESSURE GRADIENT FORCE
!-----------------------------------------------------------------------
!
      DO K=KTS,KTE
      DO I=MYIS_P4,MYIE_P4

        CNE(I,K,J2_P1)=(DFDZ(I+IHE(J+1),K,J1_P2)+DFDZ(I,K,J1_P1))*2.    &
                      *(APEL(I+IHE(J+1),K,J1_P2)-APEL(I,K,J1_P1))
        CSE(I,K,J3_P2)=(DFDZ(I+IHE(J+2),K,J1_P1)+DFDZ(I,K,J1_P2))*2.    &
                      *(APEL(I+IHE(J+2),K,J1_P1)-APEL(I,K,J1_P2))
!
        PNE(I,K,J2_P1)=(FIM(I+IHE(J+1),K,J+2)-FIM(I,K,J+1))             &
                      *(DWDT(I+IHE(J+1),K,J+2)+DWDT(I,K,J+1))
        PSE(I,K,J3_P2)=(FIM(I+IHE(J+2),K,J+1)-FIM(I,K,J+2))             &
                      *(DWDT(I+IHE(J+2),K,J+1)+DWDT(I,K,J+2))
!
        ADPDNE=DPDE(I+IHE(J+1),K,J1_P2)+DPDE(I,K,J1_P1)
        ADPDSE=DPDE(I+IHE(J+2),K,J1_P1)+DPDE(I,K,J1_P2)
        PCNE(I,K,J2_P1)=CNE(I,K,J2_P1)*ADPDNE
        PPNE(I,K,J2_P1)=PNE(I,K,J2_P1)*ADPDNE
        PCSE(I,K,J3_P2)=CSE(I,K,J3_P2)*ADPDSE
        PPSE(I,K,J3_P2)=PSE(I,K,J3_P2)*ADPDSE
      ENDDO
      ENDDO
!
!----------------------------------------------------------------------
!***  CONTINUITY EQUATION MODIFICATION
!----------------------------------------------------------------------
!
      DO K=KTS,KTE
      DO I=MYIS_P1,MYIE_P1
        PCXC(I,K,J4_P1)=(PNE(I+IVW(J+1),K,J2_P1)                       &
                        +CNE(I+IVW(J+1),K,J2_P1)                       &
                        +PSE(I+IVW(J+1),K,J3_P1)                       &
                        +CSE(I+IVW(J+1),K,J3_P1)                       &
                        -PNE(I,K,J2_00)                                &
                        -CNE(I,K,J2_00)                                &
                        -PSE(I,K,J3_P2)                                &
                        -CSE(I,K,J3_P2))                               &
                        *VBM3(I,J+1)*VTM(I,K,J+1)
      ENDDO
      ENDDO
!
!----------------------------------------------------------------------
!
      DO K=KTS,KTE
      DO I=MYIS1,MYIE1
        DIVL(I,K)=(DETA1(K)*PDTOP/(SLP_STD-PT)                         &
                  +DETA2(K)*(1.-PDTOP/(SLP_STD-PT)))*WPDAR(I,J)        &
                 *(PCXC(I+IHE(J),K,J4_00)-PCXC(I,K,J4_P1)              &
                  +PCXC(I+IHW(J),K,J4_00)-PCXC(I,K,J4_M1))
      ENDDO
      ENDDO
!
!----------------------------------------------------------------------
!
      IF(.NOT.FIRST)THEN     ! Skip at timestep 0
!
!----------------------------------------------------------------------
!***  LAT & LONG PRESSURE FORCE COMPONENTS
!----------------------------------------------------------------------
!
        DO K=KTS,KTE
        DO I=MYIS_P2,MYIE1_P2
          DPNEK=PNE(I+IVW(J+1),K,J2_P1)+PNE(I,K,J2_00)
          DPSEK=PSE(I+IVW(J+1),K,J3_P1)+PSE(I,K,J3_P2)
          PEW=DPNEK+DPSEK
          PNS=DPNEK-DPSEK
!
          ADPDX=DPDE(I+IVW(J+1),K,J1_P1)+DPDE(I+IVE(J+1),K,J1_P1)
          ADPDY=DPDE(I,K,J1_00)+DPDE(I,K,J1_P2)
          DCNEK=CNE(I+IVW(J+1),K,J2_P1)+CNE(I,K,J2_00)
          DCSEK=CSE(I+IVW(J+1),K,J3_P1)+CSE(I,K,J3_P2)
          PCEW=(DCNEK+DCSEK)*ADPDX
          PCNS=(DCNEK-DCSEK)*ADPDY
!
          DPFNEK=((PPNE(I+IVW(J+1),K,J2_P1)+PPNE(I,K,J2_00))           &
                 +(PCNE(I+IVW(J+1),K,J2_P1)+PCNE(I,K,J2_00)))*2.
          DPFSEK=((PPSE(I+IVW(J+1),K,J3_P1)+PPSE(I,K,J3_P2))           &
                 +(PCSE(I+IVW(J+1),K,J3_P1)+PCSE(I,K,J3_P2)))*2.
          DPFEW=DPFNEK+DPFSEK
          DPFNS=DPFNEK-DPFSEK
!
!----------------------------------------------------------------------
!***  UPDATE U AND V FOR PRESSURE GRADIENT FORCE
!----------------------------------------------------------------------
!
          VM=VTM(I,K,J+1)*VBM2(I,J+1)
          U(I,K,J+1)=(((DPFEW+PCEW)/ADPDX+PEW)*CPGFU(I,J+1))*VM        &
                    +U(I,K,J+1) 
          V(I,K,J+1)=(((DPFNS+PCNS)/ADPDY+PNS)*CPGFV       )*VM        &
                    +V(I,K,J+1)
        ENDDO
        ENDDO
!----------------------------------------------------------------------
!
      ENDIF    !End of IF block executed for FIRST equal to .FALSE.
!
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!
      IF(.NOT.LAST_TIME)THEN    !Do not execute block at last timestep
!
!----------------------------------------------------------------------
        DO K=KTS,KTE
        DO I=MYIS_P5,MYIE_P5
          UDY(I,K,J4_P1)=DY*U(I,K,J+1)
          VDX(I,K,J4_P1)=DX(I,J+1)*V(I,K,J+1)
        ENDDO
        ENDDO
!
!----------------------------------------------------------------------
!***  LAT & LON FLUXES & OMEGA-ALPHA COMPONENTS
!----------------------------------------------------------------------
!
        DO K=KTS,KTE
        DO I=MYIS_P3,MYIE_P3
          ADPDX=DPDE(I+IVW(J),K,J1_00)+DPDE(I+IVE(J),K,J1_00)
          DCNEK=CNE(I+IVW(J),K,J2_00)+CNE(I,K,J2_M1)
          DCSEK=CSE(I+IVW(J),K,J3_00)+CSE(I,K,J3_P1)
          TEW(I,K)=UDY(I,K,J4_00)*((DCNEK+DCSEK)*ADPDX)
          FEW(I,K,J)=UDY(I,K,J4_00)*ADPDX
!
          ADPDY=DPDE(I,K,J1_P2)+DPDE(I,K,J1_00)
          DCNEK=CNE(I+IVW(J+1),K,J2_P1)+CNE(I,K,J2_00)
          DCSEK=CSE(I+IVW(J+1),K,J3_P1)+CSE(I,K,J3_P2)
          TNS(I,K,J4_P1)=VDX(I,K,J4_P1)*((DCNEK-DCSEK)*ADPDY)
          FNS(I,K,J+1)=VDX(I,K,J4_P1)*ADPDY
!
!----------------------------------------------------------------------
!***  DIAGONAL FLUXES AND DIAGONALLY AVERAGED WIND
!----------------------------------------------------------------------
!
          ADPDNE=DPDE(I+IHE(J),K,J1_P1)+DPDE(I,K,J1_00)
          PVNEK=(UDY(I+IHE(J),K,J4_00)+VDX(I+IHE(J),K,J4_00))          &
               +(UDY(I,K,J4_P1)       +VDX(I,K,J4_P1))
          TNE(I,K,J5_00)=PVNEK*PCNE(I,K,J2_00)*2.
          FNE(I,K,J)=PVNEK*ADPDNE
!
          ADPDSE=DPDE(I+IHE(J+1),K,J1_00)+DPDE(I,K,J1_P1)
          PVSEK=(UDY(I+IHE(J+1),K,J4_P1)-VDX(I+IHE(J+1),K,J4_P1))      &
               +(UDY(I,K,J4_00)         -VDX(I,K,J4_00))
          TSE(I,K,J6_P1)=PVSEK*PCSE(I,K,J3_P1)*2.
          FSE(I,K,J+1)=PVSEK*ADPDSE
        ENDDO
        ENDDO
!
!----------------------------------------------------------------------
!***  HORIZONTAL PART OF OMEGA-ALPHA & DIVERGENCE
!----------------------------------------------------------------------
!
        DO K=KTS,KTE
        DO I=MYIS1_P1,MYIE1_P1
          HM=HTM(I,K,J)*HBM2(I,J)
          OMGALF(I,K,J)=(TEW(I+IHE(J),K)+TEW(I+IHW(J),K)               &
                        +TNS(I,K,J4_P1) +TNS(I,K,J4_M1)                &
                        +TNE(I,K,J5_00) +TNE(I+IHW(J),K,J5_M1)         &
                        +TSE(I,K,J6_00)+TSE(I+IHW(J),K,J6_P1))         &
                        /DPDE(I,K,J1_00)*FCP(I,J)*HM
          T(I,K,J)=OMGALF(I,K,J)+T(I,K,J)
          EDIV=(FEW(I+IHE(J),K,J)+FNS(I,K,J+1)                         &
               +FNE(I,K,J)+FSE(I,K,J)                                  &
              -(FEW(I+IHW(J),K,J)+FNS(I,K,J-1)                         &
               +FNE(I+IHW(J),K,J-1)+FSE(I+IHW(J),K,J+1)))*FDIV(I,J)
          DIV(I,K,J)=(EDIV+DIVL(I,K))*HM
        ENDDO
        ENDDO
!----------------------------------------------------------------------
!
      ENDIF   !End block to skip execution at last timestep
!
!----------------------------------------------------------------------
!
      ENDDO main_integration
!
!----------------------------------------------------------------------
!
      END SUBROUTINE PFDHT
!
!----------------------------------------------------------------------
!**********************************************************************

      SUBROUTINE PDTEDT(                                               & 1,4
#ifdef DM_PARALLEL
                        grid,                                          &
#endif
                        NTSD,DT,PT,ETA2,RES,HYDRO                      &
                       ,HTM,HBM2,EF4T                                  &
                       ,PD,PDSL,PDSLO,T,DWDT,RTOP,OMGALF               &
                       ,PETDT,PINT,DIV,PSDT                            &
                       ,IHE,IHW,IVE,IVW,INDX3_WRK                      &                 
                       ,IDS,IDE,JDS,JDE,KDS,KDE                        &
                       ,IMS,IME,JMS,JME,KMS,KME                        &
                       ,ITS,ITE,JTS,JTE,KTS,KTE)
!**********************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .     
! SUBPROGRAM:    PDTEDT      SURFACE PRESSURE TENDENCY CALC
!   PRGRMMR: JANJIC          ORG: W/NP2      DATE: 96-07-??      
!     
! ABSTRACT:
!     PDTEDT VERTICALLY INTEGRATES THE MASS FLUX DIVERGENCE TO
!     OBTAIN THE SURFACE PRESSURE TENDENCY AND VERTICAL VELOCITY ON
!     THE LAYER INTERFACES.  THEN IT UPDATES THE HYDROSTATIC SURFACE
!     PRESSURE, THE NONHYDROSTATIC PRESSURE, AND ADDS THE LOCAL TIME
!     DERIVATIVE AND VERTICAL ADVECTION OF NONHYDROSTATIC PRESSURE
!     CONTRIBUTION TO THE OMEGA-ALPHA TERM OF THE THERMODYNAMIC
!     EQUATION.  ALSO, THE OMEGA-ALPHA TERM IS COMPUTED FOR DIAGNOSTICS.
!     
! PROGRAM HISTORY LOG:
!   87-06-??  JANJIC     - ORIGINATOR
!   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
!   96-05-??  JANJIC     - ADDED NONHYDROSTATIC EFFECTS & MERGED THE
!                          PREVIOUS SUBROUTINES PDTE & PDNEW
!   00-01-03  BLACK      - DISTRIBUTED MEMORY AND THREADS
!   01-02-23  BLACK      - CONVERTED TO WRF FORMAT
!   01-04-11  BLACK      - REWRITTEN FOR WRF CODING STANDARDS
!     
! USAGE: CALL PDTEDT FROM SUBROUTINE SOLVE_RUNSTREAM
!   INPUT ARGUMENT LIST:
!  
!   OUTPUT ARGUMENT LIST: 
!     
!   OUTPUT FILES:
!     NONE
!     
!   SUBPROGRAMS CALLED:
!  
!     UNIQUE: NONE
!  
!     LIBRARY: NONE
!  
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$  
!***********************************************************************
#ifdef DM_PARALLEL
      USE module_domain
      USE module_dm
#endif
!-----------------------------------------------------------------------
      IMPLICIT NONE
!-----------------------------------------------------------------------
#ifdef DM_PARALLEL
      INCLUDE "mpif.h"
      TYPE (DOMAIN) :: grid
#endif
!-----------------------------------------------------------------------
      LOGICAL,INTENT(IN) :: HYDRO
!
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                     &
                           ,IMS,IME,JMS,JME,KMS,KME                     &
                           ,ITS,ITE,JTS,JTE,KTS,KTE
!
      INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
!
!***  NMM_MAX_DIM is set in configure.wrf and must agree with
!***  the value of dimspec q in the Registry/Registry
!
      INTEGER,DIMENSION(-3:3,NMM_MAX_DIM,0:6),INTENT(IN) :: INDX3_WRK
!
      INTEGER,INTENT(IN) :: NTSD
!
      REAL,INTENT(IN) :: DT,EF4T,PT
!
      REAL,DIMENSION(KMS:KME),INTENT(IN) :: ETA2
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: RES,HBM2   
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DWDT,HTM
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: DIV      & 
                                                              ,OMGALF   &
                                                              ,RTOP,T
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: PD,PDSL
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: PINT
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: PETDT
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: PDSLO,PSDT
!
!-----------------------------------------------------------------------
!
!***  LOCAL VARIABLES
!
      INTEGER :: I,IHH,IHL,IX,J,JHH,JHL,JJ,JX,K,KNT,KS,NSMUD
      INTEGER :: J1_00,J1_M1,J2_00,J2_P1
      INTEGER :: MY_IS_GLB,MY_IE_GLB,MY_JS_GLB,MY_JE_GLB
!
      REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: APDT,HBMS,PPDT,PRET,TPM
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,-1:0) :: PNE
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE, 0:1) :: PSE
!
      REAL :: DWDTP,PETDTL,RHS,TPMP
!
!jm      REAL(8) :: timef
!-----------------------------------------------------------------------
!**********************************************************************
!----------------------------------------------------------------------
!
!jw   CALL ZERO(PDSLO,IMS,IME,JMS,JME)
      DO J=JMS,JME
      DO I=IMS,IME
        PDSLO(I,J)=0.
      ENDDO
      ENDDO
!
      MY_IS_GLB=ITS
      MY_IE_GLB=ITE
      MY_JS_GLB=JTS
      MY_JE_GLB=JTE
      MYIS=ITS
      MYIE=ITE
      MYJS=JTS
      MYJE=JTE
!
      MYIS1  =MAX(IDS+1,ITS)
      MYIE1  =MIN(IDE-1,ITE)
      MYJS2  =MAX(JDS+2,JTS)
      MYJE2  =MIN(JDE-2,JTE)
!
      MYIS_P1=MAX(IDS,ITS-1)
      MYIE_P1=MIN(IDE,ITE+1)
      MYIS_P2=MAX(IDS,ITS-2)
      MYIE_P2=MIN(IDE,ITE+2)
      MYJS_P2=MAX(JDS,JTS-2)
      MYJE_P2=MIN(JDE,JTE+2)
!----------------------------------------------------------------------
!***  COMPUTATION OF PRESSURE TENDENCY & PREPARATIONS
!----------------------------------------------------------------------
!
      DO J=MYJS_P2,MYJE_P2
        DO K=KTE-1,KTS,-1
        DO I=MYIS_P2,MYIE_P2
          DIV(I,K,J)=DIV(I,K+1,J)+DIV(I,K,J)
        ENDDO
        ENDDO
      ENDDO
!----------------------------------------------------------------------
      DO J=MYJS_P2,MYJE_P2
      DO I=MYIS_P2,MYIE_P2
        PSDT(I,J)=-DIV(I,KTS,J)
        APDT(I,J)=PSDT(I,J)
        PPDT(I,J)=PSDT(I,J)
        PDSLO(I,J)=PDSL(I,J)
      ENDDO
      ENDDO
!----------------------------------------------------------------------
      DO J=MYJS_P2,MYJE_P2
      DO I=MYIS_P2,MYIE_P2
        PD(I,J)=PSDT(I,J)*DT+PD(I,J)
        PRET(I,J)=PSDT(I,J)*RES(I,J)
        PDSL(I,J)=PD(I,J)*RES(I,J)
!
        PINT(I,KTE+1,J)=PT
!
        TPM(I,J)=PT+PINT(I,KTE,J)
      ENDDO
      ENDDO
!----------------------------------------------------------------------
!***  COMPUTATION OF PETDT
!----------------------------------------------------------------------
      DO J=MYJS_P2,MYJE_P2
        DO K=KTE-1,KTS,-1
        DO I=MYIS_P2,MYIE_P2
          PETDT(I,K,J)=-(PRET(I,J)*ETA2(K+1)+DIV(I,K+1,J))             &
                        *HTM(I,K,J)*HBM2(I,J)
        ENDDO
        ENDDO
      ENDDO
!----------------------------------------------------------------------
!***  KINETIC ENERGY GENERATION TERMS IN T EQUATION
!----------------------------------------------------------------------
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        DWDTP=DWDT(I,KTE,J)
        TPMP=PINT(I,KTE,J)+PINT(I,KTE-1,J)
!
        RHS=-DIV(I,KTE,J)*RTOP(I,KTE,J)*HTM(I,KTE,J)*DWDTP*EF4T
        OMGALF(I,KTE,J)=OMGALF(I,KTE,J)+RHS
        T(I,KTE,J)=RHS*HBM2(I,J)+T(I,KTE,J)
        PINT(I,KTE,J)=PRET(I,J)*(ETA2(KTE+1)+ETA2(KTE))*DWDTP*DT       &
                   +TPM(I,J)-PINT(I,KTE+1,J)
!
        TPM(I,J)=TPMP
      ENDDO
      ENDDO
!----------------------------------------------------------------------
      DO J=MYJS,MYJE
        DO K=KTE-1,KTS+1,-1
        DO I=MYIS,MYIE
          DWDTP=DWDT(I,K,J)
          TPMP=PINT(I,K,J)+PINT(I,K-1,J)
!
          RHS=-(DIV(I,K+1,J)+DIV(I,K,J))*RTOP(I,K,J)*HTM(I,K,J)*DWDTP  &
               *EF4T
          OMGALF(I,K,J)=OMGALF(I,K,J)+RHS
          T(I,K,J)=RHS*HBM2(I,J)+T(I,K,J)
          PINT(I,K,J)=PRET(I,J)*(ETA2(K+1)+ETA2(K))*DWDTP*DT           &
                     +TPM(I,J)-PINT(I,K+1,J)
!
          TPM(I,J)=TPMP
        ENDDO
        ENDDO
      ENDDO
!----------------------------------------------------------------------
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
!
        DWDTP=DWDT(I,KTS,J)
!
        RHS=-(DIV(I,KTS+1,J)+DIV(I,KTS,J))*RTOP(I,KTS,J)*HTM(I,KTS,J)  &
             *DWDTP*EF4T
        OMGALF(I,KTS,J)=OMGALF(I,KTS,J)+RHS
        T(I,KTS,J)=RHS*HBM2(I,J)+T(I,KTS,J)
        PINT(I,KTS,J)=PRET(I,J)*(ETA2(KTS+1)+ETA2(KTS))*DWDTP*DT       &
                       +TPM(I,J)-PINT(I,KTS+1,J)
      ENDDO
      ENDDO
!----------------------------------------------------------------------
!***  SMOOTHING VERTICAL VELOCITY ALONG BOUNDARIES
!----------------------------------------------------------------------
      nonhydrostatic_smoothing: IF(.NOT.HYDRO.AND.KSMUD.GT.0)THEN
!
        NSMUD=KSMUD
!
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          HBMS(I,J)=HBM2(I,J)
        ENDDO
        ENDDO
!
        JHL=LNSDT
        JHH=JDE-JHL+1
!
        DO J=JHL,JHH
          IF(J.GE.MY_JS_GLB.AND.J.LE.MY_JE_GLB)THEN
            IHL=JHL/2+1
            IHH=IDE-IHL+MOD(J,2)
!
            DO I=IHL,IHH
              IF(I.GE.MY_IS_GLB.AND.I.LE.MY_IE_GLB)THEN
                IX=I    ! -MY_IS_GLB+1
                JX=J    ! -MY_JS_GLB+1
                HBMS(IX,JX)=0.
              ENDIF
            ENDDO
!
          ENDIF
        ENDDO
!
!----------------------------------------------------------------------
!***
!***  SMOOTH THE VERTICAL VELOCITY
!***
!----------------------------------------------------------------------
!
        DO KS=1,NSMUD
!
!----------------------------------------------------------------------
!
!***  FILL SOUTHERNMOST SLABS OF THE PNE AND PSE WORKING ARRAYS
!
          JJ=MYJS2-1
          DO K=KTS,KTE-1
          DO I=MYIS_P1,MYIE1_P1
            PNE(I,K,-1)=(PETDT(I+IHE(JJ),K,JJ+1)-PETDT(I,K,JJ))        &
                        *HTM(I,K,JJ)*HTM(I+IHE(JJ),K,JJ+1)
          ENDDO
          ENDDO
!
          DO K=KTS,KTE-1
          DO I=MYIS_P1,MYIE1_P1
            PSE(I,K,0)=(PETDT(I+IHE(JJ+1),K,JJ)-PETDT(I,K,JJ+1))       &
                       *HTM(I+IHE(JJ+1),K,JJ)*HTM(I,K,JJ+1)
          ENDDO
          ENDDO
!
          KNT=0
!
!----------------------------------------------------------------------
!
!***  PROCEED NORTHWARD WITH THE SMOOTHING.
!***  PNE AT H(I,J) LIES BETWEEN (I,J) AND THE H POINT TO THE NE.
!***  PSE AT H(I,J) LIES BETWEEN (I,J) AND THE H POINT TO THE SE.
!
          DO J=MYJS2,MYJE2
!
            KNT=KNT+1
            J1_00=-MOD(KNT+1,2)
            J1_M1=-MOD(KNT,2)
            J2_P1=MOD(KNT,2)
            J2_00=MOD(KNT+1,2)
!
            DO K=KTS,KTE-1
            DO I=MYIS_P1,MYIE1_P1
              PNE(I,K,J1_00)=(PETDT(I+IHE(J),K,J+1)-PETDT(I,K,J))      &
                             *HTM(I,K+1,J)*HTM(I+IHE(J),K+1,J+1)
            ENDDO
            ENDDO
!
            DO K=KTS,KTE-1
            DO I=MYIS_P1,MYIE1_P1
              PSE(I,K,J2_P1)=(PETDT(I+IHE(J+1),K,J)-PETDT(I,K,J+1))    &
                             *HTM(I+IHE(J+1),K+1,J)*HTM(I,K+1,J+1)
            ENDDO
            ENDDO
!
            DO K=KTS,KTE-1
            DO I=MYIS1,MYIE1
              PETDTL=(PNE(I,K,J1_00)-PNE(I+IHW(J),K,J1_M1)             &
                     +PSE(I,K,J2_00)-PSE(I+IHW(J),K,J2_P1))*HBM2(I,J)
              PETDT(I,K,J)=PETDTL*HBMS(I,J)*0.125+PETDT(I,K,J) 
            ENDDO
            ENDDO
!
          ENDDO
!
#ifdef DM_PARALLEL
!jm          btim=timef()
!jw need to add the following arrays to HALO:
!jw     3D - PETDT
!jw stencil needs to go at least 2 point in the x-dir and 2 in y
!jw       CALL EXCH(PETDT,KTE-1,2,2,IMS,IME,JMS,JME,KMS,KME) 
!jw#         include "HALO_NMM_E.inc" not now -- need to move this 
!jw                                   to the mediation layer
!jm          exch_tim=exch_tim+timef()-btim

# include <HALO_NMM_E.inc>


#endif


!----------------------------------------------------------------------
!
        ENDDO  ! End of smoothing loop
!
!----------------------------------------------------------------------
      ENDIF nonhydrostatic_smoothing
!----------------------------------------------------------------------
      END SUBROUTINE PDTEDT
!----------------------------------------------------------------------
!**********************************************************************

      SUBROUTINE DDAMP(NTSD,DT,DETA1,DETA2,PDSL,PDTOP,DIV,HBM2,VTM     & 1,3
                      ,T,U,V,DDMPU,DDMPV                               &
                      ,IHE,IHW,IVE,IVW,INDX3_WRK                       &              
                      ,IDS,IDE,JDS,JDE,KDS,KDE                         &
                      ,IMS,IME,JMS,JME,KMS,KME                         &
                      ,ITS,ITE,JTS,JTE,KTS,KTE)
!**********************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .     
! SUBPROGRAM:    DDAMP       DIVERGENCE DAMPING
!   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 94-03-08       
!     
! ABSTRACT:
!     DDAMP MODIFIES THE WIND COMPONENTS SO AS TO REDUCE THE
!     HORIZONTAL DIVERGENCE.
!     
! PROGRAM HISTORY LOG:
!   87-08-??  JANJIC     - ORIGINATOR
!   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
!   95-03-28  BLACK      - ADDED EXTERNAL EDGE
!   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
!   01-03-12  BLACK      - CONVERTED TO WRF STRUCTURE
!     
! USAGE: CALL DDAMP FROM SUBROUTINE SOLVE_RUNSTREAM
!
!   INPUT ARGUMENT LIST:
!  
!   OUTPUT ARGUMENT LIST: 
!     
!   OUTPUT FILES:
!     NONE
!     
!   SUBPROGRAMS CALLED:
!  
!     UNIQUE: NONE
!  
!     LIBRARY: NONE
!  
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$  
!**********************************************************************
!----------------------------------------------------------------------
      IMPLICIT NONE
!----------------------------------------------------------------------
!
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
                           ,IMS,IME,JMS,JME,KMS,KME                    &
                           ,ITS,ITE,JTS,JTE,KTS,KTE
!
      INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
!
!***  NMM_MAX_DIM is set in configure.wrf and must agree with
!***  the value of dimspec q in the Registry/Registry
!
      INTEGER,DIMENSION(-3:3,NMM_MAX_DIM,0:6),INTENT(IN) :: INDX3_WRK
!
      INTEGER,INTENT(IN) :: NTSD
!
      REAL,INTENT(IN) :: DT,PDTOP
!
      REAL,DIMENSION(KMS:KME-1),INTENT(IN) :: DETA1,DETA2
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: DDMPU,DDMPV        &
                                                   ,HBM2,PDSL
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: VTM
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: DIV,T   &
                                                              ,U,V
!----------------------------------------------------------------------
!
!***  LOCAL VARIABLES
!
      INTEGER :: I,IER,J,J4_00,J4_M1,J4_P1,JJ,JKNT,JSTART,K,STAT
!
      REAL :: RDPDX,RDPDY
!
!***  TYPE 4 WORKING ARRAY   ! See PFDHT
!
      REAL,DIMENSION(ITS-5:ITE+5,KTS:KTE,-1:1) :: CKE,DPDE
!
!----------------------------------------------------------------------
!**********************************************************************
!----------------------------------------------------------------------
      MYJS2  =MAX(JDS+2,JTS)
      MYJE2  =MIN(JDE-2,JTE)
!
      MYIS_P1=MAX(IDS,ITS-1)
      MYIE_P1=MIN(IDE,ITE+1)
      MYIS_P2=MAX(IDS,ITS-2)
      MYIE_P2=MIN(IDE,ITE+2)
!
!----------------------------------------------------------------------
!
!***  MARCH NORTHWARD THROUGH THE SOUTHERNMOST SLABS TO BEGIN
!***  FILLING THE WORKING ARRAY NEEDED FOR AVERAGING AND
!***  DIFFERENCING IN J
!
!----------------------------------------------------------------------
      JSTART=MYJS2
!
      DO J=-1,0
        JJ=JSTART+J
!
        DO K=KTS,KTE
        DO I=MYIS_P2,MYIE_P2
          DPDE(I,K,J)=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,JJ)
          DIV(I,K,JJ)=DIV(I,K,JJ)*HBM2(I,JJ)
        ENDDO
        ENDDO
!
      ENDDO
!
      JKNT=0
!----------------------------------------------------------------------
!
      main_integration : DO J=MYJS2,MYJE2
!
!----------------------------------------------------------------------
!***
!***  SET THE 3RD INDEX OF THE WORKING ARRAYS (SEE SUBROUTINE INIT
!***                                           AND PFDHT DIAGRAMS)
!***
!***  J[TYPE]_NN WHERE "TYPE" IS THE WORKING ARRAY TYPE SEEN IN THE
!***  LOCAL DECLARATION ABOVE (DEPENDENT UPON THE J EXTENT) AND
!***  NN IS THE NUMBER OF ROWS NORTH OF THE CENTRAL ROW WHOSE J IS
!***  THE CURRENT VALUE OF THE main_integration LOOP.
!***  (P2 denotes +2, etc.)
!***
      JKNT=JKNT+1
!
      J4_P1=INDX3_WRK(1,JKNT,4)
      J4_00=INDX3_WRK(0,JKNT,4)
      J4_M1=INDX3_WRK(-1,JKNT,4)
!
!----------------------------------------------------------------------
      DO K=KTS,KTE
      DO I=MYIS_P1,MYIE_P1
        DPDE(I,K,J4_P1)=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J+1)
        DIV(I,K,J+1)=DIV(I,K,J+1)*HBM2(I,J+1)
      ENDDO
      ENDDO
!
      DO K=KTS,KTE
      DO I=MYIS1_P1,MYIE1_P1
        RDPDX=VTM(I,K,J)/(DPDE(I+IVW(J),K,J4_00)                       &
                         +DPDE(I+IVE(J),K,J4_00))
        U(I,K,J)=U(I,K,J)+(DIV(I+IVE(J),K,J)-DIV(I+IVW(J),K,J))        &
                          *RDPDX*DDMPU(I,J)
!
        RDPDY=VTM(I,K,J)/(DPDE(I,K,J4_M1)+DPDE(I,K,J4_P1))
        V(I,K,J)=V(I,K,J)+(DIV(I,K,J+1)-DIV(I,K,J-1))                  &
                          *RDPDY*DDMPV(I,J)
      ENDDO
      ENDDO
!
!----------------------------------------------------------------------
!
      ENDDO main_integration
!
!----------------------------------------------------------------------
      END SUBROUTINE DDAMP
!----------------------------------------------------------------------
      END MODULE MODULE_IGWAVE_ADJUST
#endif