!--------------------------------------------------------------------------------------------------
!  DOXYGEN
!> @brief COMPUTES THE VERTICAL AND HORIZONTAL ADVECTION OF DZ/DT
!> @details COMPUTES THE VERTICAL AND HORIZONTAL ADVECTION OF DZ/DT
!> @author ORIGINATOR - JANJIC 
!> @date ??-??-?? \n
!> @author LUCCI 
!> @date 18-03-20 \n
!> @version V1.1.0
!> @details MODERNIZATION OF THE CODE, INCLUDING:
!!                      * F77 TO F90/F95
!!                      * INDENTATION & UNIFORMIZATION CODE
!!                      * REPLACEMENT OF COMMONS BLOCK FOR MODULES
!!                      * DOCUMENTATION WITH DOXYGEN
!!                      * OPENMP FUNCTIONALITY 
!<
!> @details <b>Use Module:</b>
!! @arg @c CONTIN
!! @arg @c CTLBLK
!! @arg @c DYNAM
!! @arg @c F77KINDS
!! @arg @c LOOPS
!! @arg @c MASKS
!! @arg @c METRCS
!! @arg @c MPPSTAFF
!! @arg @c  NHYDRO
!! @arg @c PARMETA
!! @arg @c PVRBLS
!! @arg @c VRBLS
!> @details <b>Driver:</b> 
!<
!> @details <b>Calls:</b>
!! @arg @c AVRH
!! @arg @c BOCOHMPI
!! @arg @c CORNERHM2
!--------------------------------------------------------------------------------------------------
    SUBROUTINE EPS
!--------------------------------------------------------------------------------------------------
! SUBROUTINE EPS
!
! SUBPROGRAM: EPS - EPS COMPUTES THE VERTICAL AND HORIZONTAL ADVECTION OF DZ/DT
! PROGRAMMER: JANJIC
! ORG: W/NP22
! DATE: ??-??-??
!
! ABSTRACT:
! EPS COMPUTES THE VERTICAL AND HORIZONTAL ADVECTION OF DZ/DT
!
! PROGRAM HISTORY LOG:
! ??-??-??  JANJIC     - ORIGINATOR
! 18-01-15  LUCCI      - MODERNIZATION OF THE CODE, INCLUDING:
!                        * F77 TO F90/F95
!                        * INDENTATION & UNIFORMIZATION CODE
!                        * REPLACEMENT OF COMMONS BLOCK FOR MODULES
!                        * DOCUMENTATION WITH DOXYGEN
!                        * OPENMP FUNCTIONALITY
!
! INPUT ARGUMENT LIST:                                                                                                                                            
! NONE                                       
!  
! OUTPUT ARGUMENT LIST:
! NONE 
!
! INPUT/OUTPUT ARGUMENT LIST:
! NONE                                                 
!        
! USE MODULES: CONTIN
!              CTLBLK
!              DYNAM
!              F77KINDS
!              LOOPS
!              MASKS
!              METRCS
!              MPPSTAFF
!              NHYDRO
!              PARMETA
!              PVRBLS
!              VRBLS
!
! DRIVER     : -----
!
! CALLS      : AVRH
!              BOCOHMPI
!              CORNERHM2
!-------------------------------------------------------------------------------------------------- 
    USE CONTIN
    USE CTLBLK
    USE DYNAM
    USE F77KINDS
    USE LOOPS
    USE MASKS
    USE METRCS
    USE MPPSTAFF
    USE NHYDRO
    USE PARMETA
    USE PVRBLS
    USE VRBLS
!
    INCLUDE "mpif.h"
!--------------------------------------
! LOCAL VARIABLES REQUIRING DECLARATION 
!--------------------------------------
    INTEGER(KIND=I4)    , PARAMETER :: NTSHY  =    3
!
    REAL   (KIND=R4)    , PARAMETER :: EPSQ   =    1.E-12
    REAL   (KIND=R4)    , PARAMETER :: FF1    =    0.52500
    REAL   (KIND=R4)    , PARAMETER :: FF2    =   -0.64813
    REAL   (KIND=R4)    , PARAMETER :: FF3    =    0.24520
    REAL   (KIND=R4)    , PARAMETER :: FF4    =   -0.12189
    REAL   (KIND=R4)    , PARAMETER :: HBM2IJ =    1. 
!
    REAL   (KIND=R4)    , PARAMETER :: EPSFC  =    1. / 1.05
    REAL   (KIND=R4)    , PARAMETER :: EPSN   = -EPSFC
    REAL   (KIND=R4)    , PARAMETER :: EPSP   =  EPSFC
    REAL   (KIND=R4)    , PARAMETER :: ZERO   =    1.E-06
    REAL   (KIND=R4)    , PARAMETER :: G      =    9.8
    REAL   (KIND=R4)    , PARAMETER :: CAPA   =  287.04 / CP
    REAL   (KIND=R4)    , PARAMETER :: GMA    = -287.04 * (1. - CAPA) / 2. 
    REAL   (KIND=R4)    , PARAMETER :: EPSVW  = 9999.
!
    INTEGER(KIND=I4)                                                                            ::&
    & LMP     , LAP     , LLAP    , KP 
!
    REAL   (KIND=R4)    , DIMENSION(LM+1)                                                       ::&
    & PONE    , PSTR    , PNP1    , COFF    , DFRC    , CHI     , CHIN    , WRES    , CLIM
!
    REAL   (KIND=R4)    , DIMENSION(LM)                                                         ::&
    & RDPP    , C0      , B1      , B2      , B3      , ETADTL  , AFR
!
    REAL   (KIND=R4)    , DIMENSION(LM+1)                                                       ::&
    & W3      , W4      , DWL
!
    INTEGER(KIND=I4)    , DIMENSION(LM)                                                         ::&
    & LA
!
    REAL   (KIND=R4)    , DIMENSION(0:IM+1, 0:JM+1, LM)                                         ::&
    & AFP     , AFQ     , W1      , DWST    , DVOL
!
    REAL   (KIND=R4)    , DIMENSION(0:IM+1, 0:JM+1)                                             ::&
    & DARE    , EMH     , ANE     , ASE
!
    INTEGER(KIND=I4)    , DIMENSION(0:IM+1, 0:JM+1, LM)                                         ::&
    & IFPA    , IFQA    , IFPF    , IFQF    , JFPA    , JFQA    , JFPF    , JFQF
!
    REAL   (KIND=R4)    , DIMENSION(2, LM)                                                      ::&
    & GSUMS   , XSUMS
!      
    REAL   (KIND=R4)                                                                            ::&
    & SUMPW   , SUMNW   , AFRP    , ARR     , D2PQW   , DETAL   , DELP    , DPSTR   , DPTL    ,   &
    & DPTU    , PDP     , PP1     , RDPUP   , RDPDN   , DPPL    , FCT     , PSTRUP  , PSTRDN  ,   &
    & GDT     , GDT2    , DWDTMX  , DWDTMN 
!
    REAL   (KIND=R4)    , DIMENSION(0:IM+1, 0:JM+1)                                             ::&
    & HM      , DPDE    , GDPD    , RGDP    , AGDP    , UT      , VT      , WT      , UTLD    ,   &
    & VTLD    , FXE     , FYE     , FXB     , FYB     , WXE     , WYE     , WXB     , WYB     ,   &
    & ADW 
!                              
    REAL   (KIND=R4)                                                                            ::&
    & FAC     , WADXE   , WADYE   , WADXB   , WADYB   , DWP     , RFACW   , RR      , W00     ,   &
    & W4P     , WP      , WP0     , A       , B       , C
!
    LOGICAL(KIND=L4)                                                                            ::&
    & TOP     , BOT
!
    IF (NTSD <= NTSHY) THEN 
!
!$omp parallel do
!
        DO J=0,JM+1
            DO I=0,IM+1 
                PINT(I,J,1) = PT
            END DO
        END DO
!
        DO L=1,LM
!
!$omp parallel do
!
            DO J=0,JM+1
                DO I=0,IM+1 
                     DWDT(I, J, L  ) = 1.
                    PDWDT(I, J, L  ) = 1.
                     PINT(I, J, L+1) = PDSL(I,J) * DETA(L) + PINT(I,J,L)
                END DO
            END DO
!
        END DO
!
        RETURN
!
    END IF
!
    ADDT = DT
    RDT  =  1. / ADDT
!--------------
! TIME TENDENCY
!--------------
!
!$omp parallel do
!
    DO L=1,LM
!
        DO J=1,JM
            DO I=1,IM 
                DWDT(I,J,L) = (W(I,J,L) - DWDT(I,J,L)) * HTM(I,J,L) * RDT   
            END DO
        END DO
!
    END DO
!-------------------
! VERTICAL ADVECTION 
!-------------------
!
!$omp parallel do
!$omp private (AFR      , AFRP    , ARR     , BOT     , D2PQW   , DETAL   , DWL     , DWP     ,   &
!$omp          ETADTL   , LA      , LAP     , LLAP    , LMP     , RFACW   , RR      , SUMNW   ,   &
!$omp          SUMPW    , TOP     , W00     , W3      , W4      , W4P     , WP      , WP0)
!
    DO J=1,JM
        DO I=1,IM
            LMP = LMH(I,J)
!
            DO L=1,LMP
                W3(L) =  W(I,J,L)
                W4(L) = W3(L)
            END DO
!
            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
            END DO
!
            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
                END IF
!
                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.           
                END IF
!
            END DO
!
            IF (TOP) THEN
                IF (LA(2) < 0) THEN
                    DWL(1) = -DWL(2) * DETA(2) / DETA(1)
                END IF
            END IF
!
            IF (BOT) THEN
                IF (LA(LMP-1) > 0) THEN
                    DWL(LMP) = -DWL(LMP-1) * DETA(LMP-1) / DETA(LMP)
                END IF
            END IF
!
            DO L=1,LMP
                DETAL = DETA(L)
                DWP   =  DWL(L) * DETAL
!
                IF (DWP > 0.) THEN
                    SUMPW = SUMPW + DWP
                ELSE
                    SUMNW = SUMNW + DWP
                END IF
!
            END DO
!-------------------------------
! FIRST MOMENT CONSERVING FACTOR
!-------------------------------
            IF (SUMPW > 1.E-9) THEN
                RFACW = -SUMNW / SUMPW
            ELSE
                RFACW = 1.
            END IF
!
            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
                END DO 
            ELSE
                DO L=1,LMP
                    DWP = DWL(L)
                    IF (DWP >= 0.)  DWP = DWP * RFACW
                    W4(L) = W3(L) + DWP
                END DO
            END IF
!--------------------
! ANTI-FILTERING STEP
!--------------------
        SUMPW = 0.
        SUMNW = 0.
!------------------------
! ANTI-FILTERING LIMITERS
!------------------------
        DO 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.
            END IF
!
            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
            END IF
!
        END DO
!
        DWL(1)   = 0.
!
        DWL(LMP) = 0.
!-------------------------------
! FIRST MOMENT CONSERVING FACTOR
!-------------------------------
        IF (SUMPW > 1.E-9) THEN
            RFACW = -SUMNW / SUMPW
        ELSE
            RFACW = 1.
        END IF
!
        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
            END DO
        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
                DWP         = (W3(LLAP) - W3(L)) * RR
            END DO
        END IF
!
        END DO
    END DO
!--------------------------
! END OF VERTICAL ADVECTION
!--------------------------

!---------------------
! HORIZONTAL ADVECTION
!--------------------
    DO L=1,LM
!
        HM(:,:) = HTM(0:IM+1, 0:JM+1, L)
!
        DPDE(:,:) = PDSL(0:IM+1, 0:JM+1) * DETA(L)
        GDPD(:,:) = DPDE(:,:) * SQH(:,:)
!
        DO I=1,IM
            DO J=1,JM
                RGDP(I,J) = 1. / GDPD(I,J)
            END DO
        END DO
!      
        CALL BOCOHMPI(RGDP,1)
!    
        UT(:,:) = U(0:IM+1, 0:JM+1, L)
        VT(:,:) = V(0:IM+1, 0:JM+1, L)
!    
        WT(:,:) = W(0:IM+1, 0:JM+1, L) 
!
        CALL AVRH(GDPD, AGDP, IM, JM)
        CALL CORNERHM2(AGDP, GDPD)
!
        DO IT=1,2
!-------
! FLUXES
!-------
            DO J=0,JM
                DO I=0,IM
                    UTLD(I,J) =  Q11(I,J) *   UT(I,J) + Q12(I,J) * VT(I,J)
                    VTLD(I,J) =  Q12(I,J) *   UT(I,J) + Q22(I,J) * VT(I,J)
                     FXE(I,J) = AGDP(I,J) * UTLD(I,J)  
                     FYE(I,J) = AGDP(I,J) * VTLD(I,J)  
                END DO
            END DO
!    
            DO J=1,JM
                DO I=0,IM
                    FXB(I,J) =  FXE(I,J) + FYE(I,J) + FXE(I,J-1) + FYE(I,J-1)
                END DO
            END DO
!
            DO J=0,JM
                DO I=1,IM
                    FYB(I,J) = -FXE(I,J) + FYE(I,J) - FXE(I-1,J) + FYE(I-1,J)
                END DO
            END DO     
!---------------------
!  VERTICAL VELOCITY W
!---------------------
            DO J=0,JM
                DO I=0,IM
	            WXE(I,J) = FXE(I,J) * (WT(I+1,J  ) - WT(I,J+1)) 
                    WYE(I,J) = FYE(I,J) * (WT(I+1,J+1) - WT(I,J  )) 
                END DO
            END DO
!      
            DO J=1,JM
                DO I=0,IM
                    WXB(I,J) = FXB(I,J) * (WT(I+1,J) - WT(I,J)) 
                END DO
            END DO
!      
            DO J=0,JM
                DO I=1,IM  
                    WYB(I,J) = FYB(I,J) * (WT(I,J+1) - WT(I,J)) 
                END DO
            END DO
!      
            DO J=1,JM
                DO I=1,IM 
                    FAC   = FADT       * RGDP(I,J)
!
                    WADXE = WXE(I,J-1) + WXE(I-1,J  )
                    WADYE = WYE(I,J  ) + WYE(I-1,J-1)
                    WADXB = WXB(I,J  ) + WXB(I-1,J  )
                    WADYB = WYB(I,J  ) + WYB(I  ,J-1)
!
                    ADW(I,J) = FAC * HM(I,J) * (WADXE + WADYE + WADXB + WADYB)
                END DO 
            END DO
!
            IF (IT == 1) THEN 
!  
	        DO J=1,JM
                    DO I=1,IM
                        WT(I,J) = WT(I,J) + ADW(I,J) * 0.707 
                    END DO
                END DO
!
                CALL BOCOHMPI(WT,1)
!
            END IF
!
        END DO
!
	DO J=1,JM
            DO I=1,IM
                W(I,J,L) = W(I,J,L) + ADW(I,J)
            END DO
        END DO
!
        DO J=1,JM
            DO I=1,IM 

                DWDT(I,J,L) = ADW(I,J) + DWDT(I,J,L)
        
            END DO 
        END DO
!
    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 L=1,LM
!
        DO J=1,JM
            DO I=1,IM
!
                HM(I,J) =  HTM(I,J,L)
!       
                DWDTT   = DWDT(I,J,L) * HM(I,J)
!
                DWDTMX  = AMAX1(DWDTT, DWDTMX)
                DWDTMN  = AMIN1(DWDTT, DWDTMN)
!
                IF (DWDTT < EPSN) THEN
                    DWDTT = EPSN
                    KN    = KN + 1
                END IF
!
                IF (DWDTT > EPSP) THEN
                    DWDTT = EPSP
                    KP    = KP + 1
                END IF
!
                DWDT(I,J,L) = (DWDTT / G + 1.) * (1. - WP) + PDWDT(I,J,L) * WP
!
            END DO
        END DO
!
    END DO
!------------------------------- 
! VERTICAL SOUND WAVES TREATMENT
!-------------------------------
    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      ,   &
!$omp          DPSTR    , DPTL    , DPTU    , DWDTT   , HBM2IJ  , IMAX    , INOT    , ITER    ,   &
!$omp          ITMX     , JMAX    , L       , LMP     , PDP     , PNP1    , PONE    , PP1     ,   &
!$omp          PSTR     , RDP     , RESDL   , WIL     , WRES)
!
    DO J=1,JM
        DO I=1,IM
!
            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.
                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)))      
            END DO
!
            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
            END DO
!
            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)
            END DO
!
            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)
            END DO
!
            DPTU = 0.
!
            FCT  = 0.5 / CP
!
            DO L=1,LMP
!
                DPTL          = CHI(L+1) + PSTR(L+1) - PINT(I,J,L+1)  ! DRAGAN-REMOVED HBM2IJ = HBM2
                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) ! DRAGAN-REMOVED HBM2IJ = HBM2   
                DWDT(I,J,L)   = (DELP - 1.) * HBM2IJ + 1.
!
                DPTU = DPTL
!
                IF (DWDT(I,J,L) < 0) THEN
                    PRINT*, 'DWDT<0', NTSD, I, J, L, DWDT(I,J,L), MYPE
                END IF
!
            END DO
!
            DO L=LMP+1,LM
                PINT(I,J,L+1) = PDSL(I,J) * DETA(L) + PINT(I,J,L)
            END DO
!
        END DO
    END DO
!----------------- 
! END DO MAIN_LOOP
!----------------- 
    CALL BOCOHMPI(W   , LM+1)
    CALL BOCOHMPI(DWDT, LM  )
    CALL BOCOHMPI(PINT, LM+1)
!
    RETURN
!
    END
