C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 
                            SUBROUTINE ADVEPS 
C     ****************************************************************** 
C$$$  SUBPROGRAM DOCUMENTATION BLOCK 
C                .      .    . 
C SUBPROGRAM:    ADV2        HORIZ/VERT ADVECTION OF VAPOR AND CLOUD 
C   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 96-07-19 
C 
C ABSTRACT: 
C     ADV2 CALCULATES THE CONTRIBUTION OF THE HORIZONTAL AND VERTICAL
C     ADVECTION TO THE TENDENCIES OF WATER SUBSTANCE AND TKE AND THEN 
C     UPDATES THOSE VARIABLES.  AN ANTI-FILTERING TECHNIQUE IS USED. 
C 
C PROGRAM HISTORY LOG: 
C   96-07-19  JANJIC  - ORIGINATOR 
C   96-08-    BLACK   - CRAY DIRECTIVES, CRAY OPTIMIZATION 
C   99-07-    JANJIC  - ADDED HORIZONTAL ADVECTION OF TKE 
C   99-07-    JANJIC  - ADDED VERTICAL ADVECTION 
C   00-01-06  BLACK   - DISTRIBUTED MEMORY AND THREADS
C 
C USAGE: CALL ADV2 FROM MAIN PROGRAM 
C   INPUT ARGUMENT LIST: 
C       NONE 
C 
C   OUTPUT ARGUMENT LIST 
C       NONE 
C 
C   OUTPUT FILES: 
C       NONE 
C   SUBPROGRAMS CALLED: 
C 
C     UNIQUE: NONE 
C 
C     LIBRARY: NONE 
C 
C   COMMON BLOCKS: CTLBLK 
C                  LOOPS 
C                  MASKS 
C                  DYNAM 
C                  CONTIN 
C                  VRBLS 
C                  PVRBLS 
C                  CLDWTR 
C                  INDX 
C 
C ATTRIBUTES: 
C   LANGUAGE: FORTRAN 90 
C   MACHINE : IBM SP
C$$$ 
C*********************************************************************** 
      INCLUDE "parmeta" 
      INCLUDE "mpp.h"
      INCLUDE "mpif.h"
C----------------------------------------------------------------------- 
                             P A R A M E T E R 
     &(IM1=IM-1,JAM=6+2*(JM-10),JAMD=(JAM*2-10)*3,LP1=LM+1) 
C----------------------------------------------------------------------- 
                             P A R A M E T E R
     &(EPSQ=1.E-12,EPSQ2=0.2,CLIMIT=1.E-20
     &,FF1=0.52500,FF2=-0.64813,FF3=0.24520,FF4=-0.12189) 
C-----------------------------------------------------------------------
                             P A R A M E T E R
     &(NTSHY=3,EPSFC=1./1.05,EPSN=-EPSFC,EPSP=EPSFC,ZERO=1.E-06
     &,G=9.806,CP=1004.6,CAPA=287.04/CP,GMA=-287.04*(1.-CAPA)/2.
     &,cappa=287.04/1004.6)
C----------------------------------------------------------------------- 
                             L O G I C A L
     & RUN,FIRST,RESTRT,SIGMA,TOP,BOT,nonh 
C----------------------------------------------------------------------- 
      INCLUDE "CTLBLK.comm" 
C----------------------------------------------------------------------- 
      INCLUDE "LOOPS.comm" 
C----------------------------------------------------------------------- 
      INCLUDE "MASKS.comm" 
C----------------------------------------------------------------------- 
      INCLUDE "DYNAM.comm" 
C----------------------------------------------------------------------- 
      INCLUDE "CONTIN.comm" 
C----------------------------------------------------------------------- 
      INCLUDE "VRBLS.comm" 
C----------------------------------------------------------------------- 
      INCLUDE "PVRBLS.comm" 
C-----------------------------------------------------------------------
      INCLUDE "NHYDRO.comm"
C----------------------------------------------------------------------- 
      INCLUDE "CLDWTR.comm" 
C----------------------------------------------------------------------- 
      INCLUDE "INDX.comm" 
C
C-------------LOCAL VARIABLES REQUIRING DECLARATION--------------------- 
C
C-------------LOCAL VARIABLES REQUIRING DECLARATION---------------------
                             R E A L
     & c0(lm),b1(lm),b2(lm),b3(lm)
     &,PONE(LM+1),PSTR(LM+1),PNP1(LM+1),COFF(LM+1)
     &,CHI(LM+1)
C-----------------------------------------------------------------------
      real v3(lm),v4(lm)
                             D I M E N S I O N
     & Q3(LM),W3(LM),E3(LM),Q4(LM),W4(LM),E4(LM) 
     &,petdtl(LM),DQL(LM),DWL(LM),DEL(LM),dvl(lm),LA(LM),AFR(LM)  
C 
                             D I M E N S I O N
     & IFPA(IDIM1:IDIM2,JDIM1:JDIM2,LM),IFQA(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,IFPF(IDIM1:IDIM2,JDIM1:JDIM2,LM),IFQF(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,JFPA(IDIM1:IDIM2,JDIM1:JDIM2,LM),JFQA(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,JFPF(IDIM1:IDIM2,JDIM1:JDIM2,LM),JFQF(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,AFP (IDIM1:IDIM2,JDIM1:JDIM2,LM),AFQ (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,Q1  (IDIM1:IDIM2,JDIM1:JDIM2,LM),DQST(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,W1  (IDIM1:IDIM2,JDIM1:JDIM2,LM),DWST(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,E1  (IDIM1:IDIM2,JDIM1:JDIM2,LM),DEST(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,v1  (idim1:idim2,jdim1:jdim2,lm),dvst(idim1:idim2,jdim1:jdim2,lm)
     &,E2  (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,DARE(IDIM1:IDIM2,JDIM1:JDIM2),   DVOL(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,EMH (IDIM1:IDIM2,JDIM1:JDIM2)
     &,ANE (IDIM1:IDIM2,JDIM1:JDIM2)   ,ASE (IDIM1:IDIM2,JDIM1:JDIM2)
C----------------------------------------------------------------------- 
                             R E A L
     & GSUMS(8,LM),XSUMS(8,LM)
C-----------------------------------------------------------------------
C*********************************************************************** 
C-----------------------------------------------------------------------
      IF(NTSD.LE.NTSHY.OR.HYDRO)THEN
C***
!$omp parallel do
        DO J=MYJS_P2,MYJE_P2
        DO I=MYIS_P1,MYIE_P1
          PINT(I,J,1)=PT
        ENDDO
        ENDDO
C
        DO L=1,LM
!$omp parallel do
          DO J=MYJS_P2,MYJE_P2
          DO I=MYIS_P1,MYIE_P1
            DWDT(I,J,L)=1.
            PDWDT(I,J,L)=1.
            pint(i,j,l+1)=deta1(l)*pdtop+deta2(l)*pdsl(i,j)+pint(i,j,l)
          ENDDO
          ENDDO
        ENDDO
C***
        go to 7777
C***
      ENDIF
C-----------------------------------------------------------------------
      ADDT=DT
      RDT=1./ADDT
C-----------------TIME TENDENCY-----------------------------------------
!$omp parallel do
      DO L=1,LM
        DO J=MYJS_P1,MYJE_P1
        DO I=MYIS_P1,MYIE_P1
          DWDT(I,J,L)=(W(I,J,L)-DWDT(I,J,L))*HTM(I,J,L)*HBM2(I,J)*RDT
        ENDDO
        ENDDO
      ENDDO
C-----------------------------------------------------------------------
 7777 ADDT=FLOAT(IDTAD)*DT 
C----------------------------------------------------------------------- 
C-----------------VERTICAL ADVECTION------------------------------------ 
C----------------------------------------------------------------------- 
!$omp parallel do
!$omp& private (afr,afrp,arr,bot,d2pqe,d2pqq,d2pqw,del,dep,detap,
!$omp&          dpdn,dpup,dql,dqp,dwl,dwp,e00,e3,e4,e4p,ep,ep0,
!$omp&          hbm2ij,la,lap,llap,lmp,petdtl,q00,q3,q4,q4p,qp,qp0,
!$omp&          rface,rfacq,rfacw,rr,sumne,sumnq,sumnw,
!$omp&          sumpe,sumpq,sumpw,top,w00,w3,w4,w4p,wp,wp0)
C----------------------------------------------------------------------- 
      DO 200 J=MYJS2_P2,MYJE2_P2
      DO 200 I=MYIS1_P2,MYIE1_P2
C----------------------------------------------------------------------- 
      LMP=LMH(I,J) 
C----------------------------------------------------------------------- 
      E3(1)=Q2(I,J,1)*0.5 
      DO L=2,LMP 
        E3(L)=MAX((Q2(I,J,L-1)+Q2(I,J,L))*0.5,EPSQ2) 
      ENDDO 
C
      DO L=1,LMP 
        Q3(L)=MAX(Q(I,J,L),EPSQ) 
        W3(L)=MAX(CWM(I,J,L),CLIMIT) 
        v3(l)=w(i,j,l)
        E4(L)=E3(L) 
        Q4(L)=Q3(L) 
        W4(L)=W3(L) 
        v4(l)=v3(l)
      ENDDO 
C 
      petdtl(1)=petdt(I,J,1)*0.5 
C
      DO L=2,LMP-1 
        petdtl(L)=(petdt(I,J,L-1)+petdt(I,J,L))*0.5 
      ENDDO 
C
      petdtl(LMP)=petdt(I,J,LMP-1)*0.5 
C----------------------------------------------------------------------- 
      SUMPQ=0. 
      SUMNQ=0. 
      SUMPW=0. 
      sumpv=0.
      SUMNW=0. 
      SUMPE=0. 
      SUMNE=0. 
C 
      DO 25 L=1,LMP 
      RR=petdtl(L)*(-ADDT) 
      ARR=ABS(RR) 
C
      IF(ARR.GT.0.)THEN
        LAP=RR/ARR
      ELSE
        LAP=0
      ENDIF
C 
      LA(L)=LAP 
      LLAP=L+LAP 
C 
      TOP=.FALSE. 
      BOT=.FALSE. 
C 
      IF(LLAP.GT.0.AND.LLAP.LT.LMH(I,J)+1.AND.LAP.NE.0)THEN 
        rr=arr/abs((aeta1(llap)-aeta1(l))*pdtop
     &            +(aeta2(llap)-aeta2(l))*pdsl(i,j))
        AFR(L)=(((FF4*RR+FF3)*RR+FF2)*RR+FF1)*RR 
        DQP=(Q3(LLAP)-Q3(L))*RR 
        DWP=(W3(LLAP)-W3(L))*RR 
        DEP=(E3(LLAP)-E3(L))*RR 
        dvp=(v3(llap)-v3(l))*rr
        DQL(L)=DQP 
        DWL(L)=DWP            
        DEL(L)=DEP 
        dvl(l)=dvp
      ELSE 
        TOP=LLAP.EQ.0 
        BOT=LLAP.EQ.LMH(I,J)+1 
C 
        RR=0. 
        AFR(L)=0. 
        DQL(L)=0. 
        DWL(L)=0.            
        DEL(L)=0. 
        dvl(l)=0.
      ENDIF 
   25 CONTINUE
C 
      IF(TOP)THEN 
        IF(LA(2).LT.0)THEN 
          rfc=(deta1(2)*pdtop+deta2(2)*pdsl(i,j))
     &       /(deta1(1)*pdtop+deta2(1)*pdsl(i,j))
          dql(1)=-dql(2)*rfc
          dwl(1)=-dwl(2)*rfc
          del(1)=-del(2)*rfc 
          dvl(1)=-dvl(2)*rfc
        ENDIF 
      ENDIF 
C
      IF(BOT)THEN 
        IF(LA(LMP-1).GT.0)THEN 
          rfc=(deta1(lmp-1)*pdtop+deta2(lmp-1)*pdsl(i,j))
     &       /(deta1(lmp)*pdtop+deta2(lmp)*pdsl(i,j))
          dql(lmp)=-dql(lmp-1)*rfc
          dwl(lmp)=-dwl(lmp-1)*rfc
          del(lmp)=-del(lmp-1)*rfc
          dvl(lmp)=-dvl(lmp-1)*rfc
        ENDIF 
      ENDIF 
C 
      DO L=1,LMP 
        detap=(deta1(l)*pdtop+deta2(l)*pdsl(i,j)) 
        DQP=DQL(L)*detap 
        DWP=DWL(L)*detap 
        DEP=DEL(L)*detap 
        dvp=dvl(l)*detap
C
        IF(DQP.GT.0.)THEN 
          SUMPQ=SUMPQ+DQP 
        ELSE 
          SUMNQ=SUMNQ+DQP 
        ENDIF 
C
        IF(DWP.GT.0.)THEN 
          SUMPW=SUMPW+DWP 
        ELSE 
          SUMNW=SUMNW+DWP 
        ENDIF 
C
        IF(DEP.GT.0.)THEN 
          SUMPE=SUMPE+DEP 
        ELSE 
          SUMNE=SUMNE+DEP 
        ENDIF 
!
        if(dvp.gt.0.) then
          sumpv=sumpv+dvp
        else
          sumnv=sumnv+dvp
        endif
C
      ENDDO 
C
C--------------FIRST MOMENT CONSERVING FACTOR--------------------------- 
C
      IF(SUMPQ.GT.1.E-9)THEN 
        RFACQ=-SUMNQ/SUMPQ 
      ELSE 
        RFACQ=1. 
      ENDIF 
C
      IF(SUMPW.GT.1.E-9)THEN 
        RFACW=-SUMNW/SUMPW 
      ELSE 
        RFACW=1. 
      ENDIF 
C
      IF(SUMPE.GT.1.E-9)THEN 
        RFACE=-SUMNE/SUMPE 
      ELSE 
        RFACE=1. 
      ENDIF 
C
      if(sumpv.gt.1.e-9)then
        rfacv=-sumnv/sumpv
      else
        rfacv=1.
      endif
C
      IF(RFACQ.LT.0.9.OR.RFACQ.GT.1.1) RFACQ=1. 
      IF(RFACW.LT.0.9.OR.RFACW.GT.1.1) RFACW=1. 
      IF(RFACE.LT.0.9.OR.RFACE.GT.1.1) RFACE=1. 
      if(rfacv.lt.0.9.or.rfacv.gt.1.1) rfacv=1.
C
C--------------IMPOSE CONSERVATION ON ADVECTION------------------------- 
C
      IF(RFACQ.LT.1.)THEN 
        DO L=1,LMP 
          DQP=DQL(L) 
          IF(DQP.LT.0.) DQP=DQP/RFACQ 
          Q4(L)=Q3(L)+DQP 
        ENDDO 
      ELSE 
        DO L=1,LMP 
          DQP=DQL(L) 
          IF(DQP.GE.0.) DQP=DQP*RFACQ 
          Q4(L)=Q3(L)+DQP 
        ENDDO 
      ENDIF 
C----------------------------------------------------------------------- 
      IF(RFACW.LT.1.)THEN 
        DO L=1,LMP 
          DWP=DWL(L) 
          IF(DWP.LT.0.) DWP=DWP/RFACW 
          W4(L)=W3(L)+DWP 
        ENDDO 
      ELSE 
        DO L=1,LMP 
          DWP=DWL(L) 
          IF(DWP.GE.0.) DWP=DWP*RFACW 
          W4(L)=W3(L)+DWP 
        ENDDO 
      ENDIF 
C----------------------------------------------------------------------- 
      IF(RFACE.LT.1.)THEN 
        DO L=1,LMP 
          DEP=DEL(L) 
          IF(DEP.LT.0.) DEP=DEP/RFACE 
          E4(L)=E3(L)+DEP 
        ENDDO 
      ELSE 
        DO L=1,LMP 
          DEP=DEL(L) 
          IF(DEP.GE.0.) DEP=DEP*RFACE 
          E4(L)=E3(L)+DEP 
        ENDDO 
      ENDIF 
C-----------------------------------------------------------------------
      if(rfacv.lt.1.)then
        do l=1,lmp
          dvp=dvl(l)
          if(dvp.lt.0.) dvp=dvp/rfacv
          v4(l)=v3(l)+dvp
        enddo
      else
        do l=1,lmp
          dvp=dvl(l)
          if(dvp.ge.0.) dvp=dvp*rfacv
          v4(l)=v3(l)+dvp
        enddo
      endif
C--------------ANTI-FILTERING STEP-------------------------------------- 
      SUMPQ=0. 
      SUMNQ=0. 
      SUMPW=0. 
      sumpv=0.
      SUMNW=0. 
      SUMPE=0. 
      SUMNE=0.
      sumnv=0. 
C--------------ANTI-FILTERING LIMITERS---------------------------------- 
      DO 50 L=2,LMP-1 
      detap=deta1(l)*pdtop+deta2(l)*pdsl(i,j) 
C 
      Q4P=Q4(L) 
      W4P=W4(L) 
      E4P=E4(L) 
      v4p=v4(l)
C 
      LAP=LA(L) 
C 
      IF(LAP.NE.0)THEN 
        dpdn=(aeta1(l+lap)-aeta1(l))*pdtop
     &      +(aeta2(l+lap)-aeta2(l))*pdsl(i,j)
        dpup=(aeta1(l)-aeta1(l-lap))*pdtop
     &      +(aeta2(l)-aeta2(l-lap))*pdsl(i,j)
!
        afrp=2.*afr(l)*dpdn*dpdn/(dpdn+dpup)
        d2pqq=((q4(l+lap)-q4p)/dpdn 
     1        -(q4p-q4(l-lap))/dpup)*afrp 
        d2pqw=((w4(l+lap)-w4p)/dpdn 
     1        -(w4p-w4(l-lap))/dpup)*afrp 
        d2pqe=((e4(l+lap)-e4p)/dpdn 
     1        -(e4p-e4(l-lap))/dpup)*afrp 
        d2pqv=((v4(l+lap)-v4p)/dpdn
     1        -(v4p-v4(l-lap))/dpup)*afrp
      ELSE 
        D2PQQ=0. 
        D2PQW=0. 
        D2PQE=0. 
        d2pqv=0.
      ENDIF 
C 
      QP=Q4P-D2PQQ 
      WP=W4P-D2PQW 
      EP=E4P-D2PQE 
      vp=v4p-d2pqv
C 
      Q00=Q3(L) 
      QP0=Q3(L+LAP) 
C 
      W00=W3(L) 
      WP0=W3(L+LAP) 
C 
      E00=E3(L) 
      EP0=E3(L+LAP) 
C
      v00=v3(l)
      vp0=v3(l+lap)
C 
      IF(LAP.NE.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)) 
        vp=max(vp,min(v00,vp0))
        vp=min(vp,max(v00,vp0))
      ENDIF 
C 
      DQP=QP-Q4P 
      DWP=WP-W4P 
      DEP=EP-E4P 
      dvp=vp-v4p
C 
      DQL(L)=DQP 
      DWL(L)=DWP 
      DEL(L)=DEP 
      dvl(l)=dvp
C 
      dqp=dqp*detap 
      dwp=dwp*detap 
      dep=dep*detap 
      dvp=dvp*detap
C 
      IF(DQP.GT.0.)THEN 
        SUMPQ=SUMPQ+DQP 
      ELSE 
        SUMNQ=SUMNQ+DQP 
      ENDIF 
C 
      IF(DWP.GT.0.)THEN 
        SUMPW=SUMPW+DWP 
      ELSE 
        SUMNW=SUMNW+DWP 
      ENDIF 
C 
      IF(DEP.GT.0.)THEN 
        SUMPE=SUMPE+DEP 
      ELSE 
        SUMNE=SUMNE+DEP 
      ENDIF 
C
      if(dvp.gt.0.)then
        sumpv=sumpv+dvp
      else
        sumnv=sumnv+dvp
      endif
C 
   50 CONTINUE
C----------------------------------------------------------------------- 
      DQL(1)=0. 
      DWL(1)=0. 
      DEL(1)=0. 
      dvl(1)=0.
C 
      DQL(LMP)=0. 
      DWL(LMP)=0. 
      DEL(LMP)=0. 
      dvl(lmp)=0.
C--------------FIRST MOMENT CONSERVING FACTOR--------------------------- 
      IF(SUMPQ.GT.1.E-9)THEN 
        RFACQ=-SUMNQ/SUMPQ 
      ELSE 
        RFACQ=1. 
      ENDIF 
C
      IF(SUMPW.GT.1.E-9)THEN 
        RFACW=-SUMNW/SUMPW 
      ELSE 
        RFACW=1. 
      ENDIF 
C
      IF(SUMPE.GT.1.E-9)THEN 
        RFACE=-SUMNE/SUMPE 
      ELSE 
        RFACE=1. 
      ENDIF 
C
      if(sumpv.gt.1.e-9)then
        rfacv=-sumnv/sumpv
      else
        rfacv=1.
      endif
C
      IF(RFACQ.LT.0.9.OR.RFACQ.GT.1.1) RFACQ=1. 
      IF(RFACW.LT.0.9.OR.RFACW.GT.1.1) RFACW=1. 
      IF(RFACE.LT.0.9.OR.RFACE.GT.1.1) RFACE=1. 
      if(rfacv.lt.0.9.or.rfacv.gt.1.1) rfacv=1.
C
C--------------IMPOSE CONSERVATION ON ANTI-FILTERING-------------------- 
C
      IF(RFACQ.LT.1.)THEN 
        DO L=1,LMP 
          DQP=DQL(L) 
          IF(DQP.GE.0.) DQP=DQP*RFACQ 
          Q(I,J,L)=Q4(L)+DQP 
        ENDDO 
      ELSE 
        DO L=1,LMP 
          DQP=DQL(L) 
          IF(DQP.LT.0.) DQP=DQP/RFACQ 
          Q(I,J,L)=Q4(L)+DQP 
        ENDDO 
      ENDIF 
C----------------------------------------------------------------------- 
      IF(RFACW.LT.1.)THEN 
        DO L=1,LMP 
          DWP=DWL(L) 
          IF(DWP.GE.0.) DWP=DWP*RFACW 
          CWM(I,J,L)=W4(L)+DWP 
        ENDDO 
      ELSE 
        DO L=1,LMP 
          DWP=DWL(L) 
          IF(DWP.LT.0.) DWP=DWP/RFACW 
          CWM(I,J,L)=W4(L)+DWP 
        ENDDO 
      ENDIF 
C----------------------------------------------------------------------- 
      IF(RFACE.LT.1.)THEN 
        DO L=1,LMP 
          DEP=DEL(L) 
          IF(DEP.GE.0.) DEP=DEP*RFACE 
          E3(L)=E4(L)+DEP 
        ENDDO 
      ELSE 
        DO L=1,LMP 
          DEP=DEL(L) 
          IF(DEP.LT.0.) DEP=DEP/RFACE 
          E3(L)=E4(L)+DEP 
        ENDDO 
      ENDIF 
C
      HBM2IJ=HBM2(I,J)
      Q2(I,J,1)=MAX(E3(1)+E3(1)-EPSQ2,EPSQ2)*HBM2IJ 
     1         +Q2(I,J,1)*(1.-HBM2IJ)
      DO L=2,LMP-1 
        Q2(I,J,L)=MAX(E3(L)+E3(L)-Q2(I,J,L-1),EPSQ2)*HBM2IJ
     1           +Q2(I,J,L)*(1.-HBM2IJ)
      ENDDO 
!
      if(rfac1.lt.1.)then
        do l=1,lmp
          dvp=dvl(l)
          if(dvp.ge.0.)dvp=dvp*rfacv
          dwdt(i,j,l)=dwdt(i,j,l)-(v4(l)+dwp-v3(l))*rdt
        enddo
      else
        do l=1,lmp
          dvp=dvl(l)
          if(dvp.lt.0.)dvp=dvp/rfacv
          dwdt(i,j,l)=dwdt(i,j,l)-(v4(l)+dvp-v3(l))*rdt
        ENDDO
      ENDIF
C----------------------------------------------------------------------- 
  200 CONTINUE
C--------------END OF VERTICAL ADVECTION-------------------------------- 
C*********************************************************************** 
      ENH=ADDT/(08.*DY) 
C----------------------------------------------------------------------- 
!$omp parallel do
      DO J=MYJS_P2,MYJE_P2
      DO I=MYIS_P1,MYIE_P1
        EMH (I,J)=ADDT/(08.*DX(I,J)) 
        DARE(I,J)=HBM2(I,J)*DX(I,J)*DY 
        E1(I,J,1)=MAX(Q2(I,J,1)*0.5,EPSQ2)
        E2(I,J,1)=E1(I,J,1)
      ENDDO
      ENDDO
C
C----------------------------------------------------------------------- 
C--------------HORIZONTAL ADVECTION------------------------------------- 
C----------------------------------------------------------------------- 
C
C*********************************************************************** 
!$omp parallel do
!$omp& private (destij,dqstij,dvolp,dwstij,hbm2ij,hm,
!$omp&          jfp,jfq,pp,qp,sumne,sumnq,sumnw,sumpe,sumpq,sumpw,
!$omp&          tta,ttb)
      DO 300 L=1,LM
C*********************************************************************** 
      DO J=MYJS_P2,MYJE_P2
      DO I=MYIS_P1,MYIE_P1
        dvol(i,j,l)=dare(i,j)*(deta1(l)*pdtop+deta2(l)*pdsl(i,j)) 
        HM=HTM(I,J,L) 
        Q  (I,J,L)=MAX(Q  (I,J,L),EPSQ)*HM 
        CWM(I,J,L)=MAX(CWM(I,J,L),CLIMIT)*HM 
        Q1 (I,J,L)=Q  (I,J,L) 
        W1 (I,J,L)=CWM(I,J,L) 
        v1 (i,j,l)=w  (i,j,l)*hm
      ENDDO 
      ENDDO 
C
      IF(L.GT.1)THEN 
        DO J=MYJS_P2,MYJE_P2
        DO I=MYIS_P1,MYIE_P1
          E1(I,J,L)=MAX((Q2(I,J,L-1)+Q2(I,J,L))*0.5,EPSQ2) 
          E2(I,J,L)=E1(I,J,L) 
        ENDDO 
        ENDDO 
      ENDIF 
C----------------------------------------------------------------------- 
      DO 225 J=MYJS2_P1,MYJE2_P1
      DO 225 I=MYIS1_P1,MYIE1_P1
C
      HBM2IJ=HBM2(I,J)
      DVOLP=DVOL(I,J,L)*HBM2IJ 
      TTA=(U(I,J-1,L)+U(I+IHW(J),J,L)+U(I+IHE(J),J,L)+U(I,J+1,L)) 
     2   *HBM2IJ*EMH(I,J) 
      TTB=(V(I,J-1,L)+V(I+IHW(J),J,L)+V(I+IHE(J),J,L)+V(I,J+1,L)) 
     2   *HBM2IJ*ENH 
C 
      PP=-TTA-TTB 
      QP= TTA-TTB 
C 
      JFP=INT(SIGN(1.,PP)) 
      JFQ=INT(SIGN(1.,QP)) 
C 
      IFPA(I,J,L)=IHE(J)+I+( JFP-1  )/2 
      IFQA(I,J,L)=IHE(J)+I+(-JFQ-1  )/2 
C 
      JFPA(I,J,L)=J+JFP 
      JFQA(I,J,L)=J+JFQ 
C 
      IFPF(I,J,L)=IHE(J)+I+(-JFP-1)/2 
      IFQF(I,J,L)=IHE(J)+I+( JFQ-1)/2 
C 
      JFPF(I,J,L)=J-JFP 
      JFQF(I,J,L)=J-JFQ 
C 
      PP=ABS(PP)*HTM(I,J,L)*HTM(IFPA(I,J,L),JFPA(I,J,L),L) 
      QP=ABS(QP)*HTM(I,J,L)*HTM(IFQA(I,J,L),JFQA(I,J,L),L) 
C 
      AFP (I,J,L)=(((FF4*PP+FF3)*PP+FF2)*PP+FF1)*PP 
      AFQ (I,J,L)=(((FF4*QP+FF3)*QP+FF2)*QP+FF1)*QP 
C 
      DQSTIJ=(Q  (IFPA(I,J,L),JFPA(I,J,L),L)-Q  (I,J,L))*PP 
     2      +(Q  (IFQA(I,J,L),JFQA(I,J,L),L)-Q  (I,J,L))*QP 
      DWSTIJ=(CWM(IFPA(I,J,L),JFPA(I,J,L),L)-CWM(I,J,L))*PP 
     2      +(CWM(IFQA(I,J,L),JFQA(I,J,L),L)-CWM(I,J,L))*QP 
      DESTIJ=(E1 (IFPA(I,J,L),JFPA(I,J,L),L)-E1 (I,J,L))*PP 
     2      +(E1 (IFQA(I,J,L),JFQA(I,J,L),L)-E1 (I,J,L))*QP 
      dvstij=(w  (ifpa(i,j,l),jfpa(i,j,l),l)-w  (i,j,l))*pp
     2      +(w  (ifqa(i,j,l),jfqa(i,j,l),l)-w  (i,j,l))*qp
C 
      q1(i,j,l)=q(i,j,l)+dqstij
      w1(i,j,l)=cwm(i,j,l)+dwstij
      e2(i,j,l)=e1(i,j,l)+destij
      v1(i,j,l)=w(i,j,l)+dvstij
!
  225 continue
!-----------------------------------------------------------------------
  300 CONTINUE          ! End of Vertical Loop
C--------------ANTI-FILTERING STEP-------------------------------------- 
      do 400 l=1,lm
!-----------------------------------------------------------------------
      SUMPQ=0. 
      SUMNQ=0. 
      SUMPW=0. 
      sumpv=0.
      SUMNW=0. 
      SUMPE=0. 
      SUMNE=0. 
      sumnv=0.
C--------------ANTI-FILTERING LIMITERS---------------------------------- 
      DO 350 J=MYJS2,MYJE2
      DO 350 I=MYIS1,MYIE1
C
      DVOLP=DVOL(I,J,L)
      Q1IJ =Q1(I,J,L)
      W1IJ =W1(I,J,L)
      E2IJ =E2(I,J,L)
      v1ij =v1(i,j,l)
!
      D2PQQ=((Q1(IFPA(I,J,L),JFPA(I,J,L),L)-Q1IJ               )
     2      -(Q1IJ                   -Q1(IFPF(I,J,L),JFPF(I,J,L),L))
     3      *htm(ifpf(i,j,l),jfpf(i,j,l),l))*afp(i,j,l)
     4     +((Q1(IFQA(I,J,L),JFQA(I,J,L),L)-Q1IJ               )
     5      -(Q1IJ                   -Q1(IFQF(I,J,L),JFQF(I,J,L),L))
     6      *htm(ifqf(i,j,l),jfqf(i,j,l),l))*afq(i,j,l)
C 
      D2PQW=((W1(IFPA(I,J,L),JFPA(I,J,L),L)-W1IJ               )
     2      -(W1IJ                   -W1(IFPF(I,J,L),JFPF(I,J,L),L))
     3      *htm(ifpf(i,j,l),jfpf(i,j,l),l))*afp(i,j,l)
     4     +((W1(IFQA(I,J,L),JFQA(I,J,L),L)-W1IJ               )
     5      -(W1IJ                   -W1(IFQF(I,J,L),JFQF(I,J,L),L))
     6      *htm(ifqf(i,j,l),jfqf(i,j,l),l))*afq(i,j,l)
C 
      D2PQE=((E2(IFPA(I,J,L),JFPA(I,J,L),L)-E2IJ               )
     2      -(E2IJ                   -E2(IFPF(I,J,L),JFPF(I,J,L),L))
     3      *htm(ifpf(i,j,l),jfpf(i,j,l),l))*afp(i,j,l)
     4     +((E2(IFQA(I,J,L),JFQA(I,J,L),L)-E2IJ               )
     5      -(E2IJ                   -E2(IFQF(I,J,L),JFQF(I,J,L),L))
     6      *htm(ifqf(i,j,l),jfqf(i,j,l),l))*afq(i,j,l)
C
      d2pqv=((v1(ifpa(i,j,l),jfpa(i,j,l),l)-v1ij               )
     2      -(v1ij                   -v1(ifpf(i,j,l),jfpf(i,j,l),l))
     3      *htm(ifpf(i,j,l),jfpf(i,j,l),l))*afp(i,j,l)
     4     +((v1(ifqa(i,j,l),jfqa(i,j,l),l)-v1ij               )
     5      -(v1ij                   -v1(ifqf(i,j,l),jfqf(i,j,l),l))
     6      *htm(ifqf(i,j,l),jfqf(i,j,l),l))*afq(i,j,l)
C 
      QSTIJ=Q1IJ-D2PQQ 
      WSTIJ=W1IJ-D2PQW 
      ESTIJ=E2IJ-D2PQE 
      vstij=v1ij-d2pqv
C 
      Q00=Q  (I          ,J          ,L) 
      QP0=Q  (IFPA(I,J,L),JFPA(I,J,L),L) 
      Q0Q=Q  (IFQA(I,J,L),JFQA(I,J,L),L) 
C 
      W00=CWM(I          ,J          ,L) 
      WP0=CWM(IFPA(I,J,L),JFPA(I,J,L),L) 
      W0Q=CWM(IFQA(I,J,L),JFQA(I,J,L),L) 
C 
      E00=E1 (I          ,J          ,L) 
      EP0=E1 (IFPA(I,J,L),JFPA(I,J,L),L) 
      E0Q=E1 (IFQA(I,J,L),JFQA(I,J,L),L) 
C
      w00=w  (i          ,j          ,l)
      wp0=w  (ifpa(i,j,l),jfpa(i,j,l),l)
      w0q=w  (ifqa(i,j,l),jfqa(i,j,l),l)
C 
      QSTIJ=AMAX1(QSTIJ,AMIN1(Q00,QP0,Q0Q)) 
      QSTIJ=AMIN1(QSTIJ,AMAX1(Q00,QP0,Q0Q)) 
      WSTIJ=AMAX1(WSTIJ,AMIN1(W00,WP0,W0Q)) 
      WSTIJ=AMIN1(WSTIJ,AMAX1(W00,WP0,W0Q)) 
      ESTIJ=AMAX1(ESTIJ,AMIN1(E00,EP0,E0Q)) 
      ESTIJ=AMIN1(ESTIJ,AMAX1(E00,EP0,E0Q)) 
      vstij=amax1(vstij,amin1(v00,vp0,v0q))
      vstij=amin1(vstij,amax1(v00,vp0,v0q))
C 
      dqstij=qstij-q(i,j,l) 
      dwstij=wstij-cwm(i,j,l) 
      destij=estij-e1(i,j,l)
      dvstij=vstij-w(i,j,l)
C 
      DQST(I,J,L)=DQSTIJ 
      DWST(I,J,L)=DWSTIJ 
      DEST(I,J,L)=DESTIJ 
      dvst(i,j,l)=dvstij
C 
      DQSTIJ=DQSTIJ*DVOLP 
      DWSTIJ=DWSTIJ*DVOLP 
      DESTIJ=DESTIJ*DVOLP 
      dvstij=dvstij*dvolp
C 
      IF(DQSTIJ.GT.0.)THEN
        SUMPQ=SUMPQ+DQSTIJ
      ELSE
        SUMNQ=SUMNQ+DQSTIJ
      ENDIF
C
      IF(DWSTIJ.GT.0.)THEN
        SUMPW=SUMPW+DWSTIJ
      ELSE
        SUMNW=SUMNW+DWSTIJ
      ENDIF
C
      IF(DESTIJ.GT.0.)THEN
        SUMPE=SUMPE+DESTIJ
      ELSE
        SUMNE=SUMNE+DESTIJ
      ENDIF
C
      if(dvstij.gt.0.)then
        sumpv=sumpv+dvstij
      else
        sumnv=sumnv+dvstij
      endif
C
  350 CONTINUE
C-----------------------------------------------------------------------
      XSUMS(1,L)=SUMPQ
      XSUMS(2,L)=SUMNQ
      XSUMS(3,L)=SUMPW
      XSUMS(4,L)=SUMNW
      XSUMS(5,L)=SUMPE
      XSUMS(6,L)=SUMNE
      xsums(7,l)=sumpv
      xsums(8,l)=sumnv
C
  400 CONTINUE            ! End of Vertical Loop
C-----------------------------------------------------------------------
C
C***  GLOBAL REDUCTION
C
      call mpi_allreduce(xsums,gsums,8*lm,mpi_real,mpi_sum,
     1                     mpi_comm_comp,irecv)
      IF(IRECV.NE.0)THEN
         PRINT*, ' RETURN CODE FROM 2nd ALLREDUCE IN ADV2 = ',IRECV
         STOP
      ENDIF
C
C***  END OF GLOBAL REDUCTION
C
c     if(mype.eq.0)then
c       do l=1,lm
c         read(93)gsums(1,l)
c         read(93)gsums(2,l)
c         read(93)gsums(3,l)
c         read(93)gsums(4,l)
c         read(93)gsums(5,l)
c         read(93)gsums(6,l)
c       enddo
c     endif
c     call mpi_bcast(gsums,6*lm,mpi_real,0,MPI_COMM_COMP,irtn)
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
      DO 425 L=1,LM
C
      SUMPQ=GSUMS(1,L)
      SUMNQ=GSUMS(2,L)
      SUMPW=GSUMS(3,L)
      SUMNW=GSUMS(4,L)
      SUMPE=GSUMS(5,L)
      SUMNE=GSUMS(6,L)
      sumpv=gsums(7,l)
      sumnv=gsums(8,l)
C
C--------------FIRST MOMENT CONSERVING FACTOR--------------------------- 
C
      IF(SUMPQ.GT.1.)THEN
        RFACQ=-SUMNQ/SUMPQ
      ELSE
        RFACQ=1.
      ENDIF
C
      IF(SUMPW.GT.1.)THEN
        RFACW=-SUMNW/SUMPW
      ELSE
        RFACW=1.
      ENDIF
C
      IF(SUMPE.GT.1.)THEN
        RFACE=-SUMNE/SUMPE
      ELSE
        RFACE=1.
      ENDIF
C
      if(sumpv.gt.1.)then
        rfacv=-sumnv/sumpv
      else
        rfacv=1.
      endif
C
      IF(RFACQ.LT.0.9.OR.RFACQ.GT.1.1)RFACQ=1. 
      IF(RFACW.LT.0.9.OR.RFACW.GT.1.1)RFACW=1. 
      IF(RFACE.LT.0.9.OR.RFACE.GT.1.1)RFACE=1. 
      if(rfacv.lt.0.9.or.rfacv.gt.1.1)rfacv=1.
C
C--------------IMPOSE CONSERVATION ON ANTI-FILTERING-------------------- 
C
      IF(RFACQ.LT.1.)THEN
!$omp parallel do private (dqstij,rfqij)
        DO J=MYJS2,MYJE2
        DO I=MYIS1,MYIE1
          DQSTIJ=DQST(I,J,L)
          RFQIJ=HBM2(I,J)*(RFACQ-1.)+1.
          IF(DQSTIJ.GE.0.)DQSTIJ=DQSTIJ*RFQIJ
          q(i,j,l)=q(i,j,l)+dqstij 
        ENDDO
        ENDDO
      ELSE
!$omp parallel do private (dqstij,rfqij)
        DO J=MYJS2,MYJE2
        DO I=MYIS1,MYIE1
          DQSTIJ=DQST(I,J,L)
          RFQIJ=HBM2(I,J)*(RFACQ-1.)+1.
          IF(DQSTIJ.LT.0.)DQSTIJ=DQSTIJ/RFQIJ
          q(i,j,l)=q(i,j,l)+dqstij 
        ENDDO
        ENDDO
      ENDIF
C----------------------------------------------------------------------- 
      IF(RFACW.LT.1.)THEN
!$omp parallel do private (dwstij,rfwij)
        DO J=MYJS2,MYJE2
        DO I=MYIS1,MYIE1
          DWSTIJ=DWST(I,J,L)
          RFWIJ=HBM2(I,J)*(RFACW-1.)+1.
          IF(DWSTIJ.GE.0.)DWSTIJ=DWSTIJ*RFWIJ
          cwm(i,j,l)=cwm(i,j,l)+dwstij
        ENDDO
        ENDDO
      ELSE
!$omp parallel do private (dwstij,rfwij)
        DO J=MYJS2,MYJE2
        DO I=MYIS1,MYIE1
          DWSTIJ=DWST(I,J,L)
          RFWIJ=HBM2(I,J)*(RFACW-1.)+1.
          IF(DWSTIJ.LT.0.)DWSTIJ=DWSTIJ/RFWIJ
          cwm(i,j,l)=cwm(i,j,l)+dwstij
        ENDDO
        ENDDO
      ENDIF
C----------------------------------------------------------------------- 
      IF(RFACE.LT.1.)THEN
!$omp parallel do private (destij,rfeij)
        DO J=MYJS2,MYJE2
        DO I=MYIS1,MYIE1
          DESTIJ=DEST(I,J,L)
          RFEIJ=HBM2(I,J)*(RFACE-1.)+1.
          IF(DESTIJ.GE.0.)DESTIJ=DESTIJ*RFEIJ
          e1(i,j,l)=e1(i,j,l)+destij
        ENDDO
        ENDDO
      ELSE
!$omp parallel do private (destij,rfeij)
        DO J=MYJS2,MYJE2
        DO I=MYIS1,MYIE1
          DESTIJ=DEST(I,J,L)
          RFEIJ=HBM2(I,J)*(RFACE-1.)+1.
          IF(DESTIJ.LT.0.)DESTIJ=DESTIJ/RFEIJ
          e1(i,j,l)=e1(i,j,l)+destij
        ENDDO
        ENDDO
      ENDIF
C----------------------------------------------------------------------- 
!$omp parallel do private (hm)
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        HM=HTM(I,J,L) 
        Q  (I,J,L)=AMAX1(Q  (I,J,L),EPSQ)*HM 
        CWM(I,J,L)=AMAX1(CWM(I,J,L),CLIMIT)*HM 
      ENDDO 
      ENDDO 
C***
      IF(L.EQ.LM)GO TO 425
C***
      IF(L.EQ.1)THEN 
!$omp parallel do 
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          Q2(I,J,L)=MAX(E1(I,J,L)+E1(I,J,L)-EPSQ2,EPSQ2)*HTM(I,J,1) 
        ENDDO 
        ENDDO 
      ELSE 
!$omp parallel do 
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          IF(L.LT.LMH(I,J))THEN 
            Q2(I,J,L)=MAX(E1(I,J,L)+E1(I,J,L)-Q2(I,J,L-1),EPSQ2)
     1                   *HTM(I,J,L) 
          ELSE
            Q2(I,J,L)=Q2(I,J,L-1)
          ENDIF 
        ENDDO 
        ENDDO 
      ENDIF 
!-----------------------------------------------------------------------
      if(rfacv.lt.1.)then
        do j=myjs2_p2,myje2_p2
        do i=myis1_p2,myie1_p2
          dvstij=dvst(i,j,l)
          rfvij=hbm2(i,j)*(rfacv-1.)+1.
          if(dvstij.ge.0.)dvstij=dvstij*rfvij
          dwdt(i,j,l)=(dwdt(i,j,l)-dvstij*rdt)
     1                *hbm3(i,j)
        enddo
        enddo
      else
        do j=myjs2_p2,myje2_p2
        do i=myis1_p2,myie1_p2
          dvstij=dvst(i,j,l)
          rfvij=hbm2(i,j)*(rfacv-1.)+1.
          if(dvstij.lt.0.)dvstij=dvstij/rfvij
          dwdt(i,j,l)=(dwdt(i,j,l)-dvstij*rdt)
     1                *hbm3(i,j)
        enddo
        enddo
      endif

C*********************************************************************** 
  425 CONTINUE
C----------------------------------------------------------------------- 
      if(ntsd.le.ntshy.or.hydro)  return
!-----------------------------------------------------------------------
C
C--------------RESTRICTING THE ACCELERATION ALONG THE BOUNDARIES--------
C
c     JHL=07
      JHL=0
C***
      IF(JHL.GT.1)THEN
        JHH=JM-JHL+1
C
        IHL=JHL/2+1
        IHH=IM-IHL+MOD(J,2)
C-----------------------------------------------------------------------
!$omp parallel do private (ix,jx)
        DO 450 L=1,LM
C
        DO J=1,JHL
          IF(J.GE.MY_JS_GLB-JBPAD2.AND.J.LE.MY_JE_GLB+JTPAD2)THEN
            JX=J-MY_JS_GLB+1
            DO I=1,IM
              IF(I.GE.MY_IS_GLB-ILPAD2.AND.I.LE.MY_IE_GLB+IRPAD2)THEN
                IX=I-MY_IS_GLB+1
                DWDT(IX,JX,L)=MAX(DWDT(IX,JX,L)*HBM3(IX,JX),-0.03)
                DWDT(IX,JX,L)=MIN(DWDT(IX,JX,L)*HBM3(IX,JX), 0.03)
              ENDIF
            ENDDO
          ENDIF
        ENDDO
C
        DO J=JHH,JM
          IF(J.GE.MY_JS_GLB-JBPAD2.AND.J.LE.MY_JE_GLB+JTPAD2)THEN
            JX=J-MY_JS_GLB+1
            DO I=1,IM
              IF(I.GE.MY_IS_GLB-ILPAD2.AND.I.LE.MY_IE_GLB+IRPAD2)THEN
                IX=I-MY_IS_GLB+1
                DWDT(IX,JX,L)=MAX(DWDT(IX,JX,L)*HBM3(IX,JX),-0.03)
                DWDT(IX,JX,L)=MIN(DWDT(IX,JX,L)*HBM3(IX,JX), 0.03)
              ENDIF
            ENDDO
          ENDIF
        ENDDO
C
        DO J=1,JM
          IF(J.GE.MY_JS_GLB-JBPAD2.AND.J.LE.MY_JE_GLB+JTPAD2)THEN
            JX=J-MY_JS_GLB+1
            DO I=1,IHL
              IF(I.GE.MY_IS_GLB-ILPAD2.AND.I.LE.MY_IE_GLB+IRPAD2)THEN
                IX=I-MY_IS_GLB+1
                DWDT(IX,JX,L)=MAX(DWDT(IX,JX,L)*HBM3(IX,JX),-0.03)
                DWDT(IX,JX,L)=MIN(DWDT(IX,JX,L)*HBM3(IX,JX), 0.03)
              ENDIF
            ENDDO
          ENDIF
        ENDDO
C
        DO J=1,JM
          IF(J.GE.MY_JS_GLB-JBPAD2.AND.J.LE.MY_JE_GLB+JTPAD2)THEN
            JX=J-MY_JS_GLB+1
            DO I=IHH,IM
              IF(I.GE.MY_IS_GLB-ILPAD2.AND.I.LE.MY_IE_GLB+IRPAD2)THEN
                IX=I-MY_IS_GLB+1
                DWDT(IX,JX,L)=MAX(DWDT(IX,JX,L)*HBM3(IX,JX),-0.03)
                DWDT(IX,JX,L)=MIN(DWDT(IX,JX,L)*HBM3(IX,JX), 0.03)
              ENDIF
            ENDDO
          ENDIF
        ENDDO
C
  450   CONTINUE
      ENDIF
C------------------------------------------------------------------------
      WA=0.075
C
!$omp parallel do private (ane,ase)
      DO 500 L=1,LM
C
      DO J=MYJS1_P1,MYJE2_P1
      DO I=MYIS_P1,MYIE1_P1
        ANE(I,J)=(DWDT(I+IHE(J),J+1,L)-DWDT(I,J,L))
     &           *HTM(I,J,L)*HTM(I+IHE(J),J+1,L)
      ENDDO
      ENDDO
C
      DO J=MYJS2_P1,MYJE1_P1
      DO I=MYIS_P1,MYIE1_P1
        ASE(I,J)=(DWDT(I+IHE(J),J-1,L)-DWDT(I,J,L))
     &           *HTM(I+IHE(J),J-1,L)*HTM(I,J,L)
      ENDDO
      ENDDO
C
      DO J=MYJS2,MYJE2
      DO I=MYIS1,MYIE1
        DWDT(I,J,L)=(ANE(I,J)-ANE(I+IHW(J),J-1)
     1              +ASE(I,J)-ASE(I+IHW(J),J+1))*WA*HBM2(I,J)
     2              +DWDT(I,J,L)
      ENDDO
      ENDDO
C-----------------------------------------------------------------------
  500 CONTINUE
C-----------------------------------------------------------------------
      WP=0.075
C
      KN=0
      KP=0
      DWDTMX=0.
      DWDTMN=0.
C-----------------------------------------------------------------------
!$omp parallel do
!$omp& private (dwdtmn,dwdtmx,dwdtt,hm,kn,kp)
      DO 525 L=1,LM
C-----------------------------------------------------------------------
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        HM=HTM(I,J,L)*HBM2(I,J)
        DWDTT=DWDT(I,J,L)*HM
C
        DWDTMX=AMAX1(DWDTT,DWDTMX)
        DWDTMN=AMIN1(DWDTT,DWDTMN)
C
        IF(DWDTT.LT.EPSN)THEN
          DWDTT=EPSN
          KN=KN+1
        ENDIF
C
        IF(DWDTT.GT.EPSP)THEN
          DWDTT=EPSP
          KP=KP+1
        ENDIF
C
        DWDT(I,J,L)=(DWDTT/G+1.)*(1.-WP)+PDWDT(I,J,L)*WP
      ENDDO
      ENDDO
C-----------------------------------------------------------------------
  525 CONTINUE
C-----------------------------------------------------------------------
      GDT=G*DT
      GDT2=GDT*GDT
      WGHT=0.35
      FCC=GMA/GDT2
C-----------------------------------------------------------------------
!$omp parallel do
!$omp& private (c0,b1,b2,b3,chi,coff,dpstr,dp,tfc,fcc,dfrc,rdpdn,
!$omp&          pnp1,pone,pp1,pstr,rdpdup,tmp,dptu,dptl,dwdtt,wil,delp,
!$omp&          hbm2ij)
C-----------------------------------------------------------------------
      DO 600 J=MYJS2,MYJE2
      DO 600 I=MYIS1,MYIE1
C-----------------------------------------------------------------------
      LMP=LMH(I,J)
      PDP=PDSL(I,J)
!-----------------------------------------------------------------------
      pone(1)=pt
      pstr(1)=pt
      pnp1(1)=pt
      chi(1)=0.
!
          do l=2,lmp+1
      chi(l)=0.
      pone(l)=pint(i,j,l)
      dpstr=dwdt(i,j,l-1)*(deta1(l-1)*pdtop+deta2(l-1)*pdp)
      pstr(l)=pstr(l-1)+dpstr
      pp1=pnp1(l-1)+dpstr
      dp=(pp1-pone(l))*wght
      pnp1(l)=pone(l)+dp
          enddo
!
          do l=2,lmp+1
      tfc=q(i,j,l-1)*0.608-cwm(i,j,l-1)+1.
      fcc=-r*tfc*(1.-cappa*tfc)*0.5/gdt2
      coff(l-1)=t(i,j,l-1)*(deta1(l-1)*pdtop+deta2(l-1)*pdp)*fcc
     &         /((pnp1(l-1)+pnp1(l))*0.5)**2
          enddo
!
          do l=2,lmp
      dfrc=((pstr(l-1)+pstr(l)-pone(l-1)-pone(l))*coff(l-1)
     &     +(pstr(l)+pstr(l+1)-pone(l)-pone(l+1))*coff(l))*0.5
!
      rdpdn=1./(deta1(l)*pdtop+deta2(l)*pdp)
      rdpup=1./(deta1(l-1)*pdtop+deta2(l-1)*pdp)
!
      b1(l)=(coff(l-1)*0.5+rdpup)
      b2(l)=(coff(l-1)+coff(l))*0.5-(rdpup+rdpdn)
!zj      b2(l)=min(b2(l),-0.001)
      b3(l)=(coff(l)*0.5+rdpdn)
      c0(l)=-dfrc
          enddo
!
      b1(2)=0.
      b2(lmp)=b2(lmp)+b3(lmp)
!---------------elimination---------------------------------------------
          do l=3,lmp
      tmp=-b1(l)/b2(l-1)
      b2(l)=b3(l-1)*tmp+b2(l)
      c0(l)=c0(l-1)*tmp+c0(l)
          enddo
!
      chi(1)=0.
!---------------back substitution--------------------------------------
      chi(lmp)=c0(lmp)/b2(lmp)
      chi(lmp+1)=chi(lmp)
!
          do l=lmp-1,2,-1
      chi(l)=(-b3(l)*chi(l+1)+c0(l))/b2(l)
          enddo
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      DO L=1,LMP+1
        if(abs(chi(l)).gt.4000.)print*,ntsd,i,j,l,chi(l),fis(i,j)
        PNP1(L)=CHI(L)+PSTR(L)
        PINT(I,J,L)=PNP1(L)
      ENDDO
C
      DPTU=0.
C
      HBM2IJ=HBM2(I,J)
C
      DO L=1,LMP
        DPTL=PNP1(L+1)-PONE(L+1)
        T(I,J,L)=(DPTU+DPTL)*RTOP(I,J,L)*0.5/CP*HBM2IJ+T(I,J,L)
        delp=(pnp1(l+1)-pnp1(l))/(deta1(l)*pdtop+deta2(l)*pdp)
!zj        w(i,j,l)=(w(i,j,l)+(delp-dwdt(i,j,l))*gdt)*hbm2ij
!zj        dwdt(i,j,l)=delp
        DWDTT=DWDT(I,J,L)
        WIL=W(I,J,L)+(DELP-DWDTT)*GDT
        DWDT(I,J,L)=DWDTT+(WIL-W(I,J,L))/GDT*HBM2IJ
        W(I,J,L)=WIL*HBM2IJ
C
        DPTU=DPTL
      ENDDO
C-----------------------------------------------------------------------
  600 CONTINUE
C-----------------------------------------------------------------------
                             RETURN 
                             END 
