!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

!==============================================================================|
   SUBROUTINE ADVECTION_EDGE_GCY(XFLUX,YFLUX)
!==============================================================================|
!   Calculate the Advection and Diffusion Terms of 3D Velocity Field           |
!   These Terms will be vertically integrated to form the Mean Terms in        |
!   the Gx and Gy Terms of the External Mode Equation                          |
!   Ghost cell boundary conditions are used here                               |
!==============================================================================|

   USE MOD_UTILS
   USE ALL_VARS
   USE BCS
   USE MOD_SPHERICAL
   USE MOD_NORTHPOLE
   USE MOD_WD

   IMPLICIT NONE
   REAL(SP), INTENT(OUT), DIMENSION(0:NT,KB) :: XFLUX,YFLUX
   REAL(SP) :: DIJ
   REAL(SP) :: COFA1,COFA2,COFA3,COFA4,COFA5,COFA6,COFA7,COFA8
   REAL(SP) :: XADV,YADV,TXXIJ,TYYIJ,TXYIJ,UN
   REAL(SP) :: VISCOF,VISCOF1,VISCOF2,TEMP,TPA,TPB
   REAL(SP) :: XIJA,YIJA,XIJB,YIJB,UIJ,VIJ
   REAL(SP) :: FACT,FM1
   INTEGER  :: I,IA,IB,J1,J2,K1,K2,K3,K4,K5,K6,K,II,J,I1,I2
   REAL(SP) :: ISWETTMP

#  if defined (SPHERICAL)
   REAL(DP) :: XTMP,XTMP1  !@AJ, LIMITATION@20240515
#  endif

#  if defined (LIMITED_NO)
   REAL(SP) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
   REAL(SP) :: UK1(KB),UK2(KB),UK3(KB),UK4(KB),UK5(KB),UK6(KB), &
               VK1(KB),VK2(KB),VK3(KB),VK4(KB),VK5(KB),VK6(KB)
#  else
   REAL(SP),ALLOCATABLE,DIMENSION(:) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
   REAL(SP),ALLOCATABLE,DIMENSION(:) :: UALFA,VALFA
   REAL(SP) :: UALFA_TMP,VALFA_TMP
   INTEGER :: ERROR
   REAL(SP) :: EPS
   REAL(SP) :: XAB, YAB   ! Lu Wang, LIMITATION@20240328

   REAL(SP) :: UK1,UK2,UK3,UK4,UK5,UK6,VK1,VK2,VK3,VK4,VK5,VK6
#  endif
!------------------------------------------------------------------------------|

   if(dbg_set(dbg_sbr)) write(ipt,*) "Start: advection_edge_gcy.F"


   SELECT CASE(HORIZONTAL_MIXING_TYPE)
   CASE ('closure')
      FACT = 1.0_SP
      FM1  = 0.0_SP
   CASE('constant')
      FACT = 0.0_SP
      FM1  = 1.0_SP
   CASE DEFAULT
      CALL FATAL_ERROR("UNKNOW HORIZONTAL MIXING TYPE:",&
           & TRIM(HORIZONTAL_MIXING_TYPE) )
   END SELECT

!
!--Initialize Variables--------------------------------------------------------|
!
   XFLUX = 0.0_SP
   YFLUX = 0.0_SP

!
!--Loop Over Edges and Accumulate Fluxes-For Each Element----------------------|
!
#  if !defined (LIMITED_NO)
   ALLOCATE(UIJ1(NE),VIJ1(NE),UIJ2(NE),VIJ2(NE),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays UIJ1,VIJ1,UIJ2 and VIJ2 can not be allocated.")
   
   ALLOCATE(UALFA(0:NT),VALFA(0:NT),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays UALFA,VALFA can not be allocated.")
   
   ALLOCATE(FXX(NE),FYY(NE),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays FXX,FYY can not be allocated.")

   DO K=1,KBM1
     UIJ1=0.0_SP;VIJ1=0.0_SP;UIJ2=0.0_SP;VIJ2=0.0_SP
     UALFA=1.0_SP;VALFA=1.0_SP
     FXX=0.0_SP;FYY=0.0_SP
     
     DO I=1,NE
       IA=IEC(I,1)
       IB=IEC(I,2)
#      if defined (WET_DRY)
#      if !defined (SEMI_IMPLICIT)
       IF(ISWETCT(IA)*ISWETC(IA) == 1 .OR. ISWETCT(IB)*ISWETC(IB) == 1)THEN
#      else
       IF(ISWETCT(IA) == 1 .OR. ISWETCT(IB) == 1)THEN
#      endif
#      endif
       J1=IENODE(I,1)
       J2=IENODE(I,2)

       K1=NBE(IA,1)
       K2=NBE(IA,2)
       K3=NBE(IA,3)
       K4=NBE(IB,1)
       K5=NBE(IB,2)
       K6=NBE(IB,3)
#      if defined (SPHERICAL)
       XIJA=DLTXNE(I,1)
       YIJA=DLTYNE(I,1)
       XIJB=DLTXNE(I,2)
       YIJB=DLTYNE(I,2)
!----> Lu Wang, LIMITATION@20240328
       XTMP  = XC(IA)*TPI - XC(IB)*TPI
       XTMP1 = XC(IA) - XC(IB)
       IF(XTMP1 >  180.0_SP)THEN
         XTMP = -360.0_SP*TPI+XTMP
       ELSE IF(XTMP1 < -180.0_SP)THEN
         XTMP =  360.0_SP*TPI+XTMP
       END IF
       XAB = XTMP * COS(DEG2RAD * (YC(IA)+YC(IB))*0.5_SP)
       YAB = (YC(IA) - YC(IB)) * TPI
!<----
#      else
       XIJA=XIJC(I)-XC(IA)
       YIJA=YIJC(I)-YC(IA)
       XIJB=XIJC(I)-XC(IB)
       YIJB=YIJC(I)-YC(IB)
!----> Lu Wang, LIMITATION@20240328
       XAB = XC(IA) - XC(IB)
       YAB = YC(IA) - YC(IB)
!<----
#      endif

       UK1 = U(K1,K)
       UK2 = U(K2,K)
       UK3 = U(K3,K)
       UK4 = U(K4,K)
       UK5 = U(K5,K)
       UK6 = U(K6,K)
       VK1 = V(K1,K)
       VK2 = V(K2,K)
       VK3 = V(K3,K)
       VK4 = V(K4,K)
       VK5 = V(K5,K)
       VK6 = V(K6,K)

       IF(K1 == 0) CALL GHOSTUV3(IA,1,UK1,VK1,K)
       IF(K2 == 0) CALL GHOSTUV3(IA,2,UK2,VK2,K)
       IF(K3 == 0) CALL GHOSTUV3(IA,3,UK3,VK3,K)
       IF(K4 == 0) CALL GHOSTUV3(IB,1,UK4,VK4,K)
       IF(K5 == 0) CALL GHOSTUV3(IB,2,UK5,VK5,K)
       IF(K6 == 0) CALL GHOSTUV3(IB,3,UK6,VK6,K)

       DIJ= 0.5_SP*(DT(J1)*DZ(J1,K)+DT(J2)*DZ(J2,K))

       !!FORM THE LEFT FLUX
       COFA1=A1U(IA,1)*U(IA,K)+A1U(IA,2)*UK1+A1U(IA,3)*UK2+A1U(IA,4)*UK3
       COFA2=A2U(IA,1)*U(IA,K)+A2U(IA,2)*UK1+A2U(IA,3)*UK2+A2U(IA,4)*UK3
       COFA5=A1U(IA,1)*V(IA,K)+A1U(IA,2)*VK1+A1U(IA,3)*VK2+A1U(IA,4)*VK3
       COFA6=A2U(IA,1)*V(IA,K)+A2U(IA,2)*VK1+A2U(IA,3)*VK2+A2U(IA,4)*VK3
       UIJ1(I)=COFA1*XIJA+COFA2*YIJA
       VIJ1(I)=COFA5*XIJA+COFA6*YIJA
!----> Lu Wang, LIMITATION@20240328
!       UALFA_TMP=ABS(U(IA,K)-U(IB,K))/ABS(UIJ1(I)+EPSILON(EPS))
!       VALFA_TMP=ABS(V(IA,K)-V(IB,K))/ABS(VIJ1(I)+EPSILON(EPS))
       UALFA_TMP=ABS(U(IA,K)-U(IB,K))/(ABS(COFA1*XAB + COFA2*YAB) + EPSILON(EPS))
       VALFA_TMP=ABS(V(IA,K)-V(IB,K))/(ABS(COFA5*XAB + COFA6*YAB) + EPSILON(EPS))
!<----
       IF(UALFA_TMP > 1)UALFA_TMP = 1.0_SP
       IF(VALFA_TMP > 1)VALFA_TMP = 1.0_SP
       UALFA(IA)=MIN(UALFA(IA),UALFA_TMP)
       VALFA(IA)=MIN(VALFA(IA),VALFA_TMP)

       !!FORM THE RIGHT FLUX
       COFA3=A1U(IB,1)*U(IB,K)+A1U(IB,2)*UK4+A1U(IB,3)*UK5+A1U(IB,4)*UK6
       COFA4=A2U(IB,1)*U(IB,K)+A2U(IB,2)*UK4+A2U(IB,3)*UK5+A2U(IB,4)*UK6
       COFA7=A1U(IB,1)*V(IB,K)+A1U(IB,2)*VK4+A1U(IB,3)*VK5+A1U(IB,4)*VK6
       COFA8=A2U(IB,1)*V(IB,K)+A2U(IB,2)*VK4+A2U(IB,3)*VK5+A2U(IB,4)*VK6
       UIJ2(I)=COFA3*XIJB+COFA4*YIJB
       VIJ2(I)=COFA7*XIJB+COFA8*YIJB
!----> Lu Wang, LIMITATION@20240328
!       UALFA_TMP=ABS(U(IA,K)-U(IB,K))/ABS(UIJ2(I)+EPSILON(EPS))
!       VALFA_TMP=ABS(V(IA,K)-V(IB,K))/ABS(VIJ2(I)+EPSILON(EPS))
       UALFA_TMP=ABS(U(IA,K)-U(IB,K))/(ABS(COFA3*XAB + COFA4*YAB) + EPSILON(EPS))
       VALFA_TMP=ABS(V(IA,K)-V(IB,K))/(ABS(COFA7*XAB + COFA8*YAB) + EPSILON(EPS))
!<----
       IF(UALFA_TMP > 1)UALFA_TMP = 1.0_SP
       IF(VALFA_TMP > 1)VALFA_TMP = 1.0_SP
       UALFA(IB)=MIN(UALFA(IB),UALFA_TMP)
       VALFA(IB)=MIN(VALFA(IB),VALFA_TMP)

       VISCOF1=ART(IA)*SQRT(COFA1**2+COFA6**2+0.5_SP*(COFA2+COFA5)**2)
       VISCOF2=ART(IB)*SQRT(COFA3**2+COFA8**2+0.5_SP*(COFA4+COFA7)**2)

!       VISCOF = HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2)/HPRNU + FM1)
!       VISCOF = HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2) + FM1)/HPRNU
       ! David moved HPRNU and added HVC
       VISCOF=(FACT*0.5_SP*(VISCOF1*CC_HVC(IA)+VISCOF2*CC_HVC(IB)) + FM1*0.5_SP*(CC_HVC(IA)+CC_HVC(IB)))/HPRNU

       TXXIJ=(COFA1+COFA3)*VISCOF
       TYYIJ=(COFA6+COFA8)*VISCOF
       TXYIJ=0.5_SP*(COFA2+COFA4+COFA5+COFA7)*VISCOF
       FXX(I)=DIJ*(TXXIJ*DLTYC(I)-TXYIJ*DLTXC(I))
       FYY(I)=DIJ*(TXYIJ*DLTYC(I)-TYYIJ*DLTXC(I))

#    if defined (PLBC)
       TXYIJ=0.5_SP*(COFA5+COFA7)*VISCOF
       FXX(I)=DIJ*(TXXIJ*DLTYC(I)-0.0_SP*DLTXC(I))
       FYY(I)=DIJ*(TXYIJ*DLTYC(I)-0.0_SP*DLTXC(I))
#    endif

#      if defined (WET_DRY)
       END IF
#      endif
     END DO

     DO I=1,NE
       IA=IEC(I,1)
       IB=IEC(I,2)
#      if defined (WET_DRY)
#      if !defined (SEMI_IMPLICIT)
       IF(ISWETCT(IA)*ISWETC(IA) == 1 .OR. ISWETCT(IB)*ISWETC(IB) == 1)THEN
#      else
       IF(ISWETCT(IA) == 1 .OR. ISWETCT(IB) == 1)THEN
#      endif
#      endif
       J1=IENODE(I,1)
       J2=IENODE(I,2)

       DIJ= 0.5_SP*(DT(J1)*DZ(J1,K)+DT(J2)*DZ(J2,K))
       UIJ1(I)=U(IA,K)+UALFA(IA)*UIJ1(I)
       VIJ1(I)=V(IA,K)+VALFA(IA)*VIJ1(I)
       UIJ2(I)=U(IB,K)+UALFA(IB)*UIJ2(I)
       VIJ2(I)=V(IB,K)+VALFA(IB)*VIJ2(I)

#      if defined (LIMITED_1)
       IF(UIJ1(I) > MAX(U(IA,K),U(IB,K)) .OR. UIJ1(I) < MIN(U(IA,K),U(IB,K)) .OR. &
          UIJ2(I) > MAX(U(IA,K),U(IB,K)) .OR. UIJ2(I) < MIN(U(IA,K),U(IB,K)))THEN
         UIJ1(I)=U(IA,K)
         UIJ2(I)=U(IB,K)
       END IF

       IF(VIJ1(I) > MAX(V(IA,K),V(IB,K)) .OR. VIJ1(I) < MIN(V(IA,K),V(IB,K)) .OR. &
          VIJ2(I) > MAX(V(IA,K),V(IB,K)) .OR. VIJ2(I) < MIN(V(IA,K),V(IB,K)))THEN
         VIJ1(I)=V(IA,K)
         VIJ2(I)=V(IB,K)
       END IF
#      endif       

       !!COMPUTE THE NORMAL VELOCITY ACROSS THE EDGE
       UIJ=0.5_SP*(UIJ1(I)+UIJ2(I))
       VIJ=0.5_SP*(VIJ1(I)+VIJ2(I))
       UN=VIJ*DLTXC(I) - UIJ*DLTYC(I)

#    if defined (PLBC)
       UN=0.0_SP*DLTXC(I) - UIJ*DLTYC(I)
#    endif

       !!UPWIND THE ADVECTIVE FLUX
       XADV=DIJ*UN*((1.0_SP-SIGN(1.0_SP,UN))*UIJ2(I)+(1.0_SP+SIGN(1.0_SP,UN))*UIJ1(I))*0.5_SP
       YADV=DIJ*UN*((1.0_SP-SIGN(1.0_SP,UN))*VIJ2(I)+(1.0_SP+SIGN(1.0_SP,UN))*VIJ1(I))*0.5_SP

       !!COMPUTE BOUNDARY FLUX AUGMENTERS
       TPA = FLOAT(1-ISBC(I))*EPOR(IA)
       TPB = FLOAT(1-ISBC(I))*EPOR(IB)

       !!ACCUMULATE THE FLUX
       XFLUX(IA,K)=XFLUX(IA,K)+XADV*TPA+FXX(I)*TPA
       YFLUX(IA,K)=YFLUX(IA,K)+YADV*TPA+FYY(I)*TPA
       XFLUX(IB,K)=XFLUX(IB,K)-XADV*TPB-FXX(I)*TPB
       YFLUX(IB,K)=YFLUX(IB,K)-YADV*TPB-FYY(I)*TPB

#      if defined (WET_DRY)
       END IF
#      endif
     END DO
   END DO

   DEALLOCATE(UIJ1,VIJ1,UIJ2,VIJ2,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for UIJ1,VIJ1,UIJ2 and VIJ2.")
   DEALLOCATE(UALFA,VALFA,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for UALFA,VALFA.")
   DEALLOCATE(FXX,FYY,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for FXX,FYY.")

#  else

   DO I=1,NE
     IA=IEC(I,1)
     IB=IEC(I,2)
#  if defined (WET_DRY)
#   if !defined (SEMI_IMPLICIT)
    IF(ISWETCT(IA)*ISWETC(IA) == 1 .OR. ISWETCT(IB)*ISWETC(IB) == 1)THEN
#   else
    IF(ISWETCT(IA) == 1 .OR. ISWETCT(IB) == 1)THEN
#   endif
#  endif
     J1=IENODE(I,1)
     J2=IENODE(I,2)

     K1=NBE(IA,1)
     K2=NBE(IA,2)
     K3=NBE(IA,3)
     K4=NBE(IB,1)
     K5=NBE(IB,2)
     K6=NBE(IB,3)
#    if defined (SPHERICAL)
     XIJA=DLTXNE(I,1)
     YIJA=DLTYNE(I,1)
     XIJB=DLTXNE(I,2)
     YIJB=DLTYNE(I,2)
#    else
     XIJA=XIJC(I)-XC(IA)
     YIJA=YIJC(I)-YC(IA)
     XIJB=XIJC(I)-XC(IB)
     YIJB=YIJC(I)-YC(IB)
#    endif

     UK1 = U(K1,:)
     UK2 = U(K2,:)
     UK3 = U(K3,:)
     UK4 = U(K4,:)
     UK5 = U(K5,:)
     UK6 = U(K6,:)
     VK1 = V(K1,:)
     VK2 = V(K2,:)
     VK3 = V(K3,:)
     VK4 = V(K4,:)
     VK5 = V(K5,:)
     VK6 = V(K6,:)

     IF(K1 == 0) CALL GHOSTUV3(IA,1,UK1,VK1)
     IF(K2 == 0) CALL GHOSTUV3(IA,2,UK2,VK2)
     IF(K3 == 0) CALL GHOSTUV3(IA,3,UK3,VK3)
     IF(K4 == 0) CALL GHOSTUV3(IB,1,UK4,VK4)
     IF(K5 == 0) CALL GHOSTUV3(IB,2,UK5,VK5)
     IF(K6 == 0) CALL GHOSTUV3(IB,3,UK6,VK6)

     DO K=1,KBM1

       DIJ= 0.5_SP*(DT(J1)*DZ(J1,K)+DT(J2)*DZ(J2,K))

       !!FORM THE LEFT FLUX
       COFA1=A1U(IA,1)*U(IA,K)+A1U(IA,2)*UK1(K)+A1U(IA,3)*UK2(K)+A1U(IA,4)*UK3(K)
       COFA2=A2U(IA,1)*U(IA,K)+A2U(IA,2)*UK1(K)+A2U(IA,3)*UK2(K)+A2U(IA,4)*UK3(K)
       COFA5=A1U(IA,1)*V(IA,K)+A1U(IA,2)*VK1(K)+A1U(IA,3)*VK2(K)+A1U(IA,4)*VK3(K)
       COFA6=A2U(IA,1)*V(IA,K)+A2U(IA,2)*VK1(K)+A2U(IA,3)*VK2(K)+A2U(IA,4)*VK3(K)
       UIJ1=U(IA,K)+COFA1*XIJA+COFA2*YIJA
       VIJ1=V(IA,K)+COFA5*XIJA+COFA6*YIJA

       !!FORM THE RIGHT FLUX
       COFA3=A1U(IB,1)*U(IB,K)+A1U(IB,2)*UK4(K)+A1U(IB,3)*UK5(K)+A1U(IB,4)*UK6(K)
       COFA4=A2U(IB,1)*U(IB,K)+A2U(IB,2)*UK4(K)+A2U(IB,3)*UK5(K)+A2U(IB,4)*UK6(K)
       COFA7=A1U(IB,1)*V(IB,K)+A1U(IB,2)*VK4(K)+A1U(IB,3)*VK5(K)+A1U(IB,4)*VK6(K)
       COFA8=A2U(IB,1)*V(IB,K)+A2U(IB,2)*VK4(K)+A2U(IB,3)*VK5(K)+A2U(IB,4)*VK6(K)
       UIJ2=U(IB,K)+COFA3*XIJB+COFA4*YIJB
       VIJ2=V(IB,K)+COFA7*XIJB+COFA8*YIJB

       !!COMPUTE THE NORMAL VELOCITY ACROSS THE EDGE
       UIJ=0.5_SP*(UIJ1+UIJ2)
       VIJ=0.5_SP*(VIJ1+VIJ2)
       UN=VIJ*DLTXC(I) - UIJ*DLTYC(I)

       VISCOF1=ART(IA)*SQRT(COFA1**2+COFA6**2+0.5_SP*(COFA2+COFA5)**2)
       VISCOF2=ART(IB)*SQRT(COFA3**2+COFA8**2+0.5_SP*(COFA4+COFA7)**2)

!       VISCOF = HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2)/HPRNU + FM1)
!       VISCOF = HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2) + FM1)
       ! David moved HPRNU and added HVC
       VISCOF=(FACT*0.5_SP*(VISCOF1*CC_HVC(IA)+VISCOF2*CC_HVC(IB)) + FM1*0.5_SP*(CC_HVC(IA)+CC_HVC(IB)))/HPRNU

       TXXIJ=(COFA1+COFA3)*VISCOF
       TYYIJ=(COFA6+COFA8)*VISCOF
       TXYIJ=0.5_SP*(COFA2+COFA4+COFA5+COFA7)*VISCOF
       FXX=DIJ*(TXXIJ*DLTYC(I)-TXYIJ*DLTXC(I))
       FYY=DIJ*(TXYIJ*DLTYC(I)-TYYIJ*DLTXC(I))

#    if defined (PLBC)
       TXYIJ=0.5_SP*(COFA5+COFA7)*VISCOF
       FXX=DIJ*(TXXIJ*DLTYC(I)-0.0_SP*DLTXC(I))
       FYY=DIJ*(TXYIJ*DLTYC(I)-0.0_SP*DLTXC(I))
#    endif

       !!UPWIND THE ADVECTIVE FLUX
       XADV=DIJ*UN*((1.0_SP-SIGN(1.0_SP,UN))*UIJ2+(1.0_SP+SIGN(1.0_SP,UN))*UIJ1)*0.5_SP
       YADV=DIJ*UN*((1.0_SP-SIGN(1.0_SP,UN))*VIJ2+(1.0_SP+SIGN(1.0_SP,UN))*VIJ1)*0.5_SP


       !!COMPUTE BOUNDARY FLUX AUGMENTERS
       TPA = FLOAT(1-ISBC(I))*EPOR(IA)
       TPB = FLOAT(1-ISBC(I))*EPOR(IB)


       !!ACCUMULATE THE FLUX
       XFLUX(IA,K)=XFLUX(IA,K)+XADV*TPA+FXX*TPA
       YFLUX(IA,K)=YFLUX(IA,K)+YADV*TPA+FYY*TPA
       XFLUX(IB,K)=XFLUX(IB,K)-XADV*TPB-FXX*TPB
       YFLUX(IB,K)=YFLUX(IB,K)-YADV*TPB-FYY*TPB


     END DO
#  if defined (WET_DRY)
    END IF
#  endif
   END DO
#  endif

#  if defined (SPHERICAL)
   CALL ADVECTION_EDGE_XY(XFLUX,YFLUX)
#  endif  

#  if defined (WET_DRY)
   DO I=1,N
#    if !defined (SEMI_IMPLICIT)
     ISWETTMP = ISWETCT(I)*ISWETC(I)
#    else
     ISWETTMP = ISWETCT(I)
#    endif
     DO K=1,KBM1
       XFLUX(I,K) = XFLUX(I,K)*ISWETTMP
       YFLUX(I,K) = YFLUX(I,K)*ISWETTMP
     END DO
   END DO
#  endif       	 


!
!--Boundary Conditions on Flux-------------------------------------------------|
!
   DO I=1,N
     IF(ISBCE(I) == 2)THEN
       DO K=1,KBM1
         XFLUX(I,K)=0.0_SP
         YFLUX(I,K)=0.0_SP
       END DO
     END IF
   END DO

!
!--Adjust Flux for Fresh Water Inflow------------------------------------------|
!

   IF(NUMQBC > 0) THEN
     IF(RIVER_INFLOW_LOCATION == 'node') THEN
       DO II=1,NUMQBC
         J=INODEQ(II)
         I1=NBVE(J,1)
         I2=NBVE(J,NTVE(J))
         DO K=1,KBM1
           VLCTYQ(II)=QDIS(II)/QAREA(II)
!           TEMP=0.5_SP*QDIS(II)*VQDIST(II,K)*VLCTYQ(II)
           TEMP=0.5_SP*QDIS(II)*VQDIST(II,K)*VQDIST(II,K)*VLCTYQ(II)/DZ(J,K)
!           XFLUX(I1,K)=XFLUX(I1,K)-TEMP/DZ(J,K)*COS(ANGLEQ(II))
!           XFLUX(I2,K)=XFLUX(I2,K)-TEMP/DZ(J,K)*COS(ANGLEQ(II))
!           YFLUX(I1,K)=YFLUX(I1,K)-TEMP/DZ(J,K)*SIN(ANGLEQ(II))
!           YFLUX(I2,K)=YFLUX(I2,K)-TEMP/DZ(J,K)*SIN(ANGLEQ(II))
           XFLUX(I1,K)=XFLUX(I1,K)-TEMP*COS(ANGLEQ(II))
           XFLUX(I2,K)=XFLUX(I2,K)-TEMP*COS(ANGLEQ(II))
           YFLUX(I1,K)=YFLUX(I1,K)-TEMP*SIN(ANGLEQ(II))
           YFLUX(I2,K)=YFLUX(I2,K)-TEMP*SIN(ANGLEQ(II))
         END DO
       END DO
     ELSE IF(RIVER_INFLOW_LOCATION == 'edge') THEN
       DO II=1,NUMQBC
         I1=ICELLQ(II)
         DO K=1,KBM1
           VLCTYQ(II)=QDIS(II)/QAREA(II)
!           TEMP=QDIS(II)*VQDIST(II,K)*VLCTYQ(II)
           TEMP=QDIS(II)*VQDIST(II,K)*VQDIST(II,K)*VLCTYQ(II)/DZ1(I1,K)
!           XFLUX(I1,K)=XFLUX(I1,K)-TEMP/DZ1(I1,K)*COS(ANGLEQ(II))
!           YFLUX(I1,K)=YFLUX(I1,K)-TEMP/DZ1(I1,K)*SIN(ANGLEQ(II))
           XFLUX(I1,K)=XFLUX(I1,K)-TEMP*COS(ANGLEQ(II))
           YFLUX(I1,K)=YFLUX(I1,K)-TEMP*SIN(ANGLEQ(II))
         END DO
       END DO
     END IF
   END IF

   if(dbg_set(dbg_sbr)) write(ipt,*) "End: advection_edge_gcy.F"

   RETURN
   END SUBROUTINE ADVECTION_EDGE_GCY
!==============================================================================|