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

!==============================================================================|
   SUBROUTINE advection_edge_gcn(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                          |
!==============================================================================|

   USE all_vars
   USE mod_utils
   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_TMP
   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 (THIN_DAM)
   REAL(SP) :: A1UIA1,A1UIA2,A1UIA3,A1UIA4,A2UIA1,A2UIA2,A2UIA3,A2UIA4
   REAL(SP) :: A1UIB1,A1UIB2,A1UIB3,A1UIB4,A2UIB1,A2UIB2,A2UIB3,A2UIB4
   INTEGER  :: J11,J12,J21,J22,E1,E2,ISBCE1,ISBC_TMP,IB_TMP
   LOGICAL  :: ISMATCH
!#  endif

   REAL(SP) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
!------------------------------------------------------------------------------|
   
   if(dbg_set(dbg_sbr)) write(ipt,*) "Start: advection_edge_gcn.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----------------------|
!

   DO i=1,ne
     ia=iec(i,1)
     ib=iec(i,2)

    IF(iswetct(ia)*iswetc(ia) == 1 .OR. iswetct(ib)*iswetc(ib) == 1)THEN
     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)
!#    if defined (THIN_DAM)
!     IF(IB==0.AND.E_DAM_MATCH(IA)/=0)XIJB=DLTXNE_DAM_MATCH(I)
!     IF(IB==0.AND.E_DAM_MATCH(IA)/=0)YIJB=DLTYNE_DAM_MATCH(I)
!#    endif
!#    else
!     XIJA=XIJC(I)-XC(IA)
!     YIJA=YIJC(I)-YC(IA)
!     XIJB=XIJC(I)-XC(IB)
!     YIJB=YIJC(I)-YC(IB)
!#    if defined (THIN_DAM)
!     IF(IB==0.AND.E_DAM_MATCH(IA)/=0)XIJB=XIJC(I)-XC(E_DAM_MATCH(IA))
!     IF(IB==0.AND.E_DAM_MATCH(IA)/=0)YIJB=YIJC(I)-YC(E_DAM_MATCH(IA))
!#    endif
!#    endif

     DO k=1,kbm1

     xija=xijc(i)-xc(ia)
     yija=yijc(i)-yc(ia)
     xijb=xijc(i)-xc(ib)
     yijb=yijc(i)-yc(ib)

       ib_tmp = ib
       dij= 0.5_sp*(dt(j1)*dz(j1,k)+dt(j2)*dz(j2,k))
!--------------Used for Dam Model by Jadon----------------------
       a1uia1 = a1u(ia,1)
       a1uia2 = a1u(ia,2)
       a1uia3 = a1u(ia,3)
       a1uia4 = a1u(ia,4)
       a2uia1 = a2u(ia,1)
       a2uia2 = a2u(ia,2)
       a2uia3 = a2u(ia,3)
       a2uia4 = a2u(ia,4)
       
       a1uib1 = a1u(ib_tmp,1)
       a1uib2 = a1u(ib_tmp,2)
       a1uib3 = a1u(ib_tmp,3)
       a1uib4 = a1u(ib_tmp,4)
       a2uib1 = a2u(ib_tmp,1)
       a2uib2 = a2u(ib_tmp,2)
       a2uib3 = a2u(ib_tmp,3)
       a2uib4 = a2u(ib_tmp,4)
!---------------------------------------------------------------
       !!FORM THE LEFT FLUX
       cofa1=a1uia1*u(ia,k)+a1uia2*u(k1,k)+a1uia3*u(k2,k)+a1uia4*u(k3,k)
       cofa2=a2uia1*u(ia,k)+a2uia2*u(k1,k)+a2uia3*u(k2,k)+a2uia4*u(k3,k)
       cofa5=a1uia1*v(ia,k)+a1uia2*v(k1,k)+a1uia3*v(k2,k)+a1uia4*v(k3,k)
       cofa6=a2uia1*v(ia,k)+a2uia2*v(k1,k)+a2uia3*v(k2,k)+a2uia4*v(k3,k)
!       COFA1=A1U(IA,1)*U(IA,K)+A1U(IA,2)*U(K1,K)+A1U(IA,3)*U(K2,K)+A1U(IA,4)*U(K3,K)
!       COFA2=A2U(IA,1)*U(IA,K)+A2U(IA,2)*U(K1,K)+A2U(IA,3)*U(K2,K)+A2U(IA,4)*U(K3,K)
!       COFA5=A1U(IA,1)*V(IA,K)+A1U(IA,2)*V(K1,K)+A1U(IA,3)*V(K2,K)+A1U(IA,4)*V(K3,K)
!       COFA6=A2U(IA,1)*V(IA,K)+A2U(IA,2)*V(K1,K)+A2U(IA,3)*V(K2,K)+A2U(IA,4)*V(K3,K)
       uij1=u(ia,k)+cofa1*xija+cofa2*yija
       vij1=v(ia,k)+cofa5*xija+cofa6*yija

       !!FORM THE RIGHT FLUX
       cofa3=a1uib1*u(ib_tmp,k)+a1uib2*u(k4,k)+a1uib3*u(k5,k)+a1uib4*u(k6,k)
       cofa4=a2uib1*u(ib_tmp,k)+a2uib2*u(k4,k)+a2uib3*u(k5,k)+a2uib4*u(k6,k)
       cofa7=a1uib1*v(ib_tmp,k)+a1uib2*v(k4,k)+a1uib3*v(k5,k)+a1uib4*v(k6,k)
       cofa8=a2uib1*v(ib_tmp,k)+a2uib2*v(k4,k)+a2uib3*v(k5,k)+a2uib4*v(k6,k)
!       COFA3=A1U(IB,1)*U(IB,K)+A1U(IB,2)*U(K4,K)+A1U(IB,3)*U(K5,K)+A1U(IB,4)*U(K6,K)
!       COFA4=A2U(IB,1)*U(IB,K)+A2U(IB,2)*U(K4,K)+A2U(IB,3)*U(K5,K)+A2U(IB,4)*U(K6,K)
!       COFA7=A1U(IB,1)*V(IB,K)+A1U(IB,2)*V(K4,K)+A1U(IB,3)*V(K5,K)+A1U(IB,4)*V(K6,K)
!       COFA8=A2U(IB,1)*V(IB,K)+A2U(IB,2)*V(K4,K)+A2U(IB,3)*V(K5,K)+A2U(IB,4)*V(K6,K)
       uij2=u(ib_tmp,k)+cofa3*xijb+cofa4*yijb
       vij2=v(ib_tmp,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_tmp=vij*dltxc(i) - uij*dltyc(i)

       viscof1=art(ia)*sqrt(cofa1**2+cofa6**2+0.5_sp*(cofa2+cofa5)**2)
       viscof2=art(ib_tmp)*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_tmp)) + fm1*0.5_sp*(cc_hvc(ia)+cc_hvc(ib_tmp)))/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))

       !!UPWIND THE ADVECTIVE FLUX
       xadv=dij*un_tmp*((1.0_sp-sign(1.0_sp,un_tmp))*uij2+(1.0_sp+sign(1.0_sp,un_tmp))*uij1)*0.5_sp
       yadv=dij*un_tmp*((1.0_sp-sign(1.0_sp,un_tmp))*vij2+(1.0_sp+sign(1.0_sp,un_tmp))*vij1)*0.5_sp
        

       !!COMPUTE BOUNDARY FLUX AUGMENTERS   

       isbc_tmp = isbc(i)
       tpa = float(1-isbc_tmp)*epor(ia)
       tpb = float(1-isbc_tmp)*epor(ib_tmp)
       !!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

        xflux(ia,k)=xflux(ia,k)+xadv*tpa+(fxx+3.0_sp*fxx*float(isbc_tmp))*epor(ia)
        yflux(ia,k)=yflux(ia,k)+yadv*tpa+(fyy+3.0_sp*fyy*float(isbc_tmp))*epor(ia)
        xflux(ib,k)=xflux(ib,k)-xadv*tpb-(fxx+3.0_sp*fxx*float(isbc_tmp))*epor(ib)
        yflux(ib,k)=yflux(ib,k)-yadv*tpb-(fyy+3.0_sp*fyy*float(isbc_tmp))*epor(ib)

!        XFLUX(IA,K)=XFLUX(IA,K)+XADV*TPA+(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IA)
!        YFLUX(IA,K)=YFLUX(IA,K)+YADV*TPA+(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IA)
!        XFLUX(IB,K)=XFLUX(IB,K)-XADV*TPB-(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IB)
!        YFLUX(IB,K)=YFLUX(IB,K)-YADV*TPB-(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IB)

     END DO

    END IF
   END DO



   DO i=1,n
     iswettmp = iswetct(i)*iswetc(i)
     DO k=1,kbm1
       xflux(i,k) = xflux(i,k)*iswettmp
       yflux(i,k) = yflux(i,k)*iswettmp
     END DO
   END DO


!
!--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

!
!--Adjust Flux for Open Boundary Inflow/Outflow------------------------------------------|
!

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

   RETURN
   END SUBROUTINE advection_edge_gcn
!==============================================================================|