advave_edge_gcy.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	advave_edge_gcy (XFLUX, YFLUX)

Function/Subroutine Documentation

advave_edge_gcy()

subroutine advave_edge_gcy	(	real(sp), dimension(0:nt)	XFLUX,
		real(sp), dimension(0:nt)	YFLUX
	)

Definition at line 47 of file advave_edge_gcy.f90.

!==============================================================================|

   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
   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
   REAL(SP) :: VISCOF,VISCOF1,VISCOF2,TEMP
   REAL(SP) :: XFLUX(0:NT),YFLUX(0:NT)
   REAL(SP) :: FACT,FM1,ISWETTMP
   REAL(SP) :: UAK1,UAK2,UAK3,VAK1,VAK2,VAK3


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

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

     uak1 = ua(k1)
     uak2 = ua(k2)
     uak3 = ua(k3)
     vak1 = va(k1)
     vak2 = va(k2)
     vak3 = va(k3)

     IF(k1 == 0) CALL ghostuv2(ia,1,uak1,vak1)
     IF(k2 == 0) CALL ghostuv2(ia,2,uak2,vak2)
     IF(k3 == 0) CALL ghostuv2(ia,3,uak3,vak3)

     cofa1=a1u(ia,1)*ua(ia)+a1u(ia,2)*uak1+a1u(ia,3)*uak2+a1u(ia,4)*uak3
     cofa2=a2u(ia,1)*ua(ia)+a2u(ia,2)*uak1+a2u(ia,3)*uak2+a2u(ia,4)*uak3
     cofa5=a1u(ia,1)*va(ia)+a1u(ia,2)*vak1+a1u(ia,3)*vak2+a1u(ia,4)*vak3
     cofa6=a2u(ia,1)*va(ia)+a2u(ia,2)*vak1+a2u(ia,3)*vak2+a2u(ia,4)*vak3

     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)

     uak1 = ua(k1)
     uak2 = ua(k2)
     uak3 = ua(k3)
     vak1 = va(k1)
     vak2 = va(k2)
     vak3 = va(k3)

     IF(k1 == 0) CALL ghostuv2(ib,1,uak1,vak1)
     IF(k2 == 0) CALL ghostuv2(ib,2,uak2,vak2)
     IF(k3 == 0) CALL ghostuv2(ib,3,uak3,vak3)

     cofa3=a1u(ib,1)*ua(ib)+a1u(ib,2)*uak1+a1u(ib,3)*uak2+a1u(ib,4)*uak3
     cofa4=a2u(ib,1)*ua(ib)+a2u(ib,2)*uak1+a2u(ib,3)*uak2+a2u(ib,4)*uak3
     cofa7=a1u(ib,1)*va(ib)+a1u(ib,2)*vak1+a1u(ib,3)*vak2+a1u(ib,4)*vak3
     cofa8=a2u(ib,1)*va(ib)+a2u(ib,2)*vak1+a2u(ib,3)*vak2+a2u(ib,4)*vak3

     xij=xijc(i)-xc(ib)
     yij=yijc(i)-yc(ib)
     uij2=ua(ib)+cofa3*xij+cofa4*yij
     vij2=va(ib)+cofa7*xij+cofa8*yij

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

!lwu

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

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

!lwu

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

!    ACCUMULATE FLUX
!     XFLUX(IA)=XFLUX(IA)+XADV
!     YFLUX(IA)=YFLUX(IA)+YADV
!     XFLUX(IB)=XFLUX(IB)-XADV
!     YFLUX(IB)=YFLUX(IB)-YADV

     xflux(ia)=xflux(ia)+(xadv+fxx*epor(ia))*(1.0_sp-isbc(i))*iucp(ia)
     yflux(ia)=yflux(ia)+(yadv+fyy*epor(ia))*(1.0_sp-isbc(i))*iucp(ia)
     xflux(ib)=xflux(ib)-(xadv+fxx*epor(ib))*(1.0_sp-isbc(i))*iucp(ib)
     yflux(ib)=yflux(ib)-(yadv+fyy*epor(ib))*(1.0_sp-isbc(i))*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


   
   RETURN

Here is the call graph for this function: