C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
                             SUBROUTINE VTADV
C     ******************************************************************
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    VTADV       VERTICAL ADVECTION
C   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 93-11-17
C     
C ABSTRACT:
C     VTADV CALCULATES THE CONTRIBUTION OF THE VERTICAL ADVECTION
C     TO THE TENDENCIES OF TEMPERATURE, SPECIFIC HUMIDITY, WIND
C     COMPONENTS, AND TURBULENT KINETIC ENERGY AND THEN UPDATES THOSE
C     VARIABLES.  FOR ALL VARIABLES EXCEPT SPECIFIC HUMIDITY A
C     SIMPLE CENTERED DIFFERENCE SCHEME IN SPACE IS USED IN
C     CONJUNCTION WITH THE PURE EULER-BACKWARD TIME SCHEME.  
C     A PIECEWISE LINEAR SCHEME IS USED TO CALCULATE THE VERTICAL
C     ADVECTION OF SPECIFIC HUMIDITY SO THAT NO FALSE MAXIMA OR
C     MINIMA ARE PRODUCED.
C     
C     PIECEWISE LINEAR SCHEME (MESINGER AND JOVIC, NCEP OFFICE NOTE
C     #439) HERE USED FOR ALL VARIABLES, AVOIDING FALSE ADVECTION
C     FROM BELOW GROUND, AND AT THE SAME TIME MAKING THE ETA
C     VERY NEARLY A FINITE VOLUME MODEL
C
C PROGRAM HISTORY LOG:
C   87-06-??  JANJIC     - ORIGINATOR
C   90-??-??  MESINGER   - INSERTED PIECEWISE LINEAR SCHEME FOR
C                          SPECIFIC HUMIDITY
C   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
C   95-11-20  ABELES     - PARALLEL OPTIMIZATION
C   96-03-29  BLACK      - ADDED EXTERNAL EDGE; REMOVED SCRCH COMMON
C   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
C   00-02-04  BLACK      - ADDED CLOUD WATER/ICE
C   01-12-11  BLACK      - SMOOTHING FOR CFL VIOLATION
C
C   This is the code as e-mailed April 10, 07, but with bugs corrected:
C      1)  FFD1, twice, 1 is removed,
C      2)  F4D, in loop 1320, the second pair changed into F4Q;
C      3)  ..TEND1.. 2nd time in loop 1330 repl. by ..TEND2..
C   and also:  
C      a) dead pieces of the code - abandoned finite difference vertical
C   advection of T and u,v (twice "GO TO ..." segments) - deleted;
C      b) Switching to calls every time step commented out, as this
C   takes non-negligible time, and in the Eta horizontal advection
C   is called every second time step also (consistency)
C     
C USAGE: CALL VTADV FROM MAIN PROGRAM EBU
C   INPUT ARGUMENT LIST:
C       NONE     
C  
C   OUTPUT ARGUMENT LIST: 
C     NONE
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C  
C     UNIQUE: NONE
C  
C     LIBRARY: NONE
C  
C   COMMON BLOCKS: CTLBLK
C                  MASKS
C                  DYNAM
C                  VRBLS
C                  CONTIN
C                  PVRBLS
C                  CLDWTR
C                  INDX
C   
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE : IBM SP
C$$$  
C***********************************************************************
                             P A R A M E T E R
     & (EDQMX=2.E-5,EDQMN=-2.E-5,EPSQ=1.E-12,EPSQ2=0.2,KSMUD=0)
C
                             P A R A M E T E R
     & (CFL_MAX=0.97)
C-----------------------------------------------------------------------
      INCLUDE "parmeta"
      INCLUDE "mpp.h"
#include "sp.h"
C-----------------------------------------------------------------------
                             P A R A M E T E R
     & (IMJM=IM*JM-JM/2,JAM=6+2*(JM-10)
     &, LM1=LM-1,LM2=LM-2,LP1=LM+1)
C-----------------------------------------------------------------------
                             L O G I C A L
     & RUN,FIRST,RESTRT,SIGMA,NOSLA,NOSLAW
C----------------------------------------------------------------------
      INCLUDE "CTLBLK.comm"
C-----------------------------------------------------------------------
      INCLUDE "MASKS.comm"
C-----------------------------------------------------------------------
      INCLUDE "DYNAM.comm"
C-----------------------------------------------------------------------
      INCLUDE "VRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "CONTIN.comm"
C-----------------------------------------------------------------------
      INCLUDE "PVRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "CLDWTR.comm"
C-----------------------------------------------------------------------
      INCLUDE "INDX.comm"
C-----------------------------------------------------------------------
      INCLUDE "NHYDRO.comm"
C-----------------------------------------------------------------------
                             D I M E N S I O N
     & WFA   ( LM1),WFB   ( LM1)
C
                             D I M E N S I O N
     & ETADTL(IDIM1:IDIM2,JDIM1:JDIM2)
     &,                                TQ2B  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,DQTI  (IDIM1:IDIM2,JDIM1:JDIM2),DQBI  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,RPDX  (IDIM1:IDIM2,JDIM1:JDIM2),RPDY  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,QDEDB (IDIM1:IDIM2,JDIM1:JDIM2),QDEUB (IDIM1:IDIM2,JDIM1:JDIM2)
     &,                                EDBD  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,                                DQDEB (IDIM1:IDIM2,JDIM1:JDIM2)
C
     &,DUTI  (IDIM1:IDIM2,JDIM1:JDIM2),DUBI  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,DVTI  (IDIM1:IDIM2,JDIM1:JDIM2),DVBI  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,UDEDB (IDIM1:IDIM2,JDIM1:JDIM2),UDEUB (IDIM1:IDIM2,JDIM1:JDIM2)
     &,VDEDB (IDIM1:IDIM2,JDIM1:JDIM2),VDEUB (IDIM1:IDIM2,JDIM1:JDIM2)
     &,EDBDU (IDIM1:IDIM2,JDIM1:JDIM2),EDBDV (IDIM1:IDIM2,JDIM1:JDIM2)
     &,DUDEB (IDIM1:IDIM2,JDIM1:JDIM2),DVDEB (IDIM1:IDIM2,JDIM1:JDIM2)

                             D I M E N S I O N
     & FNE (IDIM1:IDIM2,JDIM1:JDIM2),FSE (IDIM1:IDIM2,JDIM1:JDIM2)
C
                             D I M E N S I O N
     & SAM   (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,QBI   (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,Q2ST  (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,ARRAY1(IDIM1:IDIM2,JDIM1:JDIM2,LM1)
     &,ARRAY2(IDIM1:IDIM2,JDIM1:JDIM2,LM1)
C
     &,SAMU   (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,UBI    (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,SAMV   (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,VBI    (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,ARRAYU1(IDIM1:IDIM2,JDIM1:JDIM2,LM1)
     &,ARRAYU2(IDIM1:IDIM2,JDIM1:JDIM2,LM1)
     &,ARRAYV1(IDIM1:IDIM2,JDIM1:JDIM2,LM1)
     &,ARRAYV2(IDIM1:IDIM2,JDIM1:JDIM2,LM1)
C-----------------------------------------------------------------------
                             R E A L
     &,ALLOCATABLE,DIMENSION(:,:,:) :: S
C-----------------------------------------------------------------------
                             R E A L
     & VAD_TEND1(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,VAD_TEND2(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,VAD_TNDX1(LM),VAD_TNDX2(LM)
!
                             I N T E G E R
     & LBOT_CFL_T(IDIM1:IDIM2,JDIM1:JDIM2)
     &,LTOP_CFL_T(IDIM1:IDIM2,JDIM1:JDIM2)
     &,LBOT_CFL_U(IDIM1:IDIM2,JDIM1:JDIM2)
     &,LTOP_CFL_U(IDIM1:IDIM2,JDIM1:JDIM2)
     &,LBOT_CFL_V(IDIM1:IDIM2,JDIM1:JDIM2)
     &,LTOP_CFL_V(IDIM1:IDIM2,JDIM1:JDIM2)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      DTAD=IDTAD*DT
Cpazi
Cpazi (constants change to call VTADV every time step, commented out)
C      DTAD=DTAD*0.5
C      F4D=F4D*0.5
C      F4Q=F4Q*0.5
C      DO L=1,LM
C      F4Q2(L)=F4Q2(L)*0.5
C      END DO
Cpazi
C-----------------------------------------------------------------------
C--------------DEFINE ADDED UPSTREAM ADVECTION CONSTANTS----------------
C-----------------------------------------------------------------------
                             DO 25 L=1,LM1
      WFA(L)=DETA(L  )/(DETA(L)+DETA(L+1))
      WFB(L)=DETA(L+1)/(DETA(L)+DETA(L+1))
   25 CONTINUE
C--------------NO MOISTURE SLOPE ADJUSTMENT IF NOT WANTED---------------
      NOSLA=.FALSE.
      NOSLAW=.FALSE.
C       IF FALSE, NUMBER OF MOISTURE SLOPE ADJUSTMENT PASSES
      NMSAP=3
      NMSAPW=3
C--------------SMOOTHING VERTICAL VELOCITY AT H POINTS------------------
      IF(KSMUD.GT.0)THEN
!$omp parallel do 
                             DO 90 L=1,LM1
          DO 50 J=MYJS_P4,MYJE_P4
          DO 50 I=MYIS_P4,MYIE_P4
      ETADT(I,J,L)=ETADT(I,J,L)*HBM2(I,J)
   50 CONTINUE
C-----------------------------------------------------------------------
      NSMUD=KSMUD
C***
C***  THE FNE, FSE, ETADTL, AND ETADT ARRAYS
C***  ARE ON OR ASSOCIATED WITH H POINTS
C***
                             DO 90 KS=1,NSMUD
          DO 80 J=MYJS_P3,MYJE1_P3
          DO 80 I=MYIS_P3,MYIE_P3
      FNE(I,J)=(ETADT(I+IHE(J),J+1,L)-ETADT(I,J,L))
     1         *HTM(I,J,L+1)*HTM(I+IHE(J),J+1,L+1)
   80 CONTINUE
          DO 82 J=MYJS1_P3,MYJE_P3
          DO 82 I=MYIS_P3,MYIE_P3
      FSE(I,J)=(ETADT(I+IHE(J),J-1,L)-ETADT(I,J,L))
     1         *HTM(I+IHE(J),J-1,L+1)*HTM(I,J,L+1)
   82 CONTINUE
          DO 84 J=MYJS2_P1,MYJE2_P1
          DO 84 I=MYIS_P1,MYIE_P1
      ETADTL(I,J)=(FNE(I,J)-FNE(I+IHW(J),J-1)
     1            +FSE(I,J)-FSE(I+IHW(J),J+1))*HBM2(I,J)
   84 CONTINUE
          DO 86 J=MYJS2_P1,MYJE2_P1
          DO 86 I=MYIS_P1,MYIE_P1
      ETADT(I,J,L)=ETADTL(I,J)*0.125+ETADT(I,J,L)
   86 CONTINUE
   90 CONTINUE
C-----------------------------------------------------------------------
      ENDIF
!-----------------------------------------------------------------------
!
!***  IF THE CFL CRITERION IS VIOLATED THEN LOCATE VERTICAL LIMITS
!***  BETWEEN WHICH TO SMOOTH THE TENDENCIES
!
!-----------------------------------------------------------------------
!$omp parallel do
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        LTOP_CFL_T(I,J)=0
        LBOT_CFL_T(I,J)=0
        LTOP_CFL_U(I,J)=0
        LBOT_CFL_U(I,J)=0
        LTOP_CFL_V(I,J)=0
        LBOT_CFL_V(I,J)=0
      ENDDO
      ENDDO
!
      DO L=1,LM1
!
!$omp parallel do private(cfl)
        DO J=MYJS2,MYJE2
        DO I=MYIS,MYIE
!
!***  MASS POINTS
!
          CFL=ETADT(I,J,L)*DTAD*HBM2(I,J)/(0.5*(DETA(L)+DETA(L+1)))
          IF(ABS(CFL).GT.CFL_MAX)THEN
            IF(LTOP_CFL_T(I,J).EQ.0)LTOP_CFL_T(I,J)=MAX(L,2)
            IF(LBOT_CFL_T(I,J).LT.L)LBOT_CFL_T(I,J)=MAX(L,2)
          ENDIF
!
!***  U COMPONENT
!
          CFL=(ETADT(I+IVW(J),J,L)+ETADT(I+IVE(J),J,L))*DTAD*VBM2(I,J)
     1       /(DETA(L)+DETA(L+1))
!
          IF(ABS(CFL).GT.CFL_MAX)THEN
            IF(LTOP_CFL_U(I,J).EQ.0)LTOP_CFL_U(I,J)=MAX(L,2)
            IF(LBOT_CFL_U(I,J).LT.L)LBOT_CFL_U(I,J)=MAX(L,2)
          ENDIF
!
!***  V COMPONENT
!
          CFL=(ETADT(I,J-1,L)+ETADT(I,J+1,L))*DTAD*VBM2(I,J)
     1       /(DETA(L)+DETA(L+1))
!
          IF(ABS(CFL).GT.CFL_MAX)THEN
            IF(LTOP_CFL_V(I,J).EQ.0)LTOP_CFL_V(I,J)=MAX(L,2)
            IF(LBOT_CFL_V(I,J).LT.L)LBOT_CFL_V(I,J)=MAX(L,2)
          ENDIF
!
        ENDDO
        ENDDO
!
      ENDDO
!
C-----------------------------------------------------------------------
C--------------PIECEWISE LINEAR UPSTREAM VERTICAL ADVECTION ------------
C-------------------------- OF Q AND CLOUD -----------------------------
C-----------------------------------------------------------------------
      ALLOCATE(S(IDIM1:IDIM2,JDIM1:JDIM2,LM),STAT=I)
C-----------------------------------------------------------------------
C       INTIALIZE Q AT THE BOTTOM INTERFACE AND THE SLOPE ADJUSTMENT
C       MASK (SAM=1 FOR SA PERMITTED, 0 FOR NOT PERMITTED)
C-----------------------------------------------------------------------
C
C***  LOOP OVER S VARIABLES (Q,CWM,T)
C
      DO 400 NS=1,3
C-----------------------------------------------------------------------
      IF(NS.EQ.1)THEN
!$omp parallel do
        DO L=1,LM
          DO J=JDIM1,JDIM2
          DO I=IDIM1,IDIM2
            S(I,J,L)=Q(I,J,L)
          ENDDO
          ENDDO
        ENDDO
      ELSE IF (NS.EQ.2) THEN
!$omp parallel do
        DO L=1,LM
          DO J=JDIM1,JDIM2
          DO I=IDIM1,IDIM2
            S(I,J,L)=CWM(I,J,L)
          ENDDO
          ENDDO
        ENDDO
      ELSE IF (NS.EQ.3) THEN
!$omp parallel do
        DO L=1,LM
          DO J=JDIM1,JDIM2
          DO I=IDIM1,IDIM2
            S(I,J,L)=T(I,J,L)
          ENDDO
          ENDDO
        ENDDO
      ELSE
        WRITE(0,*)'ERROR in VTADV. Will stop now.'
        STOP
      ENDIF
C-----------------------------------------------------------------------
C
!$omp parallel do
                             DO 175 L=1,LM
          DO 175 J=MYJS2,MYJE2
          DO 175 I=MYIS,MYIE
      QBI(I,J,L)=S(I,J,L)
      SAM(I,J,L)=1.
  175 CONTINUE
      IF(NOSLA) GO TO 290
C--------------THE SLOPE ADJUSTMENT CODE--------------------------------
C       NO SLOPE PERMITTED AT THE TOP AND AT THE BOTTOM LAYER
C-----------------------------------------------------------------------
!$omp parallel do
          DO 190 J=MYJS2,MYJE2
          DO 190 I=MYIS,MYIE
      SAM(I,J, 1)=0.
      SAM(I,J,LM)=0.
  190 CONTINUE
C
!$omp parallel do
          DO 200 L=1,LM1
      DO 200 J=MYJS2,MYJE2
      DO 200 I=MYIS,MYIE
      SAM(I,J,L)=SAM(I,J,L)*HTM(I,J,L+1)
  200 CONTINUE
C-----------------------------------------------------------------------
C       NOW, SEARCH FOR THE MAXIMA AND MINIMA OF Q (AT THE FIRST
C       PASS) AND FOR LAYERS WHICH HAD OVERADJUSTED (AT SUBSEQUENT
C       PASSES) DUE TO ROUND-OFF ERRORS
C-----------------------------------------------------------------------
!$omp parallel do private(dqbi,dqti,extrem)
                             DO 220 L=2,LM1
      DO 220 J=MYJS2,MYJE2
      DO 220 I=MYIS,MYIE
        DQTI(I,J)=S(I,J,L)-S(I,J,L-1)
        DQBI(I,J)=S(I,J,L+1)-S(I,J,L)
        EXTREM=DQTI(I,J)*DQBI(I,J)
        IF(EXTREM.LE.0.)SAM(I,J,L)=0.
  220 CONTINUE
C
!$omp parallel do 
      DO 230 L=2,LM1
      DO 230 J=MYJS2,MYJE2
      DO 230 I=MYIS,MYIE
        ARRAY1(I,J,L)=WFA(L-1)*(1.-SAM(I,J,L-1))+WFB(L-1)
        ARRAY2(I,J,L)=WFA(L)+WFB(L)*(1.-SAM(I,J,L+1))
  230 CONTINUE
                             DO 260 MSA=1,NMSAP
C-----------------------------------------------------------------------
C       CALCULATE DQ AT INTERFACES AND ADJUST THE SLOPES WHERE
C       AND TO THE EXTENT PERMITTED OBSERVING THE MONOTONICITY
C       CONDITION (E.G. VAN LEER, J. COMP. PHYS. 1977, 276-299)
C-----------------------------------------------------------------------
!$omp parallel do
          DO 240 J=MYJS2,MYJE2
          DO 240 I=MYIS,MYIE
      DQBI(I,J)=2.*S(I,J,2)-QBI(I,J,2) -QBI(I,J,1)
  240 CONTINUE
C
                             DO 250 L=2,LM1
!$omp parallel do private(asbik,astik,dqtik)
          DO 250 J=MYJS2,MYJE2
          DO 250 I=MYIS,MYIE
      DQTIK  =DQBI(I,J)
      ASTIK  =ARRAY1(I,J,L)*DQTIK
      DQBI(I,J)=2.*S(I,J,L+1)-QBI(I,J,L+1)-QBI(I,J,L)
      ASBIK  =ARRAY2(I,J,L)*DQBI(I,J)
      QBI(I,J,L)=QBI(I,J,L)
     1         +(ASTIK-SIGN(1.,ASTIK)
     2         *DIM(ABS(ASTIK),ABS(ASBIK)))*SAM(I,J,L)
  250 CONTINUE
  260                        CONTINUE
C-----------------------------------------------------------------------
C       SLOPE ADJUSTMENT OF THE LAYERS ABOVE THAT NEXT TO THE SURFACE
C       IS DONE; NOW ADJUST THE LOWERMOST LAYER
C-----------------------------------------------------------------------
                             DO 270 L=9,LM1
!$omp parallel do 
          DO 270 J=MYJS2,MYJE2
          DO 270 I=MYIS,MYIE
      IF(HTM(I,J,L+1).EQ.0.)QBI(I,J,L)=2.*S(I,J,L)-QBI(I,J,L-1)
  270 CONTINUE
C
!$omp parallel do
          DO 280 J=MYJS2,MYJE2
          DO 280 I=MYIS,MYIE
      QBI(I,J,LM)=2.*S(I,J,LM)-QBI(I,J,LM1)
  280 CONTINUE
C-----------------------------------------------------------------------
C--------------END OF THE SLOPE ADJUSTMENT CODE-------------------------
C-----------------------------------------------------------------------
  290 CONTINUE
!$omp parallel do
          DO 300 J=MYJS2,MYJE2
          DO 300 I=MYIS,MYIE
      QDEDB(I,J)=0.
      QDEUB(I,J)=0.
      DQDEB(I,J)=2.*(QBI(I,J,1)-S(I,J,1))*RDETA(1)
      EDBD (I,J)=0.
  300 CONTINUE
C
                             DO 320 L=1,LM1
!$omp parallel do private(dqdek,edbfk,edtdk,qdedtk,qdeutk,sedbk)
          DO 320 J=MYJS2,MYJE2
          DO 320 I=MYIS,MYIE
      QDEDTK  =QDEDB(I,J)
      QDEUTK  =QDEUB(I,J)
      SEDBK   =SIGN(1.,ETADT(I,J,L))
      DQDEK   =DQDEB(I,J)
      DQDEB(I,J)=2.*(QBI(I,J,L+1)-S(I,J,L+1))*RDETA(L+1)
      EDBFK   =ETADT(I,J,L)*F4D
      QDEDB(I,J)=(1.+SEDBK)*(QBI(I,J,L)+DQDEK*EDBFK)*(-EDBFK)
      QDEUB(I,J)=(1.-SEDBK)*(2.*S(I,J,L+1)-QBI(I,J,L+1)
     1           +DQDEB(I,J)*EDBFK)*EDBFK
      EDTDK   =EDBD(I,J)
      EDBD (I,J)=ETADT(I,J,L)*(-F4Q)
      S(I,J,L)=S(I,J,L)+(QDEDTK-QDEUTK-QDEDB(I,J)+QDEUB(I,J)
     1                     +S(I,J,L)*(EDBD(I,J)-EDTDK))*RDETA(L)
  320 CONTINUE
C
!$omp parallel do 
          DO 330 J=MYJS2,MYJE2
          DO 330 I=MYIS,MYIE
      S(I,J,LM)=S(I,J,LM)+(QDEDB(I,J)-QDEUB(I,J)
     1                    +S(I,J,LM)*(-EDBD(I,J)))*RDETA(LM)
  330 CONTINUE
C-------NEGATIVE MOISTURE MAY OCCUR DUE TO VIOLATION OF THE CFL---------
                             DO 350 L=1,LM1
!$omp parallel do
          DO 350 J=MYJS2,MYJE2
          DO 350 I=MYIS,MYIE
      IF(S(I,J,L).LT.EPSQ)THEN
        DQBI(I,J)=S(I,J,L)
        S(I,J,L)=EPSQ
        S(I,J,L+1)=S(I,J,L+1)+DETA(L)*RDETA(L+1)*DQBI(I,J)
      ENDIF
  350 CONTINUE
C
!$omp parallel do
          DO 360 J=MYJS2,MYJE2
          DO 360 I=MYIS,MYIE
      IF(S(I,J,LM).LT.EPSQ)S(I,J,LM)=EPSQ
  360 CONTINUE
C-----------------------------------------------------------------------
      IF(NS.EQ.1)THEN
!$omp parallel do
        DO L=1,LM
          DO J=JDIM1,JDIM2
          DO I=IDIM1,IDIM2
            Q(I,J,L)=S(I,J,L)
          ENDDO
          ENDDO
        ENDDO
      ELSE IF (NS.EQ.2) THEN
!$omp parallel do
        DO L=1,LM
          DO J=JDIM1,JDIM2
          DO I=IDIM1,IDIM2
            CWM(I,J,L)=S(I,J,L)
          ENDDO
          ENDDO
        ENDDO
      ELSE IF (NS.EQ.3) THEN
!$omp parallel do
        DO L=1,LM
          DO J=JDIM1,JDIM2
          DO I=IDIM1,IDIM2
            T(I,J,L)=S(I,J,L)
          ENDDO
          ENDDO
        ENDDO
      ELSE
        WRITE(0,*)'ERROR in VTADV. Will stop now.'
        STOP
      ENDIF
C-----------------------------------------------------------------------
  400 CONTINUE
C
      DEALLOCATE(S,STAT=IER)
C-----------------------------------------------------------------------
C--------------VERTICAL (MATSUNO) ADVECTION OF Q2-----------------------
C-----------------------------------------------------------------------
!$omp parallel do
          DO 420 J=MYJS2,MYJE2
          DO 420 I=MYIS,MYIE
      TQ2B(I,J)=Q2(I,J,1)*ETADT(I,J,1)*F4Q2(1)
  420 CONTINUE
C
                             DO 425 L=1,LM2
!$omp parallel do private(tq2ak)
          DO 425 J=MYJS2,MYJE2
          DO 425 I=MYIS,MYIE
      TQ2AK=(Q2(I,J,L+1)-Q2(I,J,L))*(ETADT(I,J,L)+ETADT(I,J,L+1))
     1                              *F4Q2(L+1)
      Q2ST(I,J,L)=TQ2AK+TQ2B(I,J)+Q2(I,J,L)
      TQ2B(I,J)=TQ2AK
  425 CONTINUE
C
!$omp parallel do private(tq2ak)
          DO 440 J=MYJS2,MYJE2
          DO 440 I=MYIS,MYIE
      TQ2AK=(Q2(I,J,LM)-Q2(I,J,LM1))*ETADT(I,J,LM1)*F4Q2(LM)
      Q2ST(I,J,LM1)=TQ2AK+TQ2B(I,J)+Q2(I,J,LM1)
      Q2ST(I,J,LM )=Q2(I,J,LM)
  440 CONTINUE
C-----------------------------------------------------------------------
C--------------SECOND (BACKWARD) MATSUNO STEP---------------------------
C-----------------------------------------------------------------------
!$omp parallel do 
          DO 450 J=MYJS2,MYJE2
          DO 450 I=MYIS,MYIE
      TQ2B(I,J)=Q2ST(I,J,1)*ETADT(I,J,1)*F4Q2(1)
  450 CONTINUE
C
      DO L=1,LM2
!$omp parallel do private(tq2ak)
        DO J=MYJS2,MYJE2
        DO I=MYIS,MYIE
          TQ2AK  =(Q2ST(I,J,L+1)-Q2ST(I,J,L))
     1           *(ETADT(I,J,L)+ETADT(I,J,L+1))*F4Q2(L+1)
          VAD_TEND1(I,J,L)=TQ2AK+TQ2B(I,J)
          TQ2B(I,J)=TQ2AK
        ENDDO
        ENDDO
      ENDDO
!
      DO J=MYJS2,MYJE2
      DO I=MYIS,MYIE
        TQ2AK  =(Q2ST(I,J,LM)-Q2ST(I,J,LM1))*ETADT(I,J,LM1)*F4Q2(LM)
        VAD_TEND1(I,J,LM1)=TQ2AK+TQ2B(I,J)
      ENDDO
      ENDDO
!-----------------------------------------------------------------------
!
!***  IF THE CFL CRITERION IS VIOLATED THEN VERTICALLY SMOOTH
!***  THE TENDENCY
!
!-----------------------------------------------------------------------
!
!$omp parallel do
!$omp& private(cfl,lbot_cfl,lstart,lstop,ltop_cfl,vad_tend1,vad_tndx1)
      DO J=MYJS2,MYJE2
      DO I=MYIS,MYIE
!
        IF(LTOP_CFL_T(I,J).GT.0)THEN
          LSTART=LTOP_CFL_T(I,J)
          LSTOP =MIN(LBOT_CFL_T(I,J),LM-2)
!
          DO L=LSTART,LSTOP
            VAD_TNDX1(L)=(VAD_TEND1(I,J,L-1)+VAD_TEND1(I,J,L+1)
     1                   +2.*VAD_TEND1(I,J,L))*0.25
          ENDDO
          DO L=LSTART,LSTOP
            VAD_TEND1(I,J,L)=VAD_TNDX1(L)
          ENDDO
        ENDIF
!
      ENDDO
      ENDDO
C
                             DO 470 L=1,LM2
!$omp parallel do 
          DO 470 J=MYJS2,MYJE2
          DO 470 I=MYIS,MYIE
      Q2(I,J,L)=VAD_TEND1(I,J,L)+Q2(I,J,L)
      Q2(I,J,L)=AMAX1(Q2(I,J,L),EPSQ2)
  470 CONTINUE
C
!$omp parallel do
          DO 480 J=MYJS2,MYJE2
          DO 480 I=MYIS,MYIE
      Q2(I,J,LM1)=VAD_TEND1(I,J,LM1)+Q2(I,J,LM1)
      Q2(I,J,LM1)=AMAX1(Q2(I,J,LM1),EPSQ2)
  480 CONTINUE
C-----------------------------------------------------------------------
C--------------DEFINITION OF VARIABLES NEEDED AT V POINTS---------------
C-----------------------------------------------------------------------
!$omp parallel do 
                             DO 500 L=1,LM1
          DO 500 J=MYJS_P1,MYJE_P1
          DO 500 I=MYIS_P1,MYIE_P1
      ETADT(I,J,L)=ETADT(I,J,L)*PDSL(I,J)*HBM2(I,J)
  500 CONTINUE
C
!$omp parallel do
          DO 510 J=MYJS2,MYJE2
          DO 510 I=MYIS,MYIE
      RPDX(I,J)=1./(PDSL(I+IVW(J),J)+PDSL(I+IVE(J),J))
      RPDY(I,J)=1./(PDSL(I,J-1)+PDSL(I,J+1))
  510 CONTINUE

C---------- PIECEWISE LINEAR UPSTREAM VERTICAL ADVECTION OF U & V ------
C    INTIALIZE U & V AT THE BOTTOM INTERFACE AND THE SLOPE ADJUSTMENT
C    MASK (SAMU=1 AND SAMV=1 FOR SA PERMITTED, 0 FOR NOT PERMITTED)
C-----------------------------------------------------------------------
!$omp parallel do
                             DO 1175 L=1,LM
          DO 1175 J=MYJS2,MYJE2
          DO 1175 I=MYIS,MYIE
      UBI(I,J,L)=U(I,J,L)
      VBI(I,J,L)=V(I,J,L)
      SAMU(I,J,L)=1.
      SAMV(I,J,L)=1.
 1175 CONTINUE
C-----------------------------------------------------------------------
      IF(NOSLAW) GO TO 1290
C--------------THE SLOPE ADJUSTMENT CODE--------------------------------
C       NO SLOPE PERMITTED AT THE TOP AND AT THE BOTTOM LAYER
C-----------------------------------------------------------------------
!$omp parallel do
          DO 1190 J=MYJS2,MYJE2
          DO 1190 I=MYIS,MYIE
      SAMU(I,J, 1)=0.
      SAMU(I,J,LM)=0.
      SAMV(I,J, 1)=0.
      SAMV(I,J,LM)=0.   
 1190 CONTINUE
C
!$omp parallel do
          DO 1200 L=1,LM1
      DO 1200 J=MYJS2,MYJE2
      DO 1200 I=MYIS,MYIE
      SAMU(I,J,L)=SAMU(I,J,L)*VTM(I,J,L+1)
      SAMV(I,J,L)=SAMV(I,J,L)*VTM(I,J,L+1)
 1200 CONTINUE
C-----------------------------------------------------------------------
C       NOW, SEARCH FOR THE MAXIMA AND MINIMA OF U & V (AT THE FIRST
C       PASS) AND FOR LAYERS WHICH HAD OVERADJUSTED (AT SUBSEQUENT
C       PASSES) DUE TO ROUND-OFF ERRORS
C-----------------------------------------------------------------------
!$omp parallel do private(dubi, dvbi, duti, dvti, extremu, extremv)
                             DO 1220 L=2,LM1
      DO 1220 J=MYJS2,MYJE2
      DO 1220 I=MYIS,MYIE
        DUTI(I,J)=U(I,J,L)-U(I,J,L-1)
        DVTI(I,J)=V(I,J,L)-V(I,J,L-1)
        DUBI(I,J)=U(I,J,L+1)-U(I,J,L)
        DVBI(I,J)=V(I,J,L+1)-V(I,J,L)
        EXTREMU=DUTI(I,J)*DUBI(I,J)
        EXTREMV=DVTI(I,J)*DVBI(I,J)
        IF(EXTREMU.LT.0.) SAMU(I,J,L)=0.
        IF(EXTREMV.LT.0.) SAMV(I,J,L)=0.
 1220 CONTINUE
C
!$omp parallel do 
      DO 1230 L=2,LM1
      DO 1230 J=MYJS2,MYJE2
      DO 1230 I=MYIS,MYIE
        ARRAYU1(I,J,L)=WFA(L-1)*(1.-SAMU(I,J,L-1))+WFB(L-1)
        ARRAYV1(I,J,L)=WFA(L-1)*(1.-SAMV(I,J,L-1))+WFB(L-1)
        ARRAYU2(I,J,L)=WFA(L)+WFB(L)*(1.-SAMU(I,J,L+1))
        ARRAYV2(I,J,L)=WFA(L)+WFB(L)*(1.-SAMV(I,J,L+1))
 1230 CONTINUE
C-----------------------------------------------------------------------
                            DO 1260 MSA=1,NMSAPW
C-----------------------------------------------------------------------
C       CALCULATE DU & DV AT INTERFACES AND ADJUST THE SLOPES WHERE
C       AND TO THE EXTENT PERMITTED OBSERVING THE MONOTONICITY
C       CONDITION (E.G. VAN LEER, J. COMP. PHYS. 1977, 276-299,
C       scheme used here of MESINGER AND JOVIC, NCEP OFFICE NOTE #439)
C-----------------------------------------------------------------------
!$omp parallel do
          DO 1240 J=MYJS2,MYJE2
          DO 1240 I=MYIS,MYIE
      DUBI(I,J)=2.*U(I,J,2)-UBI(I,J,2) -UBI(I,J,1)
      DVBI(I,J)=2.*V(I,J,2)-VBI(I,J,2) -VBI(I,J,1)
 1240 CONTINUE
C-----------------------------------------------------------------------
!$omp parallel do private(asbiku, asbikv, astiku, astikv, dutik, dvtik)
                             DO 1250 L=2,LM1
          DO 1250 J=MYJS2,MYJE2
          DO 1250 I=MYIS,MYIE
      DUTIK  =DUBI(I,J)
      DVTIK  =DVBI(I,J)
      ASTIKU =ARRAYU1(I,J,L)*DUTIK
      ASTIKV =ARRAYV1(I,J,L)*DVTIK
      DUBI(I,J)=2.0*U(I,J,L+1)-UBI(I,J,L+1)-UBI(I,J,L)
      DVBI(I,J)=2.0*V(I,J,L+1)-VBI(I,J,L+1)-VBI(I,J,L)
      ASBIKU  =ARRAYU2(I,J,L)*DUBI(I,J)
      ASBIKV  =ARRAYV2(I,J,L)*DVBI(I,J)
      UBI(I,J,L)=UBI(I,J,L)
     1         +(ASTIKU-SIGN(1.,ASTIKU)
     2         *DIM(ABS(ASTIKU),ABS(ASBIKU)))*SAMU(I,J,L)
      VBI(I,J,L)=VBI(I,J,L)
     1         +(ASTIKV-SIGN(1.,ASTIKV)
     2         *DIM(ABS(ASTIKV),ABS(ASBIKV)))*SAMV(I,J,L)
 1250 CONTINUE
 1260                        CONTINUE
C-----------------------------------------------------------------------
C       SLOPE ADJUSTMENT OF THE LAYERS ABOVE THAT NEXT TO THE SURFACE
C       IS DONE; NOW ADJUST THE LOWERMOST LAYER
C-----------------------------------------------------------------------
                             DO 1270 L=9,LM1
!$omp parallel do 
          DO 1270 J=MYJS2,MYJE2
          DO 1270 I=MYIS,MYIE
      IF(VTM(I,J,L+1).EQ.0.) UBI(I,J,L)=2.*U(I,J,L)-UBI(I,J,L-1)
      IF(VTM(I,J,L+1).EQ.0.) VBI(I,J,L)=2.*V(I,J,L)-VBI(I,J,L-1)
 1270 CONTINUE
C
!$omp parallel do
          DO 1280 J=MYJS2,MYJE2
          DO 1280 I=MYIS,MYIE
      UBI(I,J,LM)=2.*U(I,J,LM)-UBI(I,J,LM1)
      VBI(I,J,LM)=2.*V(I,J,LM)-VBI(I,J,LM1)
 1280 CONTINUE
C-----------------------------------------------------------------------
C--------------END OF THE SLOPE ADJUSTMENT CODE-------------------------
C-----------------------------------------------------------------------
 1290 CONTINUE
C-----------------------------------------------------------------------
!$omp parallel do
          DO 1300 J=MYJS2,MYJE2
          DO 1300 I=MYIS,MYIE
      UDEDB(I,J)=0.
      VDEDB(I,J)=0.
      UDEUB(I,J)=0.
      VDEUB(I,J)=0.
      DUDEB(I,J)=2.*(UBI(I,J,1)-U(I,J,1))*RDETA(1)
      DVDEB(I,J)=2.*(VBI(I,J,1)-V(I,J,1))*RDETA(1)
      EDBDU (I,J)=0.
      EDBDV(I,J)=0.
 1300 CONTINUE
                             DO 1320 L=1,LM1
!$omp parallel do private(dqdek,edbfk,edtdk,qdedtk,qdeutk,sedbk)
          DO 1320 J=MYJS2,MYJE2
          DO 1320 I=MYIS,MYIE
      VMIJ=VTM(I,J,L)*VBM2(I,J)
      UDEDTK  =UDEDB(I,J)
      VDEDTK  =VDEDB(I,J)
      UDEUTK  =UDEUB(I,J)
      VDEUTK  =VDEUB(I,J)
      SEDBKU  =SIGN(1.,(ETADT(I+IVW(J),J,L)+ETADT(I+IVE(J),J,L)))
      SEDBKV  =SIGN(1.,(ETADT(I,J-1,L)+ETADT(I,J+1,L)))
      DUDEK   =DUDEB(I,J)
      DVDEK   =DVDEB(I,J)
      DUDEB(I,J)=2.*(UBI(I,J,L+1)-U(I,J,L+1))*RDETA(L+1)
      DVDEB(I,J)=2.*(VBI(I,J,L+1)-V(I,J,L+1))*RDETA(L+1)
C
      EDBFKU =(ETADT(I+IVW(J),J,L)+ETADT(I+IVE(J),J,L))
     1         *RPDX(I,J)*F4D
      EDBFKV =(ETADT(I,J-1,L)+ETADT(I,J+1,L))*RPDY(I,J)*F4D
C
      UDEDB(I,J)=(1.+SEDBKU)*(UBI(I,J,L)+DUDEK*EDBFKU)*(-EDBFKU)
      VDEDB(I,J)=(1.+SEDBKV)*(VBI(I,J,L)+DVDEK*EDBFKV)*(-EDBFKV)
      UDEUB(I,J)=(1.-SEDBKU)*(U(I,J,L+1)+U(I,J,L+1)-UBI(I,J,L+1)
     1           +DUDEB(I,J)*EDBFKU)*EDBFKU
      VDEUB(I,J)=(1.-SEDBKV)*(V(I,J,L+1)+V(I,J,L+1)-VBI(I,J,L+1)
     1           +DVDEB(I,J)*EDBFKV)*EDBFKV
      EDTDKU    =EDBDU(I,J)
      EDTDKV    =EDBDV(I,J)
C
      EDBDU(I,J)=(ETADT(I+IVW(J),J,L)+ETADT(I+IVE(J),J,L))
     1           *RPDX(I,J)*(-F4Q)
      EDBDV(I,J)=(ETADT(I,J-1,L)+ETADT(I,J+1,L))*RPDY(I,J)*(-F4Q)
C
      VAD_TEND1(I,J,L)=(UDEDTK-UDEUTK-UDEDB(I,J)+UDEUB(I,J)
     1                 +U(I,J,L)*(EDBDU(I,J)-EDTDKU))*RDETA(L)*VMIJ
      VAD_TEND2(I,J,L)=(VDEDTK-VDEUTK-VDEDB(I,J)+VDEUB(I,J)
     1                 +V(I,J,L)*(EDBDV(I,J)-EDTDKV))*RDETA(L)*VMIJ
 1320 CONTINUE
C-----------------------------------------------------------------------
!$omp parallel do
          DO 1330 J=3,JM-2
          DO 1330 I=2,IM-1
      VMIJ=VTM(I,J,LM)*VBM2(I,J)
      VAD_TEND1(I,J,LM)=(UDEDB(I,J)-UDEUB(I,J)
     1                  +U(I,J,LM)*(-EDBDU(I,J)))*RDETA(LM)*VMIJ
      VAD_TEND2(I,J,LM)=(VDEDB(I,J)-VDEUB(I,J)
     1                  +V(I,J,LM)*(-EDBDV(I,J)))*RDETA(LM)*VMIJ
 1330 CONTINUE
!
!-----------------------------------------------------------------------
!
!***  IF THE CFL CRITERION IS VIOLATED THEN VERTICALLY SMOOTH
!***  THE TENDENCIES
!
!-----------------------------------------------------------------------
!
!$omp parallel do
!$omp& private(lstart,lstop,vad_tndx1,vad_tndx2)
      DO J=MYJS2,MYJE2
      DO I=MYIS,MYIE
!
!***  U COMPONENT
!
        IF(LTOP_CFL_U(I,J).GT.0)THEN
          LSTART=LTOP_CFL_U(I,J)
          LSTOP =MIN(LBOT_CFL_U(I,J),LM-1)
!
          DO L=LSTART,LSTOP
            VAD_TNDX1(L)=(VAD_TEND1(I,J,L-1)+VAD_TEND1(I,J,L+1)
     1                   +2.*VAD_TEND1(I,J,L))*0.25
          ENDDO
          DO L=LSTART,LSTOP
            VAD_TEND1(I,J,L)=VAD_TNDX1(L)
          ENDDO
        ENDIF
!
!***  V COMPONENT
!
        IF(LTOP_CFL_V(I,J).GT.0)THEN
          LSTART=LTOP_CFL_V(I,J)
          LSTOP =MIN(LBOT_CFL_V(I,J),LM-1)
!
          DO L=LSTART,LSTOP
            VAD_TNDX2(L)=(VAD_TEND2(I,J,L-1)+VAD_TEND2(I,J,L+1)
     1                   +2.*VAD_TEND2(I,J,L))*0.25
          ENDDO
          DO L=LSTART,LSTOP
            VAD_TEND2(I,J,L)=VAD_TNDX2(L)
          ENDDO
        ENDIF
!
      ENDDO
      ENDDO
C
                             DO 580 L=1,LM
!$omp parallel do
          DO 580 J=MYJS2,MYJE2
          DO 580 I=MYIS,MYIE
      U(I,J,L)=U(I,J,L)+VAD_TEND1(I,J,L)
      V(I,J,L)=V(I,J,L)+VAD_TEND2(I,J,L)
  580 CONTINUE
C
C-----------------------------------------------------------------------
      IF(.NOT.HYDRO)THEN
!$omp parallel do
        DO L=1,LM1
          DO J=MYJS_P1,MYJE_P1
          DO I=MYIS_P1,MYIE_P1
            ETADT(I,J,L)=ETADT(I,J,L)/PDSL(I,J)
          ENDDO
          ENDDO
        ENDDO
      ENDIF
C
C-----------------------------------------------------------------------
Cpazi
Cpazi (Changing constats back to what they were)
C      DTAD=DTAD*2.0
C      F4D=F4D*2.0
C      F4Q=F4Q*2.0
C      DO L=1,LM
C      F4Q2(L)=F4Q2(L)*2.0
C      END DO
Cpazi
                             RETURN
                             END