mod_obcs3.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$
!/===========================================================================/

MODULE mod_obcs3

   USE all_vars
   USE mod_prec
   USE mod_obcs

   USE mod_meanflow
   USE mod_obcs2

   IMPLICIT NONE
   SAVE

!--Nonlinear Velocity Open Boundary Condition Arrays
   REAL(sp), ALLOCATABLE :: fluxobn2(:),cxobc(:),cyobc(:)
   REAL(sp), ALLOCATABLE :: fluxobc2d_x(:),     fluxobc2d_y(:)
   REAL(sp), ALLOCATABLE :: obc2d_x_tide(:),    obc2d_y_tide(:)
   REAL(sp), ALLOCATABLE :: fluxobc3d_x(:,:),   fluxobc3d_y(:,:)
   REAL(sp), ALLOCATABLE :: fluxobc3d_x_2(:,:), fluxobc3d_y_2(:,:)
   CONTAINS


!=========================================================================|
   SUBROUTINE alloc_obc3_data

   IMPLICIT NONE

   ALLOCATE(fluxobn2(0:iobcn))                 ;fluxobn2      = zero
   ALLOCATE(cxobc(0:nt))                       ;cxobc         = zero
   ALLOCATE(cyobc(0:nt))                       ;cyobc         = zero
   ALLOCATE(fluxobc2d_x(0:nmfcell_i))          ;fluxobc2d_x   = zero
   ALLOCATE(fluxobc2d_y(0:nmfcell_i))          ;fluxobc2d_y   = zero
   ALLOCATE(fluxobc3d_x(0:nmfcell_i,1:kbm1))   ;fluxobc3d_x   = zero
   ALLOCATE(fluxobc3d_y(0:nmfcell_i,1:kbm1))   ;fluxobc3d_y   = zero
   ALLOCATE(fluxobc3d_x_2(0:nmfcell_i,1:kbm1)) ;fluxobc3d_x_2 = zero
   ALLOCATE(fluxobc3d_y_2(0:nmfcell_i,1:kbm1)) ;fluxobc3d_y_2 = zero

   ALLOCATE(obc2d_x_tide(0:nmfcell),obc2d_y_tide(0:nmfcell)) 

   RETURN
   END SUBROUTINE alloc_obc3_data
!==========================================================================|

!==========================================================================|
   SUBROUTINE zero_obc3

   IMPLICIT NONE

   integer :: I

   IF (nmfcell > 0) THEN
      DO i = 0, nmfcell
         obc2d_x_tide(i) = 0.0_sp
         obc2d_y_tide(i) = 0.0_sp
      END DO
   END IF

   RETURN
   END SUBROUTINE zero_obc3
!==========================================================================|


!==============================================================================|
   SUBROUTINE setup_obc3

   IMPLICIT NONE

   INTEGER  :: I, I1, I2, J, IC, N1, N2, N3
   REAL(SP) :: DXN,DYN,DXC,DYC,CROSS

 IF (iobcn > 0) THEN
   DO i=1,iobcn
     i1 = i_obc_n(i)
     i2 = adjn_obc(i,1)
     DO j = 1, ntve(i1)
       ic = nbve(i1,j)
       n1 = nv(ic,1) ; n2 = nv(ic,2) ; n3 = nv(ic,3)
       IF( n1-i2 == 0 .OR. n2-i2 == 0 .OR. n3-i2 == 0)THEN
         dxn = vx(i2)-vx(i1) ; dyn = vy(i2)-vy(i1)
         dxc = xc(ic)-vx(i1) ; dyc = yc(ic)-vy(i1)
         cross     =  sign(1.0_sp,dxc*dyn - dyc*dxn)
         cxobc(ic) =  cross*dyn/sqrt(dxn**2 +dyn**2)
         cyobc(ic) = -cross*dxn/sqrt(dxn**2 +dyn**2)
       END IF
     END DO

     IF(nadjn_obc(i) > 1)THEN
       i2 = adjn_obc(i,2)
       DO j = 1, ntve(i1)
         ic = nbve(i1,j)
         n1 = nv(ic,1) ; n2 = nv(ic,2) ; n3 = nv(ic,3)
         IF( n1-i2 == 0 .OR. n2-i2 == 0 .OR. n3-i2 == 0)THEN
         dxn = vx(i2)-vx(i1) ; dyn = vy(i2)-vy(i1)
         dxc = xc(ic)-vx(i1) ; dyc = yc(ic)-vy(i1)
           cross     = sign(1.0_sp,dxc*dyn - dyc*dxn)
           cxobc(ic) = cross*dyn/sqrt(dxn**2 +dyn**2)
           cyobc(ic) =-cross*dxn/sqrt(dxn**2 +dyn**2)
         END IF
       END DO
     END IF
   END DO
 END IF

   RETURN
   END SUBROUTINE setup_obc3

!==============================================================================|
   SUBROUTINE flux_obn2d(KKT)

   IMPLICIT NONE

   INTEGER, INTENT(IN)  :: KKT
   INTEGER  :: I,I1,N1,N2,N3,C1,C2
   REAL(SP) :: E1_X,E1_Y,E2_X,E2_Y
   REAL(SP) :: FLUXF_OBC_1,FLUXF_OBC_2


!   IF (IOBCN > 0) THEN
!   DO I = 1, IOBCN
!     IF(NADJN_OBC(I) == 1)THEN
!       N1 = I_OBC_N(I)
!       N2 = ADJN_OBC(I,1)
!       I1 = I_OBC_NODE(N1)
!
!       E1_Y = VY(N1)-VY(N2)
!#      if defined (SPHERICAL)
!       X1_DP = VX(N2)
!       Y1_DP = VY(N2)
!       X2_DP = VX(N1)
!       Y2_DP = VY(N1)
!       CALL ARCX(X1_DP,Y1_DP,X2_DP,Y2_DP,SIDE)
!       E1_X = SIDE
!       E1_Y = TPI*E1_Y
!#      else
!       E1_X = VX(N1)-VX(N2)
!#      endif
!       
!       C1 = I_OBC_CELL2(ADJC_OBC(I,1)) 
!       FLUXF_OBC_1 = 0.5_SP*SQRT(E1_X**2+E1_Y**2)*  &
!                   (UANP(C1)*CXOBC(ADJC_OBC(I,1))+VANP(C1)*CYOBC(ADJC_OBC(I,1)))
!
!       FLUXOBN2(I) = -FLUXF_OBC_1*(H(N1)+ELT(I1)+ELP(I1)) 
!
!     ELSE
!       N1   = I_OBC_N(I)
!       N2   = ADJN_OBC(I,1)
!       N3   = ADJN_OBC(I,2)
!       I1   = I_OBC_NODE(N1)
!
!       E1_Y = VY(N1)-VY(N2)
!#      if defined (SPHERICAL)
!       X1_DP = VX(N2)
!       Y1_DP = VY(N2)
!       X2_DP = VX(N1)
!       Y2_DP = VY(N1)
!       CALL ARCX(X1_DP,Y1_DP,X2_DP,Y2_DP,SIDE)
!       E1_X = SIDE
!       E1_Y = TPI*E1_Y
!#      else
!       E1_X = VX(N1)-VX(N2)
!#      endif
!
!       E2_Y = VY(N1)-VY(N3)
!#      if defined (SPHERICAL)
!       X1_DP = VX(N3)
!       Y1_DP = VY(N3)
!       X2_DP = VX(N1)
!       Y2_DP = VY(N1)
!       CALL ARCX(X1_DP,Y1_DP,X2_DP,Y2_DP,SIDE)
!       E2_X = SIDE
!       E2_Y = TPI*E2_Y
!#      else
!       E2_X = VX(N1)-VX(N3)
!#      endif
!
!       C1 = I_OBC_CELL2(ADJC_OBC(I,1)) 
!       C2 = I_OBC_CELL2(ADJC_OBC(I,2)) 
!       FLUXF_OBC_1 = 0.5_SP*SQRT(E1_X**2+E1_Y**2)*  &
!                   (UANP(C1)*CXOBC(ADJC_OBC(I,1))+VANP(C1)*CYOBC(ADJC_OBC(I,1)))
!       FLUXF_OBC_2 = 0.5_SP*SQRT(E2_X**2+E2_Y**2)*  &
!                   (UANP(C2)*CXOBC(ADJC_OBC(I,2))+VANP(C2)*CYOBC(ADJC_OBC(I,2)))
!
!       FLUXOBN2(I) = -(FLUXF_OBC_1+FLUXF_OBC_2)*(H(N1)+ELT(I1)+ELP(I1)) 
!     END IF
!   END DO
!   END IF   

   IF (nmfcell > 0) THEN
     DO i = 1, nmfcell
        c1= i_obc_node2(node_mfcell(i,1))
        c2= i_obc_node2(node_mfcell(i,2))
        fluxobn2(c1) = fluxobn2(c1) - mfqdis(i)*rdismf(i,1)
        fluxobn2(c2) = fluxobn2(c2) - mfqdis(i)*rdismf(i,2)
     END DO
   END IF

   RETURN
   END SUBROUTINE flux_obn2d


!==============================================================================|
   SUBROUTINE flux_obc2d

   IMPLICIT NONE

   INTEGER  :: I,II,I1,I2,J,J1,J2,ITMP,JTMP
   REAL(DP)  DX12,DY12,TMP1,DTMP
   REAL(SP) :: FLUXF_OBC_1,FLUXF_OBC_2

   IF (nmfcell_i > 0) THEN
     DO i = 1, nmfcell_i
       ii= i_mfcell_n(i)
       itmp=0
       DO j=1,3
         IF(nbe(ii,j) == 0) THEN
           jtmp=j
           itmp=itmp+1
         END IF
       END DO
       IF (itmp /= 1) THEN
         print*,'something is wrong here 2'
         CALL pstop
       END IF
    
       j1=jtmp+1-int((jtmp+1)/4)*3
       j2=jtmp+2-int((jtmp+2)/4)*3
       i1=nv(ii,j1)
       i2=nv(ii,j2)
       dy12=vy(i1)-vy(i2)
       dx12=vx(i1)-vx(i2)

       dtmp = 0.5_sp*(h(i1)+h(i2)+elt(i_obc_node(i1))+elt(i_obc_node(i2)))    ! May be a problem, should be replaced by D
!       TMP1 = -(UANT(I)*cos(ANGLEMF(I))+VANT(I)*sin(ANGLEMF(I))*(SQRT(DX12**2+DY12**2))
       tmp1 = uant(i) * dy12 - vant(i) * dx12
       fluxobc2d_x(i) = dtmp * tmp1 * uant(i)
       fluxobc2d_y(i) = dtmp * tmp1 * vant(i)

     END DO
   END IF

   END SUBROUTINE flux_obc2d


!==============================================================================|
   SUBROUTINE flux_obc3d

   IMPLICIT NONE

   INTEGER  :: I,II,I1,I2,J,J1,J2,K,ITMP,JTMP
   REAL(DP)  DX12,DY12,TMP1,DTMP
   REAL(SP) :: UTMP,VTMP

   IF (nmfcell > 0) THEN

! 2-D and 3-D adjustment
      obc2d_x_tide = obc2d_x_tide/float(isplit)
      obc2d_y_tide = obc2d_y_tide/float(isplit)
                                                                                                                         
      DO i = 1, nmfcell
         ii= i_mfcell_n(i)
         itmp=0
         DO j=1,3
           IF(nbe(ii,j) == 0 .and. isonb(nv(ii,j)) /= 2) THEN
             jtmp=j
             itmp=itmp+1
           END IF
         END DO
         IF (itmp /= 1) THEN
           print*,'something is wrong here 3'
           CALL pstop
         END IF

         j1=jtmp+1-int((jtmp+1)/4)*3
         j2=jtmp+2-int((jtmp+2)/4)*3
         i1=nv(ii,j1)
         i2=nv(ii,j2)
         dtmp = 0.5_sp*(h(i1)+h(i2)+eltdt(i_obc_node(i1))+eltdt(i_obc_node(i2)))      ! May be a problem, should be replaced by D

         utmp = 0.0_sp ; vtmp = 0.0_sp
         DO k=1,kbm1
            utmp = utmp + unt(i,k)*dz1(ii,k)
            vtmp = vtmp + vnt(i,k)*dz1(ii,k)
         END DO
         utmp = utmp * dtmp
         vtmp = vtmp * dtmp
         DO k=1,kbm1
           unt(i,k) = unt(i,k) - (utmp-obc2d_x_tide(i))/dtmp
           vnt(i,k) = vnt(i,k) - (vtmp-obc2d_y_tide(i))/dtmp
         END DO
      END DO
   END IF

   IF (nmfcell_i > 0) THEN
     DO i = 1, nmfcell_i
       ii= i_mfcell_n(i)
       itmp=0
       DO j=1,3
         IF(nbe(ii,j) == 0) THEN
           jtmp=j
           itmp=itmp+1
         END IF
       END DO
       j1=jtmp+1-int((jtmp+1)/4)*3
       j2=jtmp+2-int((jtmp+2)/4)*3
       i1=nv(ii,j1)
       i2=nv(ii,j2)
       dy12=vy(i1)-vy(i2)
       dx12=vx(i1)-vx(i2)

       dtmp = 0.5_sp*(h(i1)+h(i2)+eltdt(i_obc_node(i1))+eltdt(i_obc_node(i2)))       ! May be a problem, should be replaced by D
       DO k =1, kbm1
!          TMP1 = -(UNT(I,K)*cos(ANGLEMF(I))+VNT(I,K)*sin(ANGLEMF(I))*(SQRT(DX12**2+DY12**2))
          tmp1 = unt(i,k) * dy12 - vnt(i,k) * dx12
          fluxobc3d_x(i,k) = dtmp * tmp1 * unt(i,k)
          fluxobc3d_y(i,k) = dtmp * tmp1 * vnt(i,k)
       END DO
     END DO

   END IF

   END SUBROUTINE flux_obc3d


!==============================================================================|
   SUBROUTINE flux_obc3d_2

   USE all_vars
   USE mod_obcs
   USE mod_obcs2

   IMPLICIT NONE

   INTEGER  :: I,II,I1,I2,J,J1,J2,K,ITMP,JTMP
   REAL(DP)  DX12,DY12,TMP1,DTMP
   REAL(SP) :: UTMP,VTMP

   IF (nmfcell > 0) THEN

! 2-D and 3-D adjustment
      DO i = 1, nmfcell
       ii= i_mfcell_n(i)
         utmp   = sum(unt(i,1:kbm1)*dz1(ii,1:kbm1))
         vtmp   = sum(vnt(i,1:kbm1)*dz1(ii,1:kbm1))
         unt(i,1:kbm1) = unt(i,1:kbm1) + (uant(i) - utmp) 
         vnt(i,1:kbm1) = vnt(i,1:kbm1) + (vant(i) - vtmp) 
      END DO
   END IF

   IF (nmfcell_i > 0) THEN
     DO i = 1, nmfcell_i
       ii= i_mfcell_n(i)
       itmp=0
       DO j=1,3
         IF(nbe(ii,j) == 0) THEN
           jtmp=j
           itmp=itmp+1
         END IF
       END DO
       j1=jtmp+1-int((jtmp+1)/4)*3
       j2=jtmp+2-int((jtmp+2)/4)*3
       i1=nv(ii,j1)
       i2=nv(ii,j2)
       dy12=vy(i1)-vy(i2)
       dx12=vx(i1)-vx(i2)

       dtmp = 0.5_sp*(h(i1)+h(i2)+eltdt(i_obc_node(i1))+eltdt(i_obc_node(i2)))          ! May be a problem, should be replaced by D
       DO k =1, kbm1
!          TMP1 = -(UNT(I,K)*cos(ANGLEMF(I))+VNT(I,K)*sin(ANGLEMF(I))*(SQRT(DX12**2+DY12**2))
          tmp1 = unt(i,k) * dy12 - vnt(i,k) * dx12
          fluxobc3d_x_2(i,k) = dtmp * tmp1 * unt(i,k)
          fluxobc3d_y_2(i,k) = dtmp * tmp1 * vnt(i,k)
       END DO
     END DO

   END IF

   END SUBROUTINE flux_obc3d_2

!========================================================================
END MODULE mod_obcs3