!--------------------------------------------------------------------------------------------------
!  DOXYGEN
!> @brief CALCULATES HORIZONTAL ADVECTION AND UPDATES ADVECTED VARIABLES
!> @details CALCULATES HORIZONTAL ADVECTION AND UPDATES ADVECTED VARIABLES, USING 4TH ORDER ADVECTION 
!! SCHEMES: MESINGER (1981) FOR MOMENTUM AND A-LA-RANCIC (1988) FOR TEMPERATURE AND TURBULENT 
!! KINETIC ENERGY, WITH A MODIFIED EULER-BACKWARD (OR HEUN SCHEME) FOR TIME-DIFFERENCING.
!> @author ORIGINATOR - ????? 
!> @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 
!<
!> @param[in] IM - Significado de IM
!> @param[in] JM - Significado de JM
!> @param[in] LM - Significado de LM
!> @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 PVRBLS
!! @arg @c REALPAR
!! @arg @c VRBLS
!> @details <b>Driver:</b> 
!<
!> @details <b>Calls:</b>
!! @arg @c AVRH
!! @arg @c BOCOHMPI
!! @arg @c BOCOVMPI
!! @arg @c CORNERHM2
!--------------------------------------------------------------------------------------------------
    SUBROUTINE HZADV_4TH(IM, JM, LM)
!--------------------------------------------------------------------------------------------------
! SUBROUTINE HZADV_4TH
!
! SUBPROGRAM: HZADV_4TH - CALCULATES HORIZONTAL ADVECTION AND UPDATES ADVECTED VARIABLES
! PROGRAMMER: ?????
! ORG: W/NP22
! DATE: ??-??-??
!
! ABSTRACT:
! CALCULATES HORIZONTAL ADVECTION AND UPDATES ADVECTED VARIABLES, USING 4TH ORDER ADVECTION 
! SCHEMES: MESINGER (1981) FOR MOMENTUM AND A-LA-RANCIC (1988) FOR TEMPERATURE AND TURBULENT 
! KINETIC ENERGY, WITH A MODIFIED EULER-BACKWARD (OR HEUN SCHEME) FOR TIME-DIFFERENCING.
!
! PROGRAM HISTORY LOG:
! ??-??-??  ?????   - ORIGINATOR
! 18-03-20  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:
! IM - 
! JM - 
! LM - 
!
! OUTPUT ARGUMENT LIST:
! NONE
!
! INPUT/OUTPUT ARGUMENT LIST:
! NONE
!
! USE MODULES: CONTIN
!              CTLBLK
!              DYNAM
!              F77KINDS
!              MASKS
!              METRCS
!              MPPSTAFF
!              PVRBLS
!              REALPAR
!              TPRFIL
!              VRBLS
!
! DRIVER     : -----
!
! CALLS      : AVRH
!              BOCOHMPI
!              BOCOVMPI
!              CORNERHM2
!--------------------------------------------------------------------------------------------------
    USE CONTIN
    USE CTLBLK
    USE DYNAM
    USE F77KINDS
    USE LOOPS
    USE MASKS
    USE METRCS
    USE MPPSTAFF
    USE PVRBLS
    USE REALPAR
    USE VRBLS
!
    IMPLICIT NONE
!
    INTEGER(KIND=I4)                                                      , INTENT(IN)          ::&
    & IM      , JM      , LM
!               
    REAL   (KIND=R4)    , DIMENSION(0:IM+1, 0:JM+1)                                             ::&
    & HM      , VM      , DPDE    , GDPD    , RGDP    , AGDP    ,                                 &
    & UT      , VT      , TT      , UTLD    , VTLD    , ZET     ,                                 & 
    & FXE     , FYE     , FXB     , FYB     , AZT     , VYZ     , UXZ     ,                       &
    & TXE     , TYE     , TXB     , TYB     , ADU     , ADV     , ADT     ,                       &
    & EKN     , AKN     , Q2M     , Q2L     , Q2MXE   , Q2LXE   ,                                 &
    & Q2MYE   , Q2LYE   , Q2MXB   , Q2LXB   , Q2MYB   ,                                           &
    & Q2LYB   , ADQ2M   , ADQ2L   , QXE     , QYE     , QT      ,                                 &
    & AZTH    , VXZ     , UYZ     , VYAZ    , UXAZ
! 
    REAL   (KIND=R4)    , DIMENSION(0:IM+1, 0:JM+1)                                             ::&
    & TXB2    , TYB2    , Q2MXB2  , Q2MYB2  , Q2LXB2  , Q2LYB2  , RGDPAX
!
    INTEGER(KIND=I4)                                                                            ::&
    & I       , J       , L       , IT
!
    REAL   (KIND=R4)                                                                            ::&
    & ROTU1   , ROTU2   , ROTV1   , ROTV2   , FAC     , ADXE    , ADYE    , ADXB    , ADYB    ,   &
    & ADXB2   , ADYB2   ,                                                                         &
    & FXB2    , FYB2    ,                                                                         &
    & FXEX    , FYEY    ,                                                                         &
    & FADT3
!
    LOGICAL(KIND=L4)                                                                            ::&
    & ISCORNER2
!
    REAL   (KIND=R4)    , DIMENSION(0:IM+1, 0:JM+1, LM+1)                                       ::&
    & Q2ML    , ADQ2HL
!
    FADT3 = FADT * 3.
!
    Q2ML(0:IM+1, 0:JM+1, 1) = 0.
!
    DO L=2,LM+1
         Q2ML(0:IM+1, 0:JM+1, L) = Q2(0:IM+1, 0:JM+1, L-1)
    END DO
!
    MAIN_LOOP: DO L=1,LM
!
        HM = HTM(0:IM+1, 0:JM+1, L)
        VM = VTM(0:IM+1, 0:JM+1, L)    
!
        DPDE = PDSL(0:IM+1, 0:JM+1) * DETA(L)
        GDPD = DPDE * SQH
!
        CALL AVRH(GDPD, AGDP, IM, JM)
!
        CALL CORNERHM2(AGDP, GDPD)     
!
        DO I=1,IM
            DO J=1,JM
                  RGDP(I,J) = 1. /        GDPD(I,J)
                RGDPAX(I,J) = 1. / (16. * GDPD(I,J) - AGDP(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)
        TT(:,:) = T(0:IM+1, 0:JM+1, L)
!
        Q2M = Q2ML(0:IM+1, 0:JM+1, L  )
        Q2L = Q2ML(0:IM+1, 0:JM+1, L+1)
!
        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)  
                     EKN(I,J) = (UTLD(I,J) *   UT(I,J) + VTLD(I,J) * VT(I,J)) 
                END DO
            END DO
!
            CALL AVRV(EKN,AKN,IM,JM)      
!
            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
!---------
! MOMENTUM
!---------
            DO J=1,JM
                DO I=1,IM
                    ZET(I,J) =  RGDP(I,J) * (VT(I,J-1) - VT(I-1,J) - UT(I,J) + UT(I-1,J-1))
                END DO
            END DO
!    
            CALL BOCOHMPI(ZET, 1)
!
            CALL AVRH(ZET, AZT, IM, JM) 
!     
            CALL AVRHH(ZET, AZTH, IM, JM)  
!
            CALL BOCOHMPI(AZTH,1)
!
            DO J=1,JM
                DO I=1,IM
                     UXZ(I,J) =  ZET(I,J) * (FXE(I-1,J  ) + FXE(I,J-1))
                     VXZ(I,J) =  ZET(I,J) * (FYE(I-1,J  ) + FYE(I,J-1))
                    UXAZ(I,J) = AZTH(I,J) * (FXE(I-1,J  ) + FXE(I,J-1))
                     UYZ(I,J) =  ZET(I,J) * (FXE(I-1,J-1) + FXE(I,J  ))
                     VYZ(I,J) =  ZET(I,J) * (FYE(I-1,J-1) + FYE(I,J  ))
                    VYAZ(I,J) = AZTH(I,J) * (FYE(I-1,J-1) + FYE(I,J  ))
                END DO
            END DO
!
            DO J=1,JM-1
                DO I=1,IM-1
                    ROTU1 = 4. * FYE(I,J) * AZT(I,J) - 0.5 * (VYAZ(I+1, J  ) + VYAZ(I, J+1))
                    ROTV1 = 4. * FXE(I,J) * AZT(I,J) - 0.5 * (UXAZ(I+1, J+1) + UXAZ(I, J  ))
!
                    ROTU2 = 2. * (VYZ(I  , J+1) + VYZ(I+1, J)) - (VXZ(I+1, J+1) + VXZ(I  , J))
                    ROTV2 = 2. * (UXZ(I+1, J+1) + UXZ(I  , J)) - (UYZ(I  , J+1) + UYZ(I+1, J))
!
                    ADU(I,J) = ( FADV * (ROTU1 + ROTU2) + FKIN * (AKN(I+1, J  )                   &
    &                        -                                    AKN(I  , J+1))) * VM(I,J)
!
                    ADV(I,J) = (-FADV * (ROTV1 + ROTV2) + FKIN * (AKN(I+1, J+1)                   &
    &                        -                                    AKN(I  , J  ))) * VM(I,J)
                END DO
            END DO
!-----------------------------------------
! TEMPERATURE AND TURBULENT KINETIC ENERGY
!-----------------------------------------
            DO J=0,JM
                DO I=0,IM
                      TXE(I,J) = FXE(I,J) * ( TT(I+1, J  ) -  TT(I, J+1))
                      TYE(I,J) = FYE(I,J) * ( TT(I+1, J+1) -  TT(I, J  ))
                    Q2MXE(I,J) = FXE(I,J) * (Q2M(I+1, J  ) - Q2M(I, J+1))
                    Q2LXE(I,J) = FXE(I,J) * (Q2L(I+1, J  ) - Q2L(I, J+1))
                    Q2MYE(I,J) = FYE(I,J) * (Q2M(I+1, J+1) - Q2M(I, J  ))
                    Q2LYE(I,J) = FYE(I,J) * (Q2L(I+1, J+1) - Q2L(I, J  ))
                END DO
            END DO
!
            DO J=1,JM
                DO I=0,IM
                      TXB(I,J) = FXB(I,J) * ( TT(I+1, J) -  TT(I,J))
                    Q2MXB(I,J) = FXB(I,J) * (Q2M(I+1, J) - Q2M(I,J))
                    Q2LXB(I,J) = FXB(I,J) * (Q2L(I+1, J) - Q2L(I,J))
                END DO
            END DO
!
            DO J=0,JM
                DO I=1,IM
                      TYB(I,J) = FYB(I,J) * ( TT(I, J+1) -  TT(I,J))
                    Q2MYB(I,J) = FYB(I,J) * (Q2M(I, J+1) - Q2M(I,J))
                    Q2LYB(I,J) = FYB(I,J) * (Q2L(I, J+1) - Q2L(I,J))
                END DO
            END DO
!
            DO J=1,JM
                DO I=1,IM
                    FXEX = FXE(I-1, J) + FXE(I  , J-1)
                    FYEY = FYE(I  , J) + FYE(I-1, J-1)
!
                    FXB2 = FXEX + FYEY
!
                      TXB2(I,J) = FXB2 * ( TT(I+1, J) -  TT(I-1, J))
                    Q2MXB2(I,J) = FXB2 * (Q2M(I+1, J) - Q2M(I-1, J))
                    Q2LXB2(I,J) = FXB2 * (Q2L(I+1, J) - Q2L(I-1, J))
!
                    FYB2 = -FXEX + FYEY
!
                      TYB2(I,J) = FYB2 * ( TT(I, J+1) -  TT(I, J-1))
                    Q2MYB2(I,J) = FYB2 * (Q2M(I, J+1) - Q2M(I, J-1))
                    Q2LYB2(I,J) = FYB2 * (Q2L(I, J+1) - Q2L(I, J-1))
                END DO
            END DO
!
            CALL BOCOHMPI(TXB2  , 1)
            CALL BOCOHMPI(Q2MXB2, 1)
            CALL BOCOHMPI(Q2LXB2, 1)
            CALL BOCOHMPI(TYB2  , 1)
            CALL BOCOHMPI(Q2MYB2, 1)
            CALL BOCOHMPI(Q2LYB2, 1)
!
            DO J=1,JM
                DO I=1,IM 
                    FAC  = FADT * RGDP(I,J)
                    ADXE = TXE(I, J-1) + TXE(I-1, J  )
                    ADYE = TYE(I, J  ) + TYE(I-1, J-1)
                    ADXB = TXB(I, J  ) + TXB(I-1, J  )
                    ADYB = TYB(I, J  ) + TYB(I  , J-1)
!
                    ADT(I, J) = FAC * HM(I, J) * (ADXE + ADYE + ADXB + ADYB)
!
                    ADXE = Q2MXE(I, J-1) + Q2MXE(I-1, J  )
                    ADYE = Q2MYE(I, J  ) + Q2MYE(I-1, J-1)
                    ADXB = Q2MXB(I, J  ) + Q2MXB(I-1, J  )
                    ADYB = Q2MYB(I, J  ) + Q2MYB(I  , J-1)
!
                    ADQ2M(I, J) = FAC * HM(I,J) * (ADXE + ADYE + ADXB + ADYB)
!   
                    ADXE = Q2LXE(I, J-1) + Q2LXE(I-1, J  )
                    ADYE = Q2LYE(I, J  ) + Q2LYE(I-1, J-1)
                    ADXB = Q2LXB(I, J  ) + Q2LXB(I-1, J  )
                    ADYB = Q2LYB(I, J  ) + Q2LYB(I  , J-1)
!
                    ADQ2L(I,J) = FAC * HM(I,J) * (ADXE + ADYE + ADXB + ADYB)
                END DO 
            END DO 
!-------------
! FIRST UPDATE 
!-------------
            IF (IT == 1) THEN
!
                DO J=1,JM-1
                    DO I=1,IM-1
                        UT(I,J) = UT(I,J) + ADU(I,J) * 0.75
                        VT(I,J) = VT(I,J) + ADV(I,J) * 0.75
                    END DO
                END DO
!
                DO J=1,JM
                    DO I=1,IM
                         TT(I,J) =  TT(I,J) +   ADT(I,J) * 0.75 
                        Q2M(I,J) = Q2M(I,J) + ADQ2M(I,J) * 0.75
                        Q2L(I,J) = Q2L(I,J) + ADQ2L(I,J) * 0.75
                    END DO
                END DO
!               
                CALL BOCOVMPI(UT , VT, 1)
                CALL BOCOHMPI(TT , 1)
                CALL BOCOHMPI(Q2M, 1)
                CALL BOCOHMPI(Q2L, 1)
!
            END IF
!
        END DO                         
!    
        DO J=1,JM-1
            DO I=1,IM-1
                U(I,J,L) = U(I,J,L) + ADU(I,J)
                V(I,J,L) = V(I,J,L) + ADV(I,J)
            END DO
        END DO
!
        DO J=1,JM
            DO I=1,IM
                T(I,J,L) = T(I,J,L) + ADT(I,J)
            END DO
        END DO
!
        DO J=1,JM
            DO I=1,IM
                ADQ2HL(I,J,L) = ADQ2L(I,J)
            END DO
        END DO
!
        IF (L > 1) THEN
            DO J=1,JM
                DO I=1,IM
                    Q2(I,J,L-1) = Q2(I,J,L-1) + ADQ2M(I,J) * HM(I,J) * 0.5
                END DO
            END DO
        END IF
!
    END DO MAIN_LOOP
!
    DO L=2,LM
        DO I=1,IM
            DO J=1,JM
                Q2(I,J,L-1) = Q2(I,J,L-1) + ADQ2HL(I,J,L-1) * HTM(I,J,L) * 0.5
            END DO
        END DO
    END DO
!
    DO I=1,IM
        DO J=1,JM
            DO L=1,LMH(I,J)-1
                Q2(I,J,L) = AMAX1(Q2(I,J,L), EPSQ2)
            END DO
        END DO
    END DO
!             
    CALL BOCOVMPI(U , V , LM)
    CALL BOCOHMPI(T , LM)
    CALL BOCOHMPI(Q2, LM)
!
      END SUBROUTINE HZADV_4TH
