vertvl_edge.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 SIGMA COORDINATE VERTICAL VELOCITY FOR THE 3D MODE (omega)   |
!                                                  |
!   DETERMINED FROM EQUATION:                              |
!                                              !
!   d/dt(D) + d/dx(uD) + d/dy(uD) = d/sigma(omega)                             !
!==============================================================================|

   SUBROUTINE vertvl_edge         

!------------------------------------------------------------------------------|
   USE all_vars
   USE bcs
   USE mod_wd
   USE mod_northpole
 



   IMPLICIT NONE 
   REAL(SP) :: XFLUX(MT,KBM1),WBOTTOM(MT)
   REAL(SP) :: DIJ,UIJ,VIJ,UN,EXFLUX,TMP1,DIJ1,UIJ1,VIJ1
   INTEGER  :: I,K,IA,IB,I1,I2,I3,I4,J,JJ,J1,J2
!------------------------------------------------------------------------------|
   IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: Vertvl_edge.F"
!----------------------INITIALIZE FLUX-----------------------------------------!

   xflux = 0.0_sp

!----------------------ACCUMULATE FLUX-----------------------------------------!

!!#  if !defined (1)
   DO i=1,ncv
     i1=ntrg(i)
     ia=niec(i,1)
     ib=niec(i,2)

     DO k=1,kbm1
       dij=dt1(i1)*dz1(i1,k)
       uij=u(i1,k)
       vij=v(i1,k)
       exflux=dij*(-uij*dltye(i)+vij*dltxe(i))


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


     END DO
   END DO
!!#  else
!!   DO I=1,NCV
!!     I1=NTRG(I)
!!     IA=NIEC(I,1)
!!     IB=NIEC(I,2)

!!     DO K=1,KBM1
!!#      if !defined (SEMI_IMPLICIT)
!!       DIJ=DT1(I1)*DZ1(I1,K)
!!       UIJ=US(I1,K)
!!       VIJ=VS(I1,K)
!!       EXFLUX=DIJ*(-UIJ*DLTYE(I)+VIJ*DLTXE(I))
!!#      else
!!       DIJ=DT1(I1)*DZ1(I1,K)
!!       DIJ1=D1(I1)*DZ1(I1,K)
!!       UIJ=US(I1,K)
!!       VIJ=VS(I1,K)
!!       UIJ1=UF(I1,K)
!!       VIJ1=VF(I1,K)
!!       EXFLUX=( (1.0_SP-IFCETA)*DIJ*(-UIJ*DLTYE(I)+VIJ*DLTXE(I))+IFCETA*DIJ1*(-UIJ1*DLTYE(I)+VIJ1*DLTXE(I)) )*ISWETCT(I1)*ISWETC(I1)
!!#      endif
!!       XFLUX(IA,K)=XFLUX(IA,K)-EXFLUX
!!       XFLUX(IB,K)=XFLUX(IB,K)+EXFLUX
!!     END DO
!!   END DO
!!#  endif

   
!-----------------------NULLIFY BOUNDARY FLUX----------------------------------!
! For "tide + meanflow"/"meanflow only" case, this part should be commented out;
! For "tide only" case, this part may be kept.
! However, the effect of this term is small from my experience.

      DO i=1,m
        DO k=1,kbm1
          IF(isonb(i) == 2) xflux(i,k)=0.0_sp  
        ENDDO
      ENDDO
! can be changed to (no IF statements)
!     DO I=1,IOBCN
!        DO K=1,KBM1
!           XFLUX(I_OBC_N(I),K)=0.0_SP
!        ENDDO
!     ENDDO



!-----------------------FRESH WATER INFLOW-------------------------------------!

   IF(numqbc >= 1) THEN
     IF(river_inflow_location == 'node') THEN
       DO j=1,numqbc
         jj=inodeq(j)
         DO k=1,kbm1
           xflux(jj,k)=xflux(jj,k)-qdis(j)*vqdist(j,k)    !/DZ(JJ,K)
         END DO
       END DO
     ELSE IF(river_inflow_location == 'edge') THEN
       DO j=1,numqbc
         j1=n_icellq(j,1)
         j2=n_icellq(j,2)
         DO k=1,kbm1
           xflux(j1,k)=xflux(j1,k)-qdis(j)*rdisq(j,1)*vqdist(j,k)    !/DZ1(J1,K)
           xflux(j2,k)=xflux(j2,k)-qdis(j)*rdisq(j,2)*vqdist(j,k)    !/DZ1(J2,K)
         END DO
       END DO
     END IF
   END IF


!---IF NO FRESH WATER INFLOW, OMEGA IS ZERO AT FREE SURFACE AND BOTTOM---------!

   !CLEAR OLD VALUES
   wbottom = 0.0_sp
   wts = 0.0_sp


! QXU changed sign of evap/precip
   wts(:,1) = -(qevap+qprec)*rofvros 


   ! SET BOTTOM VELOCITY
   wbottom(1:m)=  bfwdis(1:m)/art1(1:m)

!--------------------------CALCULATE OMEGA-------------------------------------!

   DO i=1,m
    IF(iswetnt(i)*iswetn(i) == 1)THEN
     DO k=1,kbm1
!       WTS(I,K+1)=WTS(I,K)+DZ(I,K)*(XFLUX(I,K)/ART1(I)+(EL(I)-ET(I))/DTI)
       wts(i,k+1)=wts(i,k)+xflux(i,k)/art1(i)+dz(i,k)*(d(i)-dt(i))/dti
     END DO
    ELSE
     DO k=1,kbm1
       wts(i,k+1)=0.0_sp
     END DO
    END IF
   END DO


!-------------------------ADJUSTMENT FOR DAM MODULE----------------------------

!# if defined (THIN_DAM)
!   DO I=1,NODE_DAM1_N
!     I1 = I_NODE_DAM1_N(I,1)
!     I2 = I_NODE_DAM1_N(I,2)
!     if(i1==nlid(2944).or.i2==nlid(2944))print*,'orginal 2944:'&
!          &,kdam(nlid(2944)),wts(nlid(2944),3)
!     if(i1==nlid(4851).or.i2==nlid(4851))print*,'orginal 4851:'&
!          &,kdam(nlid(4851)),wts(nlid(4851),3)
!
!#  if defined (1)
!    IF(ISWETNT(I1)*ISWETN(I1) == 1)THEN
!#  endif
!     DO K=1,KDAM(I1)
!       WTS(I1,K+1)=WTS(I1,K)+(XFLUX(I1,K)+XFLUX(I2,K))/   &
!                   (ART1(I1)+ART1(I2))+DZ(I1,K)*(D(I1)-DT(I1))/DTI
!     END DO
!     DO K=KDAM(I1)+1,KBM1
!       WTS(I1,K+1)=WTS(I1,K)+XFLUX(I1,K)/ART1(I1)+DZ(I1,K)*(D(I1)-DT(I1))/DTI
!     END DO
!#  if defined (1)
!    ELSE
!     DO K=1,KBM1
!       WTS(I1,K+1)=0.0_SP
!     END DO
!    END IF
!#  endif
!
!#  if defined (1)
!    IF(ISWETNT(I2)*ISWETN(I2) == 1)THEN
!#  endif
!     DO K=1,KDAM(I2)
!       WTS(I2,K+1)=WTS(I2,K)+(XFLUX(I1,K)+XFLUX(I2,K))/   &
!                   (ART1(I1)+ART1(I2))+DZ(I2,K)*(D(I2)-DT(I2))/DTI
!     END DO
!     DO K=KDAM(I2)+1,KBM1
!       WTS(I2,K+1)=WTS(I2,K)+XFLUX(I2,K)/ART1(I2)+DZ(I2,K)*(D(I2)-DT(I2))/DTI
!     END DO
!#  if defined (1)
!    ELSE
!     DO K=1,KBM1
!       WTS(I2,K+1)=0.0_SP
!     END DO
!    END IF
!#  endif
!     if(i1==nlid(2944).or.i2==nlid(2944))print*,'adjusted 2944:'&
!          &,kdam(nlid(2944)),wts(nlid(2944),3)
!     if(i1==nlid(4851).or.i2==nlid(4851))print*,'adjusted 4851:'&
!          &,kdam(nlid(4851)),wts(nlid(4851),3)
!   END DO


!   DO I=1,NODE_DAM2_N
!     I1 = I_NODE_DAM2_N(I,1)
!     I2 = I_NODE_DAM2_N(I,2)
!     I3 = I_NODE_DAM2_N(I,3)
!#  if defined (1)
!    IF(ISWETNT(I1)*ISWETN(I1) == 1)THEN
!#  endif
!     DO K=1,KDAM(I1)
!       WTS(I1,K+1)=WTS(I1,K)+(XFLUX(I1,K)+XFLUX(I2,K)+XFLUX(I3,K))/   &
!                   (ART1(I1)+ART1(I2)+ART1(I3))+DZ(I1,K)*(D(I1)-DT(I1))/DTI
!     END DO
!     DO K=KDAM(I1)+1,KBM1
!       WTS(I1,K+1)=WTS(I1,K)+XFLUX(I1,K)/ART1(I1)+DZ(I1,K)*(D(I1)-DT(I1))/DTI
!     END DO
!#  if defined (1)
!    ELSE
!     DO K=1,KBM1
!       WTS(I1,K+1)=0.0_SP
!     END DO
!    END IF
!#  endif
!
!#  if defined (1)
!    IF(ISWETNT(I2)*ISWETN(I2) == 1)THEN
!#  endif
!     DO K=1,KDAM(I2)
!       WTS(I2,K+1)=WTS(I2,K)+(XFLUX(I1,K)+XFLUX(I2,K)+XFLUX(I3,K))/   &
!                   (ART1(I1)+ART1(I2)+ART1(I3))+DZ(I2,K)*(D(I2)-DT(I2))/DTI
!     END DO
!     DO K=KDAM(I2)+1,KBM1
!       WTS(I2,K+1)=WTS(I2,K)+XFLUX(I2,K)/ART1(I2)+DZ(I2,K)*(D(I2)-DT(I2))/DTI
!     END DO
!#  if defined (1)
!    ELSE
!     DO K=1,KBM1
!       WTS(I2,K+1)=0.0_SP
!     END DO
!    END IF
!#  endif
!
!#  if defined (1)
!    IF(ISWETNT(I3)*ISWETN(I3) == 1)THEN
!#  endif
!     DO K=1,KDAM(I3)
!       WTS(I3,K+1)=WTS(I3,K)+(XFLUX(I1,K)+XFLUX(I2,K)+XFLUX(I3,K))/   &
!                   (ART1(I1)+ART1(I2)+ART1(I3))+DZ(I3,K)*(D(I3)-DT(I3))/DTI
!     END DO
!     DO K=KDAM(I3)+1,KBM1
!       WTS(I3,K+1)=WTS(I3,K)+XFLUX(I3,K)/ART1(I3)+DZ(I3,K)*(D(I3)-DT(I3))/DTI
!     END DO
!#  if defined (1)
!    ELSE
!     DO K=1,KBM1
!       WTS(I3,K+1)=0.0_SP
!     END DO
!    END IF
!#  endif
!   END DO
!
!   DO I=1,NODE_DAM3_N
!     I1 = I_NODE_DAM3_N(I,1)
!     I2 = I_NODE_DAM3_N(I,2)
!     I3 = I_NODE_DAM3_N(I,3)
!     I4 = I_NODE_DAM3_N(I,4)
!#  if defined (1)
!    IF(ISWETNT(I1)*ISWETN(I1) == 1)THEN
!#  endif
!     DO K=1,KDAM(I1)
!       WTS(I1,K+1)=WTS(I1,K)+(XFLUX(I1,K)+XFLUX(I2,K)+XFLUX(I3,K)+XFLUX(I4,K))/   &
!                   (ART1(I1)+ART1(I2)+ART1(I3)+ART1(I4))+        &
!          DZ(I1,K)*(D(I1)-DT(I1))/DTI
!     END DO
!     DO K=KDAM(I1)+1,KBM1
!       WTS(I1,K+1)=WTS(I1,K)+XFLUX(I1,K)/ART1(I1)+DZ(I1,K)*(D(I1)-DT(I1))/DTI
!     END DO
!#  if defined (1)
!    ELSE
!     DO K=1,KBM1
!       WTS(I1,K+1)=0.0_SP
!     END DO
!    END IF
!#  endif
!
!#  if defined (1)
!    IF(ISWETNT(I2)*ISWETN(I2) == 1)THEN
!#  endif
!     DO K=1,KDAM(I2)
!       WTS(I2,K+1)=WTS(I2,K)+(XFLUX(I1,K)+XFLUX(I2,K)+XFLUX(I3,K)+XFLUX(I4,K))/   &
!                   (ART1(I1)+ART1(I2)+ART1(I3)+ART1(I4))+        &
!          DZ(I2,K)*(D(I2)-DT(I2))/DTI
!     END DO
!     DO K=KDAM(I2)+1,KBM1
!       WTS(I2,K+1)=WTS(I2,K)+XFLUX(I2,K)/ART1(I2)+DZ(I2,K)*(D(I2)-DT(I2))/DTI
!     END DO
!#  if defined (1)
!    ELSE
!     DO K=1,KBM1
!       WTS(I2,K+1)=0.0_SP
!     END DO
!    END IF
!#  endif
!
!#  if defined (1)
!    IF(ISWETNT(I3)*ISWETN(I3) == 1)THEN
!#  endif
!     DO K=1,KDAM(I3)
!       WTS(I3,K+1)=WTS(I3,K)+(XFLUX(I1,K)+XFLUX(I2,K)+XFLUX(I3,K)+XFLUX(I4,K))/   &
!                   (ART1(I1)+ART1(I2)+ART1(I3)+ART1(I4))+        &
!          DZ(I3,K)*(D(I3)-DT(I3))/DTI
!     END DO
!     DO K=KDAM(I3)+1,KBM1
!       WTS(I3,K+1)=WTS(I3,K)+XFLUX(I3,K)/ART1(I3)+DZ(I3,K)*(D(I3)-DT(I3))/DTI
!     END DO
!#  if defined (1)
!    ELSE
!     DO K=1,KBM1
!       WTS(I3,K+1)=0.0_SP
!     END DO
!    END IF
!#  endif
!
!#  if defined (1)
!    IF(ISWETNT(I4)*ISWETN(I4) == 1)THEN
!#  endif
!     DO K=1,KDAM(I4)
!       WTS(I4,K+1)=WTS(I4,K)+(XFLUX(I1,K)+XFLUX(I2,K)+XFLUX(I3,K)+XFLUX(I4,K))/   &
!                   (ART1(I1)+ART1(I2)+ART1(I3)+ART1(I4))+        &
!          DZ(I4,K)*(D(I4)-DT(I4))/DTI
!     END DO
!     DO K=KDAM(I4)+1,KBM1
!       WTS(I4,K+1)=WTS(I4,K)+XFLUX(I4,K)/ART1(I4)+DZ(I4,K)*(D(I4)-DT(I4))/DTI
!     END DO
!#  if defined (1)
!    ELSE
!     DO K=1,KBM1
!       WTS(I4,K+1)=0.0_SP
!     END DO
!    END IF
!#  endif
!   END DO
!
!#  endif

!--------------------------ADJUST OMEGA----------------------------------------!
! IMPROVES MASS CONSERVATION

   DO i=1,m
     IF(abs(wts(i,kb)-wbottom(i)) > 1.0e-7_sp)THEN
       IF(isonb(i) /= 2)THEN
         tmp1=elf(i)*float(kbm1)-(wts(i,kb)-wbottom(i))*dti/dz(i,1)
         tmp1=tmp1/float(kbm1)


         dtfa(i)=tmp1+h(i)

         DO k=2,kb
           wts(i,k)=wts(i,k)-float(k-1)/float(kbm1)*(wts(i,kb)-wbottom(i))
         END DO
       END IF
     END IF
   END DO

!
!----TRANSFER OMEGA TO FACE CENTER---------------------------------------------!
!

   DO i=1,n
     DO k=1,kb
       w(i,k) = one_third*(wts(nv(i,1),k)+wts(nv(i,2),k)+wts(nv(i,3),k))
     END DO
   END DO

   DO i=1,n
     DO k=1,kb
       w(i,k) = float(iswetc(i))*w(i,k)
     END DO
   END DO

   IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: Vertvl_edge.F"

   RETURN
   END SUBROUTINE vertvl_edge
!==============================================================================|