couple2adcirc.ftn90

      MODULE Couple2Adcirc

      INTEGER   IUNIT
      INTEGER   IOSTAT, IT0, SAVITE, ILEN
      INTEGER   INERR
      INTEGER   IERR
      INTEGER   ISTAT, IF1, IL1
      CHARACTER PTYPE, PNAME *8, COMPUT *4
      CHARACTER*20 CHARS(1)
      CHARACTER*20 MSGSTR
      LOGICAL   LOPEN

      INTEGER, ALLOCATABLE :: CROSS(:)
      INTEGER, ALLOCATABLE :: BGRIDP(:)
      REAL   , ALLOCATABLE :: BSPECS(:,:,:,:)
      REAL   , ALLOCATABLE :: AC1(:,:,:), COMPDA(:,:)

      REAL, ALLOCATABLE    :: BLKND(:), BLKNDC(:)
      REAL*8, ALLOCATABLE  :: OURQT(:)

      CONTAINS



      SUBROUTINE MakeBoundariesReflective(ivert, ac2 )

      USE M_GENARR, ONLY        : spcdir
      USE SwanGridData, ONLY    : nverts
      USE SwanGridobjects, ONLY : CELLID,               &
                                  celltype,             &
                                  CELLV1,CELLV2,CELLV3, &
                                  FACEID,               &
                                  facetype,             &
                                  FACEV1,FACEV2,        &
                                  FMARKER,              &
                                  gridobject,           &
                                  verttype,             &
                                  VERTX,VERTY
      USE SWCOMM3, ONLY         : DDIR,                 &
                                  MDC,MSC
      IMPLICIT NONE

      INTRINSIC              :: ATAN
      INTRINSIC              :: COS
      INTRINSIC              :: INT
      INTRINSIC              :: MOD
      INTRINSIC              :: REAL
      INTRINSIC              :: SIN
      INTRINSIC              :: SQRT

      INTEGER                :: icell
      INTEGER                :: id
      INTEGER                :: iface
      INTEGER                :: IncD1
      INTEGER                :: IncD2
      INTEGER                :: is
      INTEGER,INTENT(IN)     :: ivert
      INTEGER                :: jc
      INTEGER                :: jf
      INTEGER                :: NumBdySegs
      INTEGER                :: NumFaces
      INTEGER                :: V1
      INTEGER                :: V2
      INTEGER                :: V3

      LOGICAL                :: IntoDomain

      REAL,INTENT(INOUT)     :: ac2(MDC,MSC,nverts)
      REAL                   :: ac2temp(MDC,MSC)
      REAL                   :: AvgDX
      REAL                   :: BdyDir(2)
      REAL                   :: BdyDirTemp
      REAL                   :: BdyDirToUse
      REAL                   :: Dir
      REAL                   :: Dist(2)
      REAL                   :: IncDir
      REAL                   :: IncCounter
      REAL                   :: Pi = 3.141592654
      REAL                   :: SubArea1
      REAL                   :: SubArea2
      REAL                   :: SubArea3
      REAL                   :: TotalArea
      REAL                   :: W1
      REAL                   :: W2
      REAL                   :: X1
      REAL                   :: X2
      REAL                   :: X3
      REAL                   :: XV
      REAL                   :: Y1
      REAL                   :: Y2
      REAL                   :: Y3
      REAL                   :: YV

      TYPE(celltype),POINTER :: cell(:)
      TYPE(facetype),POINTER :: face(:)
      TYPE(verttype),POINTER :: vert(:)

!... Initialize ac2temp.
      DO id=1,MDC
         DO is=1,MSC
            ac2temp(id,is) = 0.
         ENDDO
      ENDDO

!... Point to data.
      cell => gridobject%cell_grid
      face => gridobject%face_grid
      vert => gridobject%vert_grid

!... Find a length scale by taking the average dx between this
!... vertex and its neighbors.
      AvgDX = 0.
      NumFaces = 0
      DO jc=1,vert(ivert)%noc
         icell = vert(ivert)%cell(jc)%atti(CELLID)
         DO jf=1,cell(icell)%nof
            V1 = cell(icell)%face(jf)%atti(FACEV1)
            V2 = cell(icell)%face(jf)%atti(FACEV2)
            IF( V1==ivert .OR. V2==ivert )THEN
               X1 = vert(V1)%attr(VERTX)
               X2 = vert(V2)%attr(VERTX)
               Y1 = vert(V1)%attr(VERTY)
               Y2 = vert(V2)%attr(VERTY)
               AvgDX = AvgDX + SQRT((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1))
               NumFaces = NumFaces + 1
            ENDIF
         ENDDO
      ENDDO
      AvgDX = AvgDX / REAL(NumFaces)

!... Loop over the directions.
      DO id=1,MDC
         Dir = spcdir(id,1)

!... Determine if this direction points out of the domain.

!... Set up a fictional vertex along the direction.
         XV = vert(ivert)%attr(VERTX) + COS(Dir) * ( 0.1 * AvgDX )
         YV = vert(ivert)%attr(VERTY) + SIN(Dir) * ( 0.1 * AvgDX )

!... Loop over the adjacent cells.  Use triangle areas
!... to determine if the fictional vertex is located
!... inside an adjacent cell.
         IntoDomain = .FALSE.
         DO jc=1,vert(ivert)%noc

            icell = vert(ivert)%cell(jc)%atti(CELLID)
            V1 = cell(icell)%atti(CELLV1)
            V2 = cell(icell)%atti(CELLV2)
            V3 = cell(icell)%atti(CELLV3)

            X1 = XV
            X2 = vert(V2)%attr(VERTX)
            X3 = vert(V3)%attr(VERTX)
            Y1 = YV
            Y2 = vert(V2)%attr(VERTY)
            Y3 = vert(V3)%attr(VERTY)
            SubArea1 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

            X1 = vert(V1)%attr(VERTX)
            X2 = XV
            X3 = vert(V3)%attr(VERTX)
            Y1 = vert(V1)%attr(VERTY)
            Y2 = YV
            Y3 = vert(V3)%attr(VERTY)
            SubArea2 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

            X1 = vert(V1)%attr(VERTX)
            X2 = vert(V2)%attr(VERTX)
            X3 = XV
            Y1 = vert(V1)%attr(VERTY)
            Y2 = vert(V2)%attr(VERTY)
            Y3 = YV
            SubArea3 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

            X1 = vert(V1)%attr(VERTX)
            X2 = vert(V2)%attr(VERTX)
            X3 = vert(V3)%attr(VERTX)
            Y1 = vert(V1)%attr(VERTY)
            Y2 = vert(V2)%attr(VERTY)
            Y3 = vert(V3)%attr(VERTY)
            TotalArea = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

            IF( (SubArea1+SubArea2+SubArea3).LE.(1.001*TotalArea) )THEN

               IntoDomain = .TRUE.

            ENDIF

         ENDDO

!... If the fictional vertex does not lie inside the domain,
!... then we know that the current direction is pointing outward.
!... Reflect its wave energy back into the domain.
         IF( .NOT. IntoDomain )THEN

!... Find a representative direction for the boundary by taking
!... an average of the boundary segment(s) connected to this vertex.

!... Loop over adjacent faces to find boundary segments.
            NumBdySegs = 0
            DO jc=1,vert(ivert)%noc
               icell = vert(ivert)%cell(jc)%atti(CELLID)

               DO jf=1,cell(icell)%nof
                  iface = cell(icell)%face(jf)%atti(FACEID)
                  V1 = face(iface)%atti(FACEV1)
                  V2 = face(iface)%atti(FACEV2)

!... If this face is a boundary segment, then compute its direction.
                  IF( face(iface)%atti(FMARKER)==1 .AND. &
                     (V1==ivert .OR. V2==ivert) )THEN

                     X1 = vert(ivert)%attr(VERTX)
                     Y1 = vert(ivert)%attr(VERTY)
                     IF( V1==ivert )THEN
                        X2 = vert(V2)%attr(VERTX)
                        Y2 = vert(V2)%attr(VERTY)
                     ELSE
                        X2 = vert(V1)%attr(VERTX)
                        Y2 = vert(V1)%attr(VERTY)
                     ENDIF

                     NumBdySegs = NumBdySegs + 1
                     BdyDirTemp = ATAN( (Y2-Y1) / (X2-X1) )
                     IF( (Y2-Y1).GE.0. .AND. (X2-X1).GE.0. )THEN
                        BdyDirTemp = BdyDirTemp
                     ELSEIF( (Y2-Y1).GE.0. .AND. (X2-X1).LT.0. )THEN
                        BdyDirTemp = Pi + BdyDirTemp
                     ELSEIF( (Y2-Y1).LT.0. .AND. (X2-X1).LT.0. )THEN
                        BdyDirTemp = Pi + BdyDirTemp
                     ELSE
                        BdyDirTemp = (2.*Pi) + BdyDirTemp
                     ENDIF
                     BdyDir(NumBdySegs) = BdyDirTemp

                  ENDIF

               ENDDO

            ENDDO

!... Determine the boundary direction across which to reflect
!... the wave energy.

!... Determine the distances between the current direction
!... and the directions of the two boundary segments.
            Dist(1) = ABS(BdyDir(1)-Dir)
            IF( Dist(1).GT.Pi )THEN
               Dist(1) = ABS(BdyDir(1)-(Dir-(2.*Pi)))
            ENDIF
            Dist(2) = ABS(BdyDir(2)-Dir)
            IF( Dist(2).GT.Pi )THEN
               Dist(2) = ABS(BdyDir(2)-(Dir-(2.*Pi)))
            ENDIF

!... If the current direction is more than 90 degrees
!... from both of the boundary segments, then reflect
!... its wave energy across the average direction of the
!... boundary segments.
            IF( Dist(1).GT.(0.5*Pi) .AND. &
                Dist(2).GT.(0.5*Pi) )THEN

               BdyDirToUse = ( (BdyDir(1)-Pi) + &
                               BdyDir(2) ) / 2.0
               IF( BdyDirToUse.LT.0. )THEN
                  BdyDirToUse = BdyDirToUse + (2.*Pi)
               ENDIF

!... If the current direction is an equal distance
!... from both of the boundary segments, then reflect
!... its wave energy across the average direction of the
!... boundary segments.
            ELSEIF( Dist(1).EQ.Dist(2) )THEN

               BdyDirToUse = ( (BdyDir(1)-Pi) + &
                               BdyDir(2) ) / 2.0
               IF( BdyDirToUse.LT.0. )THEN
                  BdyDirToUse = BdyDirToUse + (2.*Pi)
               ENDIF

!... Otherwise, reflect the wave energy across
!... the closest boundary segment.
            ELSE

               IF( Dist(1).LT.Dist(2) )THEN
                  BdyDirToUse = BdyDir(1)
               ELSE
                  BdyDirToUse = BdyDir(2)
               ENDIF

            ENDIF

!... Reflect our current direction across the boundary direction.
            IncDir = 2. * BdyDirToUse - Dir
            IF( IncDir<0. )THEN
               IncDir = IncDir + (2.*Pi)
            ENDIF

!... Determine counter for which direction is IncDir.
            IncCounter = MOD( IncDir - spcdir(1,1), (2.*Pi) ) / DDIR
            IF( IncCounter<0. )THEN
               IncCounter = IncCounter + REAL(MDC)
            ENDIF

!... Determine the directions in our spectrum to which we can reflect
!... wave energy, and how much energy should go to each.
            IncD1 = 1 + INT(IncCounter)
            IncD2 = IncD1 + 1
            IF( IncD2>MDC )THEN
               IncD2 = IncD2 - MDC
            ENDIF
            W2 = IncCounter + 1. - REAL(IncD1)
            W1 = 1. - W2

!... It is sometimes possible for the incident direction to fall between
!... two directions that are not inside the active domain.  In this case,
!... one direction will be at the edge of the active domain (and near
!... the boundary), and the other direction will be just outside
!... of the domain.  We don't want to reflect energy onto that latter direction.
!... Adjust the weights to take care of this case.
            IF( IncD1==id .OR. IncD2==id )THEN
               IF( IncD1==id )THEN
                  W1 = 0.
                  W2 = 1.
               ELSE
                  W1 = 1.
                  W2 = 0.
               ENDIF
            ENDIF

!... Loop over spectrum, reflecting the wave energy as appropriate.
            DO is=1,MSC

!... Reflect components of wave energy into the domain.
               ac2temp(IncD1,is) = W1 * ac2(id,is,ivert) + ac2temp(IncD1,is)
               ac2temp(IncD2,is) = W2 * ac2(id,is,ivert) + ac2temp(IncD2,is)

            ENDDO

         ENDIF

      ENDDO

!... We must be very careful about how we add the reflected energy
!... back into the computational domain.  Our method must depend
!... on whether the direction: (1) points out of the domain,
!... (2) points into the domain through the boundary, or
!... (3) points into the domain from inside the domain.

!... First, determine if the current direction points
!... out of the domain.  Use the same method as above.
      DO id=1,MDC
         Dir = spcdir(id,1)

!... Set up a fictional vertex along the direction.
         XV = vert(ivert)%attr(VERTX) + COS(Dir) * ( 0.1 * AvgDX )
         YV = vert(ivert)%attr(VERTY) + SIN(Dir) * ( 0.1 * AvgDX )

!... Loop over the adjacent cells.  Use triangle areas
!... to determine if the fictional vertex is located
!... inside an adjacent cell.
         IntoDomain = .FALSE.
         DO jc=1,vert(ivert)%noc

            icell = vert(ivert)%cell(jc)%atti(CELLID)
            V1 = cell(icell)%atti(CELLV1)
            V2 = cell(icell)%atti(CELLV2)
            V3 = cell(icell)%atti(CELLV3)

            X1 = XV
            X2 = vert(V2)%attr(VERTX)
            X3 = vert(V3)%attr(VERTX)
            Y1 = YV
            Y2 = vert(V2)%attr(VERTY)
            Y3 = vert(V3)%attr(VERTY)
            SubArea1 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

            X1 = vert(V1)%attr(VERTX)
            X2 = XV
            X3 = vert(V3)%attr(VERTX)
            Y1 = vert(V1)%attr(VERTY)
            Y2 = YV
            Y3 = vert(V3)%attr(VERTY)
            SubArea2 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

            X1 = vert(V1)%attr(VERTX)
            X2 = vert(V2)%attr(VERTX)
            X3 = XV
            Y1 = vert(V1)%attr(VERTY)
            Y2 = vert(V2)%attr(VERTY)
            Y3 = YV
            SubArea3 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

            X1 = vert(V1)%attr(VERTX)
            X2 = vert(V2)%attr(VERTX)
            X3 = vert(V3)%attr(VERTX)
            Y1 = vert(V1)%attr(VERTY)
            Y2 = vert(V2)%attr(VERTY)
            Y3 = vert(V3)%attr(VERTY)
            TotalArea = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

            IF( (SubArea1+SubArea2+SubArea3).LE.(1.001*TotalArea) )THEN

               IntoDomain = .TRUE.

            ENDIF

         ENDDO

!... If the direction points out of the domain, then simply allow
!... its wave energy to remain constant.  Do not change its action
!... densities at all.

         IF( .NOT. IntoDomain )THEN

            CONTINUE

         ELSE

!... If the direction points into the domain, then we must determine
!... if it is coming from inside or outside of the domain.
!... Use the same method as before.
            Dir = spcdir(id,1)
            IF( Dir.LT.Pi )THEN
               Dir = Dir + Pi
            ELSE
               Dir = Dir - Pi
            ENDIF

!... Set up a fictional vertex along the direction.
            XV = vert(ivert)%attr(VERTX) + COS(Dir) * ( 0.1 * AvgDX )
            YV = vert(ivert)%attr(VERTY) + SIN(Dir) * ( 0.1 * AvgDX )

!... Loop over the adjacent cells.  Use triangle areas
!... to determine if the fictional vertex is located
!... inside an adjacent cell.
            IntoDomain = .FALSE.
            DO jc=1,vert(ivert)%noc

               icell = vert(ivert)%cell(jc)%atti(CELLID)
               V1 = cell(icell)%atti(CELLV1)
               V2 = cell(icell)%atti(CELLV2)
               V3 = cell(icell)%atti(CELLV3)

               X1 = XV
               X2 = vert(V2)%attr(VERTX)
               X3 = vert(V3)%attr(VERTX)
               Y1 = YV
               Y2 = vert(V2)%attr(VERTY)
               Y3 = vert(V3)%attr(VERTY)
               SubArea1 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

               X1 = vert(V1)%attr(VERTX)
               X2 = XV
               X3 = vert(V3)%attr(VERTX)
               Y1 = vert(V1)%attr(VERTY)
               Y2 = YV
               Y3 = vert(V3)%attr(VERTY)
               SubArea2 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

               X1 = vert(V1)%attr(VERTX)
               X2 = vert(V2)%attr(VERTX)
               X3 = XV
               Y1 = vert(V1)%attr(VERTY)
               Y2 = vert(V2)%attr(VERTY)
               Y3 = YV
               SubArea3 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

               X1 = vert(V1)%attr(VERTX)
               X2 = vert(V2)%attr(VERTX)
               X3 = vert(V3)%attr(VERTX)
               Y1 = vert(V1)%attr(VERTY)
               Y2 = vert(V2)%attr(VERTY)
               Y3 = vert(V3)%attr(VERTY)
               TotalArea = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))

               IF( (SubArea1+SubArea2+SubArea3).LE.(1.001*TotalArea) )THEN

                  IntoDomain = .TRUE.

               ENDIF

            ENDDO

!... If the direction is coming from outside of the domain,
!... then it is not receiving any natural wave energy as part
!... of the overall Swan computations.  No wave energy can propagate
!... through the boundary and into the direction.  Thus, we can
!... simply copy the reflected wave energy into these spectral bins.

            IF( .NOT. IntoDomain )THEN

               DO is=1,MSC
                  ac2(id,is,ivert) = ac2temp(id,is)
               ENDDO

!... If the direction is coming from inside the domain,
!... then it will be receiving wave energy in addition to the energy
!... that we are reflecting.  We cannot simply assign the reflected
!... energy into the spectral bins, because we must consider the
!... other wave energy.  However, we cannot simply add the reflected
!... energy, because we only want to reflect the energy once per
!... iteration.

!... So let's only reflect the energy if the wave energy
!... in this direction has just been updated.

            ELSE

               DO is=1,MSC
!                 ac2(id,is,ivert) = ac2(id,is,ivert) + ac2temp(id,is)
                  ac2(id,is,ivert) = ac2temp(id,is)
               ENDDO

            ENDIF

         ENDIF

      ENDDO

      END SUBROUTINE



      END MODULE Couple2Adcirc