!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    SUBROUTINE EPS
!     ******************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    EPS
!   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 9?-??-??

! ABSTRACT:
!     EPS COMPUTES THE VERTICAL AND HORIZONTAL ADVECTION OF DZ/DT

! PROGRAM HISTORY LOG:
!   9?-??-??  JANJIC     - ORIGINATOR
!   00-01-05  BLACK      - DISTRIBUTED MEMORY AND THREADS
!   26-11-12  RISTIC IVAN - LOT OF CHANGES AND BUG FIX

! USAGE: CALL EPS FROM MAIN PROGRAM
!   INPUT ARGUMENT LIST:
!       NONE

!   OUTPUT ARGUMENT LIST:
!     NONE

!   OUTPUT FILES:
!     NONE

!   SUBPROGRAMS CALLED:

!     UNIQUE: NONE

!     LIBRARY: NONE

!   COMMON BLOCKS: CTLBLK
!                  LOOPS
!                  MASKS
!                  DYNAM
!                  CONTIN
!                  VRBLS
!                  PVRBLS
!                  NHYDRO
!                  INDX
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$
!***********************************************************************
    INCLUDE "EXCHM.h"
    INCLUDE "parmeta.f90"
    INCLUDE "mpp.h"
    INCLUDE "mpif.h"
!-----------------------------------------------------------------------
    PARAMETER &
    (IM1=IM-1,JAM=6+2*(JM-10) &
    ,JAMD=(JAM*2-10)*3,LP1=LM+1)
!-----------------------------------------------------------------------
    PARAMETER (EPSQ=1.E-12,EPSQ2=0.2)
    PARAMETER (FF1=0.52500,FF2=-0.64813,FF3=0.24520,FF4=-0.12189)
!-----------------------------------------------------------------------
    PARAMETER &
    (NTSHY=3,EPSFC=1./1.05,EPSN=-EPSFC,EPSP=EPSFC,ZERO=1.E-06 &
    ,G=9.8,CP=1004.6,CAPA=287.04/CP,GMA=-287.04*(1.-CAPA)/2.)
!-----------------------------------------------------------------------
    LOGICAL :: &
    RUN,FIRST,RESTRT,SIGMA,TOP,BOT
!-----------------------------------------------------------------------
    INCLUDE "COMM_CTLBLK.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_LOOPS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_MASKS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_DYNAM.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_CONTIN.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_VRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_PVRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_NHYDRO.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_INDX.f90"
!-----------------------------------------------------------------------
!-------------LOCAL VARIABLES REQUIRING DECLARATION---------------------
!-----------------------------------------------------------------------
    REAL :: &
    PONE(LM+1),PSTR(LM+1),PNP1(LM+1),COFF(LM+1),DFRC(LM+1) &
    ,CHI(LM+1),CHIN(LM+1),WRES(LM+1),CLIM(LM+1)
    REAL :: &
    RDPP(LM),C0(LM),B1(LM),B2(LM),B3(LM)
!-----------------------------------------------------------------------
    REAL :: W3(LM),W4(LM)
    REAL :: ETADTL(LM),DWL(LM),AFR(LM)
    INTEGER :: LA(LM)
!-----------------------------------------------------------------------
    REAL :: &
    AFP (IDIM1:IDIM2,JDIM1:JDIM2,LM),AFQ (IDIM1:IDIM2,JDIM1:JDIM2,LM) &
    ,W1  (IDIM1:IDIM2,JDIM1:JDIM2,LM),DWST(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)
!-----------------------------------------------------------------------
    INTEGER :: &
    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)
!-----------------------------------------------------------------------
    REAL :: &
    GSUMS(2,LM),XSUMS(2,LM),SUMPW,SUMNW
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------------------------------------------------------------
    IF(NTSD <= NTSHY .OR. HYDRO)THEN
    !***
    ! omp parallel do
        DO J=MYJS_P2,MYJE_P2
            DO I=MYIS_P1,MYIE_P1
                PINT(I,J,1)=PT
            ENDDO
        ENDDO
    
        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)=PDSL(I,J)*DETA(L)+PINT(I,J,L)
                ENDDO
            ENDDO
        ENDDO
    !***
        RETURN
    !***
    ENDIF
!-----------------------------------------------------------------------
    ADDT=DT
    RDT=1./ADDT
!-----------------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

!-----------------------------------------------------------------------
!-----------------VERTICAL ADVECTION------------------------------------
!-----------------------------------------------------------------------

! omp parallel do
! omp& private (afr,afrp,arr,bot,d2pqw,detal,dwl,dwp,etadtl,
! omp&          la,lap,llap,lmp,rfacw,rr,sumnw,sumpw,top,
! omp&          w00,w3,w4,w4p,wp,wp0)
    DO 200 J=MYJS_P1,MYJE_P1
        DO 200 I=MYIS_P1,MYIE_P1
        !-----------------------------------------------------------------------
            LMP=LMH(I,J)
        !-----------------------------------------------------------------------
            DO L=1,LMP
                W3(L)=W(I,J,L)
                W4(L)=W3(L)
            ENDDO
        
            ETADTL(1)=ETADT(I,J,1)*0.5
        
            DO L=2,LMP-1
                ETADTL(L)=(ETADT(I,J,L-1)+ETADT(I,J,L))*0.5
            ENDDO
        
            ETADTL(LMP)=ETADT(I,J,LMP-1)*0.5
        !-----------------------------------------------------------------------
            SUMPW=0.
            SUMNW=0.
        
            DO L=1,LMP
                RR=ETADTL(L)*(-ADDT)
                ARR=ABS(RR)
            
                IF(ARR > 0.)THEN
                    LAP=RR/ARR
                ELSE
                    LAP=0
                ENDIF
            
                LA(L)=LAP
                LLAP=L+LAP
            
                TOP= .FALSE. 
                BOT= .FALSE. 
            
                IF(LLAP > 0 .AND. LLAP < LMH(I,J)+1 .AND. LAP /= 0)THEN
                    RR=ARR/ABS(AETA(LLAP)-AETA(L))
                    AFR(L)=(((FF4*RR+FF3)*RR+FF2)*RR+FF1)*RR
                    DWP=(W3(LLAP)-W3(L))*RR
                    DWL(L)=DWP
                ELSE
                    TOP=LLAP == 0
                    BOT=LLAP == LMH(I,J)+1
                
                    RR=0.
                    AFR(L)=0.
                    DWL(L)=0.
                ENDIF
            ENDDO
        
            IF(TOP)THEN
                IF(LA(2) < 0)THEN
                    DWL(1)=-DWL(2)*DETA(2)/DETA(1)
                ENDIF
            ENDIF
        
            IF(BOT)THEN
                IF(LA(LMP-1) > 0)THEN
                    DWL(LMP)=-DWL(LMP-1)*DETA(LMP-1)/DETA(LMP)
                ENDIF
            ENDIF
        
            DO L=1,LMP
                DETAL=DETA(L)
                DWP=DWL(L)*DETAL
                IF(DWP > 0.)THEN
                    SUMPW=SUMPW+DWP
                ELSE
                    SUMNW=SUMNW+DWP
                ENDIF
            ENDDO
        !--------------FIRST MOMENT CONSERVING FACTOR---------------------------
            IF(SUMPW > 1.E-9)THEN
                RFACW=-SUMNW/SUMPW
            ELSE
                RFACW=1.
            ENDIf
        
            IF(RFACW < 0.9 .OR. RFACW > 1.1)RFACW=1.
        !--------------IMPOSE CONSERVATION ON ADVECTION-------------------------
            IF(RFACW < 1.)THEN
                DO L=1,LMP
                    DWP=DWL(L)
                    IF(DWP < 0.)DWP=DWP/RFACW
                    W4(L)=W3(L)+DWP
                ENDDO
            ELSE
                DO L=1,LMP
                    DWP=DWL(L)
                    IF(DWP >= 0.)DWP=DWP*RFACW
                    W4(L)=W3(L)+DWP
                ENDDO
            ENDIF
        !--------------ANTI-FILTERING STEP--------------------------------------
            SUMPW=0.
            SUMNW=0.
        !--------------ANTI-FILTERING LIMITERS----------------------------------
            DO 50 L=2,LMP-1
                DETAL=DETA(L)
            
                W4P=W4(L)
            
                LAP=LA(L)
            
                IF(LAP /= 0)THEN
                    AFRP=2.*AFR(L)*(AETA(L+LAP)-AETA(L))**2 &
                    /(AETA(L+LAP)-AETA(L-LAP))
                    D2PQW=((W4(L+LAP)-W4P)/(AETA(L+LAP)-AETA(L)) &
                    -(W4P-W4(L-LAP))/(AETA(L)-AETA(L-LAP)))*AFRP
                ELSE
                    D2PQW=0.
                ENDIF
            
                WP=W4P-D2PQW
            
                W00=W3(L)
                WP0=W3(L+LAP)
            
                WP=MAX(WP,MIN(W00,WP0))
                WP=MIN(WP,MAX(W00,WP0))
            
                DWP=WP-W4P
            
                DWL(L)=DWP
            
                DWP=DWP*DETAL
            
                IF(DWP > 0.)THEN
                    SUMPW=SUMPW+DWP
                ELSE
                    SUMNW=SUMNW+DWP
                ENDIF
            
            50 END DO
        !-----------------------------------------------------------------------
            DWL(1)=0.
        
            DWL(LMP)=0.
        !--------------FIRST MOMENT CONSERVING FACTOR---------------------------
            IF(SUMPW > 1.E-9)THEN
                RFACW=-SUMNW/SUMPW
            ELSE
                RFACW=1.
            ENDIF
        
            IF(RFACW < 0.9 .OR. RFACW > 1.1)RFACW=1.
        !--------------IMPOSE CONSERVATION ON ANTI-FILTERING--------------------
            IF(RFACW < 1.)THEN
                DO L=1,LMP
                    DWP=DWL(L)
                    IF(DWP >= 0.)DWP=DWP*RFACW
                    DWDT(I,J,L)=DWDT(I,J,L)-(W4(L)+DWP-W3(L))*RDT
                ENDDO
            ELSE
                DO L=1,LMP
                    DWP=DWL(L)
                    IF(DWP < 0.)DWP=DWP/RFACW
                    DWDT(I,J,L)=DWDT(I,J,L)-(W4(L)+DWP-W3(L))*RDT
                ENDDO
            ENDIF
        !-----------------------------------------------------------------------
    200 END DO
!--------------END OF VERTICAL ADVECTION--------------------------------

!-----------------------------------------------------------------------
!***********************************************************************
    ENH=ADDT/(08.*DY)
!-----------------------------------------------------------------------
! 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
        ENDDO
    ENDDO
!-----------------------------------------------------------------------
!---------------------HORIZONTAL ADVECTION------------------------------
!-----------------------------------------------------------------------
!***********************************************************************
! omp parallel do
! omp& private (dvolp,dwstij,hm,jfp,jfq,pp,qp,
! omp&          sumnw,sumpw,tta,ttb)
    DO 300 L=1,LM
    !***********************************************************************
        DO J=MYJS_P2,MYJE_P2
            DO I=MYIS_P1,MYIE_P1
                DVOL(I,J,L)=DARE(I,J)*PDSL(I,J)*DETA(L)
                HM=HTM(I,J,L)
                W1(I,J,L)=W(I,J,L)*HM
            ENDDO
        ENDDO
    !-----------------------------------------------------------------------
        SUMPW=0.
        SUMNW=0.
    
        DO 225 J=MYJS2_P1,MYJE2_P1
            DO 225 I=MYIS1_P1,MYIE1_P1
            
                DVOLP=DVOL(I,J,L)*HBM3(I,J)
                TTA=(U(I,J-1,L)+U(I+IHW(J),J,L)+U(I+IHE(J),J,L)+U(I,J+1,L)) &
                *HBM2(I,J)*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)) &
                *HBM2(I,J)*ENH
            
                PP=-TTA-TTB
                QP= TTA-TTB
            
                JFP=INT(SIGN(1.,PP))
                JFQ=INT(SIGN(1.,QP))
            
                IFPA(I,J,L)=IHE(J)+I+( JFP-1)/2
                IFQA(I,J,L)=IHE(J)+I+(-JFQ-1)/2
            
                JFPA(I,J,L)=J+JFP
                JFQA(I,J,L)=J+JFQ
            
                IFPF(I,J,L)=IHE(J)+I+(-JFP-1)/2
                IFQF(I,J,L)=IHE(J)+I+( JFQ-1)/2
            
                JFPF(I,J,L)=J-JFP
                JFQF(I,J,L)=J-JFQ
            
                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)
            
                AFP(I,J,L)=(((FF4*PP+FF3)*PP+FF2)*PP+FF1)*PP
                AFQ(I,J,L)=(((FF4*QP+FF3)*QP+FF2)*QP+FF1)*QP
            
                DWSTIJ=(W(IFPA(I,J,L),JFPA(I,J,L),L)-W(I,J,L))*PP &
                +(W(IFQA(I,J,L),JFQA(I,J,L),L)-W(I,J,L))*QP
            
                DWST(I,J,L)=DWSTIJ
            
        225 END DO
    !***
    !***  GLOBAL SUM FOR CONSERVATION
    !***
        DO 230 J=MYJS2,MYJE2
            DO 230 I=MYIS1,MYIE1
            
                DVOLP=DVOL(I,J,L)*HBM3(I,J)
                DWSTIJ=DWST(I,J,L)*DVOLP
            
                IF(DWSTIJ > 0.)THEN
                    SUMPW=SUMPW+DWSTIJ
                ELSE
                    SUMNW=SUMNW+DWSTIJ
                ENDIF
            
        230 END DO
    
    !-----------------------------------------------------------------------
        XSUMS(1,L)=SUMPW
        XSUMS(2,L)=SUMNW
    
        CONTINUE          ! End of Vertical Loop
    300 END DO
!-----------------------------------------------------------------------

!***  GLOBAL REDUCTION

    CALL MPI_ALLREDUCE(XSUMS,GSUMS,2*LM,MPI_REAL,MPI_SUM, &
    MPI_COMM_COMP,IRECV)
    IF(IRECV /= 0)THEN
        PRINT*, ' RETURN CODE FROM 1st ALLREDUCE IN EPS = ',IRECV
        STOP
    ENDIF
!-----------------------------------------------------------------------
! omp  parallel do
! omp& private(d2pqw,dvolp,dwstij,rfacw,rfwij,sumnw,sumpw,
! omp&         w00,w0q,w1ij,wp0,wstij)
!-----------------------------------------------------------------------
    DO 400 L=1,LM
    !-----------------------------------------------------------------------
    
        SUMPW=GSUMS(1,L)
        SUMNW=GSUMS(2,L)
    
    !--------------FIRST MOMENT CONSERVING FACTOR---------------------------
        IF(SUMPW > 1.)THEN
            RFACW=-SUMNW/SUMPW
        ELSE
            RFACW=1.
        ENDIF
        IF(RFACW < 0.9 .OR. RFACW > 1.1)RFACW=1.
    !--------------IMPOSE CONSERVATION ON ADVECTION-------------------------
        IF(RFACW < 1.)THEN
            DO J=MYJS2_P1,MYJE2_P1
                DO I=MYIS1_P1,MYIE1_P1
                    DWSTIJ=DWST(I,J,L)
                    RFWIJ=HBM3(I,J)*(RFACW-1.)+1.
                    IF(DWSTIJ < 0.)DWSTIJ=DWSTIJ/RFWIJ
                    W1(I,J,L)=W(I,J,L)+DWSTIJ
                ENDDO
            ENDDO
        ELSE
            DO J=MYJS2_P1,MYJE2_P1
                DO I=MYIS1_P1,MYIE1_P1
                    DWSTIJ=DWST(I,J,L)
                    RFWIJ=HBM3(I,J)*(RFACW-1.)+1.
                    IF(DWSTIJ >= 0.)DWSTIJ=DWSTIJ*RFWIJ
                    W1(I,J,L)=W(I,J,L)+DWSTIJ
                ENDDO
            ENDDO
        ENDIF
    !--------------ANTI-FILTERING STEP--------------------------------------
        SUMPW=0.
        SUMNW=0.
    !--------------ANTI-FILTERING LIMITERS----------------------------------
    
        DO 350 J=MYJS2,MYJE2
            DO 350 I=MYIS1,MYIE1
            
                DVOLP=DVOL(I,J,L)
                W1IJ =W1(I,J,L)
            
                D2PQW=((W1(IFPA(I,J,L),JFPA(I,J,L),L)-W1IJ               ) &
                -(W1IJ                   -W1(IFPF(I,J,L),JFPF(I,J,L),L)) &
                *HTM(IFPF(I,J,L),JFPF(I,J,L),L))*AFP(I,J,L) &
                +((W1(IFQA(I,J,L),JFQA(I,J,L),L)-W1IJ               ) &
                -(W1IJ                   -W1(IFQF(I,J,L),JFQF(I,J,L),L)) &
                *HTM(IFQF(I,J,L),JFQF(I,J,L),L))*AFQ(I,J,L)
            
                WSTIJ=W1IJ-D2PQW
            
                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)
            
                WSTIJ=AMAX1(WSTIJ,AMIN1(W00,WP0,W0Q))
                WSTIJ=AMIN1(WSTIJ,AMAX1(W00,WP0,W0Q))
            
                DWSTIJ=WSTIJ-W1IJ
            
                DWST(I,J,L)=DWSTIJ
            
                DWSTIJ=DWSTIJ*DVOLP
            
                IF(DWSTIJ > 0.)THEN
                    SUMPW=SUMPW+DWSTIJ
                ELSE
                    SUMNW=SUMNW+DWSTIJ
                ENDIF
            
        350 END DO
    !-----------------------------------------------------------------------
        XSUMS(1,L)=SUMPW
        XSUMS(2,L)=SUMNW
    
        CONTINUE            ! End of Vertical Loop
    400 END DO
!-----------------------------------------------------------------------

!***  GLOBAL REDUCTION

    CALL MPI_ALLREDUCE(XSUMS,GSUMS,2*LM,MPI_REAL,MPI_SUM, &
    MPI_COMM_COMP,IRECV)
    IF(IRECV /= 0)THEN
        PRINT*, ' RETURN CODE FROM 2nd ALLREDUCE IN EPS = ',IRECV
        STOP
    ENDIF
!-----------------------------------------------------------------------

! omp  parallel do
! omp& private(dwstij,rfacw,rfwij,sumnw,sumpw)
!-----------------------------------------------------------------------
    DO 425 L=1,LM
    
        SUMPW=GSUMS(1,L)
        SUMNW=GSUMS(2,L)
    
    !--------------FIRST MOMENT CONSERVING FACTOR---------------------------
        IF(SUMPW > 1.)THEN
            RFACW=-SUMNW/SUMPW
        ELSE
            RFACW=1.
        ENDIF
        IF(RFACW < 0.9 .OR. RFACW > 1.1)RFACW=1.
    !--------------IMPOSE CONSERVATION ON ANTI-FILTERING--------------------
        IF(RFACW < 1.)THEN
            DO J=MYJS2,MYJE2
                DO I=MYIS1,MYIE1
                    DWSTIJ=DWST(I,J,L)
                    RFWIJ=HBM2(I,J)*(RFACW-1.)+1.
                    IF(DWSTIJ >= 0.)DWSTIJ=DWSTIJ*RFWIJ
                    DWDT(I,J,L)=(DWDT(I,J,L)-(W1(I,J,L)+DWSTIJ-W(I,J,L))*RDT) &
                    *HBM3(I,J)
                ENDDO
            ENDDO
        ELSE
            DO J=MYJS2,MYJE2
                DO I=MYIS1,MYIE1
                    DWSTIJ=DWST(I,J,L)
                    RFWIJ=HBM2(I,J)*(RFACW-1.)+1.
                    IF(DWSTIJ < 0.)DWSTIJ=DWSTIJ/RFWIJ
                    DWDT(I,J,L)=(DWDT(I,J,L)-(W1(I,J,L)+DWSTIJ-W(I,J,L))*RDT) &
                    *HBM3(I,J)
                ENDDO
            ENDDO
        ENDIF
    !-----------------------------------------------------------------------
    !***********************************************************************
    425 END DO
!***********************************************************************

!--------------RESTRICTING THE ACCELERATION ALONG THE BOUNDARIES--------

!     JHL=07
    JHL=0
!***
    IF(JHL > 1)THEN
        JHH=JM-JHL+1
    
        IHL=JHL/2+1
        IHH=IM-IHL+MOD(J,2)
    !-----------------------------------------------------------------------
    ! omp parallel do private (ix,jx)
        DO 450 L=1,LM
        
            DO J=1,JHL
                IF(J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2)THEN
                    JX=J-MY_JS_GLB+1
                    DO I=1,IM
                        IF(I >= MY_IS_GLB-ILPAD2 .AND. I <= 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
        
            DO J=JHH,JM
                IF(J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2)THEN
                    JX=J-MY_JS_GLB+1
                    DO I=1,IM
                        IF(I >= MY_IS_GLB-ILPAD2 .AND. I <= 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
        
            DO J=1,JM
                IF(J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2)THEN
                    JX=J-MY_JS_GLB+1
                    DO I=1,IHL
                        IF(I >= MY_IS_GLB-ILPAD2 .AND. I <= 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
        
            DO J=1,JM
                IF(J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2)THEN
                    JX=J-MY_JS_GLB+1
                    DO I=IHH,IM
                        IF(I >= MY_IS_GLB-ILPAD2 .AND. I <= 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
        
        450 END DO
    ENDIF
!------------------------------------------------------------------------
    CALL EXCH(DWDT,LM,1,1)
    WA=0.15

! omp parallel do private (ane,ase)
    DO 500 L=1,LM
    
        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
    
        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
    
        DO J=MYJS2,MYJE2
            DO I=MYIS1,MYIE1
                DWDT(I,J,L)=(ANE(I,J)-ANE(I+IHW(J),J-1) &
                +ASE(I,J)-ASE(I+IHW(J),J+1))*WA*HBM2(I,J) &
                +DWDT(I,J,L)
            ENDDO
        ENDDO
    !-----------------------------------------------------------------------
    500 END DO
!-----------------------------------------------------------------------
    WP=0.075

    KN=0
    KP=0
    DWDTMX=0.
    DWDTMN=0.
!-----------------------------------------------------------------------
! omp parallel do
! omp& private (dwdtmn,dwdtmx,dwdtt,hm,kn,kp)
    DO 525 L=1,LM
    !-----------------------------------------------------------------------
        DO J=MYJS,MYJE
            DO I=MYIS,MYIE
                HM=HTM(I,J,L)*HBM2(I,J)
                DWDTT=DWDT(I,J,L)*HM
            
                DWDTMX=AMAX1(DWDTT,DWDTMX)
                DWDTMN=AMIN1(DWDTT,DWDTMN)
            
                IF(DWDTT < EPSN)THEN
                    DWDTT=EPSN
                    KN=KN+1
                ENDIF
            
                IF(DWDTT > EPSP)THEN
                    DWDTT=EPSP
                    KP=KP+1
                ENDIF
            
                DWDT(I,J,L)=(DWDTT/G+1.)*(1.-WP)+PDWDT(I,J,L)*WP
            ENDDO
        ENDDO
    !-----------------------------------------------------------------------
    525 END DO
!-----------------------------------------------------------------------
    GDT=G*DT
    GDT2=GDT*GDT
    WGHT=0.35
    FFC=-4.*R/GDT2

!-----------------------------------------------------------------------
! omp parallel do
! omp& private (chi,chin,chmod,clim,coff,delp,dfrc,dp,dpstr,
! omp&          dptl,dptu,dwdtt,hbm2ij,imax,inot,iter,itmx,jmax,
! omp&          l,lmp,pdp,pnp1,pone,pp1,pstr,rdp,resdl,wil,wres)
!-----------------------------------------------------------------------
    DO 600 J=MYJS2,MYJE2
        DO 600 I=MYIS1,MYIE1
        !-----------------------------------------------------------------------
        !-----------------------------------------------------------------------
        
            LMP=LMH(I,J)
            PDP=PDSL(I,J)
            RPD=1./PDP
        !-----------------------------------------------------------------------
            PONE(1)=PT
            PSTR(1)=PT
            PNP1(1)=PT
            CHI(1)=0.
        
            DO L=2,LMP+1
                CHI(L)=0.
                DPPL=DETA(L-1)*PDP
                RDPP(L-1)=1./DPPL
                PONE(L)=PINT(I,J,L)
                DPSTR=DWDT(I,J,L-1)*DPPL
                PSTR(L)=PSTR(L-1)+DPSTR
                PP1=PNP1(L-1)+DPSTR
                PNP1(L)=(PP1-PONE(L))*WGHT+PONE(L)
                TFC=Q(I,J,L-1)*0.608+1.
                TTFC=-CAPA*TFC+1.
                COFF(L-1)=T(I,J,L-1)*TTFC*TFC*DPPL*FFC &
                /((PNP1(L-1)+PNP1(L))*(PNP1(L-1)+PNP1(L)))
            ENDDO
        
            PSTRUP=-(PSTR(1)+PSTR(2)-PONE(1)-PONE(2))*COFF(1)
        
            DO L=2,LMP
                RDPDN=RDPP(L)
                RDPUP=RDPP(L-1)
            
                PSTRDN=-(PSTR(L)+PSTR(L+1)-PONE(L)-PONE(L+1))*COFF(L)
            
                B1(L)=COFF(L-1)+RDPUP
                B2(L)=(COFF(L-1)+COFF(L))-(RDPUP+RDPDN)
                B3(L)=COFF(L)+RDPDN
                C0(L)=PSTRUP+PSTRDN
            
                PSTRUP=PSTRDN
            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
        !-----------------------------------------------------------------------
            HBM2IJ=HBM2(I,J)
            DPTU=0.
            FCT=0.5/CP*HBM2IJ
        
            DO L=1,LMP
                DPTL=(CHI(L+1)+PSTR(L+1)-PINT(I,J,L+1))*HBM2IJ
                PINT(I,J,L+1)=PINT(I,J,L+1)+DPTL
                T(I,J,L)=(DPTU+DPTL)*RTOP(I,J,L)*FCT+T(I,J,L)
                DELP=(PINT(I,J,L+1)-PINT(I,J,L))*RDPP(L)
                W(I,J,L)=((DELP-DWDT(I,J,L))*GDT+W(I,J,L))*HBM2IJ
                DWDT(I,J,L)=(DELP-1.)*HBM2IJ+1.
            
                DPTU=DPTL
            ENDDO
        
            DO L=LMP,LM
                PINT(I,J,L+1)=PDSL(I,J)*DETA(L)+PINT(I,J,L)
            ENDDO
        !-----------------------------------------------------------------------
    600 END DO
!-----------------------------------------------------------------------
!	write(6,*) 'for PE: ', MYPE, 'DWDT extrema: ', DWDTMN,DWDTMX
    RETURN
    END SUBROUTINE EPS
