adv_q.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 the Turbulent Kinetic Energy and Mixing Length Based  on         |
!   The Mellor-Yamada Level 2.5 Turbulent Closure Model                        |
!==============================================================================|

   SUBROUTINE adv_q(Q,QF)               

!------------------------------------------------------------------------------|
   USE mod_utils
   USE all_vars
   USE mod_par
   USE mod_wd
   USE mod_spherical
   USE mod_northpole


   IMPLICIT NONE
   REAL(SP), DIMENSION(0:MT,KB)     :: Q,QF,XFLUX
   REAL(SP), DIMENSION(0:MT)        :: PUPX,PUPY,PVPX,PVPY  
   REAL(SP), DIMENSION(0:MT)        :: PQPX,PQPY,PQPXD,PQPYD,VISCOFF
   REAL(SP), DIMENSION(3*(NT),KBM1) :: DTIJ 
   REAL(SP), DIMENSION(3*(NT),KBM1) :: UVN
   REAL(SP) :: UTMP,VTMP,SITAI,FFD,FF1 !,X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2,XI,YI
   REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2,UN
   REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF,EXFLUX,TEMP,STPOINT
   REAL(SP) :: FACT,FM1
   INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,K,JTMP,JJ,II
   REAL(SP) :: Q1MIN, Q1MAX, Q2MIN, Q2MAX

!!$#  if defined (SPHERICAL)
!!$   REAL(DP) :: TY,TXPI,TYPI
!!$   REAL(DP) :: XTMP1,XTMP
!!$   REAL(DP) :: X1_DP,Y1_DP,X2_DP,Y2_DP,XII,YII
!!$   REAL(DP) :: X11_TMP,Y11_TMP,X33_TMP,Y33_TMP
!!$   REAL(DP) :: VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP
!!$   REAL(DP) :: TXPI_TMP,TYPI_TMP
!!$#  endif

   REAL(SP) :: QMEAN1
   REAL(SP), DIMENSION(0:NT,KB)    :: UQ,VQ

   REAL(SP), ALLOCATABLE :: UQ1(:,:),VQ1(:,:)


   IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: adv_q"

!------------------------------------------------------------------------------!

   qmean1 = 1.e-8

   ALLOCATE(uq1(0,0))
   ALLOCATE(vq1(0,0))   

   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-----------------------------------------------------------!
!
   qf    = 0.0_sp
   xflux = 0.0_sp
   
   uq = 0.0_sp
   vq = 0.0_sp
   uvn = 0.0_sp
   
   
   DO k=2,kbm1
     DO i=1,nt
       uq(i,k) = (u(i,k)*dz1(i,k-1)+u(i,k-1)*dz1(i,k))/(dz1(i,k)+dz1(i,k-1))
       vq(i,k) = (v(i,k)*dz1(i,k-1)+v(i,k-1)*dz1(i,k))/(dz1(i,k)+dz1(i,k-1))
     END DO
   END DO     

!
!--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
!
   DO i=1,ncv
     i1=ntrg(i)
!     DTIJ(I)=DT1(I1)
     DO k=2,kbm1
       dtij(i,k)=dt1(i1)*dzz1(i1,k-1)
       uvn(i,k) = vq(i1,k)*dltxe(i) - uq(i1,k)*dltye(i)


     END DO
   END DO

!
!--Calculate the Advection and Horizontal Diffusion Terms----------------------!
!

   DO k=2,kbm1
     pqpx  = 0.0_sp 
     pqpy  = 0.0_sp 
     pqpxd = 0.0_sp 
     pqpyd = 0.0_sp
     DO i=1,m
       DO j=1,ntsn(i)-1
         i1=nbsn(i,j)
         i2=nbsn(i,j+1)

!         FFD=0.5_SP*(Q(I1,K)+Q(I2,K)-QMEAN1(I1,K)-QMEAN1(I2,K))
         ffd=0.5_sp*(q(i1,k)+q(i2,k)-qmean1-qmean1)
         ff1=0.5_sp*(q(i1,k)+q(i2,k))


!!$#        if defined (SPHERICAL)
!!$         XTMP  = VX(I2)*TPI-VX(I1)*TPI
!!$  XTMP1 = VX(I2)-VX(I1)
!!$  IF(XTMP1 >  180.0_SP)THEN
!!$    XTMP = -360.0_SP*TPI+XTMP
!!$  ELSE IF(XTMP1 < -180.0_SP)THEN
!!$    XTMP =  360.0_SP*TPI+XTMP
!!$  END IF  
!!$         TXPI=XTMP*COS(DEG2RAD*VY(I))
!!$         TYPI=(VY(I1)-VY(I2))*TPI
!!$
!!$         ! ERROR HERE
!!$         IF(NODE_NORTHAREA(I) == 1)THEN
!!$           VX1_TMP = REARTH * COS(VY(I1)*DEG2RAD) * COS(VX(I1)*DEG2RAD) &
!!$                     * 2._SP /(1._SP+sin(VY(I1)*DEG2RAD))
!!$           VY1_TMP = REARTH * COS(VY(I1)*DEG2RAD) * SIN(VX(I1)*DEG2RAD) &
!!$                     * 2._SP /(1._SP+sin(VY(I1)*DEG2RAD))
!!$
!!$           VX2_TMP = REARTH * COS(VY(I2)*DEG2RAD) * COS(VX(I2)*DEG2RAD) &
!!$                     * 2._SP /(1._SP+sin(VY(I2)*DEG2RAD))
!!$           VY2_TMP = REARTH * COS(VY(I2)*DEG2RAD) * SIN(VX(I2)*DEG2RAD) &
!!$                     * 2._SP /(1._SP+sin(VY(I2)*DEG2RAD))
!!$
!!$           TXPI = (VX2_TMP-VX1_TMP)/(2._SP /(1._SP+sin(VY(I)*DEG2RAD)))
!!$           TYPI = (VY1_TMP-VY2_TMP)/(2._SP /(1._SP+sin(VY(I)*DEG2RAD)))
!!$        IF(I /= NODE_NORTHPOLE)THEN
!!$      TXPI_TMP = TYPI*COS(VX(I)*DEG2RAD)-TXPI*SIN(VX(I)*DEG2RAD)
!!$      TYPI_TMP = TXPI*COS(VX(I)*DEG2RAD)+TYPI*SIN(VX(I)*DEG2RAD)
!!$      TYPI_TMP = -TYPI_TMP
!!$     
!!$      TXPI = TXPI_TMP
!!$      TYPI = TYPI_TMP
!!$    END IF  
!!$  END IF 
!!$         ! END ERROR
!!$
!!$
!!$         PQPX(I)=PQPX(I)+FF1*TYPI
!!$         PQPY(I)=PQPY(I)+FF1*TXPI
!!$         PQPXD(I)=PQPXD(I)+FFD*TYPI
!!$         PQPYD(I)=PQPYD(I)+FFD*TXPI
!!$#        else
!!$         PQPX(I)=PQPX(I)+FF1*(VY(I1)-VY(I2))
!!$         PQPY(I)=PQPY(I)+FF1*(VX(I2)-VX(I1))
!!$         PQPXD(I)=PQPXD(I)+FFD*(VY(I1)-VY(I2))
!!$         PQPYD(I)=PQPYD(I)+FFD*(VX(I2)-VX(I1))
!!$#        endif

         pqpx(i)=pqpx(i)+ff1*dltytrie(i,j)
         pqpy(i)=pqpy(i)+ff1*dltxtrie(i,j)
         pqpxd(i)=pqpxd(i)+ffd*dltytrie(i,j)
         pqpyd(i)=pqpyd(i)+ffd*dltxtrie(i,j)


       END DO
       pqpx(i)=pqpx(i)/art2(i)
       pqpy(i)=pqpy(i)/art2(i)
       pqpxd(i)=pqpxd(i)/art2(i)
       pqpyd(i)=pqpyd(i)/art2(i)
     END DO


          
     DO i=1,m 
       viscoff(i) = (viscofh(i,k)*dz(i,k-1)+viscofh(i,k-1)*dz(i,k))/  &
                    (dz(i,k)+dz(i,k-1))  
     END DO

     DO i=1,ncv_i
       ia=niec(i,1)
       ib=niec(i,2)
!!$       XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
!!$       YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
!!$#      if defined (SPHERICAL)
!!$       X1_DP=XIJE(I,1)
!!$       Y1_DP=YIJE(I,1)
!!$       X2_DP=XIJE(I,2)
!!$       Y2_DP=YIJE(I,2)
!!$       CALL ARCC(X2_DP,Y2_DP,X1_DP,Y1_DP,XII,YII)
!!$       XI=XII        
!!$       XTMP  = XI*TPI-VX(IA)*TPI
!!$       XTMP1 = XI-VX(IA)
!!$       IF(XTMP1 >  180.0_SP)THEN
!!$         XTMP = -360.0_SP*TPI+XTMP
!!$       ELSE IF(XTMP1 < -180.0_SP)THEN
!!$         XTMP =  360.0_SP*TPI+XTMP
!!$       END IF     
!!$
!!$       DXA=XTMP*COS(DEG2RAD*VY(IA))    
!!$       DYA=(YI-VY(IA))*TPI
!!$       XTMP  = XI*TPI-VX(IB)*TPI
!!$       XTMP1 = XI-VX(IB)
!!$       IF(XTMP1 >  180.0_SP)THEN
!!$         XTMP = -360.0_SP*TPI+XTMP
!!$       ELSE IF(XTMP1 < -180.0_SP)THEN
!!$         XTMP =  360.0_SP*TPI+XTMP
!!$       END IF     
!!$
!!$       DXB=XTMP*COS(DEG2RAD*VY(IB)) 
!!$       DYB=(YI-VY(IB))*TPI
!!$#      else
!!$       DXA=XI-VX(IA)
!!$       DYA=YI-VY(IA)
!!$       DXB=XI-VX(IB)
!!$       DYB=YI-VY(IB)
!!$#      endif
!!$
!!$       FIJ1=Q(IA,K)+DXA*PQPX(IA)+DYA*PQPY(IA)
!!$       FIJ2=Q(IB,K)+DXB*PQPX(IB)+DYB*PQPY(IB)

        fij1=q(ia,k)+dltxncve(i,1)*pqpx(ia)+dltyncve(i,1)*pqpy(ia)
        fij2=q(ib,k)+dltxncve(i,2)*pqpx(ib)+dltyncve(i,2)*pqpy(ib)


       q1min=minval(q(nbsn(ia,1:ntsn(ia)-1),k))
       q1min=min(q1min, q(ia,k))
       q1max=maxval(q(nbsn(ia,1:ntsn(ia)-1),k))
       q1max=max(q1max, q(ia,k))
       q2min=minval(q(nbsn(ib,1:ntsn(ib)-1),k))
       q2min=min(q2min, q(ib,k))
       q2max=maxval(q(nbsn(ib,1:ntsn(ib)-1),k))
       q2max=max(q2max, q(ib,k))
       IF(fij1 < q1min) fij1=q1min
       IF(fij1 > q1max) fij1=q1max
       IF(fij2 < q2min) fij2=q2min
       IF(fij2 > q2max) fij2=q2max
    
       un=uvn(i,k)

!       VISCOF=HORCON*(FACT*0.5_SP*(VISCOFF(IA)+VISCOFF(IB))/HPRNU + FM1)
!       VISCOF=HORCON*(FACT*0.5_SP*(VISCOFF(IA)+VISCOFF(IB)) + FM1)

       ! David moved HPRNU and added HVC
       viscof=(fact*0.5_sp*(viscoff(ia)*nn_hvc(ia)+viscoff(ib)*nn_hvc(ib)) + fm1*0.5_sp*(nn_hvc(ia)+nn_hvc(ib)))/hprnu

       txx=0.5_sp*(pqpxd(ia)+pqpxd(ib))*viscof
       tyy=0.5_sp*(pqpyd(ia)+pqpyd(ib))*viscof

       fxx=-dtij(i,k)*txx*dltye(i)
       fyy= dtij(i,k)*tyy*dltxe(i)


       exflux=-un*dtij(i,k)* &
          ((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy

       xflux(ia,k)=xflux(ia,k)+exflux
       xflux(ib,k)=xflux(ib,k)-exflux

     END DO


   END DO !!SIGMA LOOP

!
!-Accumulate Fluxes at Boundary Nodes
!

 
!--------------------------------------------------------------------
!   The central difference scheme in vertical advection
!--------------------------------------------------------------------
   DO k=2,kbm1
     DO i=1,m
       IF(iswetn(i)*iswetnt(i) == 1) THEN
         temp=wts(i,k-1)*q(i,k-1)-wts(i,k+1)*q(i,k+1)
         xflux(i,k)=xflux(i,k)+temp*art1(i)*dzz(i,k-1)/(dz(i,k-1)+dz(i,k))
       END IF
     END DO
   END DO  !! SIGMA LOOP

!
!--Update Q or QL-------------------------------------------------------------!
!

   DO i=1,m
     IF(iswetn(i)*iswetnt(i) == 1 )THEN
        DO k=2,kbm1
           qf(i,k)=(q(i,k)-xflux(i,k)/art1(i)*(dti/(dt(i)*dzz(i,k-1))))*(dt(i)/d(i))
        END DO
     ELSE
        DO k=2,kbm1
           qf(i,k)=q(i,k)
        END DO
     END IF
   END DO

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

   END SUBROUTINE adv_q
!==============================================================================|