adv_s.f90

Go to the documentation of this file.

!/===========================================================================/
! 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$
!/===========================================================================/


!==============================================================================|
!   Calculate Advection and Horizontal Diffusion Terms for Salinity            |
!==============================================================================|

SUBROUTINE adv_s               

  !------------------------------------------------------------------------------|

  USE all_vars
  USE mod_utils
  USE bcs
  USE mod_obcs
  USE mod_par
  USE mod_wd
  USE mod_spherical
  USE mod_northpole


  IMPLICIT NONE
  REAL(SP), DIMENSION(0:MT,KB)     :: XFLUX,XFLUX_ADV
  REAL(SP), DIMENSION(0:MT)        :: PUPX,PUPY,PVPX,PVPY  
  REAL(SP), DIMENSION(0:MT)        :: PSPX,PSPY,PSPXD,PSPYD,VISCOFF
  REAL(SP), DIMENSION(3*(NT),KBM1) :: DTIJ 
  REAL(SP), DIMENSION(3*(NT),KBM1) :: UVN
  REAL(SP) :: UTMP,VTMP,SITAI,FFD,FF1 !,X11,Y11,X22,Y22,X33,Y33 !,TMP1,TMP2 !,XI,YI
  REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2,UN
  REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF,EXFLUX,TEMP,STPOINT
  REAL(SP) :: FACT,FM1
  INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,K,JTMP,JJ,II
  REAL(SP) :: S1MIN, S1MAX, S2MIN, S2MAX

!!$#  if defined (SPHERICAL)
!!$  REAL(DP) :: TY,TXPI,TYPI
!!$  REAL(DP) :: XTMP1,XTMP
!!$  REAL(DP) :: X1_DP,Y1_DP,X2_DP,Y2_DP,XII,YII
!!$  REAL(DP) :: X11_TMP,Y11_TMP,X33_TMP,Y33_TMP
!!$  REAL(DP) :: VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP
!!$  REAL(DP) :: TXPI_TMP,TYPI_TMP
!!$#  endif





  IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"Start: adv_s"
  !------------------------------------------------------------------------------!


  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 Fluxes-----------------------------------------------------------!
  !
  xflux     = 0.0_sp
  xflux_adv = 0.0_sp

  !
  !--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
  !
!!#  if !defined (1)
  DO i=1,ncv
     i1=ntrg(i)
     !     DTIJ(I)=DT1(I1)
     DO k=1,kbm1
       dtij(i,k) = dt1(i1)*dz1(i1,k)
       ! USE U,V
       uvn(i,k)  = v(i1,k)*dltxe(i) - u(i1,k)*dltye(i)
     END DO
  END DO
!!#  else
!!  DO I=1,NCV
!!     I1=NTRG(I)
!!     !     DTIJ(I)=DT1(I1)
!!     DO K=1,KBM1
!!       DTIJ(I,K) = DT1(I1)*DZ1(I1,K)
!!       ! USE US,VS
!!       UVN(I,K) = VS(I1,K)*DLTXE(I) - US(I1,K)*DLTYE(I)
!!#      if defined (SEMI_IMPLICIT)
!!       DTIJ1(I,K) = D1(I1)*DZ1(I1,K)
!!       UVN1(I,K) = VF(I1,K)*DLTXE(I) - UF(I1,K)*DLTYE(I)
!!#      endif
!!     END DO
!!  END DO
!!#  endif

  !
  !--Calculate the Advection and Horizontal Diffusion Terms----------------------!
  !

  DO k=1,kbm1
     pspx  = 0.0_sp 
     pspy  = 0.0_sp 
     pspxd = 0.0_sp 
     pspyd = 0.0_sp
     DO i=1,m
        DO j=1,ntsn(i)-1
           i1=nbsn(i,j)
           i2=nbsn(i,j+1)

!J. Ge for tracer advection
           IF(backward_advection==.false.)THEN

            IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
             ffd=0.5_sp*(s1(i,k)+s1(i2,k)-smean1(i,k)-smean1(i2,k))
             ff1=0.5_sp*(s1(i,k)+s1(i2,k))
            ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
             ffd=0.5_sp*(s1(i1,k)+s1(i,k)-smean1(i1,k)-smean1(i,k))
             ff1=0.5_sp*(s1(i1,k)+s1(i,k))
            ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
             ffd=0.5_sp*(s1(i,k)+s1(i,k)-smean1(i,k)-smean1(i,k))
             ff1=0.5_sp*(s1(i,k)+s1(i,k))
            ELSE
             ffd=0.5_sp*(s1(i1,k)+s1(i2,k)-smean1(i1,k)-smean1(i2,k))
             ff1=0.5_sp*(s1(i1,k)+s1(i2,k))
            END IF

           ELSE
         
            IF(backward_step==1)THEN
              IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
               ffd=0.5_sp*((s0(i,k)+s1(i,k))*0.5+(s0(i2,k)+s1(i2,k))*0.5-smean1(i,k)-smean1(i2,k))
               ff1=0.5_sp*((s0(i,k)+s1(i,k))*0.5+(s0(i2,k)+s1(i2,k))*0.5)
              ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((s0(i1,k)+s1(i1,k))*0.5+(s0(i,k)+s1(i,k))*0.5-smean1(i1,k)-smean1(i,k))
               ff1=0.5_sp*((s0(i1,k)+s1(i1,k))*0.5+(s0(i,k)+s1(i,k))*0.5)
          ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((s0(i,k)+s1(i,k))*0.5+(s0(i,k)+s1(i,k))*0.5-smean1(i,k)-smean1(i,k))
               ff1=0.5_sp*((s0(i,k)+s1(i,k))*0.5+(s0(i,k)+s1(i,k))*0.5)
          ELSE
               ffd=0.5_sp*((s0(i1,k)+s1(i1,k))*0.5+(s0(i2,k)+s1(i2,k))*0.5-smean1(i1,k)-smean1(i2,k))
               ff1=0.5_sp*((s0(i1,k)+s1(i1,k))*0.5+(s0(i2,k)+s1(i2,k))*0.5)
          END IF 
            ELSEIF(backward_step==2)THEN
              IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
               ffd=0.5_sp*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp+(s2(i2,k)+s0(i2,k)+s1(i2,k))/3.0_sp-smean1(i,k)-smean1(i2,k))
               ff1=0.5_sp*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp+(s2(i2,k)+s0(i2,k)+s1(i2,k))/3.0_sp)
              ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((s2(i1,k)+s0(i1,k)+s1(i1,k))/3.0_sp+(s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp-smean1(i1,k)-smean1(i,k))
               ff1=0.5_sp*((s2(i1,k)+s0(i1,k)+s1(i1,k))/3.0_sp+(s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp)
          ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp+(s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp-smean1(i,k)-smean1(i,k))
               ff1=0.5_sp*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp+(s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp)
          ELSE
               ffd=0.5_sp*((s2(i1,k)+s0(i1,k)+s1(i1,k))/3.0_sp+(s2(i2,k)+s0(i2,k)+s1(i2,k))/3.0_sp-smean1(i1,k)-smean1(i2,k))
               ff1=0.5_sp*((s2(i1,k)+s0(i1,k)+s1(i1,k))/3.0_sp+(s2(i2,k)+s0(i2,k)+s1(i2,k))/3.0_sp)
          END IF 
            ENDIF

           END IF
!J. Ge for tracer advection
     

!!$#        if defined (SPHERICAL)
!!$           XTMP  = VX(I2)*TPI-VX(I1)*TPI
!!$           XTMP1 = VX(I2)-VX(I1)
!!$           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
!!$           TXPI=XTMP*COS(DEG2RAD*VY(I))
!!$           TYPI=(VY(I1)-VY(I2))*TPI
!!$           ! ERROR HERE
!!$!#    if defined (NORTHPOLE)
!!$           IF(NODE_NORTHAREA(I) == 1)THEN
!!$              VX1_TMP = REARTH * COS(VY(I1)*DEG2RAD) * COS(VX(I1)*DEG2RAD) &
!!$                   * 2._SP /(1._SP+sin(VY(I1)*DEG2RAD))
!!$              VY1_TMP = REARTH * COS(VY(I1)*DEG2RAD) * SIN(VX(I1)*DEG2RAD) &
!!$                   * 2._SP /(1._SP+sin(VY(I1)*DEG2RAD))
!!$
!!$              VX2_TMP = REARTH * COS(VY(I2)*DEG2RAD) * COS(VX(I2)*DEG2RAD) &
!!$                   * 2._SP /(1._SP+sin(VY(I2)*DEG2RAD))
!!$              VY2_TMP = REARTH * COS(VY(I2)*DEG2RAD) * SIN(VX(I2)*DEG2RAD) &
!!$                   * 2._SP /(1._SP+sin(VY(I2)*DEG2RAD))
!!$
!!$              TXPI = (VX2_TMP-VX1_TMP)/(2._SP /(1._SP+sin(VY(I)*DEG2RAD)))
!!$              TYPI = (VY1_TMP-VY2_TMP)/(2._SP /(1._SP+sin(VY(I)*DEG2RAD)))
!!$              IF(I /= NODE_NORTHPOLE)THEN
!!$                 TXPI_TMP = TYPI*COS(VX(I)*DEG2RAD)-TXPI*SIN(VX(I)*DEG2RAD)
!!$                 TYPI_TMP = TXPI*COS(VX(I)*DEG2RAD)+TYPI*SIN(VX(I)*DEG2RAD)
!!$                 TYPI_TMP = -TYPI_TMP
!!$
!!$                 TXPI = TXPI_TMP
!!$                 TYPI = TYPI_TMP
!!$              END IF
!!$           END IF
!!$# endif
!!$           ! END ERROR
!!$#        else
!!$           PSPX(I)=PSPX(I)+FF1*(VY(I1)-VY(I2))
!!$           PSPY(I)=PSPY(I)+FF1*(VX(I2)-VX(I1))
!!$           PSPXD(I)=PSPXD(I)+FFD*(VY(I1)-VY(I2))
!!$           PSPYD(I)=PSPYD(I)+FFD*(VX(I2)-VX(I1))
!!$#        endif

           
           pspx(i)=pspx(i)+ff1*dltytrie(i,j)
           pspy(i)=pspy(i)+ff1*dltxtrie(i,j)
           pspxd(i)=pspxd(i)+ffd*dltytrie(i,j)
           pspyd(i)=pspyd(i)+ffd*dltxtrie(i,j)
           

        END DO
! gather all neighboring control volumes connecting at dam node 

        pspx(i)=pspx(i)/art2(i)
        pspy(i)=pspy(i)/art2(i)
        pspxd(i)=pspxd(i)/art2(i)
        pspyd(i)=pspyd(i)/art2(i)

     END DO

     IF(k == kbm1)THEN
        DO i=1,m
           pfpxb(i) = pspx(i)
           pfpyb(i) = pspy(i)
        END DO
     END IF

     DO i = 1,m
        viscoff(i)=viscofh(i,k)
     END DO

     IF(k == kbm1) THEN
        ah_bottom(1:m) = (fact*viscoff(1:m) + fm1)*nn_hvc(1:m)
     END IF


     DO i=1,ncv_i
        ia=niec(i,1)
        ib=niec(i,2)


!!$        XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
!!$        YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
!!$#      if defined (SPHERICAL)
!!$        X1_DP=XIJE(I,1)
!!$        Y1_DP=YIJE(I,1)
!!$        X2_DP=XIJE(I,2)
!!$        Y2_DP=YIJE(I,2)
!!$        CALL ARCC(X2_DP,Y2_DP,X1_DP,Y1_DP,XII,YII)
!!$        XI=XII       
!!$        XTMP  = XI*TPI-VX(IA)*TPI
!!$        XTMP1 = XI-VX(IA)
!!$        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
!!$
!!$        DXA=XTMP*COS(DEG2RAD*VY(IA))    
!!$        DYA=(YI-VY(IA))*TPI
!!$        XTMP  = XI*TPI-VX(IB)*TPI
!!$        XTMP1 = XI-VX(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
!!$
!!$        DXB=XTMP*COS(DEG2RAD*VY(IB)) 
!!$        DYB=(YI-VY(IB))*TPI
!!$#      else
!!$        DXA=XI-VX(IA)
!!$        DYA=YI-VY(IA)
!!$        DXB=XI-VX(IB)
!!$        DYB=YI-VY(IB)
!!$#      endif

!J. Ge for tracer advection
        IF(backward_advection==.false.)THEN
          fij1=s1(ia,k)+dltxncve(i,1)*pspx(ia)+dltyncve(i,1)*pspy(ia)
          fij2=s1(ib,k)+dltxncve(i,2)*pspx(ib)+dltyncve(i,2)*pspy(ib)
        ELSE
          IF(backward_step==1)THEN
            fij1=(s0(ia,k)+s1(ia,k))*0.5+dltxncve(i,1)*pspx(ia)+dltyncve(i,1)*pspy(ia)
            fij2=(s0(ib,k)+s1(ib,k))*0.5+dltxncve(i,2)*pspx(ib)+dltyncve(i,2)*pspy(ib)
          ELSEIF(backward_step==2)THEN
            fij1=(s2(ia,k)+s0(ia,k)+s1(ia,k))/3.0_sp+dltxncve(i,1)*pspx(ia)+dltyncve(i,1)*pspy(ia)
            fij2=(s2(ib,k)+s0(ib,k)+s1(ib,k))/3.0_sp+dltxncve(i,2)*pspx(ib)+dltyncve(i,2)*pspy(ib)
          ENDIF
        ENDIF
!J. Ge for tracer advection
        s1min=minval(s1(nbsn(ia,1:ntsn(ia)-1),k))
        s1min=min(s1min, s1(ia,k))
        s1max=maxval(s1(nbsn(ia,1:ntsn(ia)-1),k))
        s1max=max(s1max, s1(ia,k))
        s2min=minval(s1(nbsn(ib,1:ntsn(ib)-1),k))
        s2min=min(s2min, s1(ib,k))
        s2max=maxval(s1(nbsn(ib,1:ntsn(ib)-1),k))
        s2max=max(s2max, s1(ib,k))
        IF(fij1 < s1min) fij1=s1min
        IF(fij1 > s1max) fij1=s1max
        IF(fij2 < s2min) fij2=s2min
        IF(fij2 > s2max) fij2=s2max

        un=uvn(i,k)

        !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
        
        ! David moved HPRNU and added HVC
        viscof=(fact*0.5_sp*(viscoff(ia)*nn_hvc(ia)+viscoff(ib)*nn_hvc(ib)) + fm1*0.5_sp*(nn_hvc(ia)+nn_hvc(ib)))
        
        txx=0.5_sp*(pspxd(ia)+pspxd(ib))*viscof
        tyy=0.5_sp*(pspyd(ia)+pspyd(ib))*viscof

        fxx=-dtij(i,k)*txx*dltye(i)
        fyy= dtij(i,k)*tyy*dltxe(i)

        exflux=-un*dtij(i,k)* &
             ((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy

        xflux(ia,k)=xflux(ia,k)+exflux
        xflux(ib,k)=xflux(ib,k)-exflux

        xflux_adv(ia,k)=xflux_adv(ia,k)+(exflux-fxx-fyy)
        xflux_adv(ib,k)=xflux_adv(ib,k)-(exflux-fxx-fyy)



     END DO



  END DO !!SIGMA LOOP

  !
  !-Accumulate Fluxes at Boundary Nodes
  !

  DO k=1,kbm1
     IF(iobcn > 0) THEN
        DO i=1,iobcn
           i1=i_obc_n(i)
           xflux_obc(i,k)=xflux_adv(i1,k)
        END DO
     END IF
  END DO

  !--Set Boundary Conditions-For Fresh Water Flux--------------------------------!
  !




  !--------------------------------------------------------------------
  !   The central difference scheme in vertical advection
  !--------------------------------------------------------------------
  DO i=1,m
     IF(iswetn(i)*iswetnt(i) == 1) THEN


        DO k=1, kbm1



!J. Ge for tracer advection
           IF(backward_advection==.false.)THEN
             IF(k == 1) THEN
              temp=-wts(i,k+1)*(s1(i,k)*dz(i,k+1)+s1(i,k+1)*dz(i,k))/   &
                   (dz(i,k)+dz(i,k+1))
             ELSE IF(k == kbm1) THEN
              temp= wts(i,k)*(s1(i,k)*dz(i,k-1)+s1(i,k-1)*dz(i,k))/(dz(i,k)+dz(i,k-1))
             ELSE
              temp= wts(i,k)*(s1(i,k)*dz(i,k-1)+s1(i,k-1)*dz(i,k))/(dz(i,k)+dz(i,k-1))-&
                   wts(i,k+1)*(s1(i,k)*dz(i,k+1)+s1(i,k+1)*dz(i,k))/(dz(i,k)+dz(i,k+1))
             END IF
           ELSE
             IF(backward_step==1)THEN
               IF(k == 1) THEN
                temp=-wts(i,k+1)*((s0(i,k)+s1(i,k))*0.5*dz(i,k+1)+(s0(i,k+1)+s1(i,k+1))*0.5*dz(i,k))/   &
                   (dz(i,k)+dz(i,k+1))
               ELSE IF(k == kbm1) THEN
                temp= wts(i,k)*((s0(i,k)+s1(i,k))*0.5*dz(i,k-1)+(s0(i,k-1)+s1(i,k-1))*0.5*dz(i,k))/(dz(i,k)+dz(i,k-1))
               ELSE
                 temp= wts(i,k)*((s0(i,k)+s1(i,k))*0.5*dz(i,k-1)+(s0(i,k-1)+s1(i,k-1))*0.5*dz(i,k))/(dz(i,k)+dz(i,k-1))-&
                   wts(i,k+1)*((s0(i,k)+s1(i,k))*0.5*dz(i,k+1)+(s0(i,k+1)+s1(i,k+1))*0.5*dz(i,k))/(dz(i,k)+dz(i,k+1))
               END IF
             ELSEIF(backward_step==2)THEN
               IF(k == 1) THEN
                temp=-wts(i,k+1)*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp*dz(i,k+1)+(s2(i,k+1)+s0(i,k+1)+s1(i,k+1))/3.0_sp*dz(i,k))/   &
                   (dz(i,k)+dz(i,k+1))
               ELSE IF(k == kbm1) THEN
                temp= wts(i,k)*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp*dz(i,k-1)+(s2(i,k-1)+s0(i,k-1)+s1(i,k-1))/3.0_sp*dz(i,k))/(dz(i,k)+dz(i,k-1))
               ELSE
                 temp= wts(i,k)*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp*dz(i,k-1)+(s2(i,k-1)+s0(i,k-1)+s1(i,k-1))/3.0_sp*dz(i,k))/(dz(i,k)+dz(i,k-1))-&
                   wts(i,k+1)*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp*dz(i,k+1)+(s2(i,k+1)+s0(i,k+1)+s1(i,k+1))/3.0_sp*dz(i,k))/(dz(i,k)+dz(i,k+1))
               END IF
             ENDIF
           ENDIF
!J. Ge for tracer advection


           IF(isonb(i) == 2) THEN
              !         XFLUX(I,K)=TEMP*ART1(I)/DZ(I,K)
              xflux(i,k)=temp*art1(i)
           ELSE
              !         XFLUX(I,K)=XFLUX(I,K)+TEMP*ART1(I)/DZ(I,K)
              xflux(i,k)=xflux(i,k)+temp*art1(i)

           END IF
        ENDDO
     END IF
  END DO  !! SIGMA LOOP

  !
  !--Set Boundary Conditions-For Fresh Water Flux--------------------------------!
  !
  IF(river_ts_setting == 'calculated') THEN
     IF(river_inflow_location == 'node') THEN
        IF(numqbc > 0) THEN
           DO j=1,numqbc
              jj=inodeq(j)
              stpoint=sdis(j)
              DO k=1,kbm1
                 !             XFLUX(JJ,K)=XFLUX(JJ,K) - QDIS(J)*VQDIST(J,K)*STPOINT/DZ(JJ,K)
                 xflux(jj,k)=xflux(jj,k) - qdis(j)*vqdist(j,k)*stpoint
              END DO
           END DO
        END IF
     ELSE IF(river_inflow_location == 'edge') THEN
        IF(numqbc > 0) THEN
           DO j=1,numqbc
              j1=n_icellq(j,1)
              j2=n_icellq(j,2)
              stpoint=sdis(j) !!ASK LIU SHOULD THIS BE STPOINT1(J1)/STPOINT2(J2)
              DO k=1,kbm1
                 !             XFLUX(J1,K)=XFLUX(J1,K)-   &
                 !                         QDIS(J)*RDISQ(J,1)*VQDIST(J,K)*STPOINT/DZ1(J1,K)
                 !             XFLUX(J2,K)=XFLUX(J2,K)-   &
                 !                         QDIS(J)*RDISQ(J,2)*VQDIST(J,K)*STPOINT/DZ1(J2,K)
                 xflux(j1,k)=xflux(j1,k)-qdis(j)*rdisq(j,1)*vqdist(j,k)*stpoint
                 xflux(j2,k)=xflux(j2,k)-qdis(j)*rdisq(j,2)*vqdist(j,k)*stpoint
              END DO
           END DO
        END IF
     END IF
  END IF
!---------------------------------------------------------------------


  ! APPLY GROUND WATER SALINITY FORCING
  IF(groundwater_on .and. groundwater_salt_on)THEN
     DO i=1,m
        xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*bfwslt(i)
     END DO
  ELSEIF(groundwater_on) THEN
     DO i=1,m
!        XFLUX(I,KBM1)=XFLUX(I,KBM1)-BFWDIS(I)*S1(I,KBM1)/DZ(I,KBM1)
!J. Ge for tracer advection
        IF(backward_advection==.false.)THEN
          xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*s1(i,kbm1)
        ELSE
          IF(backward_step==1)THEN
            xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*(s0(i,kbm1)+s1(i,kbm1))*0.5
          ELSEIF(backward_step==2)THEN
            xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*(s2(i,kbm1)+s0(i,kbm1)+s1(i,kbm1))/3.0_sp
          ENDIF
        ENDIF
!J. Ge for tracer advection
     END DO
  END IF


  !--Update Salinity-------------------------------------------------------------!
  !
  DO i=1,m
     IF(iswetn(i)*iswetnt(i) == 1 )THEN
        DO k=1,kbm1
           !       SF1(I,K)=(S1(I,K)-XFLUX(I,K)/ART1(I)*(DTI/DT(I)))*(DT(I)/DTFA(I))

           sf1(i,k)=(s1(i,k)-xflux(i,k)/art1(i)*(dti/(dt(i)*dz(i,k))))*(dt(i)/dtfa(i))

        END DO

     ELSE ! IF NOT WET S STAYS THE SAME
        DO k=1,kbm1
           sf1(i,k)=s1(i,k)
        END DO
     END IF
  END DO


  IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"End: adv_s"

END SUBROUTINE adv_s
!==============================================================================|