adv_uv_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 adv_uv_edge_gcy

!==============================================================================!
! this subroutine calculate advective, coriolis, pressure gradient, etc in     !
! x and y momentum equations except vertical diffusion terms for internal mode ! 
!==============================================================================!

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





   IMPLICIT NONE
   REAL(SP) :: XFLUX(0:NT,KB),YFLUX(0:NT,KB)
   REAL(SP) :: PSTX_TM(0:NT,KB),PSTY_TM(0:NT,KB)
   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) :: DIJ,ELIJ,TMPA,TMPB,TMP,XFLUXV,YFLUXV
   REAL(SP) :: FACT,FM1,EXFLUX,ISWETTMP
   INTEGER  :: I,IA,IB,J1,J2,K1,K2,K3,K4,K5,K6,K,II,J,I1,I2


   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: adv_uv_edge_gcy"

!------------------------------------------------------------------------------!
   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 Flux Variables------------------------------------------------!
!
   viscofm = 0.0_sp
   
   xflux  = 0.0_sp
   yflux  = 0.0_sp
   pstx_tm = 0.0_sp
   psty_tm = 0.0_sp

!
!-----Loop Over Edges and Accumulate Flux--------------------------------------!
!

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


     elij=0.5_sp*(egf(j1)+egf(j2))




     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))
       IF(iswetct(ia)*iswetc(ia) == 1 .OR. iswetct(ib)*iswetc(ib) == 1)THEN

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

       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)
!       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,k)+cofa3*xijb+cofa4*yijb
       vij2=v(ib,k)+cofa7*xijb+cofa8*yijb

       uij=0.5_sp*(uij1+uij2)
       vij=0.5_sp*(vij1+vij2)
       exflux = dij*(-uij*dltyc(i) + vij*dltxc(i))


!
!-------ADD THE VISCOUS TERM & ADVECTION TERM---------------------------------!
!

       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=FACT*0.5_SP*HORCON*(VISCOF1+VISCOF2)/HPRNU + FM1*HORCON
!       VISCOF=FACT*0.5_SP*HORCON*(VISCOF1+VISCOF2) + FM1*HORCON

       ! 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
       viscofm(ia,k) = viscofm(ia,k) + viscof
       viscofm(ib,k) = viscofm(ib,k) + viscof

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


       xadv=exflux*((1.0_sp-sign(1.0_sp,exflux))*uij2+(1.0_sp+sign(1.0_sp,exflux))*uij1)*0.5_sp
       yadv=exflux*((1.0_sp-sign(1.0_sp,exflux))*vij2+(1.0_sp+sign(1.0_sp,exflux))*vij1)*0.5_sp

       !!CALCULATE BOUNDARY FLUX AUGMENTERS
       tpa = float(1-isbc(i))*epor(ia)
       tpb = float(1-isbc(i))*epor(ib)


       !!ACCUMULATE ADVECTIVE + DIFFUSIVE + BAROTROPIC PRESSURE GRADIENT TERMS
       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 IF
! for spherical coordinator and domain across 360^o latitude
!        PSTX_TM(IA,K)=PSTX_TM(IA,K)-GRAV*DT1(IA)*ELIJ*DLTYC(I)
!        PSTY_TM(IA,K)=PSTY_TM(IA,K)+GRAV*DT1(IA)*ELIJ*DLTXC(I)
!        PSTX_TM(IB,K)=PSTX_TM(IB,K)+GRAV*DT1(IB)*ELIJ*DLTYC(I)
!        PSTY_TM(IB,K)=PSTY_TM(IB,K)-GRAV*DT1(IB)*ELIJ*DLTXC(I)
        pstx_tm(ia,k)=pstx_tm(ia,k)-grav_e(ia)*dt1(ia)*dz1(ia,k)*elij*dltyc(i)
        psty_tm(ia,k)=psty_tm(ia,k)+grav_e(ia)*dt1(ia)*dz1(ia,k)*elij*dltxc(i)
        pstx_tm(ib,k)=pstx_tm(ib,k)+grav_e(ib)*dt1(ib)*dz1(ib,k)*elij*dltyc(i)
        psty_tm(ib,k)=psty_tm(ib,k)-grav_e(ib)*dt1(ib)*dz1(ib,k)*elij*dltxc(i)

     END DO
   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
         pstx_tm(i,k)= pstx_tm(i,k)*iswettmp
         psty_tm(i,k)= psty_tm(i,k)*iswettmp
        END DO
       DO k=1,kbm1
        xflux(i,k)=xflux(i,k)+pstx_tm(i,k)
        yflux(i,k)=yflux(i,k)+psty_tm(i,k)
       END DO
      END DO

!
!-------ADD VERTICAL CONVECTIVE FLUX, CORIOLIS TERM AND BAROCLINIC PG TERM----!
!
   DO i=1,n
     IF(iswetct(i)*iswetc(i) == 1)THEN


     DO k=1,kbm1

       IF(k == 1) THEN
         xfluxv=-w(i,k+1)*(u(i,k)*dz1(i,k+1)+u(i,k+1)*dz1(i,k))/&
                 (dz1(i,k)+dz1(i,k+1))
         yfluxv=-w(i,k+1)*(v(i,k)*dz1(i,k+1)+v(i,k+1)*dz1(i,k))/&
                 (dz1(i,k)+dz1(i,k+1))
       ELSE IF(k == kbm1) THEN
         xfluxv= w(i,k)*(u(i,k)*dz1(i,k-1)+u(i,k-1)*dz1(i,k))/&
                 (dz1(i,k)+dz1(i,k-1))
         yfluxv= w(i,k)*(v(i,k)*dz1(i,k-1)+v(i,k-1)*dz1(i,k))/&
                 (dz1(i,k)+dz1(i,k-1))
      ELSE
         xfluxv= w(i,k)*(u(i,k)*dz1(i,k-1)+u(i,k-1)*dz1(i,k))/&
                 (dz1(i,k)+dz1(i,k-1))-&
                 w(i,k+1)*(u(i,k)*dz1(i,k+1)+u(i,k+1)*dz1(i,k))/&
                 (dz1(i,k)+dz1(i,k+1))
         yfluxv= w(i,k)*(v(i,k)*dz1(i,k-1)+v(i,k-1)*dz1(i,k))/&
                 (dz1(i,k)+dz1(i,k-1))-&
                 w(i,k+1)*(v(i,k)*dz1(i,k+1)+v(i,k+1)*dz1(i,k))/&
                 (dz1(i,k)+dz1(i,k+1))
       END IF
!       XFLUX(I,K)=XFLUX(I,K)+XFLUXV/DZ(K)*ART(I)&
!                 +DRHOX(I,K)-COR(I)*V(I,K)*DT1(I)*ART(I)
!       YFLUX(I,K)=YFLUX(I,K)+YFLUXV/DZ(K)*ART(I)&
!                 +DRHOY(I,K)+COR(I)*U(I,K)*DT1(I)*ART(I)
       xflux(i,k)=xflux(i,k)+xfluxv*art(i)&
                 +drhox(i,k)-cor(i)*v(i,k)*dt1(i)*dz1(i,k)*art(i)
       yflux(i,k)=yflux(i,k)+yfluxv*art(i)&
                 +drhoy(i,k)+cor(i)*u(i,k)*dt1(i)*dz1(i,k)*art(i)


     END DO
    END IF
   END DO

   
   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 AT RIVER INFLOWS
   IF(numqbc >= 1) 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
     ELSE
       print*,'RIVER_INFLOW_LOCATION NOT CORRECT'
       CALL pstop
     END IF
   END IF


   DO i=1,n
      DO k=1,kbm1
         !       UF(I,K)=U(I,K)*DT1(I,K)/D1(I,K)-DTI*XFLUX(I,K)/ART(I)/D1(I,K)
         !       VF(I,K)=V(I,K)*DT1(I,K)/D1(I,K)-DTI*YFLUX(I,K)/ART(I)/D1(I,K)
         uf(i,k)=u(i,k)*dt1(i)/d1(i)-dti*xflux(i,k)/art(i)/(d1(i)*dz1(i,k))
         vf(i,k)=v(i,k)*dt1(i)/d1(i)-dti*yflux(i,k)/art(i)/(d1(i)*dz1(i,k))
         IF(adcor_on) THEN
         ! FOR ADCOR
           ubeta(i,k)=xflux(i,k) +cor(i)*v(i,k)*dt1(i)*dz1(i,k)*art(i)*epor(i)
           vbeta(i,k)=yflux(i,k) -cor(i)*u(i,k)*dt1(i)*dz1(i,k)*art(i)*epor(i)
         ENDIF
      END DO
   END DO

   
   DO i =1,n
     IF(iswetct(i)*iswetc(i) .NE. 1)THEN
       DO k=1,kbm1
         uf(i,k)=0.0_sp
         vf(i,k)=0.0_sp
         ubeta(i,k)=0.0_sp
         vbeta(i,k)=0.0_sp
       END DO
     END IF
   END DO
   
   DO k=1,kb
    viscofm(:,k) = viscofm(:,k)/art(:)
   END DO 

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

   END SUBROUTINE adv_uv_edge_gcy
!==============================================================================!