/wrfv2_fire/dyn_nmm/module_ADVECTION.F

Large files files are truncated, but you can click here to view the full file

!----------------------------------------------------------------------
!#define BIT_FOR_BIT
!----------------------------------------------------------------------
#include "nmm_loop_basemacros.h"
#include "nmm_loop_macros.h"
!----------------------------------------------------------------------
!
!NCEP_MESO:MODEL_LAYER: HORIZONTAL AND VERTICAL ADVECTION
!
!----------------------------------------------------------------------
!
      MODULE MODULE_ADVECTION
!
!----------------------------------------------------------------------
      USE MODULE_MODEL_CONSTANTS
      USE MODULE_EXT_INTERNAL
!----------------------------------------------------------------------
#if defined(DM_PARALLEL) && !defined(STUBMPI)
      INCLUDE "mpif.h"
#endif
!----------------------------------------------------------------------
!
      REAL,PARAMETER :: FF2=-0.64813,FF3=0.24520,FF4=-0.12189
      REAL,PARAMETER :: FFC=1.533,FBC=1.-FFC
      REAL :: CONSERVE_MIN=0.9,CONSERVE_MAX=1.1
!
!----------------------------------------------------------------------
!***  CRANK-NICHOLSON OFF-CENTER WEIGHTS FOR CURRENT AND FUTURE
!***  TIME LEVELS.
!-----------------------------------------------------------------------
!
      REAL,PARAMETER :: WGT1=0.90
      REAL,PARAMETER :: WGT2=2.-WGT1
!
!***  FOR CRANK_NICHOLSON CHECK ONLY.
!
      INTEGER :: ITEST=47,JTEST=70
      REAL :: ADTP,ADUP,ADVP,TTLO,TTUP,TULO,TUUP,TVLO,TVUP
!
!----------------------------------------------------------------------
      CONTAINS
!
!*********************************************************************** 
      SUBROUTINE ADVE(NTSD,DT,DETA1,DETA2,PDTOP                         &
     &               ,CURV,F,FAD,F4D,EM_LOC,EMT_LOC,EN,ENT,DX,DY        &
     &               ,HBM2,VBM2                                         &
     &               ,T,U,V,PDSLO,TOLD,UOLD,VOLD                        &
     &               ,PETDT,UPSTRM                                      &
     &               ,FEW,FNS,FNE,FSE                                   &
     &               ,ADT,ADU,ADV                                       &
     &               ,N_IUP_H,N_IUP_V                                   &
     &               ,N_IUP_ADH,N_IUP_ADV                               &
     &               ,IUP_H,IUP_V,IUP_ADH,IUP_ADV                       &
     &               ,IHE,IHW,IVE,IVW                                   &
     &               ,IDS,IDE,JDS,JDE,KDS,KDE                           &
     &               ,IMS,IME,JMS,JME,KMS,KME                           &
     &               ,ITS,ITE,JTS,JTE,KTS,KTE)
!***********************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .     
! SUBPROGRAM:    ADVE        HORIZONTAL AND VERTICAL ADVECTION
!   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 93-10-28       
!     
! ABSTRACT:
!     ADVE CALCULATES THE CONTRIBUTION OF THE HORIZONTAL AND VERTICAL
!     ADVECTION TO THE TENDENCIES OF TEMPERATURE AND WIND AND THEN
!     UPDATES THOSE VARIABLES.
!     THE JANJIC ADVECTION SCHEME FOR THE ARAKAWA E GRID IS USED
!     FOR ALL VARIABLES INSIDE THE FIFTH ROW.  AN UPSTREAM SCHEME
!     IS USED ON ALL VARIABLES IN THE THIRD, FOURTH, AND FIFTH
!     OUTERMOST ROWS.  THE ADAMS-BASHFORTH TIME SCHEME IS USED.
!     
! PROGRAM HISTORY LOG:
!   87-06-??  JANJIC       - ORIGINATOR
!   95-03-25  BLACK        - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
!   96-03-28  BLACK        - ADDED EXTERNAL EDGE
!   98-10-30  BLACK        - MODIFIED FOR DISTRIBUTED MEMORY
!   99-07-    JANJIC       - CONVERTED TO ADAMS-BASHFORTH SCHEME
!                            COMBINING HORIZONTAL AND VERTICAL ADVECTION
!   02-02-04  BLACK        - ADDED VERTICAL CFL CHECK
!   02-02-05  BLACK        - CONVERTED TO WRF FORMAT
!   02-08-29  MICHALAKES   - CONDITIONAL COMPILATION OF MPI
!                            CONVERT TO GLOBAL INDEXING
!   02-09-06  WOLFE        - MORE CONVERSION TO GLOBAL INDEXING
!   04-05-29  JANJIC,BLACK - CRANK-NICHOLSON VERTICAL ADVECTION
!   04-11-23  BLACK        - THREADED 
!   05-12-14  BLACK        - CONVERTED FROM IKJ TO IJK
!     
! USAGE: CALL ADVE FROM SUBROUTINE SOLVE_NMM
!   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         &
                                               ,N_IUP_H,N_IUP_V         &
     &                                         ,N_IUP_ADH,N_IUP_ADV
!
      INTEGER, DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IUP_H,IUP_V     &
     &                                                 ,IUP_ADH,IUP_ADV
!
      INTEGER,INTENT(IN) :: NTSD
!
      REAL,INTENT(IN) :: DT,DY,EN,ENT,F4D,PDTOP
!
      REAL,DIMENSION(NMM_MAX_DIM),INTENT(IN) :: EM_LOC,EMT_LOC
!
      REAL,DIMENSION(KMS:KME),INTENT(IN) :: DETA1,DETA2
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CURV,DX,F,FAD,HBM2  &
     &                                             ,PDSLO,VBM2
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: ADT,ADU,ADV
!
      REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(IN) :: PETDT
!
      REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(INOUT) :: T,TOLD   &
     &                                                        ,U,UOLD   &
     &                                                        ,V,VOLD
!
      REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(OUT) :: FEW,FNE    &
     &                                                      ,FNS,FSE
!
!-----------------------------------------------------------------------
!***  LOCAL VARIABLES
!-----------------------------------------------------------------------
!
      LOGICAL :: UPSTRM
!
      INTEGER :: I,IEND,IFP,IFQ,II,IPQ,ISP,ISQ,ISTART                   &
     &          ,IUP_ADH_J,IVH,IVL                                      &
     &          ,J,J1,JA,JAK,JEND,JGLOBAL,JJ,JKNT,JP2,JSTART            &
     &          ,K,KNTI_ADH,KSTART,KSTOP                                &
     &          ,N,N_IUPH_J,N_IUPADH_J,N_IUPADV_J
!
      INTEGER :: MY_IS_GLB,MY_IE_GLB,MY_JS_GLB,MY_JE_GLB
!
      INTEGER,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: ISPA,ISQA
!
      REAL :: ADPDX,ADPDY,ARRAY3_X,CFL,CFT,CFU,CFV,CMT,CMU,CMV          &
     &       ,DTE,DTQ,F0,F1,F2,F3,FEWP,FNEP,FNSP,FPP,FSEP,HM            &
     &       ,PDOP,PDOPU,PDOPV,PP                                       &
     &       ,PVVLO,PVVLOU,PVVLOV,PVVUP,PVVUPU,PVVUPV                   &
     &       ,QP,RDP,RDPU,RDPV                                          &
     &       ,TEMPA,TEMPB,TTA,TTB,UDY                                   &
     &       ,VDX,VM,VVLO,VVLOU,VVLOV,VVUP,VVUPU,VVUPV
!
      REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: ARRAY0,ARRAY1          &
     &                                          ,ARRAY2,ARRAY3          &
     &                                          ,DPDE,RDPD,RDPDX,RDPDY  &
     &                                          ,TEW,TNE,TNS,TSE,TST    &
     &                                          ,UNE,UNED,UEW,UNS,USE   &
     &                                          ,USED,UST               &
     &                                          ,VEW,VNE,VNS,VSE        &
     &                                          ,VST
!
      REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5,KTS:KTE) :: VAD_TEND_T     &
     &                                                  ,VAD_TEND_U     &
     &                                                  ,VAD_TEND_V
!
      REAL,DIMENSION(KTS:KTE) :: CRT,CRU,CRV,DETA1_PDTOP                &
     &                          ,RCMT,RCMU,RCMV,RSTT,RSTU,RSTV          &
     &                          ,T_K,TN,U_K,UN,V_K,VN
!
!-----------------------------------------------------------------------
!***********************************************************************
!
!                         DPDE      -----  3
!                          |                      J Increasing
!                          |
!                          |                            ^
!                         FNS       -----  2            |
!                          |                            |
!                          |                            |
!                          |                            |
!                         VNS       -----  1            |
!                          |
!                          |
!                          |
!                         ADV       -----  0  ------> Current J
!                          |
!                          |
!                          |
!                         VNS       ----- -1
!                          |
!                          |
!                          |
!                         FNS       ----- -2
!                          |
!                          |
!                          |
!                         DPDE      ----- -3
!
!***********************************************************************
!-----------------------------------------------------------------------
     DO J=JTS-5,JTE+5
     DO I=ITS-5,ITE+5
       ARRAY0(I,J)=0.0
       ARRAY1(I,J)=0.0
       ARRAY2(I,J)=0.0
       ARRAY3(I,J)=0.0
       DPDE(I,J)=0.0
       RDPD(I,J)=0.0
       RDPDX(I,J)=0.0
       RDPDY(I,J)=0.0
       TEW(I,J)=0.0
       TNE(I,J)=0.0
       TNS(I,J)=0.0
       TSE(I,J)=0.0
       TST(I,J)=0.0
       UNE(I,J)=0.0
       UNED(I,J)=0.0
       UEW(I,J)=0.0
       UNS(I,J)=0.0
       USE(I,J)=0.0
       USED(I,J)=0.0
       UST(I,J)=0.0
       VEW(I,J)=0.0
       VNE(I,J)=0.0
       VNS(I,J)=0.0
       VSE(I,J)=0.0
       VST(I,J)=0.0
     ENDDO
     ENDDO
!-----------------------------------------------------------------------
!
      DTQ=DT*0.25
      DTE=DT*(0.5*0.25)
!
!-----------------------------------------------------------------------
!***
!***  PRECOMPUTE DETA1 TIMES PDTOP.
!***
!-----------------------------------------------------------------------
!
      DO K=KTS,KTE
        DETA1_PDTOP(K)=DETA1(K)*PDTOP
      ENDDO
!
!-----------------------------------------------------------------------
!***
!***  INITIALIZE SOME WORKING ARRAYS TO ZERO
!***
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!***  COMPUTE VERTICAL ADVECTION TENDENCIES USING CRANK-NICHOLSON.
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!***  FIRST THE TEMPERATURE
!-----------------------------------------------------------------------
!$omp parallel do                                                       &
!$omp& private(cft,cfu,cfv,cmt,cmu,cmv,crt,cru,crv,i,k                  &
!$omp&        ,pdop,pdopu,pdopv,pvvlo,pvvlou,pvvlov,pvvup,pvvupu,pvvupv &
!$omp&        ,rcmt,rcmu,rcmv,rdp,rdpu,rdpv,rstt,rstu,rstv,t_k,tn       &
!$omp&        ,u_k,un,v_k,vn,vvlo,vvlou,vvlov,vvup,vvupu,vvupv)
!!$omp& private(adtp,adup,advp,ttlo,ttup,tulo,tuup,tvlo,tvup)
!-----------------------------------------------------------------------
!
      main_vertical: DO J=MYJS2,MYJE2
!
!-----------------------------------------------------------------------
!
        iloop_for_t: DO I=MYIS1,MYIE1
!
!-----------------------------------------------------------------------
!***  EXTRACT T FROM THE COLUMN
!-----------------------------------------------------------------------
!
          DO K=KTS,KTE
            T_K(K)=T(I,J,K)
          ENDDO
!
!-----------------------------------------------------------------------
!
          PDOP=PDSLO(I,J)
          PVVLO=PETDT(I,J,KTE-1)*DTQ
          VVLO=PVVLO/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOP)
          CMT=-VVLO*WGT2+1.
          RCMT(KTE)=1./CMT
          CRT(KTE)=VVLO*WGT2
          RSTT(KTE)=-VVLO*WGT1*(T_K(KTE-1)-T_K(KTE))+T_K(KTE)
!
!-----------------------------------------------------------------------
!
          DO K=KTE-1,KTS+1,-1
            RDP=1./(DETA1_PDTOP(K)+DETA2(K)*PDOP)
            PVVUP=PVVLO
            PVVLO=PETDT(I,J,K-1)*DTQ
            VVUP=PVVUP*RDP
            VVLO=PVVLO*RDP
            CFT=-VVUP*WGT2*RCMT(K+1)
            CMT=-CRT(K+1)*CFT+((VVUP-VVLO)*WGT2+1.)
            RCMT(K)=1./CMT
            CRT(K)=VVLO*WGT2
            RSTT(K)=-RSTT(K+1)*CFT+T_K(K)                               &
       &            -(T_K(K)-T_K(K+1))*VVUP*WGT1                        &
       &            -(T_K(K-1)-T_K(K))*VVLO*WGT1
          ENDDO
!
!-----------------------------------------------------------------------
!
          PVVUP=PVVLO
          VVUP=PVVUP/(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOP)
          CFT=-VVUP*WGT2*RCMT(KTS+1)
          CMT=-CRT(KTS+1)*CFT+VVUP*WGT2+1.
          CRT(KTS)=0.
          RSTT(KTS)=-(T_K(KTS)-T_K(KTS+1))*VVUP*WGT1                    &
      &              -RSTT(KTS+1)*CFT+T_K(KTS)
          TN(KTS)=RSTT(KTS)/CMT
          VAD_TEND_T(I,J,KTS)=TN(KTS)-T_K(KTS)
!
          DO K=KTS+1,KTE
            TN(K)=(-CRT(K)*TN(K-1)+RSTT(K))*RCMT(K)
            VAD_TEND_T(I,J,K)=TN(K)-T_K(K)
          ENDDO
!
!-----------------------------------------------------------------------
!***  The following section is only for checking the implicit solution
!***  using back-substitution.  Remove this section otherwise.
!-----------------------------------------------------------------------
!       if(ntsd<=10.or.ntsd>=6000)then
!       IF(I==ITEST.AND.J==JTEST)THEN
!!
!         PVVLO=PETDT(I,J,KTE-1)*DT*0.25
!         VVLO=PVVLO/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOP)
!         TTLO=VVLO*(T(I,J,KTE-1)-T(I,J,KTE)                            &
!    &              +TN(KTE-1)-TN(KTE))
!         ADTP=TTLO+TN(KTE)-T(I,J,KTE)
!         WRITE(0,*)' NTSD=',NTSD,' I=',ITEST,' J=',JTEST,' K=',KTE     &
!    &,             ' ADTP=',ADTP
!         WRITE(0,*)' T=',T(I,J,KTE),' TN=',TN(KTE)                     &
!    &,               ' VAD_TEND_T=',VAD_TEND_T(I,J,KTE)
!         WRITE(0,*)' '
!!
!         DO K=KTE-1,KTS+1,-1
!           RDP=1./(DETA1_PDTOP(K)+DETA2(K)*PDOP)
!           PVVUP=PVVLO
!           PVVLO=PETDT(I,J,K-1)*DT*0.25
!           VVUP=PVVUP*RDP
!           VVLO=PVVLO*RDP
!           TTUP=VVUP*(T(I,J,K)-T(I,J,K+1)+TN(K)-TN(K+1))
!           TTLO=VVLO*(T(I,J,K-1)-T(I,J,K)+TN(K-1)-TN(K))
!           ADTP=TTLO+TTUP+TN(K)-T(I,J,K)
!           WRITE(0,*)' NTSD=',NTSD,' I=',I,' J=',J,' K=',K             &
!    &,               ' ADTP=',ADTP
!           WRITE(0,*)' T=',T(I,J,K),' TN=',TN(K)                       &
!    &,               ' VAD_TEND_T=',VAD_TEND_T(I,J,K)
!           WRITE(0,*)' '
!         ENDDO
!!
!         PVVUP=PVVLO
!         VVUP=PVVUP/(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOP)
!         TTUP=VVUP*(T(I,J,KTS)-T(I,J,KTS+1)+TN(KTS)-TN(KTS+1))
!         ADTP=TTUP+TN(KTS)-T(I,J,KTS)
!         WRITE(0,*)' NTSD=',NTSD,' I=',I,' J=',J,' K=',KTS             &
!    &,             ' ADTP=',ADTP
!         WRITE(0,*)' T=',T(I,J,KTS),' TN=',TN(KTS)                     &
!    &,             ' VAD_TEND_T=',VAD_TEND_T(I,J,KTS)
!         WRITE(0,*)' '
!       ENDIF
!       endif
!
!-----------------------------------------------------------------------
!***  End of check.
!-----------------------------------------------------------------------
!
        ENDDO iloop_for_t
!
!-----------------------------------------------------------------------
!
!***  NOW VERTICAL ADVECTION OF WIND COMPONENTS
!
!-----------------------------------------------------------------------
!
        iloop_for_uv:  DO I=MYIS1,MYIE1
!
!-----------------------------------------------------------------------
!***  EXTRACT U AND V FROM THE COLUMN
!-----------------------------------------------------------------------
!
          DO K=KTS,KTE
            U_K(K)=U(I,J,K)
            V_K(K)=V(I,J,K)
          ENDDO
!
!-----------------------------------------------------------------------
!
          PDOPU=(PDSLO(I+IVW(J),J)+PDSLO(I+IVE(J),J))*0.5
          PDOPV=(PDSLO(I,J-1)+PDSLO(I,J+1))*0.5
          PVVLOU=(PETDT(I+IVW(J),J,KTE-1)+PETDT(I+IVE(J),J,KTE-1))*DTE
          PVVLOV=(PETDT(I,J-1,KTE-1)+PETDT(I,J+1,KTE-1))*DTE
          VVLOU=PVVLOU/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOPU)
          VVLOV=PVVLOV/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOPV)
          CMU=-VVLOU*WGT2+1.
          CMV=-VVLOV*WGT2+1.
          RCMU(KTE)=1./CMU
          RCMV(KTE)=1./CMV
          CRU(KTE)=VVLOU*WGT2
          CRV(KTE)=VVLOV*WGT2
          RSTU(KTE)=-VVLOU*WGT1*(U_K(KTE-1)-U_K(KTE))+U_K(KTE)
          RSTV(KTE)=-VVLOV*WGT1*(V_K(KTE-1)-V_K(KTE))+V_K(KTE)
!
!-----------------------------------------------------------------------
!
          DO K=KTE-1,KTS+1,-1
            RDPU=1./(DETA1_PDTOP(K)+DETA2(K)*PDOPU)
            RDPV=1./(DETA1_PDTOP(K)+DETA2(K)*PDOPV)
            PVVUPU=PVVLOU
            PVVUPV=PVVLOV
            PVVLOU=(PETDT(I+IVW(J),J,K-1)+PETDT(I+IVE(J),J,K-1))*DTE
            PVVLOV=(PETDT(I,J-1,K-1)+PETDT(I,J+1,K-1))*DTE
            VVUPU=PVVUPU*RDPU
            VVUPV=PVVUPV*RDPV
            VVLOU=PVVLOU*RDPU
            VVLOV=PVVLOV*RDPV
            CFU=-VVUPU*WGT2*RCMU(K+1)
            CFV=-VVUPV*WGT2*RCMV(K+1)
            CMU=-CRU(K+1)*CFU+(VVUPU-VVLOU)*WGT2+1.
            CMV=-CRV(K+1)*CFV+(VVUPV-VVLOV)*WGT2+1.
            RCMU(K)=1./CMU
            RCMV(K)=1./CMV
            CRU(K)=VVLOU*WGT2
            CRV(K)=VVLOV*WGT2
            RSTU(K)=-RSTU(K+1)*CFU+U_K(K)                               &
     &              -(U_K(K)-U_K(K+1))*VVUPU*WGT1                       &
     &              -(U_K(K-1)-U_K(K))*VVLOU*WGT1
            RSTV(K)=-RSTV(K+1)*CFV+V_K(K)                               &
     &              -(V_K(K)-V_K(K+1))*VVUPV*WGT1                       &
     &              -(V_K(K-1)-V_K(K))*VVLOV*WGT1
          ENDDO
!
!-----------------------------------------------------------------------
!
          RDPU=1./(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOPU)
          RDPV=1./(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOPV)
          PVVUPU=PVVLOU
          PVVUPV=PVVLOV
          VVUPU=PVVUPU*RDPU
          VVUPV=PVVUPV*RDPV
          CFU=-VVUPU*WGT2*RCMU(KTS+1)
          CFV=-VVUPV*WGT2*RCMV(KTS+1)
          CMU=-CRU(KTS+1)*CFU+VVUPU*WGT2+1.
          CMV=-CRV(KTS+1)*CFV+VVUPV*WGT2+1.
          CRU(KTS)=0.
          CRV(KTS)=0.
          RSTU(KTS)=-(U_K(KTS)-U_K(KTS+1))*VVUPU*WGT1                   &
       &               -RSTU(KTS+1)*CFU+U_K(KTS)
          RSTV(KTS)=-(V_K(KTS)-V_K(KTS+1))*VVUPV*WGT1                   &
       &               -RSTV(KTS+1)*CFV+V_K(KTS)
          UN(KTS)=RSTU(KTS)/CMU
          VN(KTS)=RSTV(KTS)/CMV
          VAD_TEND_U(I,J,KTS)=UN(KTS)-U_K(KTS)
          VAD_TEND_V(I,J,KTS)=VN(KTS)-V_K(KTS)
!
          DO K=KTS+1,KTE
            UN(K)=(-CRU(K)*UN(K-1)+RSTU(K))*RCMU(K)
            VN(K)=(-CRV(K)*VN(K-1)+RSTV(K))*RCMV(K)
            VAD_TEND_U(I,J,K)=UN(K)-U_K(K)
            VAD_TEND_V(I,J,K)=VN(K)-V_K(K)
          ENDDO
!
!-----------------------------------------------------------------------
!***  The following section is only for checking the implicit solution
!***  using back-substitution.  Remove this section otherwise.
!-----------------------------------------------------------------------
!
!       if(ntsd<=10.or.ntsd>=6000)then
!       IF(I==ITEST.AND.J==JTEST)THEN
!!
!         PDOPU=(PDSLO(I+IVW(J),J)+PDSLO(I+IVE(J),J))*0.5
!         PDOPV=(PDSLO(I,J-1)+PDSLO(I,J+1))*0.5
!         PVVLOU=(PETDT(I+IVW(J),J,KTE-1)                               &
!    &           +PETDT(I+IVE(J),J,KTE-1))*DTE
!         PVVLOV=(PETDT(I,J-1,KTE-1)                                    &
!    &           +PETDT(I,J+1,KTE-1))*DTE
!         VVLOU=PVVLOU/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOPU)
!         VVLOV=PVVLOV/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOPV)
!         TULO=VVLOU*(U(I,J,KTE-1)-U(I,J,KTE)+UN(KTE-1)-UN(KTE))
!         TVLO=VVLOV*(V(I,J,KTE-1)-V(I,J,KTE)+VN(KTE-1)-VN(KTE))
!         ADUP=TULO+UN(KTE)-U(I,J,KTE)
!         ADVP=TVLO+VN(KTE)-V(I,J,KTE)
!         WRITE(0,*)' NTSD=',NTSD,' I=',I,' J=',J,' K=',KTE             &
!    &,             ' ADUP=',ADUP,' ADVP=',ADVP
!         WRITE(0,*)' U=',U(I,J,KTE),' UN=',UN(KTE)                     &
!    &,             ' VAD_TEND_U=',VAD_TEND_U(I,KTE)                    &
!    &,             ' V=',V(I,J,KTE),' VN=',VN(KTE)                     &
!    &,             ' VAD_TEND_V=',VAD_TEND_V(I,KTE)
!         WRITE(0,*)' '
!!
!         DO K=KTE-1,KTS+1,-1
!           RDPU=1./(DETA1_PDTOP(K)+DETA2(K)*PDOPU)
!           RDPV=1./(DETA1_PDTOP(K)+DETA2(K)*PDOPV)
!           PVVUPU=PVVLOU
!           PVVUPV=PVVLOV
!           PVVLOU=(PETDT(I+IVW(J),J,K-1)                               &
!    &            +PETDT(I+IVE(J),J,K-1))*DTE
!           PVVLOV=(PETDT(I,J-1,K-1)+PETDT(I,J+1,K-1))*DTE
!           VVUPU=PVVUPU*RDPU
!           VVUPV=PVVUPV*RDPV
!           VVLOU=PVVLOU*RDPU
!           VVLOV=PVVLOV*RDPV
!           TUUP=VVUPU*(U(I,J,K)-U(I,J,K+1)+UN(K)-UN(K+1))
!           TVUP=VVUPV*(V(I,J,K)-V(I,J,K+1)+VN(K)-VN(K+1))
!           TULO=VVLOU*(U(I,J,K-1)-U(I,J,K)+UN(K-1)-UN(K))
!           TVLO=VVLOV*(V(I,J,K-1)-V(I,J,K)+VN(K-1)-VN(K))
!           ADUP=TUUP+TULO+UN(K)-U(I,J,K)
!           ADVP=TVUP+TVLO+VN(K)-V(I,J,K)
!           WRITE(0,*)' NTSD=',NTSD,' I=',ITEST,' J=',JTEST,' K=',K     &
!    &,               ' ADUP=',ADUP,' ADVP=',ADVP
!           WRITE(0,*)' U=',U(I,J,K),' UN=',UN(K)                       &
!    &,               ' VAD_TEND_U=',VAD_TEND_U(I,K)                    &
!    &,               ' V=',V(I,J,K),' VN=',VN(K)                       &
!    &,               ' VAD_TEND_V=',VAD_TEND_V(I,K)
!           WRITE(0,*)' '
!         ENDDO
!!
!         PVVUPU=PVVLOU
!         PVVUPV=PVVLOV
!         VVUPU=PVVUPU/(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOPU)
!         VVUPV=PVVUPV/(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOPV)
!         TUUP=VVUPU*(U(I,J,KTS)-U(I,J,KTS+1)+UN(KTS)-UN(KTS+1))
!         TVUP=VVUPV*(V(I,J,KTS)-V(I,J,KTS+1)+VN(KTS)-VN(KTS+1))
!         ADUP=TUUP+UN(KTS)-U(I,J,KTS)
!         ADVP=TVUP+VN(KTS)-V(I,J,KTS)
!         WRITE(0,*)' NTSD=',NTSD,' I=',ITEST,' J=',JTEST,' K=',KTS     &
!    &,             ' ADUP=',ADUP,' ADVP=',ADVP
!         WRITE(0,*)' U=',U(I,J,KTS),' UN=',UN(KTS)                     &
!    &,             ' VAD_TEND_U=',VAD_TEND_U(I,KTS)                    &
!    &,             ' V=',V(I,J,KTS),' VN=',VN(KTS)                     &
!    &,             ' VAD_TEND_V=',VAD_TEND_V(I,KTS)
!         WRITE(0,*)' '
!       ENDIF
!     endif
!
!-----------------------------------------------------------------------
!***  End of check.
!-----------------------------------------------------------------------
!
        ENDDO iloop_for_uv
!
!-----------------------------------------------------------------------
!
      ENDDO main_vertical
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!***  COMPUTE HORIZONTAL ADVECTION TENDENCIES.
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!$omp parallel do                                                       &
!$omp& private(adpdx,adpdy,adt,adu,adv,array0,array1,array2,array3      &
!$omp&        ,array3_x,dpde,f0,f1,f2,f3,fewp,fnep,fnsp,fpp,fsep,hm     &
!$omp&        ,i,ifp,ifq,ii,ipq,isp,ispa,isq,isqa,iup_adh_j,j,k         &
!$omp&        ,knti_adh,n_iupadh_j,n_iupadv_j,n_iuph_j,pp,qp            &
!$omp&        ,rdpd,rdpdx,rdpdy,tew,tne,tns,tse,tst,tta,ttb             &
!$omp&        ,uew,udy,une,uned,uns,use,used,ust                        &
!$omp&        ,vdx,vew,vm,vne,vns,vse,vst)
!-----------------------------------------------------------------------
!
      main_horizontal: DO K=KTS,KTE
!
!-----------------------------------------------------------------------
!
        DO J=MYJS_P4,MYJE_P4
        DO I=MYIS_P4,MYIE_P4
          DPDE(I,J)=DETA1_PDTOP(K)+DETA2(K)*PDSLO(I,J)
          RDPD(I,J)=1./DPDE(I,J)
          TST(I,J)=T(I,J,K)*FFC+TOLD(I,J,K)*FBC
          UST(I,J)=U(I,J,K)*FFC+UOLD(I,J,K)*FBC
          VST(I,J)=V(I,J,K)*FFC+VOLD(I,J,K)*FBC
        ENDDO
        ENDDO
!
!-----------------------------------------------------------------------
!***  MASS FLUXES AND MASS POINT ADVECTION COMPONENTS
!***  THE NS AND EW FLUXES IN THE FOLLOWING LOOP ARE ON V POINTS
!***  FOR T.
!-----------------------------------------------------------------------
!
        DO J=MYJS1_P3,MYJE1_P3
        DO I=MYIS_P3,MYIE_P3
!
          ADPDX=DPDE(I+IVW(J),J)+DPDE(I+IVE(J),J)
          ADPDY=DPDE(I,J-1)+DPDE(I,J+1)
          RDPDX(I,J)=1./ADPDX
          RDPDY(I,J)=1./ADPDY
!
          UDY=U(I,J,K)*DY
          VDX=V(I,J,K)*DX(I,J)
!
          FEWP=UDY*ADPDX
          FNSP=VDX*ADPDY
!
          FEW(I,J,K)=FEWP
          FNS(I,J,K)=FNSP
!
          TEW(I,J)=FEWP*(TST(I+IVE(J),J)-TST(I+IVW(J),J))
          TNS(I,J)=FNSP*(TST(I,J+1)-TST(I,J-1))
!
          UNED(I,J)=UDY+VDX
          USED(I,J)=UDY-VDX
!
        ENDDO
        ENDDO
!
!-----------------------------------------------------------------------
!***  DIAGONAL FLUXES AND DIAGONALLY AVERAGED WIND
!***  THE NE AND SE FLUXES ARE ASSOCIATED WITH H POINTS
!***  (ACTUALLY JUST TO THE NE AND SE OF EACH H POINT).
!-----------------------------------------------------------------------
!
        DO J=MYJS1_P2,MYJE2_P2
        DO I=MYIS_P2,MYIE_P2
          FNEP=(UNED(I+IHE(J),J)+UNED(I       ,J+1))                    &
     &        *(DPDE(I       ,J)+DPDE(I+IHE(J),J+1))
          FNE(I,J,K)=FNEP
          TNE(I,J)=FNEP*(TST(I+IHE(J),J+1)-TST(I,J))
        ENDDO
        ENDDO
!
        DO J=MYJS2_P2,MYJE1_P2
        DO I=MYIS_P2,MYIE_P2
          FSEP=(USED(I+IHE(J),J)+USED(I       ,J-1))                    &
     &        *(DPDE(I       ,J)+DPDE(I+IHE(J),J-1))
          FSE(I,J,K)=FSEP
          TSE(I,J)=FSEP*(TST(I+IHE(J),J-1)-TST(I,J))
!
        ENDDO
        ENDDO
!
!-----------------------------------------------------------------------
!***  HORIZONTAL T ADVECTION TENDENCY ADT IS ON H POINTS OF COURSE.
!-----------------------------------------------------------------------
!
        DO J=MYJS5,MYJE5
        DO I=MYIS2,MYIE2
          ADT(I,J)=(TEW(I+IHW(J),J)+TEW(I+IHE(J),J)                     &
     &             +TNS(I,J-1)+TNS(I,J+1)                               &
     &             +TNE(I+IHW(J),J-1)+TNE(I,J)                          &
     &             +TSE(I,J)+TSE(I+IHW(J),J+1))                         &
     &             *RDPD(I,J)*FAD(I,J)
        ENDDO
        ENDDO
!
!
!-----------------------------------------------------------------------
!***  CALCULATION OF MOMENTUM ADVECTION COMPONENTS.
!-----------------------------------------------------------------------
!
        DO J=MYJS4_P1,MYJE4_P1
        DO I=MYIS_P1,MYIE_P1
!
!-----------------------------------------------------------------------
!***  THE NS AND EW FLUXES ARE ON H POINTS FOR U AND V.
!-----------------------------------------------------------------------
!
          UEW(I,J)=(FEW(I+IHW(J),J,K)+FEW(I+IHE(J),J,K))                &
     &            *(UST(I+IHE(J),J)-UST(I+IHW(J),J))
          UNS(I,J)=(FNS(I+IHW(J),J,K)+FNS(I+IHE(J),J,K))                &
     &            *(UST(I,J+1)-UST(I,J-1))
          VEW(I,J)=(FEW(I,J-1,K)+FEW(I,J+1,K))                          &
     &            *(VST(I+IHE(J),J)-VST(I+IHW(J),J))
          VNS(I,J)=(FNS(I,J-1,K)+FNS(I,J+1,K))                          &
     &            *(VST(I,J+1)-VST(I,J-1))
!
!-----------------------------------------------------------------------
!***  THE FOLLOWING NE AND SE FLUXES ARE TIED TO V POINTS AND ARE
!***  LOCATED JUST TO THE NE AND SE OF THE GIVEN I,J.
!-----------------------------------------------------------------------
!
          UNE(I,J)=(FNE(I+IVW(J),J,K)+FNE(I+IVE(J),J,K))                &
     &            *(UST(I+IVE(J),J+1)-UST(I,J))
          USE(I,J)=(FSE(I+IVW(J),J,K)+FSE(I+IVE(J),J,K))                &
     &            *(UST(I+IVE(J),J-1)-UST(I,J))
          VNE(I,J)=(FNE(I,J-1,K)+FNE(I,J+1,K))                          &
     &            *(VST(I+IVE(J),J+1)-VST(I,J))
          VSE(I,J)=(FSE(I,J-1,K)+FSE(I,J+1,K))                          &
     &            *(VST(I+IVE(J),J-1)-VST(I,J))
!
!-----------------------------------------------------------------------
!
        ENDDO
        ENDDO
!
!-----------------------------------------------------------------------
!***  COMPUTE THE ADVECTION TENDENCIES FOR U AND V.
!***  THE AD ARRAYS ARE ON THE VELOCITY POINTS.
!-----------------------------------------------------------------------
!
        DO J=MYJS5,MYJE5
        DO I=MYIS2,MYIE2
          ADU(I,J)=(UEW(I+IVW(J),J)+UEW(I+IVE(J),J)                     &
     &             +UNS(I,J-1)+UNS(I,J+1)                               &
     &             +UNE(I+IVW(J),J-1)+UNE(I,J)                          &
     &             +USE(I,J)+USE(I+IVW(J),J+1))                         &
     &             *RDPDX(I,J)*FAD(I+IVW(J),J)
!
          ADV(I,J)=(VEW(I+IVW(J),J)+VEW(I+IVE(J),J)                     &
     &             +VNS(I,J-1)+VNS(I,J+1)                               &
     &             +VNE(I+IVW(J),J-1)+VNE(I,J)                          &
     &             +VSE(I,J)+VSE(I+IVW(J),J+1))                         &
     &             *RDPDY(I,J)*FAD(I+IVW(J),J)
        ENDDO
        ENDDO
!
!-----------------------------------------------------------------------
!
!***  END OF JANJIC HORIZONTAL ADVECTION
!
!-----------------------------------------------------------------------
!
!***  UPSTREAM ADVECTION OF T
!
!-----------------------------------------------------------------------
!
        upstream: IF(UPSTRM)THEN
!
!-----------------------------------------------------------------------
!***
!***  COMPUTE UPSTREAM COMPUTATIONS ON THIS TASK'S ROWS.
!***
!-----------------------------------------------------------------------
!
          jloop_upstream: DO J=MYJS2,MYJE2
!
            N_IUPH_J=N_IUP_H(J)   ! See explanation in START_DOMAIN_NMM
            DO II=0,N_IUPH_J-1
!
              I=IUP_H(IMS+II,J)
              TTA=EMT_LOC(J)*(UST(I,J-1)+UST(I+IHW(J),J)                &
     &                       +UST(I+IHE(J),J)+UST(I,J+1))
              TTB=ENT       *(VST(I,J-1)+VST(I+IHW(J),J)                &
     &                       +VST(I+IHE(J),J)+VST(I,J+1))
              PP=-TTA-TTB
              QP= TTA-TTB
!
              IF(PP<0.)THEN
                ISPA(I,J)=-1
              ELSE
                ISPA(I,J)= 1
              ENDIF
!
              IF(QP<0.)THEN
                ISQA(I,J)=-1
              ELSE
                ISQA(I,J)= 1
              ENDIF
!
              PP=ABS(PP)
              QP=ABS(QP)
              ARRAY3_X=PP*QP
              ARRAY0(I,J)=ARRAY3_X-PP-QP
              ARRAY1(I,J)=PP-ARRAY3_X
              ARRAY2(I,J)=QP-ARRAY3_X
              ARRAY3(I,J)=ARRAY3_X
            ENDDO
!
!-----------------------------------------------------------------------
!
            N_IUPADH_J=N_IUP_ADH(J)
            KNTI_ADH=1
            IUP_ADH_J=IUP_ADH(IMS,J)
!
            iloop_T: DO II=0,N_IUPH_J-1
!
              I=IUP_H(IMS+II,J)
!
              ISP=ISPA(I,J)
              ISQ=ISQA(I,J)
              IFP=(ISP-1)/2
              IFQ=(-ISQ-1)/2
              IPQ=(ISP-ISQ)/2
!
!-----------------------------------------------------------------------
!
              IF(I==IUP_ADH_J)THEN  ! Upstream advection T tendencies
!
                ISP=ISPA(I,J)
                ISQ=ISQA(I,J)
                IFP=(ISP-1)/2
                IFQ=(-ISQ-1)/2
                IPQ=(ISP-ISQ)/2
!
                F0=ARRAY0(I,J)
                F1=ARRAY1(I,J)
                F2=ARRAY2(I,J)
                F3=ARRAY3(I,J)
!
                ADT(I,J)=F0*T(I,J,K)                                    &
     &                  +F1*T(I+IHE(J)+IFP,J+ISP,K)                     &
     &                  +F2*T(I+IHE(J)+IFQ,J+ISQ,K)                     &
                        +F3*T(I+IPQ,J+ISP+ISQ,K)
!
!-----------------------------------------------------------------------
!
                IF(KNTI_ADH<N_IUPADH_J)THEN
                  IUP_ADH_J=IUP_ADH(IMS+KNTI_ADH,J)
                  KNTI_ADH=KNTI_ADH+1
                ENDIF
!
              ENDIF  ! End of upstream advection T tendency IF block
!
            ENDDO iloop_T
!
!-----------------------------------------------------------------------
!
!***  UPSTREAM ADVECTION OF VELOCITY COMPONENTS
!
!-----------------------------------------------------------------------
!
            N_IUPADV_J=N_IUP_ADV(J)
!
            DO II=0,N_IUPADV_J-1
              I=IUP_ADV(IMS+II,J)
!
              TTA=EM_LOC(J)*UST(I,J)
              TTB=EN       *VST(I,J)
              PP=-TTA-TTB
              QP=TTA-TTB
!
              IF(PP<0.)THEN
                ISP=-1
              ELSE
                ISP= 1
              ENDIF
!
              IF(QP<0.)THEN
                ISQ=-1
              ELSE
                ISQ= 1
              ENDIF
!
              IFP=(ISP-1)/2
              IFQ=(-ISQ-1)/2
              IPQ=(ISP-ISQ)/2
              PP=ABS(PP)
              QP=ABS(QP)
              F3=PP*QP
              F0=F3-PP-QP
              F1=PP-F3
              F2=QP-F3
!
              ADU(I,J)=F0*U(I,J,K)                                      &
     &                +F1*U(I+IVE(J)+IFP,J+ISP,K)                       &
     &                +F2*U(I+IVE(J)+IFQ,J+ISQ,K)                       &
     &                +F3*U(I+IPQ,J+ISP+ISQ,K)
!
              ADV(I,J)=F0*V(I,J,K)                                      &
     &                +F1*V(I+IVE(J)+IFP,J+ISP,K)                       &
     &                +F2*V(I+IVE(J)+IFQ,J+ISQ,K)                       &
     &                +F3*V(I+IPQ,J+ISP+ISQ,K)
!
            ENDDO
!
          ENDDO jloop_upstream
!
!-----------------------------------------------------------------------
!
        ENDIF upstream
!
!-----------------------------------------------------------------------
!
!***  END OF HORIZONTAL ADVECTION
!
!-----------------------------------------------------------------------
!
!***  NOW SUM THE VERTICAL AND HORIZONTAL TENDENCIES,
!***  CURVATURE AND CORIOLIS TERMS.
!
!-----------------------------------------------------------------------
!
        DO J=MYJS2,MYJE2
        DO I=MYIS1,MYIE1
          HM=HBM2(I,J)
          VM=VBM2(I,J)
          ADT(I,J)=(VAD_TEND_T(I,J,K)+2.*ADT(I,J))*HM
!
          FPP=CURV(I,J)*2.*UST(I,J)+F(I,J)*2.
          ADU(I,J)=(VAD_TEND_U(I,J,K)+2.*ADU(I,J)+VST(I,J)*FPP)*VM
          ADV(I,J)=(VAD_TEND_V(I,J,K)+2.*ADV(I,J)-UST(I,J)*FPP)*VM
        ENDDO
        ENDDO
!
!-----------------------------------------------------------------------
!***  SAVE THE OLD VALUES FOR TIMESTEPPING
!-----------------------------------------------------------------------
!
        DO J=MYJS_P4,MYJE_P4
        DO I=MYIS_P4,MYIE_P4
          TOLD(I,J,K)=T(I,J,K)
          UOLD(I,J,K)=U(I,J,K)
          VOLD(I,J,K)=V(I,J,K)
        ENDDO
        ENDDO
!
!-----------------------------------------------------------------------
!***  FINALLY UPDATE THE PROGNOSTIC VARIABLES
!-----------------------------------------------------------------------
!
        DO J=MYJS2,MYJE2
        DO I=MYIS1,MYIE1
          T(I,J,K)=ADT(I,J)+T(I,J,K)
          U(I,J,K)=ADU(I,J)+U(I,J,K)
          V(I,J,K)=ADV(I,J)+V(I,J,K)
        ENDDO
        ENDDO
!
!-----------------------------------------------------------------------
!
      ENDDO main_horizontal
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE ADVE
!
!-----------------------------------------------------------------------
!
!***********************************************************************
      SUBROUTINE VAD2(NTSD,DT,IDTAD,DX,DY                               &
     &               ,AETA1,AETA2,DETA1,DETA2,PDSL,PDTOP,HBM2           &
     &               ,Q,Q2,CWM,PETDT                                    &
     &               ,N_IUP_H,N_IUP_V                                   &
     &               ,N_IUP_ADH,N_IUP_ADV                               &
     &               ,IUP_H,IUP_V,IUP_ADH,IUP_ADV                       &
     &               ,IHE,IHW,IVE,IVW                                   &
     &               ,IDS,IDE,JDS,JDE,KDS,KDE                           &
     &               ,IMS,IME,JMS,JME,KMS,KME                           &
     &               ,ITS,ITE,JTS,JTE,KTS,KTE)
!***********************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    VAD2        VERTICAL ADVECTION OF H2O SUBSTANCE AND TKE
!   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 96-07-19
!
! ABSTRACT:
!     VAD2 CALCULATES THE CONTRIBUTION OF THE VERTICAL ADVECTION
!     TO THE TENDENCIES OF WATER SUBSTANCE AND TKE AND THEN UPDATES
!     THOSE VARIABLES.  AN ANTI-FILTERING TECHNIQUE IS USED.
!
! PROGRAM HISTORY LOG:
!   96-07-19  JANJIC   - ORIGINATOR
!   98-11-02  BLACK    - MODIFIED FOR DISTRIBUTED MEMORY
!   99-03-17  TUCCILLO - INCORPORATED MPI_ALLREDUCE FOR GLOBAL SUM
!   02-02-06  BLACK    - CONVERTED TO WRF FORMAT
!   02-09-06  WOLFE    - MORE CONVERSION TO GLOBAL INDEXING
!   04-11-23  BLACK    - THREADED
!   05-12-14  BLACK    - CONVERTED FROM IKJ TO IJK
!   07-08-14  janjic   - bc & no conservation in the advection step
!
! USAGE: CALL VAD2 FROM SUBROUTINE SOLVE_NMM
!   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
      INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: N_IUP_H,N_IUP_V          &
     &                                        ,N_IUP_ADH,N_IUP_ADV
      INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IUP_H,IUP_V      &
     &                                                ,IUP_ADH,IUP_ADV
!
      INTEGER,INTENT(IN) :: IDTAD,NTSD
!
      REAL,INTENT(IN) :: DT,DY,PDTOP
!
      REAL,DIMENSION(KMS:KME),INTENT(IN) :: AETA1,AETA2,DETA1,DETA2
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: DX,HBM2,PDSL
!
      REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(IN) :: PETDT
!
      REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(INOUT) :: CWM,Q,Q2
!
!***  LOCAL VARIABLES
!----------------------------------------------------------------------
!
      REAL,PARAMETER :: FF1=0.500
!
      LOGICAL,SAVE :: TRADITIONAL=.TRUE.
!
      INTEGER :: I,IRECV,J,JFP,JFQ,K,LAP,LLAP
!
      INTEGER,DIMENSION(KTS:KTE) :: LA
!
      REAL*8 :: ADDT,AFRP,D2PQE,D2PQQ,D2PQW,DEP,DETAP,DQP               &
     &       ,DWP,E00,E4P,EP,EP0,HADDT,HBM2IJ                           &
     &       ,Q00,Q4P,QP,QP0                                            &
     &       ,rdpdn,rdpup,sfacek,sfacqk,sfacwk,RFC,RR                   &
     &       ,SUMNE,SUMNQ,SUMNW,SUMPE,SUMPQ,SUMPW                       &
     &       ,W00,W4P,WP,WP0
!
      REAL,DIMENSION(KTS:KTE) :: AFR,DEL,DQL,DWL,E3,E4,PETDTK           &
     &                          ,RFACE,RFACQ,RFACW,Q3,Q4,W3,W4
!
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------------------------------------------------------------
!
      ADDT=REAL(IDTAD)*DT
!
!-----------------------------------------------------------------------
!$omp parallel do                                                       &
!$omp& private(afr,afrp,bot,d2pqe,d2pqq,d2pqw,del,dep,detap,dpdn,dpup   &
!$omp&        ,dql,dqp,dwl,dwp,e00,e3,e4,e4p,ep,ep0,haddt,i,j,k         &
!$omp&        ,la,lap,llap,petdtk,q00,q3,q4,q4p,qp,qp0,rfacek,rfacqk    &
!$omp&        ,rfacwk,rfc,rr,sumne,sumnq,sumnw,sumpe,sumpq,sumpw,top    &
!$omp&        ,w00,w3,w4,w4p,wp,wp0)
!-----------------------------------------------------------------------
!
      main_integration: DO J=MYJS2,MYJE2
!
!-----------------------------------------------------------------------
!
        main_iloop: DO I=MYIS1_P1,MYIE1_P1
!
!-----------------------------------------------------------------------
!
          E3(KTE)=Q2(I,J,KTE)*0.5
!
          DO K=KTE-1,KTS,-1
            E3(K)=MAX((Q2(I,J,K+1)+Q2(I,J,K))*0.5,EPSQ2)
          ENDDO
!
          DO K=KTS,KTE
            Q3(K)=MAX(Q(I,J,K),EPSQ)
            W3(K)=MAX(CWM(I,J,K),CLIMIT)
            E4(K)=E3(K)
            Q4(K)=Q3(K)
            W4(K)=W3(K)
          ENDDO
!
          IF(TRADITIONAL)THEN
            PETDTK(KTE)=PETDT(I,J,KTE-1)*0.5
!
            DO K=KTE-1,KTS+1,-1
              PETDTK(K)=(PETDT(I,J,K)+PETDT(I,J,K-1))*0.5
            ENDDO
!
            PETDTK(KTS)=PETDT(I,J,KTS)*0.5
!
          ELSE
!
!-----------------------------------------------------------------------
!***    PERFORM HORIZONTAL AVERAGING OF VERTICAL VELOCITY
!-----------------------------------------------------------------------
!
            PETDTK(KTE)=(PETDT(I+IHW(J-1),J-1,KTE-1)                    &
     &                  +PETDT(I+IHE(J-1),J-1,KTE-1)                    &
     &                  +PETDT(I+IHW(J+1),J+1,KTE-1)                    &
     &                  +PETDT(I+IHE(J+1),J+1,KTE-1)                    &
     &                  +PETDT(I,J,KTE-1)*4.        )*0.0625
!
            DO K=KTE-1,KTS+1,-1
              PETDTK(K)=(PETDT(I+IHW(J-1),J-1,K-1)                      &
                        +PETDT(I+IHE(J-1),J-1,K-1)                      &
     &                  +PETDT(I+IHW(J+1),J+1,K-1)                      &
     &                  +PETDT(I+IHE(J+1),J+1,K-1)                      &
     &                  +PETDT(I+IHW(J-1),J-1,K  )                      &
     &                  +PETDT(I+IHE(J-1),J-1,K  )                      &
     &                  +PETDT(I+IHW(J+1),J+1,K  )                      &
     &                  +PETDT(I+IHE(J+1),J+1,K  )                      &
     &                  +(PETDT(I,J,K-1)+PETDT(I,J,K))*4.               &
     &                                                   )*0.0625
            ENDDO
!
            PETDTK(KTS)=(PETDT(I+IHW(J-1),J-1,KTS)                      &
     &                  +PETDT(I+IHE(J-1),J-1,KTS)                      &
     &                  +PETDT(I+IHW(J+1),J+1,KTS)                      &
     &                  +PETDT(I+IHE(J+1),J+1,KTS)                      &
     &                  +PETDT(I,J,KTS)*4.        )*0.0625

          ENDIF
!
!-----------------------------------------------------------------------
!
          HADDT=-ADDT*HBM2(I,J)
!
          DO K=KTE,KTS,-1
            RR=PETDTK(K)*HADDT
!
            IF(RR<0.)THEN
              LAP=1
            ELSE
              LAP=-1
            ENDIF
!
            LA(K)=LAP
            LLAP=K+LAP
!
            if(llap.gt.kts-1.and.llap.lt.kte+1) then ! internal and outflow pts.
              rr=abs(rr &
     &          /((aeta1(llap)-aeta1(k))*pdtop &
     &           +(aeta2(llap)-aeta2(k))*pdsl(i,j)))
              if(rr.gt.0.999) rr=0.999   
!
              AFR(K)=(((FF4*RR+FF3)*RR+FF2)*RR+FF1)*RR
              dql(k)=(q3(llap)-q3(k))*rr
              dwl(k)=(w3(llap)-w3(k))*rr
              del(k)=(e3(llap)-e3(k))*rr
            elseif(llap.eq.kts-1) then
!
!chem              rr=abs(rr &
!chem                /((1.-aeta2(kts))*pdsl(i,j)))
!chem              afr(kts)=0.
!chem              dql(kts)=(epsq  -q3(kts))*rr
!chem              dwl(kts)=(climit-w3(kts))*rr
!chem              del(kts)=(epsq2 -e3(kts))*rr
!
              rr=0.
              afr(kts)=0.
              dql(kts)=0.
              dwl(kts)=0.
              del(kts)=0.
            else
              rr=abs(rr &
                /(aeta1(kte)*pdtop))
              afr(kte)=0.
              dql(kte)=(epsq  -q3(kte))*rr
              dwl(kte)=(climit-w3(kte))*rr
              del(kte)=(epsq2 -e3(kte))*rr
            endif
          ENDDO
!
!-----------------------------------------------------------------------
!
          DO K=KTS,KTE
            Q4(K)=Q3(K)+DQL(K)
            W4(K)=W3(K)+DWL(K)
            E4(K)=E3(K)+DEL(K)
          ENDDO
!
!-----------------------------------------------------------------------
!***  ANTI-FILTERING STEP
!-----------------------------------------------------------------------
!
          SUMPQ=0.
          SUMNQ=0.
          SUMPW=0.
          SUMNW=0.
          SUMPE=0.
          SUMNE=0.
!
!***  ANTI-FILTERING LIMITERS
!
          antifilter: DO K=KTE-1,KTS+1,-1
!
            DETAP=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J)
!
            DQL(K)=0.
            DWL(K)=0.
            DEL(K)=0.
!
            Q4P=Q4(K)
            W4P=W4(K)
            E4P=E4(K)
!
            LAP=LA(K)
!
            if(lap.ne.0)then
              rdpdn=1./((aeta1(k+lap)-aeta1(k))*pdtop &
     &                 +(aeta2(k+lap)-aeta2(k))*pdsl(i,j))
              rdpup=1./((aeta1(k)-aeta1(k-lap))*pdtop &
     &                 +(aeta2(k)-aeta2(k-lap))*pdsl(i,j))
!
              afrp=afr(k)*detap
!
              d2pqq=((q4(k+lap)-q4p)*rdpdn &
     &              -(q4p-q4(k-lap))*rdpup)*afrp
              d2pqw=((w4(k+lap)-w4p)*rdpdn &
     &              -(w4p-w4(k-lap))*rdpup)*afrp
              d2pqe=((e4(k+lap)-e4p)*rdpdn &
     &              -(e4p-e4(k-lap))*rdpup)*afrp
            ELSE
              D2PQQ=0.
              D2PQW=0.
              D2PQE=0.
            ENDIF
!
            QP=Q4P-D2PQQ
            WP=W4P-D2PQW
            EP=E4P-D2PQE
!
            Q00=Q3(K)
            QP0=Q3(K+LAP)
!
            W00=W3(K)
            WP0=W3(K+LAP)
!
            E00=E3(K)
            EP0=E3(K+LAP)
!
            IF(LAP/=0)THEN
              QP=MAX(QP,MIN(Q00,QP0))
              QP=MIN(QP,MAX(Q00,QP0))
              WP=MAX(WP,MIN(W00,WP0))
              WP=MIN(WP,MAX(W00,WP0))
              EP=MAX(EP,MIN(E00,EP0))
              EP=MIN(EP,MAX(E00,EP0))
            ENDIF
!
            dqp=qp-q4p
            dwp=wp-w4p
            dep=ep-e4p
!
            DQL(K)=DQP
            DWL(K)=DWP
            DEL(K)=DEP
!
            DQP=DQP*DETAP
            DWP=DWP*DETAP
            DEP=DEP*DETAP
!
            IF(DQP>0.)THEN
              SUMPQ=SUMPQ+DQP
            ELSE
              SUMNQ=SUMNQ+DQP
            ENDIF
!
            IF(DWP>0.)THEN
              SUMPW=SUMPW+DWP
            ELSE
              SUMNW=SUMNW+DWP
            ENDIF
!
            IF(DEP>0.)THEN
              SUMPE=SUMPE+DEP
            ELSE
              SUMNE=SUMNE+DEP
            ENDIF
!
          ENDDO antifilter
!
!-----------------------------------------------------------------------
!
          DQL(KTS)=0.
          DWL(KTS)=0.
          DEL(KTS)=0.
!
          DQL(KTE)=0.
          DWL(KTE)=0.
          DEL(KTE)=0.
!
!-----------------------------------------------------------------------
!***  FIRST MOMENT CONSERVING FACTOR
!-----------------------------------------------------------------------
!
          if(sumpq*(-sumnq).gt.1.e-9) then
            sfacqk=-sumnq/sumpq
          else
            sfacqk=0.
          endif
!
          if(sumpw*(-sumnw).gt.1.e-9) then
            sfacwk=-sumnw/sumpw
          else
            sfacwk=0.
          endif
!
          if(sumpe*(-sumne).gt.1.e-9) then
            sfacek=-sumne/sumpe
          else
            sfacek=0.
          endif
!
!-----------------------------------------------------------------------
!***  IMPOSE CONSERVATION ON ANTI-FILTERING
!-----------------------------------------------------------------------
!
          DO K=KTE,KTS,-1
!
            dqp=dql(k)
            if(sfacqk.gt.0.) then
              if(sfacqk.ge.1.) then
                if(dqp.lt.0.) dqp=dqp/sfacqk
              else
                if(dqp.gt.0.) dqp=dqp*sfacqk
              endif
            else
              dqp=0.
            endif
            q  (i,j,k)=q4(k)+dqp
!
            dwp=dwl(k)
            if(sfacwk.gt.0.) then
              if(sfacwk.ge.1.) then
                if(dwp.lt.0.) dwp=dwp/sfacwk
              else
                if(dwp.gt.0.) dwp=dwp*sfacwk
              endif
            else
              dwp=0.
            endif
            cwm(i,j,k)=w4(k)+dwp
!
            dep=del(k)
            if(sfacek.gt.0.) then
              if(sfacek.ge.1.) then
                if(dep.lt.0.) dep=dep/sfacek
              else
                if(dep.gt.0.) dep=dep*sfacek
              endif
            else
              dep=0.
            endif
            e3 (    k)=e4(k)+dep
!
          ENDDO
!-----------------------------------------------------------------------
          HBM2IJ=HBM2(I,J)
          Q2(I,J,KTE)=MAX(E3(KTE)+E3(KTE)-EPSQ2,EPSQ2)*HBM2IJ           &
       &             +Q2(I,J,KTE)*(1.-HBM2IJ)
          DO K=KTE-1,KTS+1,-1
            Q2(I,J,K)=MAX(E3(K)+E3(K)-Q2(I,J,K+1),EPSQ2)*HBM2IJ         &
       &             +Q2(I,J,K)*(1.-HBM2IJ)
          ENDDO
!------------------------------------------…