C******************************************************************************
C                                                                             *
C                                                                             *
C  PROGRAM TO EXTRACT TIDAL INFORMATION FROM THE ADCIRC 2DDI TIDAL DATA BASE  *
C     FOR USE AS BOUNDARY CONDITIONS IN ANOTHER MODEL RUN 				*
C                                                                             *
C                                                                             *
C     -  DEVELOPED FOR                                                        *
C                                                                             *
C          THE U.S. ARMY ENGINEERS WATERWAYS EXPERIMENT STATION               *
C             UNDER CONTRACT WITH THE DREDGED RESEARCH PROGRAM                *
C                                                                             *
C                                                                             *
C     -  DEVELOPED BY                                                         *
C                                                                             *
C          R.A. LUETTICH, JR                                                  *
C             UNIVERSITY OF NORTH CAROLINA AT CHAPEL HILL                     *
C             INSTITUTE OF MARINE SCIENCES                                    *
C                                                                             *
C          J.J. WESTERINK                                                     *
C             DEPARTMENT OF CIVIL ENGINEERING AND GEOLOGICAL SCIENCES         *
C             UNIVERSITY OF NOTRE DAME                                        *
C                                                                             *
C                                                                             *
C     -  CODED BY                                                             *
C                                                                             *
C          R.A. LUETTICH    VERSION  1.00-1.02                                *
C              "            VERSION  1.04 for ec_95d database                 *
C          C.W. FULCHER     VERSION  1.05 bug fix (in determining the         *
C                                    ascending order of nodes the N2<N3<N1    *
C                                    and the N3<N1<N2 definitions of NO(n)    *
C                                    were switched.)  Corrected 3-22-99       *
C          R.A. LUETTICH    VERSION  1.05a dimensioned for 8 constituents     * 
C	     E.A. SPARGO	    VERISON  1.07a for ec2001 tidal database eta,u,v  *
C                                                                             *
C******************************************************************************
C                                                                             *
C                                                                             *
C     -  PARAMETERS WHICH MUST BE SET WITHIN THE MAIN CODE AND THE            *
C          SUBROUTINES TO CONTROL THE DIMENSIONING OF ARRAYS ARE AS FOLLOWS   *
C          (TO FIND THE LINES WHICH MUST BE SET, FIND "CUSER" COMMENT LINES): *
C                                                                             *
C          MNP = MAXIMUM NUMBER OF NODAL POINTS                               *
C          MNHARF = MAXIMUM NUMBER OF CONSTITUENTS IN DATA BASE               *
C                                                                             *
C******************************************************************************
C                                                                             *
C     -  STANDARD INPUT FILES ARE AS FOLLOWS:                                 *
C                                                                             *
C          UNIT 12 : FILE LISTING THE X AND Y COORDINATES OF ALL LOCATIONS FOR*
C                    INTERPOLATED OUTPUT.		                              *
C          UNIT 14 : ADCIRC FINITE ELEMENT GRID FILE                          *
C          UNIT 105 : HARMONIC CONSTITUENT ELEVATIONS AND VELOCITIES AT ALL   *
C                     NODES.  THIS FILE WAS GENERATED BY RUNNING ASTROEX.FOR  *
C                     ON THE UNIT 53 AND UNIT 54 FILES GENERATED BY ADCIRC.   *
C												      *
C     -  DESCRIPTION OF INPUT VARIABLES READ IN FROM UNIT 12 (tides.in)		*
C													*
C          NOUT = NUMBER OF OUTPUT POINTS IN THE FILE					*
C          XOUT(I),YOUT(I), I=1,NOUT  =  X AND Y COORDINATES OF EACH POINT TO *
C               INTERPOLATE OUTPUT AT							*
C                                                                             *
C     -  DESCRIPTION OF INPUT VARIABLES READ IN FROM UNIT 14 (ec2001.grd)     *
C                                                                             *
C          AGRID = ALPAHANUMERIC GRID IDENTIFICATION  (<=24 CHARACTERS)       *
C          NE,NP = NUMBER OF ELEMENTS AND NUMBER OF NODAL POINTS RESPECTIVELY *
C          JKI,X(JKI),Y(JKI),DP(JKI) , JKI=1,NP  = NODE NUMBER, X AND Y       *
C               COORDINATES, BATHYMETRIC VALUE ; NODES CAN BE INPUT IN        *
C               ANY ORDER ; IF ICS=1, THEN X,Y REPRESENT STANDARD             *
C               CARTESIAN COORDINATES SPECIFIED IN LENGTH UNITS.              *
C               IF ICS=2, THEN X,Y REPRESENT DEGREES LONGITUDE                *
C               (DEGREES EAST OF GREENWICH IS POSITIVE AND DEGREES WEST OF    *
C               GREENWICH IS NEGATIVE) AND DEGREES LATITUDE (DEGREES          *
C               NORTH OF THE EQUATOR BEING POSITIVE AND DEGREES SOUTH OF THE  *
C               EQUATOR IS NEGATIVE) RESPECTIVELY                             *
C          JKI,NHY,NM(JKI,1),NM(JKI,2),NM(JKI,3) , JKI=1,NE  = ELEMENT        *
C               NUMBER, ELEMENT TYPE, AND ELEMENT CONNECTIVITY SPECIFIED      *
C               WITH A COUNTERCLOCKWISE ORIENTATION ; NOTE THAT THE ELEMENT   *
C               TYPE IS NOT AN ACTIVE VARIABLE AND THAT ONLY 3 NODE           *
C               LINEAR TRIANGLES ARE OPERATIONAL IN THIS VERSION OF THE CODE; *
C               ELEMENTS CAN BE READ IN IN ANY ORDER                          *
C                                                                             *
C     -  DESCRIPTION OF INPUT VARIABLES READ FROM UNIT 105 (ec2001.tdb)       *
C                                                                             *
C          NHARFR = NUMBER OF CONSTITUENTS IN FILE                            *
C          J=1,NHARF                                                          *
C             FREQUENCY, NODAL FACTOR, EQUILIBRIUM ARGUEMENT, CONSTITUENT NAME*
C          NP = NUMBER OF NODES IN GRID                                       *
C          J=1,NP                                                             *
C             J,(HAMP(I),HPHASE(I) I=1,NHARFR)                                *
C               (UAMP(I),UPHASE(I) I=1,NHARFR)                                *
C               (VAMP(I),VPHASE(I) I=1,NHARFR)                                *
C                                                                             *
C******************************************************************************
C                                                                             *
C     -  STANDARD OUTPUT FILE IS AS FOLLOWS:                                  *
C                                                                             *
C          UNIT 1 : TIDAL CONSTITUENT RESULTS                                 *
C          UNIT 2 : DIAGNOSTIC INFORMATION						*
C                                                                             *
C                                                                             *
C     -  DESCRIPTION OF OUTPUT DATA WRITTEN TO UNIT 1  (tides.out)            *
C          self explanitory.  see output file.                                *
C                                                                             *
C     -  DESCRIPTION OF OUTPUT DATA WRITTEN TO UNIT 2  (tides.dia)            *
C          self explanitory.  see output file.                                *
C                                                                             *
C******************************************************************************

CUSER...
CUSER...SET PARAMETER STATEMENTS
CUSER....THE PARAMETERS HAVE TO BE RE-SET PRIOR TO EXECUTION TO INSURE
CUSER....SUFFICIENT SPACE HAS BEEN ALLOCATED FOR A SPECIFIC PROBLEM
CUSER....PARAMETER STATEMENTS MUST ALSO BE ADJUSTED IN THE SUBROUTINES
CUSER....AND FUNCTION STATEMENTS
CUSER...
        PARAMETER(MNP=400000)
        PARAMETER(MNHARF=10)
CUSER...
CUSER...END OF PARAMETER STATEMENTS
CUSER....

C...
C...DIMENSION ALL ARRAYS AND DEFINE COMMON BLOCKS
C...
      DIMENSION X(MNP),Y(MNP),XOUT(MNP),YOUT(MNP)
	DIMENSION NM1(2*MNP),NM2(2*MNP),NM3(2*MNP)
      DIMENSION FREQ(MNHARF),NFACT(MNHARF),EQARG(MNHARF),
     &          NAME(MNHARF)
      DIMENSION EAMP(3,MNHARF),EPHA(3,MNHARF),ETAMP(MNHARF,MNP),
     &          ETPHA(MNHARF,MNP)
      DIMENSION UAMP(3,MNHARF),UPHA(3,MNHARF),UTAMP(MNHARF,MNP),
     &          UTPHA(MNHARF,MNP)
      DIMENSION VAMP(3,MNHARF),VPHA(3,MNHARF),VTAMP(MNHARF,MNP),
     &          VTPHA(MNHARF,MNP)
      DIMENSION NN(3),NO(3)
      CHARACTER*10 NAME
      REAL NFACT

      OPEN(1,FILE='tides.out')
	OPEN(2,FILE='tides.dia')
	OPEN(12,FILE='tides.in')
      OPEN(14,FILE='ec2001.grd')

C...
C...DEFINE PI
C...
      PI=3.141592653589793

C...
C...OUTPUT FILE HEADER
C...
      WRITE(1,3900)
 3900 FORMAT(//,8X,'Constituent',13x,'Elevation',16x,'East Velocity',
     &          12X,'North Velocity')
      WRITE(1,3901)
 3901 FORMAT(11x,'Name/',8x,3(6x,'Amplitude     Phase'))
      WRITE(1,3902)
 3902 FORMAT(5X,'Lon',10X,'Lat',12X,'(m)',8X,'(deg)',2(9x,'(m/s)',6X,
     &                                                  '(deg)'),/)

C...
C...READ INPUT LOCATIONS FOR HARMONIC ANALYSIS OUTPUT
C...
      READ(12,*) NOUT
      DO I=1,NOUT
        READ(12,*) XOUT(I),YOUT(I)
        WRITE(2,1000) XOUT(I),YOUT(I)
 1000   FORMAT(//,' OUTPUT WILL BE GENERATED FOR THE POSITION: ',/,
     &            5X,F11.6,' E LONGITUDE , ',F11.6,' N LATITUDE',//)
	  END DO
 
C...
C...INPUT GRID INFORMATION FROM UNIT 14
C...
      READ(14,'(A24)') AGRID
      READ(14,*) NE,NP


	open(25,file='diag.out')

C...NODAL COORDINATES

      WRITE(*,*) 'READING NODAL COORDINATES......'
      WRITE(*,*) ' '
      DO I=1,NP
        READ(14,*) JKI,X(JKI),Y(JKI)
        END DO

C....CONNECTIVITY TABLE 

      WRITE(*,*) 'READING ELEMENT TABLE......'
      WRITE(*,*) ' '
      DO I=1,NE
        READ(14,*) JKI,NHY,NM1(JKI),NM2(JKI),NM3(JKI)
	  END DO

C...
C....COMPUTE ELEMENT IN WHICH EACH LOCATION LIES.
C...
      WRITE(*,*) 'COMPUTING LOCATION AND INTERPOLATING......'
	WRITE(*,*) ' '
    	DO I=1,NOUT
        AEMIN=1.0E+25
        KMIN=0
        X4=XOUT(I)
        Y4=YOUT(I)
        DO K=1,NE
          X1=X(NM1(K))
          X2=X(NM2(K))
          X3=X(NM3(K))
          Y1=Y(NM1(K))
          Y2=Y(NM2(K))
          Y3=Y(NM3(K))
 	    AREAS=(X1-X3)*(Y2-Y3)-(X3-X2)*(Y3-Y1)
          A1=(X4-X3)*(Y2-Y3)-(X3-X2)*(Y3-Y4)
          A2=(X1-X3)*(Y4-Y3)-(X3-X4)*(Y3-Y1)
          A3=(X1-X4)*(Y2-Y4)-(X4-X2)*(Y4-Y1)
          AA=ABS(A1)+ABS(A2)+ABS(A3)
          AE=ABS(AA-AREAS)/AREAS
          IF(AE.LT.AEMIN) THEN
            AEMIN=AE
            NNE=K
            N1=NM1(K)
            N2=NM2(K)
            N3=NM3(K)
            AREA=AREAS
            ENDIF
          END DO

C......PRINT WARNING IF NODE OUTSIDE THE DOMAIN

        IF(AEMIN.GT.1.0E-5) THEN                  !OUTSIDE AN ELEMENT
          WRITE(*,2000) AEMIN
          WRITE(2,2000) AEMIN
 2000     FORMAT(///,' WARNING -  SPECIFIED LOCATION DOES NOT LIE',
     &           ' WITHIN ANY ELEMENT IN THE DOMAIN.',/,' CHECK THE',
     &           ' LONGITUDE AND LATITUDE FOR THIS LOCATION',
     &           ' PROGRAM WILL ESTIMATE NEAREST ELEMENT',/,' THE',
     &           ' PROXIMITY INDEX FOR THIS LOCATION EQUALS ',E15.6)
	    ENDIF

C......COMPUTE INFORMATION REQUIRED TO INTERPOLATE AT OUTPUT LOCATION

        X1=X(N1)
        X2=X(N2)
        X3=X(N3)
        Y1=Y(N1)
        Y2=Y(N2)
        Y3=Y(N3)
        STA2=((X4-X1)*(Y3-Y1)-(Y4-Y1)*(X3-X1))/AREA
        STA3=(-(X4-X1)*(Y2-Y1)+(Y4-Y1)*(X2-X1))/AREA
        STA1=1.0-STA2-STA3


C......DETERMINE ASCENDING ORDER OF NODES

        IF((N1.LT.N2).AND.(N2.LT.N3)) THEN
          NN(1)=N1
          NN(2)=N2
          NN(3)=N3
          NO(1)=1
          NO(2)=2
          NO(3)=3
          GOTO 200
          ENDIF
        IF((N2.LT.N1).AND.(N1.LT.N3)) THEN
          NN(1)=N2
          NN(2)=N1
          NN(3)=N3
          NO(2)=1
          NO(1)=2
          NO(3)=3
          GOTO 200
          ENDIF
        IF((N1.LT.N3).AND.(N3.LT.N2)) THEN
          NN(1)=N1
          NN(2)=N3
          NN(3)=N2
          NO(1)=1
          NO(3)=2
          NO(2)=3
          GOTO 200
          ENDIF
c     Bug fix for v 1.05 on 3-22-99 starts here:
        IF((N2.LT.N3).AND.(N3.LT.N1)) THEN
          NN(1)=N2
          NN(2)=N3
          NN(3)=N1
          NO(3)=1
          NO(1)=2
          NO(2)=3
          GOTO 200
          ENDIF
        IF((N3.LT.N1).AND.(N1.LT.N2)) THEN
          NN(1)=N3
          NN(2)=N1
          NN(3)=N2
          NO(2)=1
          NO(3)=2
          NO(1)=3
          GOTO 200
          ENDIF
c     End of bug fix (the above definitions of NO(.) were switched in previous version
        IF((N3.LT.N2).AND.(N2.LT.N1)) THEN
          NN(1)=N3
          NN(2)=N2
          NN(3)=N1
          NO(3)=1
          NO(2)=2
          NO(1)=3
          GOTO 200
          ENDIF

C......WRITE THE ELEMENT CONTAINING THE OUTPUT LOCATION

  200   WRITE(*,2100) NNE,NN(1),NN(2),NN(3)
        WRITE(2,2100) NNE,NN(1),NN(2),NN(3)
 2100   FORMAT(' SPECIFIED LOCATION WAS FOUND IN ELEMENT ',I8,/,
     &      5X,' WHICH IS MADE UP OF NODES ',3I8,//)

C......READ HARMONIC CONSTITUENT INFORMATION

        OPEN(105,FILE='ec2001.tdb')

        READ(105,*) NHC
        DO J=1,NHC
          READ(105,3100) FREQ(J),NFACT(J),EQARG(J),NAME(J)
 3100     FORMAT(E21.10,F11.7,F13.8,2X,A10)
          END DO
        READ(105,*) NPP

 3200   FORMAT(//)
        DO K=1,3
          IBEG=1
          IF(K.GE.2) IBEG=NN(K-1)+1
          DO II=IBEG,NN(K)-1
            READ(105,3200)
            END DO
          READ(105,*) NODE,(EAMP(NO(K),J),EPHA(NO(K),J),J=1,NHC)
          READ(105,*)      (UAMP(NO(K),J),UPHA(NO(K),J),J=1,NHC)
          READ(105,*)      (VAMP(NO(K),J),VPHA(NO(K),J),J=1,NHC)
          IF(NODE.NE.NN(K)) THEN
            WRITE(*,*) 'ERROR FINDING NODE ',NN(K),' PROGRAM HAS FOUND',
     &                ' NODE ',NODE,' IN ITS PLACE'
            WRITE(*,*)' PROGRAM WILL NOW BE TERMINATED'
            STOP
            ENDIF
          DO J=1,NHC
            EPHA(NO(K),J)=PI*EPHA(NO(K),J)/180.
            UPHA(NO(K),J)=PI*UPHA(NO(K),J)/180.
            VPHA(NO(K),J)=PI*VPHA(NO(K),J)/180.
            END DO
         END DO

	  CLOSE(105)

C......COMPUTE HARMONIC CONSTITUENTS AT THE OUTPUT LOCATION

        DO J=1,NHC
          E1R=EAMP(1,J)*COS(EPHA(1,J))
          E1I=EAMP(1,J)*SIN(EPHA(1,J))
          E2R=EAMP(2,J)*COS(EPHA(2,J))
          E2I=EAMP(2,J)*SIN(EPHA(2,J))
          E3R=EAMP(3,J)*COS(EPHA(3,J))
          E3I=EAMP(3,J)*SIN(EPHA(3,J))
          U1R=UAMP(1,J)*COS(UPHA(1,J))
          U1I=UAMP(1,J)*SIN(UPHA(1,J))
          U2R=UAMP(2,J)*COS(UPHA(2,J))
          U2I=UAMP(2,J)*SIN(UPHA(2,J))
          U3R=UAMP(3,J)*COS(UPHA(3,J))
          U3I=UAMP(3,J)*SIN(UPHA(3,J))
          V1R=VAMP(1,J)*COS(VPHA(1,J))
          V1I=VAMP(1,J)*SIN(VPHA(1,J))
          V2R=VAMP(2,J)*COS(VPHA(2,J))
          V2I=VAMP(2,J)*SIN(VPHA(2,J))
          V3R=VAMP(3,J)*COS(VPHA(3,J))
          V3I=VAMP(3,J)*SIN(VPHA(3,J))
          ETR=E1R*STA1+E2R*STA2+E3R*STA3
          ETI=E1I*STA1+E2I*STA2+E3I*STA3
          UTR=U1R*STA1+U2R*STA2+U3R*STA3
          UTI=U1I*STA1+U2I*STA2+U3I*STA3
          VTR=V1R*STA1+V2R*STA2+V3R*STA3
          VTI=V1I*STA1+V2I*STA2+V3I*STA3
          ETAMP(J,I)=SQRT(ETR*ETR+ETI*ETI)
          UTAMP(J,I)=SQRT(UTR*UTR+UTI*UTI)
          VTAMP(J,I)=SQRT(VTR*VTR+VTI*VTI)
          IF(ETAMP(J,I).EQ.0.) THEN
             ETPHA(J,I)=0.
             ELSE
             ETPHA(J,I)=180.*ACOS(ETR/ETAMP(J,I))/PI
             IF(ETI.LT.0.) ETPHA(J,I)=360.-ETPHA(J,I)
             ENDIF
          IF(UTAMP(J,I).EQ.0.) THEN
             UTPHA(J,I)=0.
             ELSE
             UTPHA(J,I)=180.*ACOS(UTR/UTAMP(J,I))/PI
             IF(UTI.LT.0.) UTPHA(J,I)=360.-UTPHA(J,I)
             ENDIF
          IF(VTAMP(J,I).EQ.0.) THEN
             VTPHA(J,I)=0.
             ELSE
             VTPHA(J,I)=180.*ACOS(VTR/VTAMP(J,I))/PI
             IF(VTI.LT.0.) VTPHA(J,I)=360.-VTPHA(J,I)
             ENDIF
          END DO

        END DO

C......WRITE THE RESULTS

      DO J=1,NHC
	  WRITE(1,3999) NAME(J)
 3999   FORMAT(//,1X,A10,/)
	  DO I=1,NOUT
          WRITE(1,4000) XOUT(I),YOUT(I),ETAMP(J,I),ETPHA(J,I),
     &           		UTAMP(J,I),UTPHA(J,I),VTAMP(J,I),VTPHA(J,I)
 4000     FORMAT(1X,2(F11.6,2X),3(E12.5,2X,F8.3,3X))
          END DO
	  END DO


      WRITE(*,*) ' '                   
      WRITE(*,4100)
 4100 FORMAT('**** RESULTS HAVE BEEN STORED IN FILE: TIDES.OUT ****')

      CLOSE(14)
      CLOSE(1)

      STOP
      END
