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

!==============================================================================|
!   CALCULATE CONVECTION AND DIFFUSION FLUXES FOR EXTERNAL MODE                !
!==============================================================================|
   SUBROUTINE advave_edge_gcn(XFLUX,YFLUX)
!==============================================================================|

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





   IMPLICIT NONE
   INTEGER  :: I,J,K,IA,IB,J1,J2,K1,K2,K3,I1,I2,II
   REAL(SP) :: DIJ,ELIJ,XIJ,YIJ,UIJ,VIJ
   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
   REAL(SP) :: XFLUX(0:NT),YFLUX(0:NT)
   REAL(SP) :: FACT,FM1,ISWETTMP
   INTEGER  :: STAT_MF    
   REAL(SP) :: TPA,TPB


   REAL(SP) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY

!#  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) :: BTPS
   REAL(SP) :: U_TMP,V_TMP,UAC_TMP,VAC_TMP,WUSURF_TMP,WVSURF_TMP,WUBOT_TMP,WVBOT_TMP,UAF_TMP,VAF_TMP 

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: advave_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 FLUXES------------------------------------!
!
   xflux = 0.0_sp
   yflux = 0.0_sp
   pstx  = 0.0_sp
   psty  = 0.0_sp


!
!-------------------------ACCUMULATE FLUX OVER ELEMENT EDGES-------------------!
!
   
   DO i=1,ne
     ia=iec(i,1)
     ib=iec(i,2)

     j1=ienode(i,1)
     j2=ienode(i,2)

     dij=0.5_sp*(d(j1)+d(j2))
     elij=0.5_sp*(el(j1)+el(j2))




     IF(iswetce(ia)*iswetc(ia) == 1 .OR. iswetce(ib)*iswetc(ib) == 1)THEN
!    FLUX FROM LEFT
     k1=nbe(ia,1)
     k2=nbe(ia,2)
     k3=nbe(ia,3)

     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)
!---------------------------------------------------------------
     cofa1=a1uia1*ua(ia)+a1uia2*ua(k1)+a1uia3*ua(k2)+a1uia4*ua(k3)
     cofa2=a2uia1*ua(ia)+a2uia2*ua(k1)+a2uia3*ua(k2)+a2uia4*ua(k3)
     cofa5=a1uia1*va(ia)+a1uia2*va(k1)+a1uia3*va(k2)+a1uia4*va(k3)
     cofa6=a2uia1*va(ia)+a2uia2*va(k1)+a2uia3*va(k2)+a2uia4*va(k3)

!     COFA1=A1U(IA,1)*UA(IA)+A1U(IA,2)*UA(K1)+A1U(IA,3)*UA(K2)+A1U(IA,4)*UA(K3)
!     COFA2=A2U(IA,1)*UA(IA)+A2U(IA,2)*UA(K1)+A2U(IA,3)*UA(K2)+A2U(IA,4)*UA(K3)
!     COFA5=A1U(IA,1)*VA(IA)+A1U(IA,2)*VA(K1)+A1U(IA,3)*VA(K2)+A1U(IA,4)*VA(K3)
!     COFA6=A2U(IA,1)*VA(IA)+A2U(IA,2)*VA(K1)+A2U(IA,3)*VA(K2)+A2U(IA,4)*VA(K3)
     
     xij=xijc(i)-xc(ia)
     yij=yijc(i)-yc(ia)
     uij1=ua(ia)+cofa1*xij+cofa2*yij
     vij1=va(ia)+cofa5*xij+cofa6*yij

!    FLUX FROM RIGHT
     k1=nbe(ib,1)
     k2=nbe(ib,2)
     k3=nbe(ib,3)
     ib_tmp = ib

     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)

     cofa3=a1uib1*ua(ib_tmp)+a1uib2*ua(k1)+a1uib3*ua(k2)+a1uib4*ua(k3)
     cofa4=a2uib1*ua(ib_tmp)+a2uib2*ua(k1)+a2uib3*ua(k2)+a2uib4*ua(k3)
     cofa7=a1uib1*va(ib_tmp)+a1uib2*va(k1)+a1uib3*va(k2)+a1uib4*va(k3)
     cofa8=a2uib1*va(ib_tmp)+a2uib2*va(k1)+a2uib3*va(k2)+a2uib4*va(k3)

!     COFA3=A1U(IB,1)*UA(IB)+A1U(IB,2)*UA(K1)+A1U(IB,3)*UA(K2)+A1U(IB,4)*UA(K3)
!     COFA4=A2U(IB,1)*UA(IB)+A2U(IB,2)*UA(K1)+A2U(IB,3)*UA(K2)+A2U(IB,4)*UA(K3)
!     COFA7=A1U(IB,1)*VA(IB)+A1U(IB,2)*VA(K1)+A1U(IB,3)*VA(K2)+A1U(IB,4)*VA(K3)
!     COFA8=A2U(IB,1)*VA(IB)+A2U(IB,2)*VA(K1)+A2U(IB,3)*VA(K2)+A2U(IB,4)*VA(K3)
     
     xij=xijc(i)-xc(ib_tmp)
     yij=yijc(i)-yc(ib_tmp)
     uij2=ua(ib_tmp)+cofa3*xij+cofa4*yij
     vij2=va(ib_tmp)+cofa7*xij+cofa8*yij

!    NORMAL VELOCITY
     uij=0.5_sp*(uij1+uij2)
     vij=0.5_sp*(vij1+vij2)
     un_tmp=-uij*dltyc(i) + vij*dltxc(i)

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

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

!    ADD CONVECTIVE AND VISCOUS FLUXES
     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

!    ACCUMULATE FLUX
     isbc_tmp = isbc(i)
     xflux(ia)=xflux(ia)+(xadv+fxx*epor(ia))*(1.0_sp-isbc_tmp)*iucp(ia)
     yflux(ia)=yflux(ia)+(yadv+fyy*epor(ia))*(1.0_sp-isbc_tmp)*iucp(ia)
     xflux(ib)=xflux(ib)-(xadv+fxx*epor(ib))*(1.0_sp-isbc_tmp)*iucp(ib)
     yflux(ib)=yflux(ib)-(yadv+fyy*epor(ib))*(1.0_sp-isbc_tmp)*iucp(ib)

     END IF



!    ACCUMULATE BAROTROPIC FLUX
!for spherical coordinator and domain across 360^o latitude         
     pstx(ia)=pstx(ia)-grav_e(ia)*d1(ia)*elij*dltyc(i)
     psty(ia)=psty(ia)+grav_e(ia)*d1(ia)*elij*dltxc(i)
     pstx(ib)=pstx(ib)+grav_e(ib)*d1(ib)*elij*dltyc(i)
     psty(ib)=psty(ib)-grav_e(ib)*d1(ib)*elij*dltxc(i)

   END DO

!#  if !defined (SEMI_IMPLICIT)

   DO i = 1,n
     iswettmp = iswetce(i)*iswetc(i)
     xflux(i) = xflux(i)*iswettmp
     yflux(i) = yflux(i)*iswettmp
   END DO


!
!-------------------------SET BOUNDARY VALUES----------------------------------!
!

!  MODIFY BOUNDARY FLUX
      DO i=1,n
        IF(isbce(i) == 2) THEN
          xflux(i)=(xflux(i)+fluxobn(i)*ua(i))*iucp(i)
          yflux(i)=(yflux(i)+fluxobn(i)*va(i))*iucp(i)
        ENDIF
      END DO


!  ADJUST FLUX FOR RIVER INFLOW
   IF(numqbc > 0) THEN
     IF(river_inflow_location == 'node')THEN
       DO k=1,numqbc
         j=inodeq(k)
         i1=nbve(j,1)
         i2=nbve(j,ntve(j))
         vlctyq(k)=qdis(k)/qarea(k)
         xflux(i1)=xflux(i1)-0.5_sp*qdis(k)*vlctyq(k)*cos(angleq(k))
         yflux(i1)=yflux(i1)-0.5_sp*qdis(k)*vlctyq(k)*sin(angleq(k))
         xflux(i2)=xflux(i2)-0.5_sp*qdis(k)*vlctyq(k)*cos(angleq(k))
         yflux(i2)=yflux(i2)-0.5_sp*qdis(k)*vlctyq(k)*sin(angleq(k))
       END DO
     ELSE IF(river_inflow_location == 'edge') THEN
       DO k=1,numqbc
         i1=icellq(k)
         vlctyq(k)=qdis(k)/qarea(k)
         temp=qdis(k)*vlctyq(k)
         xflux(i1)=xflux(i1)-temp*cos(angleq(k))
         yflux(i1)=yflux(i1)-temp*sin(angleq(k))
       END DO
     END IF
   END IF

!  ADJUST FLUX FOR OPEN BOUNDARY MEAN FLOW



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

   END SUBROUTINE advave_edge_gcn
!==============================================================================|