advection_edge_gcy.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_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

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

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

   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)
     xija=xijc(i)-xc(ia)
     yija=yijc(i)-yc(ia)
     xijb=xijc(i)-xc(ib)
     yijb=yijc(i)-yc(ib)

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


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

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

   RETURN
   END SUBROUTINE advection_edge_gcy
!==============================================================================|