adv_uv_edge_gcn.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	adv_uv_edge_gcn

Function/Subroutine Documentation

adv_uv_edge_gcn()

subroutine adv_uv_edge_gcn

(

)

Definition at line 43 of file adv_uv_edge_gcn.f90.

!==============================================================================!
! 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 bcs
   USE mod_utils
   USE mod_spherical
   USE mod_northpole
   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
   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


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

   IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: adv_uv_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 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)
!#    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
       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

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

       ib_tmp = ib
!----------------------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)
!---------------------------------------------------------------------

!!       if(k==1)then
!!       if((egid(ia)==6794.and.egid(ib_tmp)==6793).or.(egid(ia)=&
!!            &=6793.and.egid(ib_tmp)==6794))then
!!        print*,'k1,k2,k3:',egid(k1),egid(k2),egid(k3)
!!        print*,'k4,k5,k6:',egid(k4),egid(k5),egid(k6)
!!        print*,'ia_a1u:',egid(ia),A1UIA1,A1UIA2,A1UIA3,A1UIA4
!!        print*,'ia_a2u:',egid(ia),A2UIA1,A2UIA2,A2UIA3,A2UIA4
!!        print*,'ib_a1u:',egid(ib_tmp),A1UIB1,A1UIB2,A1UIB3,A1UIB4
!!        print*,'ib_a2u:',egid(ib_tmp),A2UIB1,A2UIB2,A2UIB3,A2UIB4
!!       endif
!!       if((egid(ia)==3350.and.egid(ib_tmp)==3351).or.(egid(ia)=&
!!            &=3351.and.egid(ib_tmp)==3350))then
!!        print*,'k1,k2,k3:',egid(k1),egid(k2),egid(k3)
!!        print*,'k4,k5,k6:',egid(k4),egid(k5),egid(k6)
!!        print*,'ia_a1u:',egid(ia),A1UIA1,A1UIA2,A1UIA3,A1UIA4
!!        print*,'ia_a2u:',egid(ia),A2UIA1,A2UIA2,A2UIA3,A2UIA4
!!        print*,'ib_a1u:',egid(ib_tmp),A1UIB1,A1UIB2,A1UIB3,A1UIB4
!!        print*,'ib_a2u:',egid(ib_tmp),A2UIB1,A2UIB2,A2UIB3,A2UIB4
!!       endif
!!       endif

       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

       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

       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_tmp)*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_tmp)) + fm1*0.5_sp*(cc_hvc(ia)+cc_hvc(ib_tmp)))/hprnu
       viscofm(ia,k) = viscofm(ia,k) + viscof
       viscofm(ib_tmp,k) = viscofm(ib_tmp,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

       isbc_tmp = isbc(i)
       tpa = float(1-isbc_tmp)*epor(ia)
       tpb = float(1-isbc_tmp)*epor(ib_tmp)

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

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

   !ADJUST FLUX AT OPEN BOUNDARY MEAN FLOW

   
   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
         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_gcn"

Here is the call graph for this function:

Here is the caller graph for this function: