adv_t.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 Temperature         |
!==============================================================================|

SUBROUTINE adv_t               

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

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

  IMPLICIT NONE
  REAL(SP), DIMENSION(0:MT,KB)      :: XFLUX,XFLUX_ADV,RF
  REAL(SP), DIMENSION(0:MT)         :: PUPX,PUPY,PVPX,PVPY  
  REAL(SP), DIMENSION(0:MT)         :: PTPX,PTPY,PTPXD,PTPYD,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,TTIME,ZDEP
  REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF,EXFLUX,TEMP,STPOINT,STPOINT1,STPOINT2
  REAL(SP) :: FACT,FM1
  INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,K,JTMP,JJ,II
  REAL(SP) :: T1MIN, T1MAX, T2MIN, T2MAX

!!$#  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_t"
  !------------------------------------------------------------------------------!

  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

  !
  !--Add the Shortwave Radiation Body Force--------------------------------------!
  !
  rf = 0.0_sp
  IF(heating_on) THEN
     SELECT CASE(heating_type)
     CASE('body')

        ! REMOVED START TIME FOR BODY HEATING
        !     TTIME=THOUR
        !     if(ttime < thour_hs) then
        !       RF=0.0_SP
        !     else

        DO  k=1,kbm1
           DO  i=1,m
              !         IF(TTIME < THOUR_HS)THEN
              !           RF(I,K)=0.0_SP
              !         ELSE
              ! REMOVED DEPTH CHECK FROM ECOM-SI, FVCOM DOES NOT NEED IT!
              !          ZEDP = 0.0_SP
              !           IF(DT(I) > 0.0_SP) THEN
              zdep=0.5_sp*(z(i,k)+z(i,k+1))*dt(i)
              !           END IF
              rf(i,k)=-swrad(i)*((rheat/zeta1)*exp(zdep/zeta1) &
                   +((1-rheat)/zeta2)*exp(zdep/zeta2))*dt(i)
              !         END IF
           END DO
        END DO
        !     endif

     CASE('flux')
        rf = 0.0_sp
     CASE DEFAULT
        CALL fatal_error('The surface heating type is set incorrectly:',&
             & trim(heating_type))
     END SELECT
  END IF

!--ADJUST VOLUME'S HEAT CONTENT FOR EVAPORATION AND PRECIPITATION ------------------!
  IF (precipitation_on) THEN
     rf(:,1)=rf(:,1)+rofvros*(qevap+qprec)*t1(:,1)
  END IF

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

  DO k=1,kbm1
     ptpx  = 0.0_sp
     ptpy  = 0.0_sp
     ptpxd = 0.0_sp
     ptpyd = 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*(t1(i,k)+t1(i2,k)-tmean1(i,k)-tmean1(i2,k))
            ff1=0.5_sp*(t1(i,k)+t1(i2,k))
           ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
            ffd=0.5_sp*(t1(i1,k)+t1(i,k)-tmean1(i1,k)-tmean1(i,k))
            ff1=0.5_sp*(t1(i1,k)+t1(i,k))
       ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
            ffd=0.5_sp*(t1(i,k)+t1(i,k)-tmean1(i,k)-tmean1(i,k))
            ff1=0.5_sp*(t1(i,k)+t1(i,k))
       ELSE
            ffd=0.5_sp*(t1(i1,k)+t1(i2,k)-tmean1(i1,k)-tmean1(i2,k))
            ff1=0.5_sp*(t1(i1,k)+t1(i2,k))
       END IF 
         ELSE
           IF(backward_step==1)THEN
              IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
               ffd=0.5_sp*((t0(i,k)+t1(i,k))*0.5+(t0(i2,k)+t1(i2,k))*0.5-tmean1(i,k)-tmean1(i2,k))
               ff1=0.5_sp*((t0(i,k)+t1(i,k))*0.5+(t0(i2,k)+t1(i2,k))*0.5)
              ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((t0(i1,k)+t1(i1,k))*0.5+(t0(i,k)+t1(i,k))*0.5-tmean1(i1,k)-tmean1(i,k))
               ff1=0.5_sp*((t0(i1,k)+t1(i1,k))*0.5+(t0(i,k)+t1(i,k))*0.5)
          ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((t0(i,k)+t1(i,k))*0.5+(t0(i,k)+t1(i,k))*0.5-tmean1(i,k)-tmean1(i,k))
               ff1=0.5_sp*((t0(i,k)+t1(i,k))*0.5+(t0(i,k)+t1(i,k))*0.5)
          ELSE
               ffd=0.5_sp*((t0(i1,k)+t1(i1,k))*0.5+(t0(i2,k)+t1(i2,k))*0.5-tmean1(i1,k)-tmean1(i2,k))
               ff1=0.5_sp*((t0(i1,k)+t1(i1,k))*0.5+(t0(i2,k)+t1(i2,k))*0.5)
          END IF 
           ELSEIF(backward_step==2)THEN
              IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
               ffd=0.5_sp*((t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp+(t2(i2,k)+t0(i2,k)+t1(i2,k))/3.0_sp-tmean1(i,k)-tmean1(i2,k))
               ff1=0.5_sp*((t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp+(t2(i2,k)+t0(i2,k)+t1(i2,k))/3.0_sp)
              ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((t2(i1,k)+t0(i1,k)+t1(i1,k))/3.0_sp+(t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp-tmean1(i1,k)-tmean1(i,k))
               ff1=0.5_sp*((t2(i1,k)+t0(i1,k)+t1(i1,k))/3.0_sp+(t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp)
          ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp+(t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp-tmean1(i,k)-tmean1(i,k))
               ff1=0.5_sp*((t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp+(t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp)
          ELSE
               ffd=0.5_sp*((t2(i1,k)+t0(i1,k)+t1(i1,k))/3.0_sp+(t2(i2,k)+t0(i2,k)+t1(i2,k))/3.0_sp-tmean1(i1,k)-tmean1(i2,k))
               ff1=0.5_sp*((t2(i1,k)+t0(i1,k)+t1(i1,k))/3.0_sp+(t2(i2,k)+t0(i2,k)+t1(i2,k))/3.0_sp)
          END IF 
            ENDIF
        ENDIF
!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
!!$           ! END ERROR
!!$           
!!$           PTPX(I)=PTPX(I)+FF1*TYPI
!!$           PTPY(I)=PTPY(I)+FF1*TXPI
!!$           PTPXD(I)=PTPXD(I)+FFD*TYPI
!!$           PTPYD(I)=PTPYD(I)+FFD*TXPI
!!$#        else
!!$           PTPX(I)=PTPX(I)+FF1*(VY(I1)-VY(I2))
!!$           PTPY(I)=PTPY(I)+FF1*(VX(I2)-VX(I1))
!!$           PTPXD(I)=PTPXD(I)+FFD*(VY(I1)-VY(I2))
!!$           PTPYD(I)=PTPYD(I)+FFD*(VX(I2)-VX(I1))
!!$#        endif

           
           ptpx(i)=ptpx(i)+ff1*dltytrie(i,j)
           ptpy(i)=ptpy(i)+ff1*dltxtrie(i,j)
           ptpxd(i)=ptpxd(i)+ffd*dltytrie(i,j)
           ptpyd(i)=ptpyd(i)+ffd*dltxtrie(i,j)

        END DO
! gather all neighboring control volumes connecting at dam node 
        ptpx(i)=ptpx(i)/art2(i)
        ptpy(i)=ptpy(i)/art2(i)
        ptpxd(i)=ptpxd(i)/art2(i)
        ptpyd(i)=ptpyd(i)/art2(i)

     END DO

     IF(k == kbm1)THEN
        DO i=1,m
           pfpxb(i) = ptpx(i)
           pfpyb(i) = ptpy(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
!!$        FIJ1=T1(IA,K)+DXA*PTPX(IA)+DYA*PTPY(IA)
!!$        FIJ2=T1(IB,K)+DXB*PTPX(IB)+DYB*PTPY(IB)

!J. Ge for tracer advection
        IF(backward_advection==.false.)THEN
          fij1=t1(ia,k)+dltxncve(i,1)*ptpx(ia)+dltyncve(i,1)*ptpy(ia)
          fij2=t1(ib,k)+dltxncve(i,2)*ptpx(ib)+dltyncve(i,2)*ptpy(ib)
        ELSE
          IF(backward_step==1)THEN
            fij1=(t0(ia,k)+t1(ia,k))*0.5+dltxncve(i,1)*ptpx(ia)+dltyncve(i,1)*ptpy(ia)
            fij2=(t0(ib,k)+t1(ib,k))*0.5+dltxncve(i,2)*ptpx(ib)+dltyncve(i,2)*ptpy(ib)
          ELSEIF(backward_step==2)THEN
            fij1=(t2(ia,k)+t0(ia,k)+t1(ia,k))/3.0_sp+dltxncve(i,1)*ptpx(ia)+dltyncve(i,1)*ptpy(ia)
            fij2=(t2(ia,k)+t0(ib,k)+t1(ib,k))/3.0_sp+dltxncve(i,2)*ptpx(ib)+dltyncve(i,2)*ptpy(ib)
          ENDIF
        ENDIF
!J. Ge for tracer advection

        t1min=minval(t1(nbsn(ia,1:ntsn(ia)-1),k))
        t1min=min(t1min, t1(ia,k))
        t1max=maxval(t1(nbsn(ia,1:ntsn(ia)-1),k))
        t1max=max(t1max, t1(ia,k))
        t2min=minval(t1(nbsn(ib,1:ntsn(ib)-1),k))
        t2min=min(t2min, t1(ib,k))
        t2max=maxval(t1(nbsn(ib,1:ntsn(ib)-1),k))
        t2max=max(t2max, t1(ib,k))
        IF(fij1 < t1min) fij1=t1min
        IF(fij1 > t1max) fij1=t1max
        IF(fij2 < t2min) fij2=t2min
        IF(fij2 > t2max) fij2=t2max

        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*(ptpxd(ia)+ptpxd(ib))*viscof
        tyy=0.5_sp*(ptpyd(ia)+ptpyd(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 !! K LOOP


  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)*(t1(i,k)*dz(i,k+1)+t1(i,k+1)*dz(i,k))/   &
                   (dz(i,k)+dz(i,k+1))
             ELSE IF(k == kbm1) THEN
              temp= wts(i,k)*(t1(i,k)*dz(i,k-1)+t1(i,k-1)*dz(i,k))/(dz(i,k)+dz(i,k-1))
             ELSE
              temp= wts(i,k)*(t1(i,k)*dz(i,k-1)+t1(i,k-1)*dz(i,k))/(dz(i,k)+dz(i,k-1))-&
                   wts(i,k+1)*(t1(i,k)*dz(i,k+1)+t1(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)*((t0(i,k)+t1(i,k))*0.5*dz(i,k+1)+(t0(i,k+1)+t1(i,k+1))*0.5*dz(i,k))/   &
                   (dz(i,k)+dz(i,k+1))
               ELSE IF(k == kbm1) THEN
                temp= wts(i,k)*((t0(i,k)+t1(i,k))*0.5*dz(i,k-1)+(t0(i,k-1)+t1(i,k-1))*0.5*dz(i,k))/(dz(i,k)+dz(i,k-1))
               ELSE
                 temp= wts(i,k)*((t0(i,k)+t1(i,k))*0.5*dz(i,k-1)+(t0(i,k-1)+t1(i,k-1))*0.5*dz(i,k))/(dz(i,k)+dz(i,k-1))-&
                   wts(i,k+1)*((t0(i,k)+t1(i,k))*0.5*dz(i,k+1)+(t0(i,k+1)+t1(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)*((t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp*dz(i,k+1)+(t2(i,k+1)+t0(i,k+1)+t1(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)*((t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp*dz(i,k-1)+(t2(i,k-1)+t0(i,k-1)+t1(i,k-1))/3.0_sp*dz(i,k))/(dz(i,k)+dz(i,k-1))
               ELSE
                 temp= wts(i,k)*((t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp*dz(i,k-1)+(t2(i,k-1)+t0(i,k-1)+t1(i,k-1))/3.0_sp*dz(i,k))/(dz(i,k)+dz(i,k-1))-&
                   wts(i,k+1)*((t2(i,k)+t0(i,k)+t1(i,k))/3.0_sp*dz(i,k+1)+(t2(i,k+1)+t0(i,k+1)+t1(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  !!K 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=tdis(j)  
              DO k=1,kbm1
              !   STPOINT    = T1(JJ,K)
                 !             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=tdis(j) !!ASK LIU SHOULD THIS BE STPOINT1(J1)/STPOINT2(J2)
              DO k=1,kbm1
                 !STPOINT1 = T1(J1,K)
                 !STPOINT2 = T1(J2,K)
                 !             XFLUX(J1,K)=XFLUX(J1,K)-  &
                 !                         QDIS(J)*RDISQ(J,1)*VQDIST(J,K)*STPOINT1/DZ1(J1,K)
                 !             XFLUX(J2,K)=XFLUX(J2,K)-  &
                 !                         QDIS(J)*RDISQ(J,2)*VQDIST(J,K)*STPOINT2/DZ1(J2,K)
                 xflux(j1,k)=xflux(j1,k)-qdis(j)*rdisq(j,1)*vqdist(j,k)*stpoint   !1
                 xflux(j2,k)=xflux(j2,k)-qdis(j)*rdisq(j,2)*vqdist(j,k)*stpoint   !2
              END DO
           END DO
        END IF
     END IF
  END IF
  !---------------------------------------------------------------------


  ! APPLY GROUND WATER TEMPERATURE FORCING
  IF(groundwater_on .and. groundwater_temp_on)THEN
     DO i=1,m
        xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*bfwtmp(i)
     END DO
  ELSEIF(groundwater_on) THEN
     DO i=1,m
 !J. Ge for tracer advection
        IF(backward_advection==.false.)THEN
          xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*t1(i,kbm1)
        ELSE
          IF(backward_step==1)THEN
            xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*(t0(i,kbm1)+t1(i,kbm1))*0.5
          ELSEIF(backward_step==2)THEN
            xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*(t2(i,kbm1)+t0(i,kbm1)+t1(i,kbm1))/3.0_sp
          ENDIF
        ENDIF
!J. Ge for tracer advection
     END DO
  END IF


  !
  !--Update Temperature----------------------------------------------------------!
  !


  DO i=1,m
     IF(iswetn(i)*iswetnt(i) == 1 )THEN
        DO k=1,kbm1
           xflux(i,k) = xflux(i,k) - rf(i,k)*art1(i)    !/DZ(I,K)
           !       TF1(I,K)=(T1(I,K)-XFLUX(I,K)/ART1(I)*(DTI/DT(I)))*(DT(I)/DTFA(I))

           tf1(i,k)=(t1(i,k)-xflux(i,k)/art1(i)*(dti/(dt(i)*dz(i,k))))*(dt(i)/dtfa(i))

        END DO
     ELSE
        DO k=1,kbm1
           tf1(i,k)=t1(i,k)
        END DO
     END IF
  END DO



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

END SUBROUTINE adv_t
!==============================================================================|