extel_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 FLUXES OF FREE SURFACE ELEVATION (CONTINUITY) EQUATION            |
!==============================================================================|
   SUBROUTINE extel_edge(K)       
!==============================================================================|
   USE all_vars
   USE bcs
   USE mod_obcs



!  ggao 0903/2007



   IMPLICIT NONE
   INTEGER, INTENT(IN) :: K
   REAL(SP) :: XFLUX(0:MT)
   REAL(SP) :: DIJ,UIJ,VIJ,DTK,UN,EXFLUX
   INTEGER  :: I,J,I1,IA,IB,JJ,J1,J2


!==============================================================================|
   IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: extel_edge.F",k
!----------INITIALIZE FLUX ARRAY ----------------------------------------------!

   xflux = 0.0_sp

!---------ACCUMULATE FLUX BY LOOPING OVER CONTROL VOLUME HALF EDGES------------!

   DO i=1,ncv
     i1  = ntrg(i)
     ia  = niec(i,1)
     ib  = niec(i,2)
     dij = d1(i1)

     uij = ua(i1)
     vij = va(i1)
     exflux = dij*(-uij*dltye(i) + vij*dltxe(i))  


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


   END DO

!   write(ipt,*) "after control volume flux",maxval(xflux),minval(xflux)
!   write(ipt,*) "after control volume flux",maxloc(xflux),minloc(xflux)


!--ADD EVAPORATION AND PRECIPITATION TERMS-------------------------------------!
   IF (precipitation_on) THEN

!qxu   XFLUX = XFLUX+(QEVAP2-QPREC2)*ROFVROS*ART1 
!qxu---the evap is negative for evaporating in ocean
      xflux = xflux-(qevap2+qprec2)*rofvros*art1 
   END IF
    
!--ADD GROUND WATER TERM-------------------------------------------------------!

   IF(groundwater_on) THEN
      xflux = xflux - bfwdis2
   END IF

!--SAVE ACCUMULATED FLUX ON OPEN BOUNDARY NODES AND ZERO OUT OPEN BOUNDARY FLUX!

   IF(iobcn > 0) THEN  
     DO i=1,iobcn
       xflux_obcn(i)=xflux(i_obc_n(i))
       xflux(i_obc_n(i)) = 0.0_sp
     END DO
   END IF


!---------ADJUST FLUX FOR FRESH WATER DISCHARGE--------------------------------!

   IF(numqbc >= 1) THEN   
     IF(river_inflow_location == 'node') THEN
       DO j=1,numqbc
         jj=inodeq(j)
         xflux(jj)=xflux(jj)-qdis(j)
       END DO
     ELSE IF(river_inflow_location == 'edge') THEN
       DO j=1,numqbc
         j1=n_icellq(j,1)
         j2=n_icellq(j,2)
         xflux(j1)=xflux(j1)-qdis(j)*rdisq(j,1)
         xflux(j2)=xflux(j2)-qdis(j)*rdisq(j,2)
       END DO
     END IF
   END IF


!----------PERFORM UPDATE ON ELF-----------------------------------------------!

   dtk = alpha_rk(k)*dte
   elf = elrk - dtk*xflux/art1
!!#  if defined (THIN_DAM)
!!   DO I=1,NODE_DAM1_N
!!     JN=I_NODE_DAM1_N(I,1)
!!     ELF(JN)=ELRK(JN)-DTK*(XFLUX(JN)+XFLUX(I_NODE_DAM1_N(I,2)))&
!!          &/(ART1(JN)+ART1(I_NODE_DAM1_N(I,2)))
!!     ELF(I_NODE_DAM1_N(I,2))=ELF(JN)
!!   END DO
!!   DO I=1,NODE_DAM2_N
!!     JN=I_NODE_DAM2_N(I,1)
!!     ELF(JN)=ELRK(JN)-DTK*(XFLUX(JN)+XFLUX(I_NODE_DAM2_N(I,2))&
!!          &+XFLUX(I_NODE_DAM2_N(I,3)) )&
!!          &/(ART1(JN)+ART1(I_NODE_DAM2_N(I,2))+ART1(I_NODE_DAM2_N(I,3)))
!!     ELF(I_NODE_DAM2_N(I,2))=ELF(JN)
!!     ELF(I_NODE_DAM2_N(I,3))=ELF(JN)
!!   END DO
!!   DO I=1,NODE_DAM3_N
!!     JN=I_NODE_DAM3_N(I,1)
!!     ELF(JN)=ELRK(JN)-DTK*(XFLUX(JN)+XFLUX(I_NODE_DAM3_N(I,2))&
!!          &+XFLUX(I_NODE_DAM3_N(I,3))+XFLUX(I_NODE_DAM3_N(I,4)) )&
!!          &/(ART1(JN)+ART1(I_NODE_DAM3_N(I,2))+ART1(I_NODE_DAM3_N(I&
!!          &,3))+ART1(I_NODE_DAM3_N(I,4)))
!!     ELF(I_NODE_DAM3_N(I,2))=ELF(JN)
!!     ELF(I_NODE_DAM3_N(I,3))=ELF(JN)
!!     ELF(I_NODE_DAM3_N(I,4))=ELF(JN)
!!   END DO
!!#  endif

!
!--STORE VARIABLES FOR MOMENTUM BALANCE CHECK----------------------------------|
!
  
   IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "END: extel_edge.F"

   RETURN
   END SUBROUTINE extel_edge
!==============================================================================|