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

   SUBROUTINE shape_coef_gcy

!----------------------------------------------------------------------!
!  This subroutine is used to calculate the coefficient for a linear   !
!  function on the x-y plane, i.e.:                                    !
!                                                                      !
!                     r(x,y;phai)=phai_c+cofa1*x+cofa2*y               !
!                                                                      !
!  This subroutine is used for ghost cell boundary condition cases     !
!----------------------------------------------------------------------!

   USE all_vars
   USE mod_utils
   IMPLICIT NONE
   REAL(DP) X1,X2,X3,Y1,Y2,Y3,DELT,AI1,AI2,AI3,BI1,BI2,BI3,CI1,CI2,CI3
   REAL(DP) DELTX,DELTY,TEMP1,ANG1,ANG2,B1,B2,ANGLE
   INTEGER  I,II,J,JJ,J1,J2
   REAL(SP) AA1,AA2,BB1,BB2,CC1,CC2,XTMP1,YTMP1,XTMP2,YTMP2
!
!---------------interior and boundary cells------------------------------------!
!

   IF(dbg_set(dbg_log)) THEN
      WRITE(ipt,*) "!"
      WRITE(ipt,*) "!           SETTING UP LINEAR INTEROPOLATION COEFFICIENTS"
   END IF

   DO i = 1, n
     IF(isbce(i) == 0) THEN
       y1 = yc(nbe(i,1))-yc(i)
       y2 = yc(nbe(i,2))-yc(i)
       y3 = yc(nbe(i,3))-yc(i)
       x1 = xc(nbe(i,1))-xc(i)
       x2 = xc(nbe(i,2))-xc(i)
       x3 = xc(nbe(i,3))-xc(i)
     ELSE IF(isbce(i) == 1) THEN
       DO j = 1, 3
         IF(nbe(i,j) == 0) jj = j
       END DO
       j1 = jj+1-int((jj+1)/4)*3
       j2 = jj+2-int((jj+2)/4)*3



       IF(jj == 1)THEN
         x1 = xtmp1-xc(i)
         y1 = ytmp1-yc(i)
         x2 = xc(nbe(i,j1))-xc(i)
         y2 = yc(nbe(i,j1))-yc(i)
         x3 = xc(nbe(i,j2))-xc(i)
         y3 = yc(nbe(i,j2))-yc(i)
       ELSE IF(jj == 2)THEN
         x1 = xc(nbe(i,j2))-xc(i)
         y1 = yc(nbe(i,j2))-yc(i)
         x2 = xtmp1-xc(i)
         y2 = ytmp1-yc(i)
         x3 = xc(nbe(i,j1))-xc(i)
         y3 = yc(nbe(i,j1))-yc(i)
       ELSE
         x1 = xc(nbe(i,j1))-xc(i)
         y1 = yc(nbe(i,j1))-yc(i)
         x2 = xc(nbe(i,j2))-xc(i)
         y2 = yc(nbe(i,j2))-yc(i)
         x3 = xtmp1-xc(i)
         y3 = ytmp1-yc(i)
       END IF
!     else if(isbce(i) == 2) then
!       do j=1,4
!         a1u(i,j)=0.0_SP
!         a2u(i,j)=0.0_SP
!       end do
     ELSE IF(isbce(i) == 2) THEN
       DO j = 1, 3
         IF(nbe(i,j) == 0) jj = j
       END DO
       j1 = jj+1-int((jj+1)/4)*3
       j2 = jj+2-int((jj+2)/4)*3

       deltx = vx(nv(i,j1))-vx(nv(i,j2))
       delty = vy(nv(i,j1))-vy(nv(i,j2))

       alpha(i) = atan2(delty,deltx)
       alpha(i) = alpha(i)-3.1415926_sp/2.0_sp

       aa1 = -delty
       bb1 = deltx
       cc1 = -aa1*vx(nv(i,j1))-bb1*vy(nv(i,j1))

       aa2 = bb1
       bb2 = -aa1
       cc2 = -aa2*xc(i)-bb2*yc(i)

       xtmp1 = -(cc1*bb2-cc2*bb1)/(aa1*bb2-aa2*bb1)
       ytmp1 = -(cc1*aa2-cc2*aa1)/(bb1*aa2-bb2*aa1)

       IF(jj == 1)THEN
         x1 = (xtmp1-xc(i))*2.0_sp
         y1 = (ytmp1-yc(i))*2.0_sp
         x2 = xc(nbe(i,j1))-xc(i)
         y2 = yc(nbe(i,j1))-yc(i)
         x3 = xc(nbe(i,j2))-xc(i)
         y3 = yc(nbe(i,j2))-yc(i)
       ELSE IF(jj == 2)THEN
         x1 = xc(nbe(i,j2))-xc(i)
         y1 = yc(nbe(i,j2))-yc(i)
         x2 = (xtmp1-xc(i))*2.0_sp
         y2 = (ytmp1-yc(i))*2.0_sp
         x3 = xc(nbe(i,j1))-xc(i)
         y3 = yc(nbe(i,j1))-yc(i)
       ELSE
         x1 = xc(nbe(i,j1))-xc(i)
         y1 = yc(nbe(i,j1))-yc(i)
         x2 = xc(nbe(i,j2))-xc(i)
         y2 = yc(nbe(i,j2))-yc(i)
         x3 = (xtmp1-xc(i))*2.0_sp
         y3 = (ytmp1-yc(i))*2.0_sp
       END IF
     ELSE IF(isbce(i) == 3)THEN
       DO j = 1, 3
         IF(nbe(i,j) /= 0) jj = j
       END DO
       j1 = jj+1-int((jj+1)/4)*3
       j2 = jj+2-int((jj+2)/4)*3



       IF(jj == 1)THEN
         x1 = xc(nbe(i,jj))-xc(i)
         y1 = yc(nbe(i,jj))-yc(i)
         x2 = xtmp1-xc(i)
         y2 = ytmp1-yc(i)
         x3 = xtmp2-xc(i)
         y3 = ytmp2-yc(i)
       ELSE IF(jj == 2)THEN
         x1 = xtmp2-xc(i)
         y1 = ytmp2-yc(i)
         x2 = xc(nbe(i,jj))-xc(i)
         y2 = yc(nbe(i,jj))-yc(i)
         x3 = xtmp1-xc(i)
         y3 = ytmp1-yc(i)
       ELSE
         x1 = xtmp1-xc(i)
         y1 = ytmp1-yc(i)
         x2 = xtmp2-xc(i)
         y2 = ytmp2-yc(i)
         x3 = xc(nbe(i,jj))-xc(i)
         y3 = yc(nbe(i,jj))-yc(i)
       END IF
     END IF

     x1 = x1/1000.0_sp
     x2 = x2/1000.0_sp
     x3 = x3/1000.0_sp
     y1 = y1/1000.0_sp
     y2 = y2/1000.0_sp
     y3 = y3/1000.0_sp

     delt = (x1*y2-x2*y1)**2+(x1*y3-x3*y1)**2+(x2*y3-x3*y2)**2
     delt = delt*1000._sp

     a1u(i,1) = (y1+y2+y3)*(x1*y1+x2*y2+x3*y3)- &
                (x1+x2+x3)*(y1**2+y2**2+y3**2)
     a1u(i,1) = a1u(i,1)/delt
     a1u(i,2) = (y1**2+y2**2+y3**2)*x1-(x1*y1+x2*y2+x3*y3)*y1
     a1u(i,2) = a1u(i,2)/delt
     a1u(i,3) = (y1**2+y2**2+y3**2)*x2-(x1*y1+x2*y2+x3*y3)*y2
     a1u(i,3) = a1u(i,3)/delt
     a1u(i,4) = (y1**2+y2**2+y3**2)*x3-(x1*y1+x2*y2+x3*y3)*y3
     a1u(i,4) = a1u(i,4)/delt

     a2u(i,1) = (x1+x2+x3)*(x1*x1+x2*x2+x3*x3)- &
                (y1+y2+y3)*(x1**2+x2**2+x3**2)
     a2u(i,1) = a2u(i,1)/delt
     a2u(i,2) = (x1**2+x2**2+x3**2)*y1-(x1*y1+x2*y2+x3*y3)*x1
     a2u(i,2) = a2u(i,2)/delt
     a2u(i,3) = (x1**2+x2**2+x3**2)*y2-(x1*y1+x2*y2+x3*y3)*x2
     a2u(i,3) = a2u(i,3)/delt
     a2u(i,4) = (x1**2+x2**2+x3**2)*y3-(x1*y1+x2*y2+x3*y3)*x3
     a2u(i,4) = a2u(i,4)/delt
!     end if

     x1 = vx(nv(i,1))-xc(i)
     x2 = vx(nv(i,2))-xc(i)
     x3 = vx(nv(i,3))-xc(i)
     y1 = vy(nv(i,1))-yc(i)
     y2 = vy(nv(i,2))-yc(i)
     y3 = vy(nv(i,3))-yc(i)


     ai1 = y2-y3
     ai2 = y3-y1
     ai3 = y1-y2
     bi1 = x3-x2
     bi2 = x1-x3
     bi3 = x2-x1
     ci1 = x2*y3-x3*y2
     ci2 = x3*y1-x1*y3
     ci3 = x1*y2-x2*y1

     aw0(i,1) = -ci1/2./art(i)
     aw0(i,2) = -ci2/2./art(i)
     aw0(i,3) = -ci3/2./art(i)
     awx(i,1) = -ai1/2./art(i)
     awx(i,2) = -ai2/2./art(i)
     awx(i,3) = -ai3/2./art(i)
     awy(i,1) = -bi1/2./art(i)
     awy(i,2) = -bi2/2./art(i)
     awy(i,3) = -bi3/2./art(i)
   END DO


   ang1=359.9_sp/180.0_sp*3.1415926_sp
   ang2=-0.1_sp/180.0_sp*3.1415926_sp

   do i=1,m
     if((isonb(i).eq.1).and.(ntve(i).gt.2)) then
       angle=alpha(nbve(i,ntve(i)))-alpha(nbve(i,1))
       if(angle.gt.ang1) then
         angle=100000.0_sp
       else if(angle.gt.3.1415926_sp) then
         angle=angle-2.0_sp*3.1415926_sp
       else if(angle.lt.-3.1415926_sp) then
         angle=angle+2.0_sp*3.1415926_sp
       else if(angle.lt.ang2) then
         angle=100000.0_sp
       end if
       do j=2,ntve(i)-1
         ii=nbve(i,j)
         if(isbce(ii).ne.1) then
           alpha(ii)=alpha(nbve(i,1))+ &
                     angle/float(ntve(i)-1)*float(j-1)
         end if
       end do
     end if
   end do
!end>

   IF(dbg_set(dbg_log)) THEN
       WRITE(ipt,*) "!"
       WRITE(ipt,*) "!  INTERP COEFFICIENTS   :    COMPLETE"
    END IF

   RETURN
   END SUBROUTINE shape_coef_gcy