cccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM truxton1

cccccccccccccccccccccccccccccccccccccccccccccc
c
c 15 dec 98, drl:
c  Reversed the definition of (uerr, verr) at the point of its creation:  
c     now = observed velocity - model velocity 
c     was = model velocity    - observed velocity
c  This change made so that the 'error' is now the unexplained observation.
c  In the absence of any model result, the error is equal to the observation.
c  These arrays created solely for use in subroutine truxout; 
c    todays change has no internal effect.
c
c dec 15 1998
c FEW, DJM modified getobs to screen out obs not in mesh
c
cccccccccccccccccccccccccccccccccccccccccccccc

C
c june 10 '97, drl: eliminated the .bel input to subr. BuildW.
c   Uses boundary node list LRVS instead, assuming the nodes are 
c   numbered sequentially, counterclockwise.  
c   Subr. BNodeOrder is intended to guarantee this ordering, 
c   after reading the .inp file; but it is presently null and 
c   writes a message to that effect.
c
c june 9 1997, drl: added TruxOut and new arrays for same:
C   complex ZCamp(nbcdim,nfdim) 
C   real Uerr(nobsdim), Verr(nobsdim)
C   real Upartial(nobsdim, nfdim), Vpartial(nobsdim,nfdim) 
C
c june 5 1997, drl:
c     f5xcm3SEI_3 with I/O streamlined and some minor cleanup.
c     include file changed to TRUXTON1ARRAYS.DIM, throughout.
c       Fundy5 include plus these extras:
c         NBCdim  >= number of boundary nodes.
c         NOBSdim >= number of velocity observations
c         NFdim   >= number of tidal constituents to be used            
c         Additional parameters derived from them:
c            NOBS2dim  = 2*NOBSdim
c            NBC2dim   = 2*NBCdim
c            NBC2NFDIM = NBC2dim*NFdim
c
cccccccccccccccccccccccccccccccccccccccc
c
c      PROGRAM f5xcm3SEI_3
c
c may 23 '97 drl: added covariance calculation to f5xcm3SEI
c
c    1)  Add Cov(i,j) to main: 
c           this is the real covariance of the real and -im parts of 
c           the bc's in j-space (svd or wls) or k-space (fourier)
c        Cov is computed in invertWLS3 and invert3, along with SE as
c           before. 
c        Cov is normalized by mean squared velocity
c        SE is normalized by rms velocity, i.e. se in the
c           detiding ms. SE(i)=[sqrt{Cov(i,i)}]
c
c    2)  Add CCov(i,j) to main: this is the complex covariance of the
c           bc's, derived by combining the proper parts of Cov
c        Do this by sending Cov to subr. CmplxCov; get back CCov
c
c    3)  Convert Subr.'s SVD_se and Fourier_se to operate on Ccov
c
c    4)  Convert the j-space se's to the correct form taking account 
c           of covariance among k-space results.
c
cccccccccccccccccccccccccccccccccccccccc
c
c      PROGRAM f5xcm3SEI
c
c may 13 '97 drl: added f5xcm3SE features to f5xcm3I
c
c        f5xcm3SEI calculates std error and coef. of variation 
c        Basic info computed from svd vectors in invert2 is returned to main
c        For SVD or Fourier inverse:
c            New subr's fourier_se and svd_se compute and report std errors.
c            Reusults written to file se.out.  
c            In Fourier case, results also converted from k- to j-space.
c        For WLS: above is invalid; se.out not written.
c
c        Subr. invert2 same as invert1 but passes se array back to main
c        New array: real se(nbc2NFdim)
c
cccccccccccccccccccccccccccccccccccccccc
c
c may 7 '97 drl: f5xcm3I.f implements inference of selected frequencies. 
c        Basic data input is at subr. GetFreqI which replaces subr. GetFreq.
c        Call GetFreqI(NF, NF1, NF2, Freq, Infer, Amprat, NFDIM)
c             NF = total # of frequencies = NF1 + NF2
c             NF1= # of primary frequencies
c             NF2= # of frequencies inferred from the primary frequencies
c                        (nested inference not allowed)
c             Freq(i) = frequency of ith constituent. Primaries: i = 1,     NF1
c                                                     Inferred:  i = NF1+1, NF
c             Infer(i) = index of frequency to be used for inference for
c                        frequency i.  Infer(i)=i for i=1, NF1
c                                      Infer(i).le.NF1 for i=NF1+1, NF
c             Amprat(i) = amplitude ratio for inference of constituent i.
c                         Amprat(i) = cmplx(1., 0.) for i=1, NF1
c                         Amprat(i) = data          for i=NF1+1, NF
c             NFDIM.ge.NF
c
c New declarations: 
c      Integer Infer(NFDIM), NF1, NF2
c      Complex Amprat(NFDIM), ZZ1, ZZ2, ZZ3, ZZ4, ZZ5, one
c 
cccccccccccccccccccccccccccccccccccccccc
c
c may 7 '97 drl: added changes to subr. invert1.  It now calculates the 
c        Standard error of parameter estimates and writes them to invert.out.  
c        No change outside of invert1.
c
cccccccccccccccccccccccccccccccccccccccc
c
c feb 21 '97 drl
c from feb 10 97 version
c repaired to account for negative and positive parts of the spectrum
c
c User is asked for Nmode, the half-wdith of the spatial spectrum. 
c The spatial spectrum used will be 
c    -Nmode <= k <= Nmode
c Its length (width) is Lspec = 2*Nmode +1. 
c 
c In the case where 2*Nmode = NBC, the spectrum is truncated 
c to length Lspec = 2*Nmode by ignoring the top mode k=Nmode.  That mode 
c is indistinguishable from the bottom mode k=-Nmode, and therefore not needed.
c Under these conditions the spectrum is complete.
c
c The FP columns will be in the natural order of ascending wavenumber, 
c beginning with k=-Nmode and working up.
c  
cccccccccccccccccccccccccccccccccccccccc
c
c feb 10 '97 drl
c uses fourier basis around-bank for unknown elevations at each frequency
c
c   - Adds the COMPLEX fourier projection matrix FP:
c        {zeta} = [FP]{Zeta}
c     zeta is elevation, Zeta is its spatial fourier transform around-bank.
c     The around-bank coordinate is node number, not distance.
c     There are Nmode fourier modes; [FP] is dimensioned (Nbc, Nmode)
c
c   - Nmode is queried in the input stream.  Nmode = 0 indicates that this
c     feature is not to be used.
c
c   - Mode 1 is spatially constant (wavenumber 0)
c
c from feb 10'97 version of f5xcm1.f
c 
cccccccccccccccccccccccccccccccccccccccc
c
c feb 8 '97 drl: 
c
c   - Adds general least squares weighting matrix Wls, 
c       which is passed to subr. InvertWLS1.
c   - Wls: positive definite symmetric, size (2*nbc*nfreq, 2*nbc*nfreq)
c   - Call to subr. BuildW to build Wls
c   - BuildW in turn calls BuildStS
c   - BuildStS reads boundary element info for assembly of S, then S*S, 
c       for a single frequency
c   - BuildW then makes W as block diagonal matrix 
c       comprising w0[I] + w1{StS] for each frequency.
c   - The weights are queried in the input stream.  
c        They regulate the relative importance of the WLS error criterion:
c          w0 is the weight assigned to ||elevation||
c          w1 is the weight assigned to bndry integral of ||boundary slope||
c           1 is the weight assigned to ||errors in velocity||
c   - Note the dimensions:
c          w0: [velocity^2/elevation^2]
c          w1: [velocity^2/elevation^2]*ds       
c                             
c  also added feb 8 '97: 
C
c   - Rcond added to input stream, passed to both InvertWLS1 and Invert1
c
cccccccccccccccccccccccccccccccccccccccc
C
C f5xcm.f
C     allows for multiple frequencies 
C based on f5xc.f 
C    finds the inverse Type I problem by SVD
C    assumes a periodic solution is sampled in the time domain
C Jan 20 '97, drl
C Jan 24 '97, cgh
C Jan 25 '97, cgh  dowvert flag for VELOCITY5xc
C                  GETFREQ for getting frequencies from user.
C
cccccccccccccccccccccccccccccccccccccccc
c  
c Original text from FUNDY5 follows:
c      
C-----------------------------------------------------------------------
C Harmonic solution of 3-D linearized shallow water wave equations.
C
C Method: Finite Element Solution - Integral lumping on linear elements
C                   2-D triangular horizontal elements
C                   1-D vertical elements
C
C Forcings: Barotropic Tides, Boundary Transport, Wind, and/or 
C                    Baroclinic.
C
C Units:  MKS
C
C Coordinate System: Horizontal dimensions - X(I),Y(I)
C                    Vertical dimension - Z(I,J) is positive upward
C                                         Z(I,J)=-HDOWN(I) at the bottom
C                                         Z(I,J)=0.0 at J=NNV
C-----------------------------------------------------------------------
C Arrays and Constants:
C
C   Index range: I=1,NN, J=1,NNV, L=1,NE, K=1,NBE, M=1,NGEO
C 
C   1-D Vertical Arrays (Z) are dimensioned by vertical node number (J).
C          J=1 at the bottom of the water column
C          J=NNV at the top of the water column
C   2-D Horizontal Arrays (X,Y) are dimensioned either by horizontal
C          node number (I) or horizontal element number (L). 
C   3-D Arrays are dimensioned by horizontal node number, vertical node
C          number (I,J).
C
C   Loaded/Computed in Subroutine READINPF5
C          NN - number of horizontal nodes
C          NE - the number of horizontal 2-D triangular elements
C          X(I),Y(I) - cartesian nodal coordinates
C          IN(3,L) - LABLIB STANDARD triangular element incidence list
C          HDOWN(I) - Bathymetric depth at nodes
C          SUMX(I) - cosine value for boundary nodes where V.n is 
C                          specified
C          SUMY(I) - sine value for boundary nodes where V.n is 
C                          specified
C          DSBDRY(I) - surface integral information
C          INBE(4,K) - LABLIB STANDARD boundary element incidence list
C          NBE - The number of boundary elements
C          INGEO(2,M) - LABLIB STANDARD geostrophic boundary element
C                          incidence list
C          NGEO - The number of geostrophic boundary elements
C          NORM(I) - boundary node numbers where normal velocity is set
C                          I =1,NVN       are V.n .NE. 0.0
C                          I =KST,NNDIM   are V.n=0.0
C          LRVS(I) - node numbers of set depth and set velocity nodes
C                          I=1,NR         are set depth nodes
C                          I=NR+1,NR+NPCN are V=0.0
C                          I=LST,NNDIM    are V .NE. 0.0
C          COR - Coriolis Parameter
C          OMEGACB - Radian frequency (rad/sec)
C   Loaded/Computed in Subroutine VERTGRID5
C          NNV - number of vertical nodes
C          AK(I) - Linear bottom stress coefficient
C          ENZ(I,J) - Vertical viscosity
C          Z(I,J) - Vertical coordinate of 3-D nodes
C          RHOX(I,J) - x-component of baroclinic forcing
C          RHOY(I,J) - y-component of baroclinic forcing
C          HRBARELX(L) - HDOWN(L)*(Vertical integral of RHOX(L,J) 
C                          over entire water column)
C          HRBARELX(L) - HDOWN(L)*(Vertical integral of RHOY(L,J)
C                          over entire water column)
C   Loaded/Computed in Subroutine ATMOS
C          ATMX(I) - x-component of atmospheric forcing (H(I)*PSIX(I))
C          ATMY(I) - y-component of atmospheric forcing (H(I)*PSIY(I))
C   Computed in Subroutine ELEVATION5
C          RHS(I) - the free surface elevation at nodes
C   Computed in Subroutine BUILDPPXPPY
C          IQ(I),JQ(I*NFTR) - SPRSPAK STANDARD pointers for PPX and PPY
C          PPX(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix
C          PPY(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix 
C   Computed in Subroutine VELOCITY5
C          UVERT(I,J),VVERT(I,J),WVERT(I,J) - x,y,and z components of 
C                          velocity
C          UBAR(I) - Vertical average of UVERT(I,J)
C          VBAR(I) - Vertical average of VVERT(I,J)
C-----------------------------------------------------------------------
C  Include File:  The file 'FUNDY5xcARRAYS.DIM' is included by each 
C                  routine which requires local variable assignment.
C                  The dimensioning parameters specified must be at
C                  least as large as the following solution dependent
C                  quantities.
C
C                  NNDIM  .GTE. NN
C                  NEDIM  .GTE. NE
C                  NNVDIM .GTE. NNV
C                  NBWE   .GTE. NBW
C                  NFTRDIM .GTE. NFTR
C
C                  NN - Number of horizontal nodes
C                  NE - Number of horizontal elements
C                  NNV - Number of vertical nodes under each horizontal
C                          node
C                  NBW - Nonsymmetric bandwidth for elevation matrix
C                  NFTR - SPRSPAK sparsity factor
C-----------------------------------------------------------------------
C History:  Developed by Daniel R. Lynch and Christopher E. Naimie
C           Dartmouth College
C
C FUNDY5: FROM FUNDY4, 3/25/93
C   FROM FUNDY3, 8/14/90 (BEGAN AS FUNDY3A, 6/28/90)
C     FROM: FUNDY3, 6/22/90
C       FROM: FUNDY2, 3/22/90
C         FROM: FUNDY1, 3/2/90
C           FROM: FUNDY, THE 9/27/89 MARKUP FOR FIRST BAROCLINIC TRY
C             FROM: CISCO.PAPERA.CASE3A (3-D IJNMF PART I)
C               FROM: WAVETLS (SEPTEMBER 10,1984) (2-D VERSION)
C
C-----------------------------------------------------------------------
C
C
C Include file for parameter statements
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Arrays for TRUXTON
C
      Real Xobs(NOBSdim), Yobs(NOBSdim)
      Real Tobs(NOBSdim), Tobsgmt0(NOBSdim)
      Real Ubarobs(NOBSdim), Vbarobs(NOBSdim)
      Real Zmaxobs(NOBSdim), Zminobs(NOBSdim), Zbotobs(NOBSdim)
      Real Els(NOBS2dim), Esamp(NOBS2dim)
      Real  Als(NOBS2dim,NBC2NFdim), Bls(NOBS2dim) 
      Complex Ubarsamp(NOBSdim), Vbarsamp(NOBSdim)
      Real Zls(NBC2NFdim)
      Real Zlsone(NBC2dim)
      Real UBARls(NOBSdim),VBARls(NOBSdim)
      Complex QG(NBWE,NNDIM)
      complex eye, ph
C
c june 9 1997, drl: new arrays for TruxOut:
C
      complex ZCamp(nbcdim,nfdim) 
      real Uerr(nobsdim), Verr(nobsdim)
      real Upartial(nobsdim, nfdim), Vpartial(nobsdim,nfdim) 
c      
c may 23 '97 drl:
c  
      real Cov(NBC2NFdim,NBC2NFdim)
      complex CCov(NBCdim*NFdim,NBCdim*NFdim)
c      
c 13 May '97 drl: array for std error coming from subr. invert2  
c
      real se(NBC2NFdim)    
c
c may 7 '97 drl for f5xcm3I:
c
      Integer Infer(NFDIM), NF1, NF2
      Complex Amprat(NFDIM), ZZ1, ZZ2, ZZ3, ZZ4, ZZ5, one
c
c 8 feb '97 drl: general square matrix Wls for the least squares problem: 
c
      real Wls(NBC2NFdim,NBC2NFdim)
c
c 10 feb '97 drl: complex matrix of fourier projections FP:
c
      complex FP(nbcdim, nbcdim), ffact, ffact1
c
c f5xcm 
c
      REAL    FREQ(NFDIM)
      INTEGER IFREQ
      LOGICAL DOWVERT
C
C-----------------------------------------------------------------------
C
C End of TRUXTON1 declarations; proceed with FUNDY5 declarations
C-----------------------------------------------------------------------
CEN
CEN Arrays filled by Subroutine F5READ
CEN
CEN      INTEGER INBE(2,NBEDIM),IN(3,NEDIM)
CEN      INTEGER NORM(NNDIM),LRVS(NNDIM)
CEN      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),SUMX(NNDIM),
CEN     &     SUMY(NNDIM),DSBDRY(NNDIM)
C
C Arrays filled by Subroutine INPUTF6
C
      INTEGER IN(3,NEDIM),INBE(4,NBEDIM),INGEO(2,NBEDIM),
     &NORM(NNDIM),LRVS(NNDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),SUMX(NNDIM),SUMY(NNDIM),
     &DSBDRY(NNDIM),ZL(NLEVDIM)
      CHARACTER*72 FILINQ,MESHNAME
C
C Arrays filled by Subroutine VERTGRID5
C
      REAL AK(NNDIM)
      REAL Z(NNDIM,NNVDIM),ENZ(NNDIM,NNVDIM)
      COMPLEX HRBARELX(NEDIM),HRBARELY(NEDIM)
      COMPLEX RHOX(NNDIM,NNVDIM),RHOY(NNDIM,NNVDIM)
C
C Arrays filled by Subroutine ATMOS
C
      COMPLEX ATMX(NNDIM),ATMY(NNDIM)
C
C Elevation solution array computed in Subroutine ELEVATION
C
      COMPLEX RHS(NNDIM)
C
C Sparse matrix arrays filled by Subroutine BUILDPPXPPY
C
      INTEGER IQ(NNDIM),JQ(NFTRDIM*NNDIM)
      REAL PPX(NFTRDIM*NNDIM),PPY(NFTRDIM*NNDIM)
C
C Velocity solution arrays computed by Subroutine VELOCITY
C
      COMPLEX UBAR(NNDIM),VBAR(NNDIM) 
      COMPLEX UVERT(NNDIM,NNVDIM),VVERT(NNDIM,NNVDIM),
     &        WVERT(NNDIM,NNVDIM)
C
C-----------------------------------------------------------------------
C
C End of Declarations; Read and Initialize
C End of Declarations; Read and Initialize
C
C-----------------------------------------------------------------------
CEN
CEN Load Fundy5 input deck information
CEN
CEN      CALL READINPF5(NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
CEN     &INBE,NBE,NORM,NVN,KST,LRVS,NR,NPCN,LST,COR,OMEGACB)
C-----------------------------------------------------------------------
C INPUT SECTION
C
C Objectives: Read .inq file, input mesh geometry data, open and begin
C writing to echo file, and set all simulation control parameters.
C
C Enter the name of the input file (".inq") of the simulation
C
      WRITE(*,'(A)')'Enter the name of the input file (".inq")'
      READ (*,'(A)')FILINQ
C
C Load run information from .inq and associated input files
C
      CALL INPUTF6(FILINQ,
     &DEG,NNV,
     &MESHNAME,NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &INBE,NBE,INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,NLEV,ZL)
C
C Standard physical constants (REAL)
C
CEN      G=9.806
      OMEGACB=0.7292E-4
CEN
CEN Call SUBROUTINE SETF6 to enable the user to change standard 
CEN physical constants and mass variable character strings
CEN
CEN      CALL SETF6(G,OMEGACB)
C
C Set coriolis parameter and number of elements in the vertical
C
      COR=2.0*OMEGACB*SIN(DEG/180.*ACOS(-1.))
CEN
CEN Call SUBROUTINE BCF6 to set the tidal frequency being solved for
CEN (note: UBAR and VBAR are used as dummy variables here) 
CEN
CEN      CALL BCF6(FREQF6,UBAR,VBAR,0)
      FREQF6=0.0
C
C Write settings to echo file 
C
      CALL ECHOSETF6(G,OMEGACB,DEG,COR,FREQF6)
C
      WRITE(2,'(/A)')'INPUT SECTION COMPLETED'
C-----------------------------------------------------------------------
C
      NBC=NR
      NBC2=2*NBC
      eye=cmplx(0.0, 1.0)
      one=cmplx(1.0,0.0)
c      
c Initialize some arrays
c
      do i=1, NBC2NFdim
      se(i)=0.
      do j=1, NBC2NFdim
      Cov(i,j)=0.
      enddo
      enddo
c
      do i=1, NOBS2dim
      do j=i, NBC2NFdim
      als(i,j)=0.0
      enddo
      enddo      
C
C Get the frequencies to be used in inversion
C
      Call GetFreqI(NF, NF1, NF2, Freq, Infer, Amprat, NFDIM)
      write(2,*) ' '
      write(2,*)'END OF BASIC FUNDY5 INPUT ECHO'
      write(2,20)
      write(2,*)'INVERSE SPECTRUM'
      write(2,*)' '
      write(2,*)'Number of constituents         ', NF      
      write(2,*)'Number of primary constituents ',NF1
      write(2,*)'          i, FREQ(i)'
      do i = 1,NF1
	 write(2,*)i,FREQ(i)
	 if(infer(i).ne.i.or.Amprat(i).ne.one) then
	   write(2,*)'error in inference data; quitting'     
	   write(2,*)'I, FREQ(i), Infer(i), Amprat(i)'
	   write(2,*)i,FREQ(i), Infer(i), Amprat(i)
	   close(2)
	   stop
	 endif
      enddo      
c      
      write(2,*)'Number of inferred constituents', NF2
      if(nf2.gt.0) then
      write(2,*)'          I, FREQ(i),           Infer(i), Amprat(i)'
      do i = NF1+1, NF
	 write(2,*)i,FREQ(i), Infer(i), Amprat(i)
	 if(infer(i).lt.1.or.infer(i).gt.nf1) then	 
	   write(2,*)'error in inference data; quitting'
	   close(2)
	   stop
	 endif
      enddo
      endif
c 
C Get the observatons
c     
      write(2,20)       

cDJM

C     Call GETOBS(NOBS,Xobs,Yobs,Tobs,Tobsgmt0,IYEAR,
C    $Ubarobs,Vbarobs,Zmaxobs,Zminobs,Zbotobs,NOBSDIM)

      Call GETOBS(NOBS,Xobs,Yobs,Tobs,Tobsgmt0,IYEAR,
     $Ubarobs,Vbarobs,Zmaxobs,Zminobs,Zbotobs,NOBSDIM,
     $nndim,nedim,in,x,y,nn,ne)

cDJM

      If(NOBS.gt.NOBSdim) then
         WRITE(2,*) Nobs, Nobsdim
         WRITE(2,*)'EXECUTION TERMINATED - NOBS.gt.NOBSdim'
         close(2)
         stop
      endif
      NOBS2=2*NOBS
c
c Check for only Dirichlet BC's; if not, stop.
c
      Call BnodeCheck(NBC, NBCdim, NGEO, NVN, 
     &  KST, NNDIM, NPCN, LST)
C
C Load vertical structure and baroclinic forcing
C
CEN      CALL VERTGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,X,Y,IN,
CEN     &NNV,Z,AK,ENZ,RHOX,RHOY,HRBARELX,HRBARELY,HDOWN)
      CALL VERTGRIDF6(NN,NE,X,Y,HDOWN,IN,NLEV,ZL,NNV,
     &Z,AK,ENZ,RHOX,RHOY,HRBARELX,HRBARELY)
C
C Check for oversize of NNV
C
      IF (NNV.GT.NNVDIM) THEN
      WRITE(2,*)'NNV =', NNV
      WRITE(2,*)'EXECUTION TERMINATED - NNVDIM IS LESS THAN NNV'
      close(2)
      STOP
      ENDIF
      WRITE(2,110)NNV
C
C Initialize atmospheric forcing then load nonzero values
C Truxton1 does not call atmos; requires no wind stress.
c
      DO 10 I=1,NN
         ATMX(I)=cmplx(0.,0.)
         ATMY(I)=cmplx(0.,0.)
 10   CONTINUE
C      CALL ATMOS(ATMX,ATMY,NNDIM)
c
C Compute IQ,JQ,PPX,PPY
C
      CALL BUILDPPXPPY5(NN,NE,X,Y,IN,
     &IQ,JQ,PPX,PPY)
c
c Determine if Fourier Inverse is to be used; if so, 
c   get Nmode, Lspec, and build FP matrix.
c     Nmode is the Fourier half-width
c     Lspec=2*Nmode+1 is the length of the Fourier spectrum
c     Lspec=0 indicates Fourier is inactive
c
      Nmode=0
      Lspec=0
      write(*,*) 'Input Nmode, the number of Fourier modes to be used',
     & ' around the boundary'
      write(*,*) ' (If Nmode=0 this option will be circumvented.)'
      read(*,*) Nmode
      if (2*Nmode.gt.Nbc) then
        write(2,*) '2*Nmode.gt.Nbc; stopping'
        close(2)
        stop
      endif
      if (Nmode.gt.0) then  
        Lspec=min(2*Nmode+1,Nbc)    
        call BuildFP1(FP, Nbc, Nmode,Lspec)
      endif 
c 
C-----------------------------------------------------------------------
c 
c Build the Least Squares System
c Build the Least Squares System
c
C-----------------------------------------------------------------------
c  
c Loop over frequencies 
c Loop over frequencies
c
      write(2,20)
      write(2,*) 'BUILDING THE LEAST SQUARES SYSTEM' 
      write(2,*) ' '     
      do 2000 ifreq = 1, NF    
      write(2,*)'Frequency ',ifreq,' = ',freq(ifreq)
      omega = freq(ifreq)
C      
C Elevaton matrix for this frequency:
C   Elevation5x brings back LU decompositon of QG; no RHS info for TYPE I BCs
C  (Elevation5x zeroes out QG and RHS before assembling) 
C
      CALL ELEVATION5X(QG,NHBW,NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &INGEO,NGEO,NORM,NVN,LRVS,NR,LST,COR,OMEGA,
     &NNV,Z,AK,ENZ,RHOX,RHOY,HRBARELX,HRBARELY,ATMX,ATMY,
     &RHS)
C
C BEGIN basic loop over Type I BC's for unit repsonses
C If Lspec.gt.0 then the loop is over fourier modes of the bc's
c
      limit=nbc
      if(Lspec.gt.0)limit=Lspec
      Do 1000 IBC=1,limit
c      
c Set the RHS's and solve for elevation
c  
      Do I=1,NN
         RHS(I) = cmplx(0.0,0.0)
      Enddo
      if(Lspec.gt.0) then 
        do i=1, nbc
          rhs(lrvs(i))=fp(i,ibc)
        enddo
      else     
        RHS(LRVS(IBC))= cmplx(1.0,0.0)
      endif
      CALL CSOLVET(2,QG,RHS,NN,NHBW,NNDIM,NBWE)
C
C Velocity Solution
C
      DOWVERT = .FALSE.
      CALL VELOCITY5xc(NN,X,Y,HDOWN,SUMX,SUMY,
     &INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,COR,OMEGA,
     &NNV,Z,AK,ENZ,RHOX,RHOY,ATMX,ATMY,RHS,IQ,JQ,PPX,PPY,
     &UVERT,VVERT,WVERT,UBAR,VBAR,DOWVERT)
C
C Call Sampler to project the FEM sol'n onto the observational grid     
C
      Call Sampler(IBC, Xobs, Yobs, Nobs, Zmaxobs, Zminobs,
     & Ubarsamp, Vbarsamp, 
     & X, Y, Z, IN, Ubar, Vbar, UVERT, VVERT, NN, NE, NNV)
C      
C  Load the phased velocity samples into the Least Squares matrix ALS
C
      joffset = (Infer(IFREQ)-1)*2*nbc
      if(Lspec.gt.0) joffset = (Infer(IFREQ)-1)*2*Lspec
      DO I=1, NOBS
         ph=cexp(eye*omega*Tobs(I))
         ALS(2*I-1, joffset + 2*IBC-1)= ALS(2*I-1, joffset + 2*IBC-1) 
     &                        +  real( ph*Ubarsamp(I)*Amprat(ifreq) )
         ALS(2*I-1, joffset + 2*IBC  )= ALS(2*I-1, joffset + 2*IBC  )
     &                        + aimag( ph*Ubarsamp(I)*Amprat(ifreq) )
         ALS(2*I  , joffset + 2*IBC-1)= ALS(2*I  , joffset + 2*IBC-1)
     &                        +  real( ph*Vbarsamp(I)*Amprat(ifreq) )
         ALS(2*I  , joffset + 2*IBC  )= ALS(2*I  , joffset + 2*IBC  )
     &                        + aimag( ph*Vbarsamp(I)*Amprat(ifreq) )
      Enddo
C
C END the loop over the Type I BC's
C END the loop over the Type I BC's
C
 1000 Continue
C
C END loop over frequencies; the Least Squares Matrix ALS is now built
C END loop over frequencies; the Least Squares Matrix ALS is now built
C
 2000 Continue
C
C  Build the Least Squares right-hand side BLS
C
      DO I=1, NOBS
      BLS(2*I-1)=Ubarobs(I)
      BLS(2*I  )=Vbarobs(I)
      Enddo            
c 
C-----------------------------------------------------------------------
c 
c The Least Squares System is Built; Invert It.
c The Least Squares System is Built; Invert It.
c
C-----------------------------------------------------------------------
C      
      write(2,20)
      write(2,*) 'INVERTING THE SYSTEM'
      write(2,*)' '
      nrow=nobs2
      ncol=nbc2*nf1
c
      if(Lspec.gt.0) then
         write(2,*) 'Fourier modes active; Nmode=', Nmode
         Iknow=0
         ncol=Lspec*2*nf1
      else
         write(*,*) 'enter 1 for Generalized Weighted Least Squares,',
     &     ' 0 for regular SVD'
         read(*,*) Iknow
      endif
c
      if(Iknow.eq.1) then
         write(*,*) 'enter Rcond, w0, w1'
         read(*,*) Rcond, w0, w1
         write(2,*) 'Weighted Least Squares: '
         write(2,*) 'Rcond, w0, w1:', Rcond, w0, w1 
         call BuildW(Wls,nobs,nbc,nf1,w0,w1,LRVS,x,y)
         call InvertWLS3(nrow,ncol,als,zls,bls,els,
     &                            wls,se,Cov,rcond,irank,RMSV)            
      else 
         write(*,*) 'enter Rcond'
         read(*,*) Rcond           
         write(2,*) 'Solution by SVD with Rcond=',Rcond     
         call Invert3(nrow,ncol,als,zls,bls,els,se,Cov,
     &                                   rcond,irank,RMSV)       
      endif
c 
c Invert3 and InvertWLS3 return these arrays:
c
c Als, Bls unchanged; 
c ZLS is the solution vector
c    Als*Zls=Bls  
c Els is the least square residual error vector
c    Els=Als*Zls-Bls i.e. (model - observation)
c se is the normalized standard error for Zls.  
c    se is SE/(RMS Velocity error)
c Cov is the covariance matrix, normalized by the mean square velocity error.
c    se(i)=sqrt(cov(i,i))
C
c RMSV is the RMS velocity error.
C
c All arrays are real i.e. they treat x- and y- errors and Re and -Im parts of
c   Zls as separate real variables.
c
C-----------------------------------------------------------------------
c 
c The Least Squares solution is finished; process it     
c The Least Squares solution is finished; process it 
c
C-----------------------------------------------------------------------
c      
c  Get the complex covariance matrix CCov from its real parts Cov
c    
      call CmplxCov(CCov, Cov, ncol)
c
c Recover the inferred constitutents;
c Also fill in se for the inferred constituents
c 23 may: leaving the inferred se feature but note there is no 
c    parallel calculation for Cov for inferred constituents.
c 
      if(nf1.lt.nf)then
        Limit=NBC
        if(lspec.gt.0) Limit=lspec      
        do K=nf1+1, nf
          Koffset =       (K-1)*Limit*2
          KIoffset=(Infer(K)-1)*Limit*2
          fact=cabs(amprat(k))
          do i=1, Limit   
            ZZ1=cmplx(zls(KIoffset+2*i-1), -zls(KIoffset+2*i))
            ZZ2=amprat(k)*zz1
            Zls(Koffset+2*i-1) = real(zz2)
            Zls(Koffset+2*i  ) =-aimag(zz2)
            se(Koffset + 2*i-1)=se(KIoffset +2*i-1)*fact
            se(Koffset + 2*i  )=se(KIoffset +2*i  )*fact
          enddo
        enddo
      endif
c      
c Interpret and write the std errors
c
      write(2,20)
      if(lspec.gt.0)then
      call fourier_se(se,Ccov,nbc,nf,zls,lspec,fp,
     &                    freq,rcond,irank,ncol,RMSV)
      else 
      call svd_se(se,Ccov,nbc,nf,zls,
     &                    freq,rcond,irank,ncol,RMSV)
      endif
C 
C Recover the LS velocities ( = Az) from Az-b=e
C
      Do I=1,NOBS
      UBARls(I) = Bls(2*I-1)+Els(2*I-1)
      VBARls(I) = Bls(2*I  )+Els(2*I)
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 15 dec 98 drl:
c sign of (uerr, verr) reversed so error=observation-model.
c this is the most useful definition for truxout purposes.
c
c      Uerr(I)   = els(2*I-1)
c      Verr(I)   = els(2*I)
      Uerr(I)   = -els(2*I-1)
      Verr(I)   = -els(2*I)
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
      Enddo
c 
c We now have in addition the real arrays 
c    (UBarls, Vbarls) = the modeled velocity at the observation points; 
c    (Uerr,   Verr)   = the modeled error at the observation points;
c And the complex array of Covariances CCov 
c    Like Cov, these are normalized by the mean squared velocity error. 
C 
C-----------------------------------------------------------------------
c
c Solve the forward problems using the optimal BC's
c Solve the forward problems using the optimal BC's
c
C-----------------------------------------------------------------------
c
C Solve the composite system using Zls; call a standard Fundy5 Output
C   subroutine for each frequency.
C
      write(2,20)
      write(2,*) 'RECONSTRUCTION OF DOMAIN SOLUTION WITH INVERSE BCS'
      write(2,*) ' '
C      
      Do 3000 IFREQ = 1, NF
	 write(2,*)'   frequency ',ifreq,' = ',freq(ifreq)
	 omega = freq(ifreq)
C
C Reconstruct QC at this freq.
      CALL ELEVATION5X(QG,NHBW,NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &INGEO,NGEO,NORM,NVN,LRVS,NR,LST,COR,OMEGA,
     &NNV,Z,AK,ENZ,RHOX,RHOY,HRBARELX,HRBARELY,ATMX,ATMY,
     &RHS)
C
C Assign RHS's at this freq.
      Do I=1,NN
          RHS(I) = Cmplx(0.0,0.0)
      Enddo
      If(Lspec.gt.0) then
          joffset = (IFREQ-1)*2*Lspec
          do IBC=1,nbc
          do j=1,Lspec
            rhs(lrvs(ibc)) = rhs(lrvs(ibc)) + 
     &        FP(ibc,j)*Cmplx( Zls(joffset+2*j-1), -Zls(joffset+2*j))
          enddo
          enddo
      else
          joffset = (IFREQ-1)*2*NBC
          Do IBC=1,NBC
            RHS(LRVS(IBC))= 
     $	      Cmplx( Zls(joffset+2*IBC-1), -Zls(joffset+2*IBC))
          enddo
      endif
c
c Store the BC's for use below at Subr. TruxOut
      do ibc=1, nbc
        ZCamp(ibc,ifreq)=rhs(lrvs(ibc))
      enddo            
C
C Solve for elevation and velocity at this freq.
      CALL CSOLVET(2,QG,RHS,NN,NHBW,NNDIM,NBWE)      
      DOWVERT = .TRUE.
      CALL VELOCITY5xc(NN,X,Y,HDOWN,SUMX,SUMY,
     &INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,COR,OMEGA,
     &NNV,Z,AK,ENZ,RHOX,RHOY,ATMX,ATMY,RHS,IQ,JQ,PPX,PPY,
     &UVERT,VVERT,WVERT,UBAR,VBAR,DOWVERT)
C
C Write the results at this freq.  
CEN      CALL OUTPUT5(X,Y,Z,RHS,UBAR,VBAR,UVERT,VVERT,WVERT,OMEGA,
CEN     & NN,NNV,NNDIM,NNVDIM)      
      CALL OUTPUTT6(MESHNAME,NN,NE,X,Y,IN,NNV,Z,
     &RHS,UBAR,VBAR,UVERT,VVERT,WVERT,OMEGA,IYEAR,Tobsgmt0(1),LRVS,NR)
C     
3000  continue
c
C-----------------------------------------------------------------------
C
C Prepare arrays for output and call TruxOut
C Prepare arrays for output and call TruxOut
C
C-----------------------------------------------------------------------
C
c Build (Upartial, Vpartial)_ij = the modeled contrbutions 
c   to observation point i at frequency j
c
      number=nbc
      if (lspec.gt.0) number=lspec
      DO i=1, nobs 
      DO j=1, nf1
        k1=(j-1)*2*number+1
        k2=(j  )*2*number
        sumU=0.
        sumV=0.
        Do k=k1,k2
           sumU=sumU+ALS(2*i-1, k)*zls(k)
           sumV=sumV+ALS(2*i,   k)*zls(k)
        Enddo
        Upartial(i,j)=sumU
        Vpartial(i,j)=sumV
      Enddo
      Enddo
C      
      Call TruxOut6a(NOBS,Xobs,Yobs,Tobs,Tobsgmt0,IYEAR,
     & Ubarobs,Vbarobs,Zmaxobs,Zminobs,Zbotob,
     & Nf, Nf1, Nf2, freq, infer, amprat, 
     & NBC,ZCamp,LRVS,CCov, 
     & UBARls, VBARls, Upartial, Vpartial, 
     & Uerr, Verr, RMSV,rcond,irank,ncol)
C
C-----------------------------------------------------------------------
C
C Write successful completion message to ECHO file; Done
C Write successful completion message to ECHO file; Done
C
C-----------------------------------------------------------------------
C
      write(2,20)
      write(2,*) 'Truxton1 run successfully completed'
      write(2,20)
      CLOSE(2)
C
C Format statements
C   
   20 format(/72('*')/72('*')/)   
 110  FORMAT('Set in Subroutine VERTGRID:',
     &/5X,'The number of nodes in the vertical = ',I6)
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
      STOP
      END
C
C***********************************************************************
C***********************************************************************
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C***********************************************************************
C*****************************************************************FUNDY5
C
C***********************************************************************
C***********************************************************************
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C***********************************************************************
C*****************************************************************FUNDY5
C
C***********************************************************************
C***********************************************************************
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C***********************************************************************
C*****************************************************************FUNDY5
C
C***********************************************************************
C***********************************************************************
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE READINPF5(NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &INBE,NBE,NORM,NVN,KST,LRVS,NR,NPCN,LST,COR,OMEGA)
C-----------------------------------------------------------------------
C
C Purpose:  To read data required to solve the shallow Water wave
C    equations using FUNDY5 from a LABLIB Standard .inp file 
C    (See DRL 11/11/91 Memo for a description of a .inp file)
C
C Inputs upon execution:  .inp filename
C                         echo filename (run statistics are written to
C                                        this file)
C
C Index range: I=1,NN, J=1,NNV, L=1,NE, K=1,NBE  
C
C outputs: NN - number of horizontal nodes
C          NE - the number of horizontal 2-D triangular elements
C          X(I),Y(I) - cartesian nodal coordinates
C          IN(3,L) - LABLIB STANDARD triangular element incidence list
C          HDOWN(I) - Bathymetric depth at nodes
C          SUMX(I) - cosine value for boundary nodes where V.n is 
C                          specified
C          SUMY(I) - sine value for boundary nodes where V.n is 
C                          specified
C          DSBDRY(I) - surface integral information
C          INBE(2,K) - LABLIB STANDARD geostrophic boundary element
C                          incidence list
C          NBE - The number of geostrophic boundary elements
C          NORM(I) - boundary node numbers where normal velocity is set
C                          I =1,NVN       are V.n .NE. 0.0
C                          I =KST,NNDIM   are V.n=0.0
C          LRVS(I) - node numbers of set depth and set velocity nodes
C                          I=1,NR         are set depth nodes
C                          I=NR+1,NR+NPCN are V=0.0
C                          I=LST,NNDIM    are V .NE. 0.0
C          COR - Coriolis Parameter
C          OMEGA - Radian Frequency(rad/sec)
C
C History: Developed by Daniel R. Lynch and Christopher E. Naimie
C          Dartmouth College
C          FUNDY4 - Part of Main Routine written by Daniel R. Lynch
C          FUNDY5 - Converted to a Subroutine by Christopher. E. Naimie
C
C-----------------------------------------------------------------------
C
C INCLUDE FILE FOR PARAMETER STATEMENTS
C
      INCLUDE 'TRUXTON6ARRAYS.DIM' 
C
C GLOBAL ARRAYS
C
      integer INBE(2,NBEDIM),IN(3,NEDIM)
      integer NORM(NNDIM),LRVS(NNDIM)
      real X(NNDIM),Y(NNDIM),HDOWN(NNDIM),SUMX(NNDIM),SUMY(NNDIM),
     &     DSBDRY(NNDIM)
C
C LOCAL ASSIGNMENTS
C
      integer LRT(20)
      real CHEZY(NNDIM),VEL(NNDIM)
      dimension TITLE(20)
C
      CHARACTER *4,TITLE,START,EXIT,AST,
     * INCI,NODA,SET,SPEC,KNOW,
     * AVER,BATH,INIT,HEIG,VELO,
     * TMDT,CHEZ,MANN,CONS,CARD,UNITS,
     * ASTE,SIU,GEOS
      CHARACTER*72 INPFILE,ECHOFILE
C
      DATA START/'WAVE'/,EXIT/'XXXX'/,AST/' *'/,
     * INCI/'INCI'/,NODA/'NODA'/,SET/'SET '/,SPEC/'SPEC'/,KNOW/'KNOW'/,
     * AVER/'AVER'/,BATH/'BATH'/,INIT/'INIT'/,HEIG/'HEIG'/,VELO/'VELO'/,
     * TMDT/'T-DT'/,CHEZ/'CHEZ'/,MANN/'MANN'/,CONS/'CONS'/,
     * SIU/'  SI'/,GEOS/'GEOS'/
C
C START OF EXECUTABLE CODE
C
C Read input and echo file names
C
      WRITE(6,*) ' Enter filename for input data'
      READ(5,1000) INPFILE
      WRITE(6,*) ' Enter filename for input echo'
      READ(5,1000) ECHOFILE
      OPEN (UNIT=1,FILE=INPFILE,STATUS='OLD')
      OPEN (UNIT=2,FILE=ECHOFILE,STATUS='unknown')
      rewind(2)
C
C Write the head of the ECHO file and dimensioning parameters
C
      WRITE(2,3000)
      WRITE(2,3010)NNDIM,NEDIM,NBWE,NBEDIM,NNVDIM,NFTRDIM,RHOREF,G
C
C Read input file down to units specification
C  Halt execution if SI Units are not used
C
   10 READ(1,2000) CARD,UNITS
      IF (CARD.EQ.EXIT) STOP
      IF (CARD.NE.START) GO TO 10
      IF(UNITS.NE.SIU)THEN
      WRITE(2,*)'EXECUTION TERMINATED - THE INPUT DECK MUST BE IN SI UNI
     &TS'
      STOP
      ENDIF
C
C Read title then write title and units to ECHO file
C
      READ (1,2000) TITLE
      WRITE(2,3020)INPFILE,TITLE
C
C Continue reading INPUT file
C
      READ(1,*) NN,NE,NH,NV,NVEL
C
C     WHERE-
C     NN = THE NUMBER OF NODES
C     NE = NUMBER OF ELEMENTS
C     NH = NUMBER OF NODES AT WHICH DEPTH IS SET (ESTIMATED)
C     NVEL = NUMBER OF NODES AT WHICH VELOCITY IS SET (ESTIMATED)
C     NV = NUMBER OF NODES AT WHICH NORMAL VELOCITY IS SET (ESTIMATED)
C
C Halt execution if dimensioning parameters are to small
C
      IF(NNDIM.GE.NN .AND. NEDIM.GE.NE)GO TO 15
      WRITE(2,*)' '
      WRITE(2,*)'NN = ',NN
      WRITE(2,*)'NE = ',NE
      WRITE(2,*)'EXECUTION TERMINATED ... INSUFFICIENT DIMENSIONS'
      IF(NNDIM.LT.NN)WRITE(2,*)'     NNDIM IS LESS THAN NN'
      IF(NEDIM.LT.NE)WRITE(2,*)'     NEDIM IS LESS THAN NE'
      STOP
   15 CONTINUE
C
C     COR  =  degrees latitude of the area under study
C
      READ(1,2001) CORL
      CALL BC(OMEGA,X,Y,NNDIM,0)
      W=OMEGA
      PER=1.0E37
      IF(W.NE.0.0) PER=2*3.14159265/W/3600.
      COR=1.4584E-04*SIN(CORL*3.141593/180.)
C
C FE input and initialization
C
C  Read element incidence list
C    Keyword: INCIDENCE LISTS
C
      READ(1,2000) (TITLE(I),I=1,4)
      WRITE(2,2000) (TITLE(I),I=1,4)
      IF(TITLE(1).EQ.INCI) GO TO 74
      WRITE(2,*)'Error finding beginning of incidence list'
      WRITE(2,2000) (TITLE(I),I=1,4)
   12 WRITE(2,3030) (TITLE(I),I=1,4)
      STOP
   74 CONTINUE
      DO 75 K=1,NE
   75 READ(1,*) J,(IN(I,J),I=1,3)
      READ(1,2003) ASTE
      IF(ASTE.NE.AST)THEN
         WRITE(2,*)'Error finding end of incidence list'
         GO TO 12
      ENDIF
C
C  Compute bandwidth
C
      NB=0
      DO 52 L=1,NE
      MI=NN
      MA=0
      DO 51 I=1,3
      J=IN(I,L)
      IF(J.LT.MI) MI=J
      IF(J.GT.MA) MA=J
   51 CONTINUE
      IDIFF=MA-MI
      IF(IDIFF.GT.NB) NB=IDIFF
   52 CONTINUE
      NHBW=NB
      NBW=2*NB+1
      NDIAG=NB+1
C
C  Halt execution if dimensioning parameter is to small
C
      IF(NBW.GT.NBWE) THEN
      WRITE(2,*)' '
      WRITE(2,*)'NBW = ',NBW
      WRITE(2,*)'EXECUTION TERMINATED ... INSUFFICIENT DIMENSIONS'
      WRITE(2,*)'     NBWE IS LESS THAN NBW'
      STOP
      END IF
C
C  Read nodal coordinates and scale appropriately
C    Keyword: NODAL LOCATIONS
C
      READ(1,2004) (TITLE(I),I=1,4),XFCTR,YFCTR
      IF(TITLE(1).NE.NODA)THEN
         WRITE(2,*)'Error finding beginning of nodal locations'
         GO TO 12
      ENDIF

      IF(XFCTR.EQ.0.) XFCTR=1.0
      IF(YFCTR.EQ.0.) YFCTR=1.0
      READ(1,*)(J,X(J),Y(J),K=1,NN)
      DO 18 J=1,NN
      X(J)=X(J)*XFCTR
   18 Y(J)=Y(J)*YFCTR
      READ (1,2003) ASTE
      IF(ASTE.NE.AST)THEN
         WRITE(2,*)'Error finding end of node locations'
         GO TO 12
      ENDIF
C
C Read boundary node lists
C
C   Boundary nodes with total depth set
C     Keyword: SET DEPTH NODES
C     NR=# of set depth nodes
C     Node numbers are stored in locations 1 thru NR of LRVS
C
      NR=0
      READ(1,2000) (TITLE(I),I=1,4)
      WRITE(2,2000) (TITLE(I),I=1,4)
      IF(TITLE(1).NE.SET )THEN
         WRITE(2,*)'Error finding beginning of set depth nodes'
         GO TO 12
      ENDIF
  180 READ(1,2005) (LRT(ITT),ITT=1,13),ASTE
      IF(ASTE.EQ.AST) GO TO 200
      DO 190 I=1,13
      J=LRT(I)
      IF(J.EQ.0) GO TO 180
      IF(J.GT.0 .AND. J.LE.NN) GO TO 185
C
C   Halt execution due to nonexistance of set depth node
C
      WRITE(2,*)' '
      WRITE(2,*)'EXECUTION TERMINATED - SET DEPTH NODE ',J,
     &' DOES NOT EXIST'
      STOP
  185 NR=NR+1
      LRVS(NR)=J
  190 CONTINUE
      GO TO 180
  200 CONTINUE
C
C   Boundary nodes where normal velocity is specified
C     Keyword: SPECIFIED NORMAL
C     If LRT .GT. 0 normal velocity is set to zero
C     If LRT .LT. 0 normal velocity is non zero
C     NVV is the number of set normal velocity nodes
C     NVN is the number of these nodes with non zero normal velocity
C     Node numbers are stored in NORM(I):
C        I =   1, NVN   are V.N . NE. 0.0
C        I = KST, NNDIM are V.N   =   0.0
C
      READ(1,2000) (TITLE(I),I=1,4)
      WRITE(2,2000) (TITLE(I),I=1,4)
      IF(TITLE(1).NE. SPEC)THEN
         WRITE(2,*)'Error finding beginning of specified normal nodes'
         GO TO 12
      ENDIF
      KST=NNDIM+1
      NVN=0
  171 READ(1,2005)(LRT(ITT),ITT=1,13),ASTE
      IF(ASTE.EQ.AST) GO TO 173
      DO 172 I=1,13
      J=IABS(LRT(I))
      IF(J.EQ.0) GO TO 171
      IF(J.GT.0  .AND.  J.LE.NN)GO TO 161
C
C   Halt execution due to nonexistance of set normal velocity node
C
      WRITE(2,*)' '
      WRITE(2,*)'EXECUTION TERMINATED - KNOWN NORMAL VELOCITY NODE ',J,
     &' DOES NOT EXIST'
      STOP
C
  161 IF(LRT(I).LT.0) GO TO 162
      KST=KST-1
      NORM(KST)=J
      GO TO 172
  162 NVN=NVN+1
      NORM(NVN)=J
  172 CONTINUE
      GO TO 171
  173 CONTINUE
      NVV=NVN+(NNDIM+1-KST)
C
C
C   Physical corner nodes ... both components of velocity specified
C     Keyword: KNOWN VELOCITIES
C     If LRT .GT. 0 the velocity is set equal to zero
C     If LRT .LT. 0 the velocity is non zero
C     NPCN is the number of zero velocity nodes
C     NNZV is the number of non zero velocity nodes
C     Node numbers are stored in LRVS(I):
C        I = NR+1 thru NR+NPCN are V   =  0.0
C        I = LST  thru NNDIM   are V .NE. 0.0
C
      NNZV=0
      LST=NNDIM+1
      NPCN=0
      READ(1,2000) (TITLE(I),I=1,4)
      WRITE(2,2000) (TITLE(I),I=1,4)
      IF(TITLE(1).NE.KNOW)THEN
         WRITE(2,*)'Error finding beginning of known velocities'
         GO TO 12
      ENDIF
  174 READ(1,2005) (LRT(ITT),ITT=1,13),ASTE
      IF(ASTE.EQ.AST) GO TO 177
      DO 175 I=1,13
      J=IABS(LRT(I))
      IF(J.EQ.0) GO TO 174
      IF(J.GT.0 .AND. J.LE.NN) GO TO 176
C
C   Halt execution due to nonexistance of set velocity node
C
      WRITE(2,*)' '
      WRITE(2,*)'EXECUTION TERMINATED - SET VELOCITY NODE ',J,
     &' DOES NOT EXIST'
      STOP
C
  176 IF(LRT(I).LT.0) GO TO 178
      NPCN=NPCN+1
      IA=NPCN+NR
      LRVS(IA) = J
      GO TO 175
  178 LST=LST-1
      LRVS(LST)=J
  175 CONTINUE
      GO TO 174
  177 CONTINUE
      NNZV=NNDIM+1-LST
      III=NPCN+NNZV
      IA=NR+1
      IB=NR+NPCN
C
C   Read and compute normal directions for boundary nodes
C     Keyword: AVERAGE NORMALS
C
      READ(1,2000) (TITLE(I),I=1,4)
      WRITE(2,2000) (TITLE(I),I=1,4)
      IF(TITLE(1).NE.AVER)THEN
         WRITE(2,*)'Error finding beginning of average normals'
         GO TO 12
      ENDIF
      DO 179 I=1,NN
      DSBDRY(I)=0.0
      SUMX(I)=0.0
  179 SUMY(I)=1.0
      KK=NVV+NPCN+NNZV
      IF(KK.LE.0)THEN
         READ(1,2005) (LRT(ITT),ITT=1,13),ASTE
         IF(ASTE.NE.AST)THEN
            WRITE(2,*)'Error finding end of average normals'
            GO TO 12
         ENDIF
         GO TO 183
      ENDIF
      DO 182 I=1,KK
      READ(1,2005) (LRT(ITT),ITT=1,13),ASTE
      IF(ASTE.NE.AST) GO TO 181
C
C   Halt execution due to nonexistance of set velocity node
C
      WRITE(2,*)' '
      WRITE(2,*)'EXECUTION TERMINATED - INSUFFICIENT DATA FOR COMPUTATIO
     &N OF NORMAL DIRECTIONS'
      STOP
C
  181 J=LRT(1)
      J1=LRT(2)
      J2=LRT(3)
      DEX=X(J2)-X(J1)
      DEY=Y(J2)-Y(J1)
      DS=SQRT(DEX*DEX+DEY*DEY)
      DSBDRY(J)=DS
      SUMX(J)=-DEX/DS
      SUMY(J) = DEY/DS
  182 CONTINUE
      READ(1,2003) ASTE
      IF(ASTE.NE.AST)THEN
         WRITE(2,*)'Error finding end of average normals'
         GO TO 12
      ENDIF
  183 CONTINUE
C
C   Read geostrophic boundary elements
C     Keyword: GEOSTROPHIC BE'S 
C     NBE=number of geostrophic boundary elements
C     INBE(2,NBEDIM) is the boundary element incidence list
C     counterclockwise convention
C
      NBE=0
      READ(1,2000) (TITLE(I),I=1,4)
      WRITE(2,2000) (TITLE(I),I=1,4)
      WRITE(2,*)KK,NVV,NPCN,NNZV
      IF(TITLE(1).NE.GEOS)THEN
         WRITE(2,*)'Error finding beginning of geostrophic bes'
         GO TO 12
      ENDIF
  280 READ(1,2005) (LRT(ITT),ITT=1,13),ASTE
      IF(ASTE.EQ.AST) GO TO 281
      NBE=NBE+1
      INBE(1,NBE)=LRT(1)
      INBE(2,NBE)=LRT(2)
      GO TO 280
  281 CONTINUE
C
C   Halt execution if dimensioning parameter is to small
C
      IF(NBE.GT.NBEDIM) THEN
      WRITE(2,*)' '
      WRITE(2,*)'NBE = ',NBE
      WRITE(2,*)'EXECUTION TERMINATED ... INSUFFICIENT DIMENSIONS'
      WRITE(2,*)'     NBEDIM IS LESS THAN NBE'
      STOP
      ENDIF
C
C   Finalize
C
      KNOWNV=NPCN+NNZV
      JSET=NN-KNOWNV
      N=NN-NR
      NZEROV=NNDIM+1-KST
C
C     NN =     number of nodes
C     NR =     number of nodes with depth eliminated from consideration
C     JSET =   number of unknown velocity nodes
C     N =      number of unknown depths
C     KNOWNV = number of known velocity nodes
C     NVV =    number of nodes where normal velocity is set
C     NZEROV = number of nodes where the normal velocity is set to zero
C  
C Read bathymetric depths
C     Keyword: BATHYMETRIES
C
      READ(1,2004) (TITLE(ITT),ITT=1,4),FCTR,HMIN
      IF(TITLE(1).NE.BATH)THEN
         WRITE(2,*)'Error finding beginning of bathymetries'
         GO TO 12
      ENDIF
      IF(FCTR.EQ.0.0) FCTR=1.0
      DO 401 I=1,NN
      READ(1,2006) J,VAL,ASTE
      IF(ASTE.EQ.AST)THEN
         WRITE(2,*)ASTE
         GO TO 12
      ENDIF
      HDOWN(J)=VAL*FCTR
      IF(HDOWN(J).LT.HMIN) HDOWN(J)=HMIN
  401 CONTINUE
      READ(1,2003) ASTE
      IF(ASTE.NE.AST)THEN
         WRITE(2,*)'Error finding end of bathymetries'
         GO TO 12
      ENDIF
C
C The rest of the INPUT file is ignored by FUNDY5
C
C  INITIAL TOTAL DEPTH VALUES
C       (NOT USED)
C     KEYWORD: INITIAL HEIGHTS
      READ(1,2004) (TITLE(ITT),ITT=1,4),FCTR1
      IF(TITLE(1).NE.INIT)THEN
         WRITE(2,*)TITLE(1)
         GO TO 12
      ENDIF
  402 READ(1,2006) J,VAL,ASTE
      IF(ASTE.EQ.AST) GO TO 403
      GO TO 402
  403 CONTINUE
C
C  TOTAL DEPTH AT TIME T-DT
C       (NOT USED)
C     KEYWORD: HEIGHTS AT T-DT
C
      READ(1,2004) (TITLE(ITT),ITT=1,4),FCTR1
      IF(TITLE(1).NE.HEIG)THEN
         WRITE(2,*)TITLE(1)
         GO TO 12
      ENDIF
  404 READ(1,2006) J,VAL,ASTE
      IF(ASTE.EQ.AST) GO TO 405
      GO TO 404
  405 CONTINUE
C
C  INITIAL VELOCITIES
C     USED TO COMPUTE LINEARIZED TAU
C     KEYWORD: VELOCITIES
C
      READ(1,2004) (TITLE(I),I=1,4),XFCTR1,YFCTR1
      IF(TITLE(1).NE.VELO)THEN
         WRITE(2,*)TITLE(1)
         GO TO 12
      ENDIF
      IF(XFCTR1.EQ.0.0) XFCTR1=1.0
      IF(YFCTR1.EQ.0.0) YFCTR1=1.0
      DO 408 I=1,NN
  408 VEL(I)=0.0
  409 READ(1,2007) J,XVAL,YVAL,ASTE
      IF(ASTE.EQ.AST) GO TO 410
      VEL(J)=((XFCTR1*XVAL)**2+(YFCTR1*YVAL)**2)**.5
      GO TO 409
  410 CONTINUE
C
C  VELOCITIES AT TIME T-DT
C       (NOT USED)
C     KEYWORD: T-DT VELOCITIES
C
      READ(1,2004) (TITLE(I),I=1,4),XFCTR1,YFCTR1
      IF(TITLE(1).NE.TMDT)THEN
         WRITE(2,*)TITLE(1)
         GO TO 12
      ENDIF
  413 READ (1,2007) J,XVAL,YVAL,ASTE
      IF(ASTE.EQ.AST) GO TO 414
      GO TO 413
  414 CONTINUE
C
C  FRICTION COEFFICIENTS
C     KEYWORDS: MANNINGS N / CHEZY PARAMETERS / CONSTANT TAU
C     (NOT USED)
C
      READ(1,2004) (TITLE(I),I=1,4),FCTR1
      IFRIC=-2
      IF(TITLE(1).EQ.CONS) IFRIC=-1
      IF(TITLE(1).EQ.MANN) IFRIC=0
      IF(TITLE(1).EQ.CHEZ) IFRIC=1
      IF(IFRIC.EQ.-2)THEN
         WRITE(2,*)TITLE(1)
         GO TO 12
      ENDIF
      DO 415 I=1,NN
  415 CHEZY(I)=FCTR1
  416 READ(1,2006) J,VAL,ASTE
      IF(ASTE.EQ.AST) GO TO 417
      CHEZY(J)=VAL
      GO TO 416
  417 CONTINUE
C
C  Done Reading close input file
      close(1)

C Format Statements for reading USER INPUT
C
 1000 FORMAT(A)
C
C Format Statements for reading the INPUT file
C
 2000 FORMAT (20A4)
 2001 FORMAT (40X,F10.0/10X)
 2003 FORMAT(78X,A2)
 2004 FORMAT(4A4,4X,2G10.0)
 2005 FORMAT(13I6,A2)
 2006 FORMAT(I5,G10.0,63X,A2)
 2007 FORMAT(I5,2G10.0,53X,A2)
C
C The INPUT file has been properly read - write data to ECHO file
C
      WRITE(2,3040)PER,CORL
      WRITE(2,3050)NN,NE,NBW
      WRITE(2,3060)XFCTR,YFCTR,FCTR,HMIN
      WRITE(2,3065)NR,NVN,NNDIM+1-LST,NNDIM+1-KST,NPCN,NBE
      WRITE(2,3070)
      IF(NR.NE.NH) WRITE(2,3080) NR,NH
      IF(NR.GT.0) WRITE(2,3090) (LRVS(I),I=1,NR)
      WRITE(2,3100)
      IF(NVV.NE.NV) WRITE(2,3110) NVV,NV
      IF(NVN.GT.0) THEN
        WRITE(2,3120)
        WRITE(2,3090) (NORM(I),I=1,NVN)
      ENDIF
      IF(KST.LE.NNDIM) THEN
        WRITE(2,3130)
        WRITE(2,3090) (NORM(I),I=KST,NNDIM)
      ENDIF
      WRITE(2,3140)
      IF(III.NE.NVEL) WRITE(2,3150) III,NVEL
      IF(NPCN.GT.0) THEN
        WRITE(2,3130) 
        WRITE(2,3090) (LRVS(I),I=IA,IB)
      ENDIF
      IF(LST.LE.NNDIM) THEN
        WRITE(2,3120) 
        WRITE(2,3090) (LRVS(I),I=LST,NNDIM)
      ENDIF
      WRITE(2,3160)
      DO I=1,NBE
         WRITE(2,3170) INBE(1,I),INBE(2,I)
      ENDDO
C
C Format statements for writing ECHO file
C
 3000 format(72(1H*)/72(1H*)/2(1H*),29X,' TRUXTON1 ',29X,2(1H*)/
     &2(1H*),68X,2(1H*)/2(1H*),17X,'MULTICONSTITUENT INVERSE OF FUNDY5',
     &17X,2(1H*)/
     &2(1H*),68X,2(1H*)/
     &2(1H*),20X,'LINEARIZED HARMONIC EQUATIONS',19X,2(1H*)/
     &2(1H*),68X,2(1H*)/
     &2(1H*),19X,'WITH LINEAR TRIANGULAR ELEMENTS',18X,2(1H*)/
     &2(1H*),68X,2(1H*)/2(1H*),27X,'June 10, 1997',28X,2(1H*)/
     &72(1H*)/72(1H*))
 3010 FORMAT(//'Parameters specified in ', 
     & 'include file TRUXTON1ARRAYS.DIM:',
     &//5X,'NNDIM   ',I10,/5X,'NEDIM   ',I10,/5X,'NBWE    ',I10,
     &/5X,'NBEDIM  ',I10,/5X,'NNVDIM  ',I10,/5X,'NFTRDIM ',I10,
     &/5X,'RHOREF  ',F10.2,/5X,'G       ',F10.4)
 3020 FORMAT (///,'Input Deck:  ',A,
     &          /,'     Title:  ',20A4,//5X,'UNITS:  SI')
 3030 FORMAT (1X,10(1H*),1X,'ERROR IN DATA INPUT',1X,4A4)
 3040 FORMAT(5X,'PERIOD IN HOURS ',E10.3,/5X,'DEGREES LATITUDE',F10.2)
 3050 FORMAT(//5X,'Finite Element Data:',
     &/10X,'Total number of nodes         ',I6
     &/10X,'Total number of elements      ',I6
     &/10X,'Maximum bandwidth required    ',I6)
 3060 FORMAT(//5X,'Scale Factors:',
     &/10X,'X scaling factor      ',E11.4
     &/10X,'Y scaling factor      ',E11.4
     &/10X,'Depth scaling factor  ',E11.4
     &/10X,'Minimum depth set at  ',E11.4)
 3065 FORMAT(//5X,'Boundary condition totals:'
     &/8X,'Node totals for boundary conditions set in Subroutine BC:'
     &/10X,'Set depth nodes               ',I6
     &/10X,'Set non zero normal velocity  ',I6
     &/10X,'Set non zero velocity         ',I6
     &/8X,'Totals for boundary conditions computed by FUNDY5:'
     &/10X,'Zero normal velocity nodes    ',I6
     &/10X,'Zero velocity nodes           ',I6
     &/10x,'Geostrophic boundary elements ',i6)
 3070 FORMAT(///20X,'NODES AT WHICH TOTAL DEPTH IS SET',/20X,33(1H-))
 3080 FORMAT(20X,'NUMBER OF CONSTANT DEPTH NODES READ',I6,
     &      /20X,'DISAGREES WITH NUMBER ANTICIPATED  ',I6)
 3090 FORMAT (12I6)
 3100 FORMAT(///17X,'BOUNDARY NODES WITH NORMAL VELOCITY SET',
     &/17X,39(1H-))
 3110 FORMAT(17X,'NUMBER OF SET NORMAL VELOCITY NODES READ',I6,
     &      /17X,'DISAGREES WITH NUMBER ANTICIPATED       ',I6)
 3120 FORMAT (33X,'NONZERO')
 3130 FORMAT (34X,'ZERO')
 3140 FORMAT(///16X,'PHYSICAL CORNER NODES ... VELOCITY IS SET',
     &/16X,41(1H-))
 3150 FORMAT(16X,'NUMBER OF SET VELOCITY NODES READ',I6,
     &      /16X,'DISAGREES WITH NUMBER ANTICIPATED',I6)
 3160 FORMAT(///22X,'GEOSTROPHIC BOUNDARY ELEMENTS'/22X,29(1H-))
 3170 FORMAT (30X,2I6)
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE ELEVATION5x(QG,NHBW,NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,
     & DSBDRY,
     & INBE,NBE,NORM,NVN,LRVS,NR,LST,COR,OMEGA,
     & NNV,Z,AK,ENZ,RHOX,RHOY,HRBARELX,HRBARELY,ATMX,ATMY,
     & RHS)
C-----------------------------------------------------------------------
C
C Purpose: To solve Equation 37 in the Users Manual, which is a 2-D 
C            Helmholtz equation, for the free surface elevation.
C
C Method:  Finite Element Method on 2-D, triangular, C0 (linear)
C            elements, using nodal quadrature. 
C            (See Equations 38-42 of the Users manual)
C
C Index range: I=1,NN, J=1,NNV, L=1,NE, K=1,NBE  
C
C inputs:  NN - number of horizontal nodes
C          NE - the number of horizontal 2-D triangular elements
C          X(I),Y(I) - cartesian nodal coordinates
C          IN(3,L) - LABLIB STANDARD triangular element incidence list
C          HDOWN(I) - Bathymetric depth at nodes
C          SUMX(I) - cosine value for boundary nodes where V.n is 
C                          specified
C          SUMY(I) - sine value for boundary nodes where V.n is 
C                          specified
C          DSBDRY(I) - surface integral information
C          INBE(2,K) - LABLIB STANDARD geostrophic boundary element
C                          incidence list
C          NBE - The number of geostrophic boundary elements
C          NORM(I) - boundary node numbers where normal velocity is set
C                          I =1,NVN       are V.n .NE. 0.0
C                          I =KST,NNDIM   are V.n=0.0
C          LRVS(I) - node numbers of set depth and set velocity nodes
C                          I=1,NR         are set depth nodes
C                          I=NR+1,NR+NPCN are V=0.0
C                          I=LST,NNDIM    are V .NE. 0.0
C          COR - Coriolis Parameter
C          OMEGA - Radian Frequency  (rad/sec)
C          NNV - number of vertical nodes
C          AK(I) - Linear bottom stress coefficient
C          ENZ(I,J) - Vertical viscosity
C          Z(I,J) - Vertical coordinate of 3-D nodes
C          RHOX(I,J) - x-component of baroclinic forcing
C          RHOY(I,J) - y-component of baroclinic forcing
C          HRBARELX(L) - HDOWN(L)*(Vertical integral of RHOX(I,L) 
C                          over entire water column)
C          HRBARELX(L) - HDOWN(L)*(Vertical integral of RHOY(I,L)
C                          over entire water column)
C          ATMX(I) - x-component of atmospheric forcing (H(L)*PSIX(L))
C          ATMY(I) - y-component of atmospheric forcing (H(L)*PSIY(L))
C
C Outputs: RHS(I) - the free surface elevation at nodes
C
C History: Developed by Daniel R. Lynch and Christopher E. Naimie
C          Dartmouth College
C          FUNDY4 - Part of Main Routine written by Daniel R. Lynch
C          FUNDY5 - Converted to a Subroutine by Christopher. E. Naimie
C
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Global Arrays
C
      INTEGER IN(3,NEDIM),INBE(2,NBEDIM) 
      INTEGER LRVS(NNDIM),NORM(NNDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),AK(NNDIM),SUMX(NNDIM), 
     &     SUMY(NNDIM),DSBDRY(NNDIM)
      REAL Z(NNDIM,NNVDIM),ENZ(NNDIM,NNVDIM)
      COMPLEX ATMX(NNDIM),ATMY(NNDIM),RHS(NNDIM)
      COMPLEX RHOX(NNDIM,NNVDIM),RHOY(NNDIM,NNVDIM)
      COMPLEX HRBARELX(NEDIM),HRBARELY(NEDIM) 
C f5x: QG now global     
      COMPLEX QG(NBWE,NNDIM)
C
C Local Arrays
C
      REAL DELX(3),DELY(3)
C      COMPLEX QG(NBWE,NNDIM)
      COMPLEX TAUPRIME(NNDIM),FPRIME(NNDIM),ATMPX(NNDIM),ATMPY(NNDIM),
     &        BETAX(NNDIM),BETAY(NNDIM),RBARX(NNDIM),RBARY(NNDIM)
      COMPLEX P1P(NNVDIM),P1M(NNVDIM),
     &        P2P(NNVDIM),P2M(NNVDIM),
     &        P3P(NNVDIM),P3M(NNVDIM),
     &        RP(NNVDIM),RM(NNVDIM)
C
C Constant declarations
C
      COMPLEX A,EYE,DENOM,SUMT,SUMF,DHDX,DHDY,UU,VV,WW,
     * SUMAX,SUMAY,VNORM,F1,F2,HRBARX,HRBARY
      COMPLEX TAUP,TAUM,ALPHAP,ALPHAM,BETAP,BETAM
      COMPLEX GP,GM,PSIP,PSIM,UVP,UVM,RBP,RBM,DP,DM
C
C Constants
C
      EYE=cmplx(0.0,1.0)
      OMEGAP=OMEGA+COR
      OMEGAM=OMEGA-COR
C
C Vertical structure
C
C  Compute vertical average of nodal based baroclinic pressure gradient
C
      CALL VERTAVG(RHOX,RBARX,Z,NN,NNV,NNDIM,NNVDIM)
      CALL VERTAVG(RHOY,RBARY,Z,NN,NNV,NNDIM,NNVDIM)
C
C  Nodal loop for vertical structure
C
      DO 20 I=1,NN
         DO 10 J=1,NNV
            RP(J)=.5*(RHOX(I,J)+EYE*RHOY(I,J))
            RM(J)=.5*(RHOX(I,J)-EYE*RHOY(I,J))
 10      CONTINUE
C
C  Call VERTSTRUC5 to compute vertical structure of P1+-,P2+-,P3+-
C    (Equations 10-12,16-17 in Users Manual)
C
         CALL VERTSTRUC5(NN,NNV,I,OMEGAP,Z,ENZ,AK,RP,
     &                   P1P,P2P,P3P)
         CALL VERTSTRUC5(NN,NNV,I,OMEGAM,Z,ENZ,AK,RM,
     &                   P1M,P2M,P3M)
C
C  Compute the +,- bottom stress quantities
C    (Equations 26-28 in Users Manual with EQU's 21-22 substituted in)
C
         AKH=AK(I)/HDOWN(I)
         RBP=0.5*(RBARX(I)+EYE*RBARY(I))
         RBM=0.5*(RBARX(I)-EYE*RBARY(I))
         DP=1.0/(1.-AKH*P1P(1))
         DM=1.0/(1.-AKH*P1M(1))
C
         TAUP=EYE*OMEGAP*AKH*P1P(1)*DP
         TAUM=EYE*OMEGAM*AKH*P1M(1)*DM
         ALPHAP=AKH*P1P(1)*(1.-AKH*P2P(1))*DP - AKH*P2P(1)
         ALPHAM=AKH*P1M(1)*(1.-AKH*P2M(1))*DM - AKH*P2M(1)
         BETAP=AKH*P1P(1)*(RBP-AKH*P3P(1))*DP - AKH*P3P(1)
         BETAM=AKH*P1M(1)*(RBM-AKH*P3M(1))*DM - AKH*P3M(1)
C
C  Compute the prime quantities for this node
C    (Equations 31-33 in Users Manual)
C
         FPRIME(I)=COR-EYE*0.5*(TAUP-TAUM)
         TAUPRIME(I)=      0.5*(TAUP+TAUM)
         ATMPX(I)=ATMX(I)+ATMX(I)*0.5*(ALPHAP+ALPHAM)
     &                   +EYE*ATMY(I)*0.5*(ALPHAP-ALPHAM)
         ATMPY(I)=ATMY(I)+ATMY(I)*0.5*(ALPHAP+ALPHAM)
     &                   -EYE*ATMX(I)*0.5*(ALPHAP-ALPHAM)
C
C   BETA is computed alone here, since it is a nodal quantity.  RBAR is
C     specified as an element based quantity for the solution of the
C     Hemholtz equation for elevation. (See Equation 36 in Users Manual)
C   BETA is HDOWN*BETA for compatibility below with ATM, which is
C     defined as HDOWN*PSI 
C
         BETAX(I)=(BETAP+BETAM)*HDOWN(I)
         BETAY(I)=-EYE*(BETAP-BETAM)*HDOWN(I)
 20   CONTINUE
C
C Elevation Solution (Equation 38-39 in Users Manual)
C  The matrix equation to be built is QG*ZETA=RHS
C
C  Initialize QG and RHS
C
      DO 30 i=1,NN
         RHS(I)=cmplx(0.0,0.0)
         DO 30 J=1,NBWE
           QG(J,I)=cmplx(0.0,0.0)
 30   CONTINUE	    
C
C Compute NHBW,NBW,NDIAG
C
      NHBW=0
      do 40 L=1,NE
         L12=ABS(IN(1,L)-IN(2,L))
         L13=ABS(IN(1,L)-IN(3,L))
         L23=ABS(IN(2,L)-IN(3,L))
         NHBW=MAX0(NHBW,L12,L13,L23)
 40   CONTINUE
      NBW=2*NHBW+1
      NDIAG=NHBW+1
C
C Begin element loop to compute inner products
C
      DO 60 L=1,NE
         I1=IN(1,L)
         I2=IN(2,L)
         I3=IN(3,L)
         DELX(1)=X(I2)-X(I3)
         DELX(2)=X(I3)-X(I1)
         DELX(3)=X(I1)-X(I2)
         DELY(1)=Y(I2)-Y(I3)
         DELY(2) =Y(I3)-Y(I1)
         DELY(3)=Y(I1)-Y(I2)
         AR=0.5*(X(I1)*DELY(1) + X(I2)*DELY(2) + X(I3)*DELY(3))
         AR4=1.0/(4.0*AR)
         SUMT=cmplx(0.,0.)
         SUMF=cmplx(0.,0.)
         SUMAX=cmplx(0.,0.)
         SUMAY=cmplx(0.,0.)
C
C  Basic 2-D nodal quadrature assembly
C  The 3-D adjustments are:
C   substitute FPRIME(II) for COR
C   substitute TAUPRIME(II) for FR(II)
C   substitute (BETA+ATMP+HRBAREL) for ATM
C     This final adjustment is more subtle since BETA and ATMP are nodal
C      based while HRBAREL is element based
C
         DO 50 I=1,3
            II=IN(I,L)
            DENOM=1.0/((EYE*OMEGA + TAUPRIME(II))**2 + FPRIME(II)**2)
            SUMT=SUMT+(EYE*OMEGA + TAUPRIME(II))*G*HDOWN(II)*DENOM
            SUMF=SUMF+ FPRIME(II)*G*HDOWN(II)*DENOM
            SUMAX=SUMAX+((EYE*OMEGA+TAUPRIME(II))*(BETAX(II)+ATMPX(II)
     &           +HRBARELX(L))
     &           +FPRIME(II)*(BETAY(II)+ATMPY(II)+HRBARELY(L)))*DENOM
            SUMAY=SUMAY+((EYE*OMEGA+TAUPRIME(II))*(BETAY(II)+ATMPY(II)
     &           +HRBARELY(L))
     &           -FPRIME(II)*(BETAX(II)+ATMPX(II)+HRBARELX(L)))*DENOM
 50      CONTINUE
         DO 60 I=1,3
            II=IN(I,L)
            RHS(II)=RHS(II) + (SUMAX*DELY(I) -SUMAY*DELX(I) )/6.
            DO 60 J=1,3
               IJ = IN(J,L)
               A=0.
               IF(I.NE.J) GO TO 1
               A=EYE*OMEGA*AR/3.
 1             A=A+SUMT/3.*(DELX(I)*DELX(J) + DELY(I)*DELY(J))*AR4
               A=A-SUMF/3.*(-DELY(J)*DELX(I) + DELX(J)*DELY(I))*AR4
C
C  Global assembly
C
               ICOL=NDIAG+IJ-II
               QG(ICOL,II)=QG(ICOL,II)+A
   60 CONTINUE
C
C End of element loop
C
C Apply Boundary Conditions
C
C  Adjust for geostrophic boundary elements
C   Recall that BETAP and ATM are already multiplied by HDOWN
C
      IF(NBE.GT.0) THEN
         DO 70 K=1,NBE
            N1=INBE(1,K)
            N2=INBE(2,K)
            DX=X(N2)-X(N1)
            DY=Y(N2)-Y(N1)
            F1=1.0/FPRIME(N1)
            F2=1.0/FPRIME(N2)
C   Find domain element on which this boundary element resides
            DO 68 J=1,NE
               DO 68 JJ=1,3
                  IF(IN(JJ,J).EQ.N1)THEN
                     JJJ=JJ+1
                     IF(JJ.EQ.3)JJJ=1
                     IF(IN(JJJ,J).EQ.N2)THEN
                        NELEM=J
                        GO TO 69
                     ENDIF
                  ENDIF
 68         CONTINUE
C   Set elemental HRBAR
 69         HRBARX=HRBARELX(NELEM)
            HRBARY=HRBARELY(NELEM)
C   Make adjustments
            A=G*(HDOWN(N1)*F1/3.0+HDOWN(N2)*F2/6.0)
            QG(NDIAG,N1)=QG(NDIAG,N1)+A
            QG(NDIAG+N2-N1,N1)=QG(NDIAG+N2-N1,N1)-A
            A=G*(HDOWN(N1)*F1/6.0+HDOWN(N2)*F2/3.0)
            QG(NDIAG+N1-N2,N2)=QG(NDIAG+N1-N2,N2)+A
            QG(NDIAG,N2)=QG(NDIAG,N2)-A
            A=+(F1*(HRBARX+BETAX(N1)+ATMPX(N1))*DX
     &         +F1*(HRBARY+BETAY(N1)+ATMPY(N1))*DY)/3.
     &        +(F2*(HRBARX+BETAX(N2)+ATMPX(N2))*DX
     &         +F2*(HRBARY+BETAY(N2)+ATMPY(N2))*DY)/6.
            RHS(N1)=RHS(N1)-A
            A=+(F1*(HRBARX+BETAX(N1)+ATMPX(N1))*DX
     &         +F1*(HRBARY+BETAY(N1)+ATMPY(N1))*DY)/6.
     &        +(F2*(HRBARX+BETAX(N2)+ATMPX(N2))*DX
     &         +F2*(HRBARY+BETAY(N2)+ATMPY(N2))*DY)/3.
            RHS(N2)=RHS(N2)-A
 70      CONTINUE      
      END IF
C
C  Initialize nodal arrays used by Subroutine BC to set BC's
C   Note:  The local arrays RBARX, and RBARY are initialized and
C            reused in this section, to reduce memory requirements. 
C
      DO 80 I=1,NN
         RBARX(I)=cmplx(0.0,0.0)
         RBARY(I)=cmplx(0.0,0.0)
 80   CONTINUE
      NVIOL=0
      DEF=-2097152.0
C
C  Write line to ECHO file
C
      WRITE(2,86)
 86   FORMAT(//,'Elevation solution boundary condition diagnostics:')
C
C  Adjust for nonzero flux nodes
C
      IF (NVN.LT.1) GO TO 91
C   
      DO 85 I=1,NVN
         N=NORM(I)
         RBARX(N)=cmplx(0.0,DEF)
 85   CONTINUE
C   
      CALL BC(W,RBARX,BETAY,NNDIM,2)
C
      WRITE(2,*)'    Nodes where a required non zero normal velocity BC 
     &was not set:'
      DO 89 I=1,NVN
         N=NORM(I)
         IF(AIMAG(RBARX(N)).NE.DEF)THEN
            RHS(N)=RHS(N)-DSBDRY(N)*.5*HDOWN(N)*RBARX(N)
         ELSE
            WRITE(2,*)N
            NVIOL=NVIOL+1
         ENDIF
         RBARX(N)=cmplx(0.0,0.0)
 89   CONTINUE
C
      WRITE(2,*)'    Nodes where a non zero normal velocity BC was incor
     &rectly set:'
      DO 90 I=1,NN
         IF(REAL(RBARX(I)).EQ.0.0.AND.AIMAG(RBARX(I)).EQ.0.0)GO TO 90
         WRITE(2,*)I
         NVIOL=NVIOL+1
         RBARX(I)=cmplx(0.0,0.0)
 90   CONTINUE
C
C  Adjust for nonzero velocity nodes (corners)
C
 91   IF(NNDIM.LT.LST) GO TO 101
      WRITE(2,*)'    Nodes where a required non zero velocity BC was not
     & set:'
C   
      DO 95 I=LST,NNDIM
         N=LRVS(I)
         RBARX(N)=cmplx(0.0,DEF)
         RBARY(N)=cmplx(0.0,DEF)
 95   CONTINUE
C   
      CALL BC(W,RBARX,RBARY,NNDIM,3)
      DO 99 I=LST,NNDIM
         N=LRVS(I)
         SI=SUMX(N)
         CO=SUMY(N)
         IF(AIMAG(RBARX(N)).NE.DEF.AND.AIMAG(RBARY(N)).NE.DEF)THEN
            VNORM=RBARX(N)*CO + RBARY(N)*SI
            RHS(N)=RHS(N) -DSBDRY(N)*.5*HDOWN(N)*VNORM
         ELSE
            WRITE(2,*)N
            NVIOL=NVIOL+1
         ENDIF
         RBARX(N)=cmplx(0.0,0.0)
         RBARY(N)=cmplx(0.0,0.0)
  99  CONTINUE
C
      WRITE(2,*)'    Nodes where a non zero velocity BC was incorrectly 
     &set:'
      DO 100 I=1,NN
         IF(REAL(RBARX(I)).EQ.0.0.AND.AIMAG(RBARX(I)).EQ.0.0.AND.
     &      REAL(RBARY(I)).EQ.0.0.AND.AIMAG(RBARY(I)).EQ.0.0)GO TO 100
         WRITE(2,*)I
         NVIOL=NVIOL+1
         RBARX(I)=cmplx(0.0,0.0)
         RBARY(I)=cmplx(0.0,0.0)
 100  CONTINUE
C
C  Adjust for specified elevation nodes
C
C f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5xC
C f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5xC
C
C f5x: Here we adjust the matrix for type I BC's and call CSOLVEt(1,...).
C      We do not assign a right hand side or solve.  
C      The object is to return the LU decomposition of QG to Main.
C 
  101  Continue
C      IF(NR.EQ.0) GO TO 111
C   
C      DO 105 I=1,NN
C         RBARX(I)=cmplx(0.0,DEF)
C 105  CONTINUE
C   
C      CALL BC(W,RBARX,BETAY,NNDIM,1)
C
      DO 109 I=1,NR
         N=LRVS(I)
         DO 108 M=1,NBW
            QG(M,N)=cmplx(0.0,0.0)
 108     CONTINUE
         QG(NDIAG,N)=cmplx(1.0,0.0)
C   Apply BC from Subroutine BC if a value is given for this node
C         IF(AIMAG(RBARX(N)).NE.DEF)THEN
C            RHS(N)=RBARX(N)
C   Assign elevation=0 if Subroutine BC contained no value for this node
C         ELSE
C            RHS(N)=cmplx(0.0,0.0)
C         ENDIF
C         RBARX(N)=cmplx(0.0,DEF)
 109  CONTINUE
C
C      WRITE(2,*)'    Nodes where an elevation BC was incorrectly set:'
C      DO 110 I=1,NN
C         IF(REAL(RBARX(I)).EQ.0.0.AND.AIMAG(RBARX(I)).EQ.DEF)GO TO 110
C         WRITE(2,*)I
C         NVIOL=NVIOL+1
C 110  CONTINUE
C
C Write statement regarding total BC violations to ECHO file
C
 111  WRITE(2,*)' '
      WRITE(2,*)'    Total number of boundary condition violations:  '
     &,NVIOL
      IF(NVIOL.NE.0)THEN
      WRITE(2,*)' '
      WRITE(2,*)'EXECUTION TERMINATED DUE TO BOUNDARY CONDITION VIOLATIO
     &NS'
      STOP
      ENDIF
C
C Solve QG*ZETA=RHS for surface elevation, ZETA comes back as RHS
C
      CALL CSOLVET(1,QG,RHS,NN,NHBW,NNDIM,NBWE)
C      CALL CSOLVET(2,QG,RHS,NN,NHBW,NNDIM,NBWE)
C
C f5x: QG is now replaced with its LU decomposition; we are done
C
C f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5xC
C f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5xC
C
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE BUILDPPXPPY5(NN,NE,X,Y,IN,
     &IQ,JQ,PPX,PPY)
C-----------------------------------------------------------------------
C
C Purpose:  To compute SPRSPAK STANDARD Normalized Derivative matircies
C              for use in the elevation solution.
C
C Index range: I=1,NN, J=1,NNV, L=1,NE, K=1,NBE  
C
C inputs:  NN - number of horizontal nodes
C          NE - the number of horizontal 2-D triangular elements
C          X(I),Y(I) - cartesian nodal coordinates
C          IN(3,L) - LABLIB STANDARD triangular element incidence list
C
C outputs: IQ(I),JQ(I*NFTR) - SPRSPAK STANDARD pointers for PPX and PPY
C          PPX(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix
C          PPY(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix 
C
C History: Developed by Daniel R. Lynch and Christopher E. Naimie
C          Dartmouth College
C          FUNDY4 - Part of Main Routine written by Daniel R. Lynch
C          FUNDY5 - Converted to a Subroutine by Christopher. E. Naimie
C
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      INCLUDE 'TRUXTON6ARRAYS.DIM' 
C
C Global Arrays 
C
C   Sent to routine and returned unaltered
C
      INTEGER IN(3,NEDIM)
      REAL X(NNDIM),Y(NNDIM)
C
C   Computed in routine and returned to calling routine
C
      INTEGER IQ(NNDIM),JQ(NFTRDIM*NNDIM)
      REAL PPX(NFTRDIM*NNDIM),PPY(NFTRDIM*NNDIM)
C
C Local arrays
C
      REAL DELX(3),DELY(3)
      REAL SV(NNDIM)
C
C Prep for compact storage
C
      CALL SPRSBLD1(IN,3,3,IQ,JQ,NN,NE,NFTRDIM,1)
      NFTR=IQ(NN)/NN +1
      WRITE(2,110)NFTR
      IF(NFTRDIM.LT.NFTR)THEN
      WRITE(2,*)'EXECUTION TERMINATED - NFTRDIM IS LESS THAN NFTR'
      STOP
      END IF
C
C Initialize the sprspak arrays PPX, PPY
C
      DO 10 I=1,NFTRDIM*NNDIM
         PPX(I)=0.
         PPY(I)=0.
 10   CONTINUE
      DO 20 I=1,NN
         SV(I)=0.0
 20   CONTINUE
C
C Compute diagonal mass matrix SV(I)=<Pi> and derivative matricies in 
C   sprspak form
C
C     PPX = <Pi*dP/dxj>
C     PPY = <Pi*dP/dyj>
C
      DO 30 L=1,NE
         I1=IN(1,L)
         I2=IN(2,L)
         I3=IN(3,L)
         DELX(1)=X(I2)-X(I3)
         DELX(2)=X(I3)-X(I1)
         DELX(3)=X(I1)-X(I2)
         DELY(1)=Y(I2)-Y(I3)
         DELY(2) =Y(I3)-Y(I1)
         DELY(3)=Y(I1)-Y(I2)
         ARR=0.5*(X(I1)*DELY(1) + X(I2)*DELY(2) + X(I3)*DELY(3))
         DO 30 I=1,3
            II=IN(I,L)
            SV(II)=SV(II)+ARR/3.
            DO 30 J=1,3
               IJ=IN(J,L)
               K=KAY(II,IJ,IQ,JQ)
               PPX(K)=PPX(K)+DELY(J)/6.
               PPY(K)=PPY(K)-DELX(J)/6.
 30   CONTINUE
C
C Premulitply PPX and PPY by inverse of SV so that now,
C 
C     PPX = <Pi*dP/dxj>/<Pi>
C     PPY = <Pi*dP/dyj>/<Pi>
C
      CALL SPRSINVMLT(PPX,IQ,SV,NN)
      CALL SPRSINVMLT(PPY,IQ,SV,NN)
C
C Normalize 
C
      CALL SPRSNORM(PPX,IQ,NN)
      CALL SPRSNORM(PPY,IQ,NN)
C
C Format Statements
C
 110  FORMAT(//,'Determined during sparse matrix setup:',
     &        /,'     Sparse factor NFTR = ',I6)
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE VELOCITY5xc(NN,X,Y,HDOWN,SUMX,SUMY,
     &INBE,NBE,NORM,NVN,KST,LRVS,NR,NPCN,LST,COR,OMEGA,
     &NNV,Z,AK,ENZ,RHOX,RHOY,ATMX,ATMY,RHS,IQ,JQ,PPX,PPY,
     &UVERT,VVERT,WVERT,UBAR,VBAR,DOWVERT)
C-----------------------------------------------------------------------
C
C Purpose: To recover the 3-D velocity solution after the surface 
C            elevation solution has been computed.
C
C Method:  Finite Element Method on 2-D triangular C0 (linear) 
C            horizontal elements and 1-D linear vertical elements. 
C            using nodal quadrature. 
C
C Index range: I=1,NN, J=1,NNV, L=1,NE, K=1,NBE  
C
C inputs:  NN - number of horizontal nodes
C          X(I),Y(I) - cartesian nodal coordinates
C          HDOWN(I) - Bathymetric depth
C          SUMX(I) - cosine value for boundary nodes where V.n is 
C                          specified
C          SUMY(I) - sine value for boundary nodes where V.n is 
C                          specified
C          INBE(2,K) - LABLIB STANDARD geostrophic boundary element
C                          incidence list
C          NBE - The number of geostrophic boundary elements
C          NORM(I) - boundary node numbers where normal velocity is set
C                          I =1,NVN       are V.n .NE. 0.0
C                          I =KST,NNDIM   are V.n=0.0
C          LRVS(I) - node numbers of set depth and set velocity nodes
C                          I=1,NR         are set depth nodes
C                          I=NR+1,NR+NPCN are V=0.0
C                          I=LST,NNDIM    are V .NE. 0.0
C          COR - Coriolis Parameter
C          OMEGA - Radian Frequency (rad/sec)
C          NNV - number of vertical nodes
C          Z(I,J) - Vertical coordinate of 3-D nodes
C          AK(I) - Linear bottom stress coefficient
C          ENZ(I,J) - Vertical viscosity
C          RHOX(I,J) - X component of baroclinic forcing
C          RHOY(I,J) - Y component of baroclinic forcing
C          ATMX(I) - X component of atmospheric forcing (H*PSI)
C          ATMY(I) - Y component of atmospheric forcing (H*PSI)
C          RHS(I) - the free surface elevation at nodes
C          IQ(I),JQ(I*NFTR) - SPRSPAK STANDARD pointers for PPX and PPY
C          PPX(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix
C          PPY(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix 
C
C Outputs: UVERT(I,J),VVERT(I,J),WVERT(I,J) - x,y,and z components of 
C                          velocity
C          UBAR(I) - Vertical average of UVERT(I,J)
C          VBAR(I) - Vertical average of VVERT(I,J)
C
C History: Developed by Daniel R. Lynch and Christopher E. Naimie
C          Dartmouth College
C          FUNDY4 - Part of Main Routine written by Daniel R. Lynch
C          FUNDY5 - Converted to a Subroutine by Christopher. E. Naimie
C
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Global Arrays
C
      INTEGER INBE(2,NBEDIM)
      INTEGER IQ(NNDIM),JQ(NFTRDIM*NNDIM)
      REAL PPX(NFTRDIM*NNDIM),PPY(NFTRDIM*NNDIM)
      INTEGER LRVS(NNDIM),NORM(NNDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),AK(NNDIM),
     &SUMX(NNDIM),SUMY(NNDIM)
      REAL Z(NNDIM,NNVDIM),ENZ(NNDIM,NNVDIM)
      COMPLEX ATMX(NNDIM),ATMY(NNDIM),UBAR(NNDIM),VBAR(NNDIM),
     &        RHS(NNDIM)
      COMPLEX RHOX(NNDIM,NNVDIM),RHOY(NNDIM,NNVDIM),UVERT(NNDIM,NNVDIM),
     *        VVERT(NNDIM,NNVDIM),WVERT(NNDIM,NNVDIM)
C
C Local Arrays
C
      REAL DZZDX(NNVDIM),DZZDY(NNVDIM)
      COMPLEX DZDX(NNDIM),DZDY(NNDIM),VN(NNDIM)
      COMPLEX P1P(NNVDIM),P1M(NNVDIM),P2P(NNVDIM),P2M(NNVDIM),
     &        P3P(NNVDIM),P3M(NNVDIM),RP(NNVDIM),RM(NNVDIM),
     &        DUZDX(NNVDIM),DUZDY(NNVDIM),DVZDX(NNVDIM),DVZDY(NNVDIM),
     &        ACONT(NNVDIM),BCONT(NNVDIM),CCONT(NNVDIM),DCONT(NNVDIM)
C
C Constant declarations
C
      COMPLEX EYE,DENOM,UU,VV,WW,DZTOP,DU,DIVTOP,DIVBOT,DZBOT,DV
      COMPLEX GP,GM,PSIP,PSIM,UVP,UVM,RBP,RBM,DP,DM

      LOGICAL DOWVERT
C
C Constants
C
      EYE=cmplx(0.0,1.0)
      OMEGAP=OMEGA+COR
      OMEGAM=OMEGA-COR
C
C Compute the nodal gravity terms
C
      CALL SPRSMLTRC(PPX,IQ,JQ,RHS,DZDX,NN)
      CALL SPRSMLTRC(PPY,IQ,JQ,RHS,DZDY,NN)
C
C DZDX,DZDY are the nodal values of GRAD ZETA
C
C Begin nodal loop for velocity
C
      DO 20 I=1,NN
C
C  Compute G+, G-, PSI+, PSI- at this node
C    (Equations 13-14 in Users Manual)
C
         DENOM=1.0/HDOWN(I)
         GP=-G*0.5*(DZDX(I)+EYE*DZDY(I))
         GM=-G*0.5*(DZDX(I)-EYE*DZDY(I))
         PSIP= 0.5*(ATMX(I)+EYE*ATMY(I))*DENOM
         PSIM= 0.5*(ATMX(I)-EYE*ATMY(I))*DENOM
C
C  Compute vertical profiles of RP, RM
C    (Equation 15 in Users Manual)
C
         DO 10 J=1,NNV
            RP(J)=.5*(RHOX(I,J)+EYE*RHOY(I,J))
            RM(J)=.5*(RHOX(I,J)-EYE*RHOY(I,J))
 10      CONTINUE
C
C  Call VERTSTRUC5 to compute vertical structure of P1+-,P2+-,P3+-
C    (Equations 10-12,16-17 in Users Manual)
C
         CALL VERTSTRUC5(NN,NNV,I,OMEGAP,Z,ENZ,AK,RP,
     &                   P1P,P2P,P3P)
         CALL VERTSTRUC5(NN,NNV,I,OMEGAM,Z,ENZ,AK,RM,
     &                   P1M,P2M,P3M)
         DO 20 J=1,NNV
C
C  Assemble NU+,NU-
C    (Equation 16 in Users Manual)
C
            UVP=GP*P1P(J)+PSIP*P2P(J)+P3P(J)
            UVM=GM*P1M(J)+PSIM*P2M(J)+P3M(J)
C
C  Convert to U,V
C    (Equation 9 in Users Manual)
C
            UVERT(I,J)=UVP+UVM
            VVERT(I,J)=(UVP-UVM)/EYE
 20   CONTINUE
C
C End of horizontal velocity node loop.
C UVERT, VVERT now contain the complete velocity solution.
C Compute the vertically averaged velocities
C
      CALL VERTAVG(UVERT,UBAR,Z,NN,NNV,NNDIM,NNVDIM)
      CALL VERTAVG(VVERT,VBAR,Z,NN,NNV,NNDIM,NNVDIM)
C
C Apply horizontal velocity BC's
C  definitions: CO=cosine of unit normal vector
C               SI=sine of unit normal vector
C               UU=normal velocity    = CO*U+SI*V
C               VV=tangential velocity=-SI*U+CO*V
C               U=x velocity          = CO*UU-SI*VV
C               V=y velocity          = SI*UU+CO*VV
C
C  Apply V.n=0 BC's
C
      DO 30 I=KST,NNDIM
         N=NORM(I)
         SI=SUMX(N)
         CO=SUMY(N)
         UU=0.0
         VV=-SI*UBAR(N)+CO*VBAR(N)
         UBAR(N)=CO*UU-SI*VV
         VBAR(N)=SI*UU+CO*VV
         DO 30 J=1,NNV
            VV=-SI*UVERT(N,J)+CO*VVERT(N,J)
            UVERT(N,J)=CO*UU-SI*VV
            VVERT(N,J)=SI*UU+CO*VV
 30   CONTINUE
C      
C  Initialize arrays used by Subroutine BC to set boundary conditions 
C   Note: DZDY is used as dummy arrays in this section
C
      DO 49 I=1,NN
         VN(I)=cmplx(0.0,0.0)
 49   CONTINUE
C
C  Apply V.n.NE.0  BC's
C
      IF(NVN.LT.1) GO TO 51
      CALL BC(W,VN,DZDY,NNDIM,2)
      DO 50 I=1,NVN
         N=NORM(I)
         SI=SUMX(N)
         CO=SUMY(N)
         UU=VN(N)
         VV=-SI*UBAR(N)+CO*VBAR(N)
         UBAR(N)=CO*UU-SI*VV
         VBAR(N)=SI*UU+CO*VV
 50   CONTINUE
C
C  Apply fixed velocity node BC's (corners)
C
C   V=0 nodes 
C
 51   J1=NR+1
      J2=NR+NPCN
      IF(J2.LT.J1) GO TO 60
      DO 59 I=J1,J2
         N=LRVS(I)
         UBAR(N)=(0.,0.)
         VBAR(N)=(0.,0.)
         DO 59 J=1,NNV
            UVERT(N,J)=cmplx(0.,0.)
            VVERT(N,J)=cmplx(0.,0.)
 59   CONTINUE
 60   CONTINUE
C
C   V.NE.0 nodes 
C
      IF (LST.LE.NNDIM) CALL BC(W,UBAR,VBAR,NNDIM,3)
C
C All horizontal velocity BC's have been applied
C
C Compute the vertical velocities using d(continuity)/dZ in node loop
C

      IF(DOWVERT) THEN
      DO 70 I=1,NN
         CALL SPRSCONVRC(I,PPX,PPY,IQ,JQ,UVERT,VVERT,Z,
     &   DUZDX,DUZDY,DVZDX,DVZDY,DZZDX,DZZDY,NNV,NNDIM)
C
C  Build tridiagonal system for interior nodes
C
         DO 69 J=2,NNV-1
            DZTOP= 1.0/(Z(I,J+1)- Z(I,J))
            DU=UVERT(I,J+1)-UVERT(I,J)
            DV=VVERT(I,J+1)-VVERT(I,J)
            DIVTOP=DUZDX(J+1)-DZZDX(J+1)*DU*DZTOP
     &            +DVZDY(J+1)-DZZDY(J+1)*DV*DZTOP
            DZBOT= 1.0/(Z(I,J)- Z(I,J-1))
            DU=UVERT(I,J)-UVERT(I,J-1)
            DV=VVERT(I,J)-VVERT(I,J-1)
            DIVBOT=DUZDX(J-1)-DZZDX(J-1)*DU*DZBOT
     &           + DVZDY(J-1)-DZZDY(J-1)*DV*DZBOT
            ACONT(J)=DZBOT
            CCONT(J)=DZTOP
            BCONT(J)=-ACONT(J)-CCONT(J)
            DCONT(J)=0.5*(-DIVTOP+DIVBOT)
 69      CONTINUE
C
C  Set bottom and surface BC's
C
         WVERT(I,1)=UVERT(I,1)*DZZDX(1) + VVERT(I,1)*DZZDY(1)
         ACONT(1)=cmplx(0.,0.)
         BCONT(1)=cmplx(1.,0.)
         CCONT(1)=cmplx(0.,0.)
         DCONT(1)=WVERT(I,1)
         WVERT(I,NNV)=EYE*OMEGA*RHS(I)
         ACONT(NNV)=cmplx(0.,0.)
         BCONT(NNV)=cmplx(1.,0.)
         CCONT(NNV)=cmplx(0.,0.)
         DCONT(NNV)=WVERT(I,NNV)
C
C
C  Solve the tridiagonal system
C
         CALL CTHOMAS(3,ACONT,BCONT,CCONT,DCONT,NNV)
C
C  Set DCONT to WVERT at this node
C
         DO 70 J=1,NNV
           WVERT(I,J)=DCONT(J)
 70   CONTINUE
      ELSE 
         DO 80 J=1,NNV
           WVERT(I,J)=CMPLX(0.0,0.0) 
 80   CONTINUE
	 

      ENDIF
C
C The velocity solution is complete
C
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE UNIGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,NNV,EN0,AK0,
     *      Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
C
C Uniform mesh and coefficients version of vertgrid
C
C
C Global Arrays
C
      DIMENSION Z(NNDIM,NNVDIM),AK(NNDIM),ENZ(NNDIM,NNVDIM),
     * HDOWN(NNDIM)
      COMPLEX RX(NNDIM,NNVDIM),RY(NNDIM,NNVDIM)
      COMPLEX HRBARELX(NEDIM),HRBARELY(NEDIM)
C
      DO 20 I=1,NN
         DZ=HDOWN(I)/(NNV-1.)
         ZZ=-HDOWN(I)
         AK(I)=AK0
         DO 10 J=1,NNV
            ENZ(I,J)=EN0
            Z(I,J)=ZZ
            RX(I,J)=cmplx(0.0, 0.0)
            RY(I,J)=cmplx(0.0, 0.0)
            ZZ=ZZ+DZ
 10      CONTINUE
         Z(I,NNV)=0.0
 20   CONTINUE
      DO 30 I=1,NE
         HRBARELX(I)=0.0
         HRBARELY(I)=0.0
 30   CONTINUE
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE SINEGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,NNV,EN0,AK0,DZBL,
     *      Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
C
C Sinusoidal mesh and coefficients version of vertgrid
C
C
C Global Arrays
C
      DIMENSION Z(NNDIM,NNVDIM),AK(NNDIM),ENZ(NNDIM,NNVDIM),
     * HDOWN(NNDIM)
      COMPLEX RX(NNDIM,NNVDIM),RY(NNDIM,NNVDIM)
      COMPLEX HRBARELX(NEDIM),HRBARELY(NEDIM)
C
C Set sinusoidal grid
C
      NEV=NNV-1
      PI=ACOS(-1.)
      DO 20 I=1,NN
         A=(HDOWN(i)-DZBL*NEV)/(NEV*SIN(2.*PI/NEV))
         IF(A.LT.0.)A=0.
         Z(I,1)=-HDOWN(I)
         DO 10 J=2,NNV
            EPS=(J-1.0)/(NNV-1.0)
            Z(i,j)=-HDOWN(i)*(1.-EPS) -A*SIN(2.*PI*EPS)
  10     CONTINUE
         Z(i,NNV)=0.
  20  CONTINUE
C
C Initialize other arrays
C
      DO 30 I=1,NN
         AK(I)=AK0
         DO 30 J=1,NNV
            ENZ(I,J)=EN0
            RX(I,J)=cmplx(0.0, 0.0)
            RY(I,J)=cmplx(0.0, 0.0)
   30 CONTINUE
      DO 40 I=1,NE
         HRBARELX(I)=0.0
         HRBARELY(I)=0.0
 40   CONTINUE
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE RHOGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,NNV,EN0,AK0,
     *      Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
C
C Version of vertgrid where Z,RX,RY,HRBARELX,HRBARELY are read from a 
C    .rho file
C
C Global Arrays
C
      DIMENSION Z(NNDIM,NNVDIM),AK(NNDIM),ENZ(NNDIM,NNVDIM)
      COMPLEX RX(NNDIM,NNVDIM),RY(NNDIM,NNVDIM)
      COMPLEX HRBARELX(NEDIM),HRBARELY(NEDIM)
C
C Local assignments
C
      CHARACTER*10 HEADER
      CHARACTER*72 RHOFILE
C
C Read Z,RX,RY from user specified .rho file
C Initialize AK and ENO
C
      write(*,*)'Enter the name of the .rho baroclinic forcing file'
      read(*,101)RHOFILE
      open(19,FILE=RHOFILE,STATUS='OLD')
      read(19,101)HEADER
      read(19,101)HEADER
      read(19,*)NNV
      write(*,*)'    Loading vertical grid/baroclinic forcing from:'
      write(*,*)'      ',RHOFILE
      DO 10 I=1,NN
         IF(MOD(I,500).EQ.0)write(*,*)'        Loading data for node',i
         AK(i)=AK0
         DO 10 J=1,NNV
            read(19,*)INNN,Z(I,J),XRHO,YRHO
            RX(I,J)=CMPLX(XRHO,0.0)
            RY(I,J)=CMPLX(YRHO,0.0)
            ENZ(I,J)=EN0
 10   CONTINUE
      read(19,101)HEADER
C
C Read HRBARELX,HRBARELY from specified .rho file
C
      do 20 I=1,NE
         read(19,*)INE,HRBARX,HRBARY
         HRBARELX(I)=CMPLX(HRBARX,0.0)
         HRBARELY(I)=CMPLX(HRBARY,0.0)
 20   CONTINUE
      CLOSE(19)
C
C Format statements
C
 101  FORMAT(a)
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE VERTSTRUC5(NN,NNV,I,OMEGA,Z,ENZ,AK,R,
     * VG,VW,VR)
C
C 1-D Galerkin FE
C Linear basis
C Arbitrary mesh spacing Z(I,J)
C Nodal interpolation of ENZ
C Exact integration
C All computations in original, dimensional terms
C
C
C Include file for parameter statements
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Global Arrays
C
      COMPLEX R(NNVDIM)
      DIMENSION Z(NNDIM,NNVDIM),ENZ(NNDIM,NNVDIM),AK(NNDIM)
      COMPLEX VG(NNVDIM),VW(NNVDIM),VR(NNVDIM)
C
C Local Arrays and variables
C
      COMPLEX A(NNVDIM),B(NNVDIM),C(NNVDIM)
      COMPLEX EYE
      EYE=(0.0,1.0)
C
C VG, VW, and VR are the right-sides of the tridiagonal system for
C   gravity, wind, and density gradient, respectively.
C
C Equation 1 (bottom node)
C
      J=1
      DZT=Z(I,J+1) - Z(I,J)
      DZB=0.0
      ENT=( ENZ(I,J+1) + ENZ(I,J)   )/2.
      ENB=0.0
      C(J)=EYE*OMEGA*DZT/6. -ENT/DZT
      A(J)=cmplx(0.,0.)
      B(J)=EYE*OMEGA*(DZT)/3. + ENT/DZT
      B(J)=B(J)+AK(I)
      VG(J)=(DZT)/2.
      VW(J)=cmplx(0.,0.)
      VR(J)=R(J+1)*DZT/6. + R(J)*(DZT/3.)
C
C Interior equations
C
      DO 10 J=2,NNV-1
         DZT=Z(I,J+1) - Z(I,J)
         DZB=Z(I,J)   - Z(I,J-1)
         ENT=( ENZ(I,J+1) + ENZ(I,J)   )/2.
         ENB=( ENZ(I,J)   + ENZ(I,J-1) )/2.
         C(J)=EYE*OMEGA*DZT/6. -ENT/DZT
         A(J)=EYE*OMEGA*DZB/6. -ENB/DZB
         B(J)=EYE*OMEGA*(DZT+DZB)/3. + ENT/DZT + ENB/DZB
         VG(J)=(DZT+DZB)/2.
         VW(J)=cmplx(0.,0.)
         VR(J)=R(J+1)*DZT/6. + R(J)*(DZT/3. +DZB/3.) + R(J-1)*DZB/6.
   10 CONTINUE
C
C Equation NNV (surface node)
C
      J=NNV
      DZT=0.0
      DZB=Z(I,J) - Z(I,J-1)
      ENT=0.0
      ENB=( ENZ(I,J)   + ENZ(I,J-1) )/2.
      C(J)=cmplx(0.,0.)
      A(J)=EYE*OMEGA*DZB/6. -ENB/DZB
      B(J)=EYE*OMEGA*(DZT+DZB)/3. + ENB/DZB
      VG(J)=(DZB)/2.
      VW(J)=Z(I,NNV)-Z(I,1)
C   PSI=1, THEREFORE H*PSI = HDOWN.
      VR(J)= R(J)*(DZB/3.) + R(J-1)*DZB/6.
C
C Done with setup; solve
C
      CALL CTHOMAS(3,A,B,C,VG,NNV)
      CALL CTHOMAS(2,A,B,C,VW,NNV)
      CALL CTHOMAS(2,A,B,C,VR,NNV)
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY4
C
      SUBROUTINE CTHOMAS(KKK,A,B,C,R,N)
C
C THOMAS ALGORITHM FOR TRIDIAGONAL MATRICES  -  SYSTEM FORTRAN-77
C COMPLEX VERSION
C A, B, AND C = LHS COLUMNS, DIMENSIONED (=>N) IN CALLING PROGRAM
C R = RIGHT-HAND-SIDE, DIMENSIONED (=>N) IN CALLING PROGRAM
C N = # OF EQUATIONS
C SOLUTION RETURNS IN R WHEN KKK => 2
C KKK = 1 PERFORMS LU DECOMPOSITION DESTRUCTIVELY
C KKK = 2 PERFORMS BACK SUBSTITUTION
C KKK = 3 PERFORMS OPTIONS 1 AND 2
C
      COMPLEX A(N),B(N),C(N),R(N)
C
C  STEP 1 - TRIANGULARIZE MATRIX A USING DOOLITTLE METHOD
C
      IF(KKK.EQ.2)GO TO 45
      C(1)=C(1)/B(1)
      DO 40 I=2,N
      IM=I-1
      C(I)=C(I)/(B(I)-A(I)*C(IM))
   40 CONTINUE
      IF(KKK.EQ.1)RETURN
C
C STEP 2 - MODIFY LOAD VECTOR R
C
   45 CONTINUE
      R(1)=R(1)/B(1)
      DO 41 I=2,N
      IM=I-1
      R(I)=(R(I)-A(I)*R(IM))/(B(I)-A(I)*C(IM))
   41 CONTINUE
 
C
C BACK SUBSTITUTE SOLUTION INTO R VECTOR
C
      DO 50 I=1,N-1
      J=N-I
      JP=J+1
      R(J)=R(J)-C(J)*R(JP)
   50 CONTINUE
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE CSOLVET(KKK,B,R,NEQ,IHALFB,NDIM,MDIM)
C
C   ASYMMETRIC BAND MATRIX EQUATION SOLVER
C   COMPLEX VERSION
C   DOCTORED TO IGNORE ZEROS IN LU DECOMP. STEP
C
C   KKK=1 TRIANGULARIZES THE BAND MATRIX B
C   KKK=2 SOLVES FOR RIGHT SIDE R, SOLUTION RETURNS IN R
C
      COMPLEX  B(MDIM,NDIM),R(NDIM)
      COMPLEX  PIVOT,C,SUM
      NRS=NEQ-1
      IHBP=IHALFB+1
      IF (KKK.EQ.2) GO TO 30
C
C  TRIANGULARIZE MATRIX A USING DOOLITTLE METHOD
C
      DO 10 K=1,NRS
      PIVOT=B(IHBP,K)
      KK=K+1
      KC=IHBP
      DO 21 I=KK,NEQ
      KC=KC-1
      IF(KC.LE.0) GO TO 10
      C=-B(KC,I)/PIVOT
      B(KC,I)=C
      IF(CABS(C).EQ.0.0) GO TO 21
      KI=KC+1
      LIM=KC+IHALFB
      DO 20 J=KI,LIM
      JC=IHBP+J-KC
   20 B(J,I)=B(J,I)+C*B(JC,K)
   21 CONTINUE
   10 CONTINUE
      GO TO 100
C
C  MODIFY LOAD VECTOR R
C
   30 NN=NEQ+1
      IBAND=2*IHALFB+1
      DO 40 I=2,NEQ
      JC=IHBP-I+1
      JI=1
      IF (JC.LE.0) GO TO 50
      GO TO 60
   50 JC=1
      JI=I-IHBP+1
   60 SUM=0.0
      DO 70 J=JC,IHALFB
      SUM=SUM+B(J,I)*R(JI)
   70 JI=JI+1
   40 R(I)=R(I)+SUM
C
C   BACK SOLUTION
C
      R(NEQ)=R(NEQ)/B(IHBP,NEQ)
      DO 80 IBACK=2,NEQ
      I=NN-IBACK
      JP=I
      KR=IHBP+1
      MR=MIN0(IBAND,IHALFB+IBACK)
      SUM=0.0
      DO 90 J=KR,MR
      JP=JP+1
   90 SUM=SUM+B(J,I)*R(JP)
   80 R(I)=(R(I)-SUM)/B(IHBP,I)
  100 RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      FUNCTION PHASELAGD(S)
C
C COMPUTES THE NEGATIVE ARGUMENT OF A COMPLEX NUMBER S, IN DEGREES.
C RESULT IS REPORTED IN THE RANGE 0. TO 360. 
C 
      COMPLEX S
      PI=3.1415926535 
      FACTOR=180./PI
      PHASELAGD = 0.0
      IF(CABS(S).GT.0.)  PHASELAGD = -FACTOR*ATAN2(AIMAG(S), REAL(S))
      IF(PHASELAGD.LT.0.) PHASELAGD = PHASELAGD + 360.  
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY4
C
      FUNCTION PHASELAG(S)
C
C COMPUTES THE NEGATIVE ARGUMENT OF A COMPLEX NUMBER S, IN RADIANS. 
C RESULT IS REPORTED IN THE RANGE 0. TO 2*PI. 
C 
      COMPLEX S
      PI=3.1415926535 
      PHASELAG = 0.0
      IF(CABS(S).GT.0.)  PHASELAG = -ATAN2(AIMAG(S), REAL(S))
      IF(PHASELAG.LT.0.) PHASELAG = PHASELAG + 2.*PI
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY4
C
        SUBROUTINE ELLIPSE (U,V,amaj,amin,ainc,g,apl,amn,gpl,gmin)
C
C COMPUTES TIDAL ELLIPSE PARAMETERS FROM COMPLEX VELOCITY AMPLITUDES U,V.
C
        COMPLEX U,V
        PI=3.1415926535898
        FAC=PI/180.
C
C TAKEN FROM TEMP.FOR AT IOS
C
        CX=REAL(U)
        SX=-AIMAG(U)
        CY=REAL(V)
        SY=-AIMAG(V)
	cxpsy=0.5*(cx+sy)
	cymsx=0.5*(cy-sx)
	cxmsy=0.5*(cx-sy)
	cypsx=0.5*(cy+sx)
	apl=sqrt(cxpsy**2+cymsx**2)
	amn=sqrt(cxmsy**2+cypsx**2)
	if(apl.gt.1.e-5) then
	epl=atan2(cymsx,cxpsy)/fac
	else
	epl=0.
	end if
	if(amn.gt.1.e-5) then
	emin=atan2(cypsx,cxmsy)/fac
	else
	emin=0.
	end if
	amaj=apl+amn
	amin=apl-amn
	gpl=-epl
	gmin=emin
	if(gmin-gpl.lt.0.) gmin=gmin+360.
	g=0.5*(gpl+gmin)
	ainc=0.5*(gmin-gpl)
	if(gmin.lt.0.) gmin=gmin+360.
	if(gmin.gt.360.) gmin=gmin-360.
	if(gpl.lt.0.) gpl=gpl+360.
        RETURN
        END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE DUMP5(X,Y,Z,ZETA,UBAR,VBAR,U,V,W,
     * NN,NNV,NNDIM,NNVDIM)
      DIMENSION X(NNDIM),Y(NNDIM),Z(NNDIM,NNVDIM)
      COMPLEX ZETA(NNDIM),UBAR(NNDIM),VBAR(NNDIM), 
     *     U(NNDIM,NNVDIM),V(NNDIM,NNVDIM),W(NNDIM,NNVDIM)
C        
C-----------------------------------------------------------------------
C
C A STANDARD OUTPUT ROUTINE TO WRITE ALL OUTPUT TO A USER-SPECIFIED FILE
C
C X,Y,Z ARE THE NODAL COORDINATES
C ZETA(I) IS THE SURFACE ELEVATION      
C UBAR,VBAR ARE THE VERTICALLY AVERAGED VELOCITIES
C U,V,W ARE THE VELOCITIES AT DEPTH
C
C THERE ARE NN NODES IN THE HORIZONTAL MESH
C THERE ARE NNV NODES BENEATH EACH HORIZONTAL NODE
C INDICES FOR 3-D QUANTITIES: (I,J)
C      I=HORIZONTAL NODE NUMBER
C      J=VERTICAL NODE NUMBER
C          J=1 IS THE BOTTOM
C          J=NNV IS THE TOP
C
C NNDIM, NNVDIM ARE PASSED STRICTLY FOR DIMENSIONING PURPOSES. 
C        
C-----------------------------------------------------------------------        
C 
      CHARACTER*80 FILENAME
      WRITE(6,*) ' ENTER NAME FOR DUMP FILE CONTAINING ALL OUTPUT'
      READ(5,7001)FILENAME
 7001 FORMAT(A)
      OPEN (UNIT=4,FILE=FILENAME,STATUS='NEW')
      DO 500 I=1,NN
      WRITE(4,*)I,X(I),Y(I),-Z(I,1),ZETA(I),UBAR(I),VBAR(I)
      DO 500 J=1,NNV
      WRITE(4,*)J,Z(I,J),U(I,J),V(I,J),W(I,J)
  500 CONTINUE
c drl 9 july 97
      close(4)
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE VERTAVG(F,FBAR,Z,NN,NNV,NNDIM,NNVDIM)
      COMPLEX F(NNDIM,NNVDIM), FBAR(NNDIM)
      DIMENSION Z(NNDIM,NNVDIM)
C
C COMPUTES THE VERTICAL AVERAGE OF F, PUTS IT IN FBAR.
C DOES THIS FOR ALL HORIZONTAL NODES
C
      DO 100 I=1,NN
      FBAR(I)=(0.,0.)
      DO 200 J=2,NNV
      FBAR(I)=FBAR(I)+0.5*(F(I,J) + F(I,J-1))*(Z(I,J) - Z(I,J-1))
  200 CONTINUE
      FBAR(I)=FBAR(I)/(Z(I,NNV)-Z(I,1))
  100 CONTINUE
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE BAROCLINIC5(nn,ne,nnv,x,y,in,z,rx,ry,hrbarelx,hrbarely)
C-----------------------------------------------------------------------
C purpose: This routine computes real baroclinic terms for zero 
C           frequency FUNDY5 solutions and transforms them into the 
C           complex arrays required by FUNDY5.  This routine is designed
C           to be called from Subroutine VERTGRID5.
C
C index range:  i=1,nn, j=1,nnv, k=1,ne
C
C inputs:  nn - number of horizontal nodes
C          ne - the number of elements
C          nnv - number of vertical nodes
C          x(i),y(i) - cartesian nodal coordinates
C          in - incidence list
C          z(i,j) - z coordinate of sigma mesh 3-D nodes
C
C complex outputs:
C          rx(i,j) - x-directed nodal based baroclinic forcing 
C          ry(i,j) - y-directed nodal based baroclinic forcing
C          hrbarelx(k) - element based depth*vertical integral of rx
C          hrbarely(k) - element based depth*vertical integral of ry
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           JANUARY 4, 1993
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Global variables
C
      integer in(3,nedim)
      real x(nndim),y(nndim)
      real z(nndim,nnvdim)
      complex rx(nndim,nnvdim),ry(nndim,nnvdim)
      complex hrbarelx(nedim),hrbarely(nedim)
C
C Local variables
C
      real rhox(nndim,nnvdim),rhoy(nndim,nnvdim)
      real hrbarx(nedim),hrbary(nedim)
C
C Compute real baroclinic forcing
C
      CALL RHOHRBAR(nn,ne,nnv,x,y,in,z,rhox,rhoy,hrbarx,hrbary)
C
C Transform real baroclinic forcings to complex forcings
C
      do 10 j=1,nnv
         do 10 i=1,nn
            rx(i,j)=CMPLX(rhox(i,j),0.0)
            ry(i,j)=CMPLX(rhoy(i,j),0.0)
 10   continue
      do 20 k=1,ne
         hrbarelx(k)=CMPLX(hrbarx(k),0.0)
         hrbarely(k)=CMPLX(hrbary(k),0.0)
 20   continue
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE RHOHRBAR(nn,ne,nnv,x,y,in,zs,rx,ry,hrbarx,hrbary)
C-----------------------------------------------------------------------
C purpose: This routine computes baroclinic forcing terms required by
C           FUNDY5.
C
C limitations:  This routine, and all subroutines called by it compute
C                real quantities.
C
C index range:  i=1,nn, j=1,nnv, k=1,ne, l=1,nlev
C
C inputs:  nn - number of horizontal nodes
C          ne - the number of elements
C          nnv - number of vertical nodes
C          x(i),y(i) - cartesian nodal coordinates
C          in - incidence list
C          zs(i,j) - z coordinate of sigma mesh 3-D nodes
C          nlev - number of level surfaces in .lst file
C          zl(l) - z coordinate of level mesh nodes
C          sigtl(i,l) - level mesh sigmat data
C
C outputs: rx(i,j) - x-directed nodal based baroclinic forcing 
C          ry(i,j) - y-directed nodal based baroclinic forcing
C          hrbarx(k) - depth*vertical integral of rx
C          hrbary(k) - depth*vertical integral of ry
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           JANUARY 1, 1993
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Global arrays
C
C  level mesh
      real sigtl(nndim,nlevels),zl(nlevels)
C  sigma mesh
      integer in(3,nedim)
      real x(nndim),y(nndim)
      real zs(nndim,nnvdim)
      real rx(nndim,nnvdim),ry(nndim,nnvdim)
      real hrbarx(nedim),hrbary(nedim)
C
C Local arrays
C
C  level mesh 
      dimension drdxl(nndim,nlevels),drdyl(nndim,nlevels)
      dimension rxl(nndim,nlevels),ryl(nndim,nlevels)
      dimension drdxle(nedim,nlevels),drdyle(nedim,nlevels)
      dimension rxle(nedim,nlevels),ryle(nedim,nlevels)
C  sigma mesh 
      real xe(nedim),ye(nedim)
      real zse(nedim,nnvdim)
      real rxe(nedim,nnvdim),rye(nedim,nnvdim)
C
C Local assignments
C
      character*72 filename
C
C Check oversize of arrays
C
      nflag=0
      if(nn.GT.nndim)then
      write(*,*)'The number of nodes is to large for current parameters.
     &'
      write(*,*)'  => Increase NNDIM to at least ',nn,' in the include f
     &ile TRUXTON1ARRAYS.DIM and recompile.'
      nflag=1
      endif
      if(ne.GT.nedim)then
      write(*,*)'The number of elements is to large for current paramete
     &rs.'
      write(*,*)'  => Increase NEDIM to at least ',ne,' in the include f
     &ile TRUXTON1ARRAYS.DIM and recompile.'
      nflag=1
      endif
      if(nnv.GT.nnvdim)then
      write(*,*)'The number of vertical nodes is to large for current pa
     &rameters.'
      write(*,*)'  => Increase NNVDIM to at least ',nnv,' in the include
     & file TRUXTON1ARRAYS.DIM and recompile.'
      nflag=1
      endif
      if(nflag.EQ.1)stop
C
C Inform user of values of rhoref
C
      write(*,*)' '
      write(*,*)'The reference density (MKS units) is set in the Include
     & File TRUXTON1ARRAYS.DIM:  Currently RHOREF =',rhoref
      write(*,*)'If another value is desired, change the value of RHOREF
     & in the TRUXTON1ARRAYS.DIM and recompile.'
      write(*,*)' '
C
C Load level mesh nodal sigmat data
C
      write(*,*)'Enter the name of the level surface sigmat .lst file'
      read(*,101)filename
      CALL LOADLST(nndim,nlevels,nn,nlev,filename,zl,sigtl)
C
C Compute level mesh nodal values of drdxl and drdyl
C
      CALL LST2GRADRHO(nn,ne,x,y,in,nlev,sigtl,drdxl,drdyl)
C
C Compute level mesh nodal values of RHOX and RHOY by integrating on
C  level mesh
C
      CALL LGRA2LRHO(nndim,nlevels,g,nn,nlev,zl,drdxl,drdyl,rxl,ryl)
C
C Compute sigma mesh nodal values of RHOX and RHOY by integrating on
C  level mesh
C
      CALL LRHO2SRHO(nndim,nnvdim,nlevels,g,nn,nnv,nlev,
     &zl,drdxl,drdyl,rxl,ryl,zs,rx,ry)
C
C Compute element sigma and level mesh quantities from nodal quantities
C
      CALL NOD2EL(nndim,nedim,nnvdim,nlevels,nn,ne,nnv,rhoref,
     &x,y,in,zs,nlev,sigtl,xe,ye,zse,drdxle,drdyle)
C
C Compute level mesh element values of RHOX and RHOY by integrating on
C  level mesh
C
      CALL LGRA2LRHO(nedim,nlevels,g,ne,nlev,zl,drdxle,drdyle,rxle,ryle)
C
C NOTE: RHOXL and RHOYL now contain level mesh element gradients
C
C Compute sigma mesh elemental values of RHOXL and RHOYL by integrating
C  on level mesh
C
      CALL LRHO2SRHO(nedim,nnvdim,nlevels,g,ne,nnv,nlev,
     &zl,drdxle,drdyle,rxle,ryle,zse,rxe,rye)
C
C Compute sigma mesh values of H*RBAR by integrating on sigma mesh
C
      CALL SRHO2SHRBAR(nedim,nnvdim,ne,nnv,zse,rxe,rye,hrbarx,hrbary)
C
C Format statements
C
 101  format(a)
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE LOADLST(nndim,nlevels,nn,nlev,filename,zl,sigtl)
C-----------------------------------------------------------------------
C purpose: This routine loads a level surface sigmat .lst file
C
C index range: i=1,nn, l=1,nlev
C
C inputs:  nndim - parametered size of arrays indexed with i
C          nlevels - parametered size of arrays indexed with l
C          nn - number of horizontal nodes
C          filename - name of .lst file
C
C outputs: nlev - number of level surfaces in .lst file
C          zl(l) - z coordinate of level mesh nodes
C          sigtl(i,l) - sigmat at level mesh locations
C 
C .lst file format: line 1 - meshname
C                   line 2 - header
C                   line 3 - nlev
C                   line 4 - (zl(l),l=1,nlev)
C                   next nn*nlev lines - sigt(i,l)
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           JANUARY 1, 1993
C-----------------------------------------------------------------------
C
C Global arrays
C
      dimension zl(nlevels), sigtl(nndim,nlevels)
C
C Local assignments
C
      character*72 filename,header
C
C Open file
C
      open(63,file=filename,status='old')
C
C Read headers
C
      read(63,101)header
      read(63,101)header
      read(63,*)nlev
C
C Check for oversize of arrays
C
      nflag=0
      if(nn.GT.nndim)then
      write(*,*)'The number of nodes is to large for current parameters.
     &'
      write(*,*)'  => Increase NNDIM to at least ',nn
      nflag=1
      endif
      if(nlev.GT.nlevels)then
      write(*,*)'The number of level surfaces is to large for current pa
     &rameters.'
      write(*,*)'  => Increase nlevels to at least ',nlev
      nflag=1
      endif
      if(nflag.EQ.1)stop
C
C Read zl
C
      read(63,*)(zl(l),l=1,nlev)
C
C Read sigmat data
C
      do 10 i=1,nn
         do 10 l=1,nlev
            read(63,*)sigtl(i,l)
            ii=i
 10   continue
      close(63)
C
C Message to standard out
C
      write(*,*)ii,'x',nlev,' real 3-D sigmat values loaded from ',
     &filename
C
C Format statements
C
 101  format(a)
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE LST2GRADRHO(nn,ne,x,y,in,nlev,sigtl,drdxl,drdyl)
C-----------------------------------------------------------------------
C purpose: This routine computes level mesh nodal values of:
C            drdxl(i,l) - 1/rhoref*dr/dx 
C            drdyl(i,l) - 1/rhoref*dr/dy 
C
C index range: i=1,nn, l=1,nlev
C
C
C inputs:  nn - number of horizontal nodes
C          ne - the number of elements
C          x(i),y(i) - cartesian nodal coordinates
C          in - incidence list
C          nlev - number of level surfaces in .lst file
C          sigtl(i,l) - level mesh sigmat data
C
C outputs: drdxl(i,l) - array discussed above
C          drdyl(i,l) - array discussed above
C
C history:  Christopher E. Naimie/from Dr. D.R.Lynch's gradrho1a.f
C           Dartmouth College
C           JANUARY 1, 1993
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Global arrays
C
      dimension x(nndim),y(nndim),in(3,nedim)
      dimension sigtl(nndim,nlevels),drdxl(nndim,nlevels),
     &drdyl(nndim,nlevels)
C
C Local arrays
C
      dimension rightx(nndim),righty(nndim),sv(nndim),delx(3),dely(3)
C
C Initialize gradients
C
      do 10 i=1,nn
         do 10 l=1,nlev
            drdxl(i,l)=0.
            drdyl(i,l)=0.
 10   continue
C 
C Main loop: over vertical levels to differentiate horizontally 
C
      do 40 l=1,nlev
C
C initialize the sv and rhs arrays
C
         do 20 i=1,nn
            sv(i)=0.
            rightx(i)=0.
            righty(i)=0.
 20      continue
C
C element loop to assemble rhs's and stiffness matrix
C
         do 30 k=1,ne
            i1=in(1,k)
            i2=in(2,k)
            i3=in(3,k)
            delx(1)=x(i2)-x(i3)
            delx(2)=x(i3)-x(i1)
            delx(3)=x(i1)-x(i2)
            dely(1)=y(i2)-y(i3)
            dely(2) =y(i3)-y(i1)
            dely(3)=y(i1)-y(i2)
            arr=0.5*(x(i1)*dely(1) + x(i2)*dely(2) + x(i3)*dely(3))
      rx=(sigtl(i1,l)*dely(1)+sigtl(i2,l)*dely(2)+sigtl(i3,l)*dely(3))
     &       /(2.*arr*rhoref)
      ry=-(sigtl(i1,l)*delx(1)+sigtl(i2,l)*delx(2)+sigtl(i3,l)*delx(3))
     &       /(2.*arr*rhoref)
            do 30 i=1,3
               ii=in(i,k)
               sv(ii)=sv(ii)+arr/3.
               rightx(ii)=rightx(ii)+rx*arr/3.
               righty(ii)=righty(ii)+ry*arr/3.
 30      continue
C
C Done with assembly; invert for the nodal gradients
C
         do 40 i=1,nn
            if(sv(i).eq.0.)sv(i)=1.
            drdxl(i,l)=rightx(i)/sv(i)
            drdyl(i,l)=righty(i)/sv(i)
 40   continue
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE LGRA2LRHO(ndim,nlevels,g,nn,nlev,zl,drdxl,drdyl,
     &rxl,ryl)
C-----------------------------------------------------------------------
C purpose: This routine performs the level mesh integrations:
C            rxl(i,l)=-g*(vertical integral of drdxl(i,z) from
C                          z=zl(i) to z=0) 
C            ryl(i,l)=-g*(vertical integral of drdyl(i,z) from
C                          z=zl(i) to z=0) 
C
C index range: i=1,nn, l=1,nlev
C
C
C inputs:  ndim - parametered size of arrays indexed with i
C          nlevels - parametered size of arrays indexed with l
C          g - gravitational acceleration (MKS units) = 9.806 m/sec
C          nn - number of horizontal nodes
C          nlev - number of level surfaces in .lst file
C          zl(l) - z coordinate of level mesh nodes
C          drdxl(i,l) - 1/rhoref*dr/dx 
C          drdyl(i,l) - 1/rhoref*dr/dy 
C
C outputs: rxl(i,l) - integral presented above
C          ryl(i,l) - integral presented above
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           JANUARY 1, 1993
C-----------------------------------------------------------------------
C
C Global arrays
C
      real zl(nlevels),drdxl(ndim,nlevels),drdyl(ndim,nlevels)
      real rxl(ndim,nlevels),ryl(ndim,nlevels)
C
C Compute rxl and ryl by direct integration on the level mesh
C
      do 10 i=1,nn
         rxl(i,nlev)=0.0
         ryl(i,nlev)=0.0
         do 10 l=nlev-1,1,-1
            rxl(i,l)=rxl(i,l+1)
     &               -g*.5*(drdxl(i,l+1)+drdxl(i,l))*(zl(l+1)-zl(l))
            ryl(i,l)=ryl(i,l+1)
     &               -g*.5*(drdyl(i,l+1)+drdyl(i,l))*(zl(l+1)-zl(l))
 10   continue
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE LRHO2SRHO(ndim,nnvdim,nlevels,g,nn,nnv,nlev,
     &zl,drdxl,drdyl,rxl,ryl,zs,rx,ry)
C-----------------------------------------------------------------------
C purpose: This routine computes sigma mesh values of the following
C            integrals from level mesh quantities:
C            rx(i,j)=rxl(i,n2)-g*(vertical integral of drdxl(i,z) 
C                          from z=zs(j) to z=zl(n2)) 
C            ry(i,j)=ryl(i,n2)-g*(vertical integral of drdyl(i,z) 
C                          from z=zs(j) to z=zl(n2)) 
C            n2 - level mesh node imediately above z(j)
C
C index range: i=1,nn, j=1,nnv, l=1,nlev
C
C inputs:  ndim - parametered size of arrays indexed with i
C          nnvdim - parametered size of arrays indexed with j
C          nlevels - parametered size of arrays indexed with l
C          g - gravitational acceleration (MKS units) = 9.806 m/sec
C          nn - number of horizontal nodes
C          nnv - number of vertical nodes in sigma mesh
C          nlev - number of level surfaces in .lst file
C          zl(l) - z coordinate of level mesh nodes
C          drdxl(i,l) - 1/rhoref*dr/dx 
C          drdyl(i,l) - 1/rhoref*dr/dy 
C          rxl(i,l) - integral of drdxl on level mesh
C          ryl(i,l) - integral of drdyl on level mesh
C          zs(i,j) - z coordinate of sigma mesh nodes
C
C outputs: rx(i,j) - integral presented above
C          ry(i,j) - integral presented above
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           JANUARY 1, 1993
C-----------------------------------------------------------------------
C
C global arrays
C
      real zl(nlevels),
     &     drdxl(ndim,nlevels),drdyl(ndim,nlevels),
     &     rxl(ndim,nlevels),ryl(ndim,nlevels)
      real zs(ndim,nnvdim)    
      real rx(ndim,nnvdim),ry(ndim,nnvdim)
C
C Begin node loop for interpolation
C
      do 20 i=1,nn
C
C  Loop over sigma mesh nodes and perform interpolation
C
         do 20 j=1,nnv
            CALL BASIS1D(nlev,zl,zs(i,j),b2,n2,b1,n1)
C
C  If n1=0, then zs is below the minimum zl
C
            if(n1.EQ.0)then
               rx(i,j)=rxl(i,1)
               ry(i,j)=ryl(i,1)
C
C  If n2=0, then zs is above the maximum zl
C
            else if(n2.EQ.0)then
               rx(i,j)=0.0
               ry(i,j)=0.0
            else
               drdxzij=b1*drdxl(i,n1)+b2*drdxl(i,n2)
               drdyzij=b1*drdyl(i,n1)+b2*drdyl(i,n2)
               rx(i,j)=rxl(i,n2)
     &               -g*.5*(drdxl(i,n2)+drdxzij)*(zl(n2)-zs(i,j))
               ry(i,j)=ryl(i,n2)
     &               -g*.5*(drdyl(i,n2)+drdyzij)*(zl(n2)-zs(i,j))
            endif
 20   continue
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE NOD2EL(nndim,nedim,nnvdim,nlevels,nn,ne,nnv,rhoref,
     &x,y,in,zs,nlev,sigtl,xe,ye,zse,drdxle,drdyle)
C-----------------------------------------------------------------------
C purpose: This routine computes level and mesh element values from
C            nodal quantities.
C
C index range: i=1,nn, j=1,nnv, l=1,nlev, k=1,ne
C
C inputs:  ndim - parametered size of arrays indexed with i
C          nedim - parametered size of arrays indexed with k
C          nnvdim - parametered size of arrays indexed with j
C          nlevels - parametered size of arrays indexed with l
C          nn - number of horizontal nodes
C          ne - number of horizontal elements
C          nnv - number of vertical nodes in sigma mesh
C          rhoref - reference density (kg/m**3)
C          x(i),y(i) - cartesian nodal coordinates
C          in - incidence list
C          zs(i,j) - z coordinate of sigma mesh nodes
C          nlev - number of level surfaces in .lst file
C          sigtl(i,l) - level mesh sigmat data
C
C outputs: xe(k),ye(k) - element centroid coordinates
C          zse(k,j) - z coordinate of sigma mesh element centroids
C          drdxle(k,l) - element values of 1/rhoref*dr/dx 
C          drdyle(k,l) - element values of 1/rhoref*dr/dy 
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           JANUARY 1, 1993
C-----------------------------------------------------------------------
C
C Nodal arrays sent to routine from calling routine
C
      integer in(3,nedim)
      real x(nndim),y(nndim)
      real zs(nndim,nnvdim)
      real sigtl(nndim,nlevels)
C
C Elemental arrays returned to calling routine
C
      real xe(nedim),ye(nedim)
      real zse(nedim,nnvdim)
      real drdxle(nedim,nlevels),drdyle(nedim,nlevels)
C
C Element loop to compute sigma mesh element quantities
C
      do 10 k=1,ne
         i1=in(1,k)
         i2=in(2,k)
         i3=in(3,k)
         xe(k)=(x(i1)+x(i2)+x(i3))/3.0
         ye(k)=(y(i1)+y(i2)+y(i3))/3.0
         do 10 j=1,nnv
             zse(k,j)=(zs(i1,j)+zs(i2,j)+zs(i3,j))/3.0
 10   continue
C
C Element loop to compute level mesh element gradients
C
      do 20 l=1,nlev
         do 20 k=1,ne
            i1=in(1,k)
            i2=in(2,k)
            i3=in(3,k)
            dx1=x(i2)-x(i3)
            dx2=x(i3)-x(i1)
            dx3=x(i1)-x(i2)
            dy1=y(i2)-y(i3)
            dy2=y(i3)-y(i1)
            dy3=y(i1)-y(i2)
            arr=0.5*(x(i1)*dy1 + x(i2)*dy2 + x(i3)*dy3)
          drdxle(k,l)=(sigtl(i1,l)*dy1+sigtl(i2,l)*dy2+sigtl(i3,l)*dy3)
     &                /(2.*arr*rhoref)
          drdyle(k,l)=-(sigtl(i1,l)*dx1+sigtl(i2,l)*dx2+sigtl(i3,l)*dx3)
     &                /(2.*arr*rhoref)
 20   continue
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE SRHO2SHRBAR(ndim,nnvdim,nn,nnv,zs,rx,ry,hrbarx,hrbary)
C-----------------------------------------------------------------------
C purpose: This routine performs sigma mesh integrations to determine:
C            hrbarx(i)=vertical integral of rx(i,j) 
C                        from z=z(i,1) to z=z(i,nnv)=0.0
C            hrbary(i)=vertical integral of ry(i,j) 
C                        from z=z(i,1) to z=z(i,nnv)=0.0
C
C index range: i=1,nn, j=1,nnv
C
C inputs:  ndim - parametered size of arrays indexed with i
C          nnvdim - parametered size of arrays indexed with j
C          nn - number of horizontal nodes
C          nnv - number of vertical nodes for sigma mesh
C          zs(i,j) - z coordinate of sigma mesh nodes
C          rx(i,j) - quantity to be integrated 
C          ry(i,j) - quantity to be integrated 
C
C outputs: hrbarx(i) - integral presented above
C          hrbary(i) - integral presented above
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           JANUARY 1, 1993
C-----------------------------------------------------------------------
C
C Global arrays
C
      real rx(ndim,nnvdim),ry(ndim,nnvdim)
      real zs(ndim,nnvdim)    
      real hrbarx(ndim),hrbary(ndim)
C
C   Compute H*RBAR by integrating on sigma mesh
C
      do 20 i=1,nn
         hrbx=0.0
         hrby=0.0
         do 10 j=2,nnv            
            hrbx=hrbx+0.5*(rx(i,j)+rx(i,j-1))*(zs(i,j)-zs(i,j-1))
            hrby=hrby+0.5*(ry(i,j)+ry(i,j-1))*(zs(i,j)-zs(i,j-1))
 10      continue
         hrbarx(i)=hrbx
         hrbary(i)=hrby
 20   continue
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE BASIS1D(nn,zref,zdes,b2,n2,b1,n1)
C-----------------------------------------------------------------------
C purpose: This subroutine evaluates the nodes in a reference array 
C           which bracket a desired location and determines the values 
C           of their 1-D linear basis functions at that location.
C
C index range: k=1,nn
C
C inputs:  nn - number of entries in the reference 1-D reference array
C          zref(k) - node locations in 1-D reference array
C               NOTE: zref must increase from zref(1) to zref(nn)
C                     => element 1: head=zref(2),tail=zref(1)
C                        element 2: head=zref(3),tail=zref(2)
C                        .
C                        .
C                        element nn-1: head =zref(nn),tail=zref(nn-1)
C          zdes - the z value at which the basis functions are desired
C
C outputs: n1 - reference node immediately below zdes
C          b1 - value of n1's basis function at zres
C          n2 - reference node immediately above zdes
C          b2 - value of n2's basis function at zres
C
C notes:   zdes<zref(1) => n1=0,b1=0.0,n2=1,b2=1.0 
C          zdes>zref(nn) => n1=nn,b1=1.0,n2=0,b2=0.0 
C
C history: Written by Christopher E. Naimie
C          Dartmouth College
C          26 AUGUST 1992
C          REV1 - 1 JANUARY 1993/CEN
C-----------------------------------------------------------------------
C
C ARGUMENTS
C
      dimension zref(nn)
C
C START OF EXECUTABLE CODE
C
      do 10 k=1,nn 
         if(zdes.LT.zref(k))then
C
C if zdes<zref(1), zdes is below the reference data => set b1=n1=0
C
            n2=k
            n1=k-1
            if(n2.EQ.1)then
               b1=0.0
               b2=1.0
               go to 20
            endif
            zdiff=zref(n2)-zref(n1)
            b1=(zref(n2)-zdes)/zdiff
            b2=(zdes-zref(n1))/zdiff  
            go to 20
         endif
 10   continue
C
C if zdes>zref(nn), zdes is above the reference data => set b2=n2=0
C
      b1=1.0
      n1=nn
      b2=0.0
      n2=0.0
 20   continue
C
C END OF ROUTINE
C
      return
      end
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE SPRSBLD1(IN,NDIM,NNE,IQ,JQ,NN,NE,MN,NDF)
C
C  NOVEMBER 30, 1990
C  D.LYNCH REVISION OF K.PAULSEN REVISION OF WEIPING VERSION.
C  SIMPLE ADAPTATION OF SPRSBLDW1 BY REVERSING ORDER OF IN(,) ARGUMENTS.
C
C  BUILDS SPRSPAK ARRAYS IQ,JQ FROM ELEMENT INCIDENCE LIST IN(K,L)
C     WITH L=GLOBAL ELEMENT NUMBER, K=LOCAL NODE NUMBER
C     FOR NDF DEGREES OF FREEDOM AT EACH NODE -- i.e. 
C     NDF EQUATIONS AND NDF UNKNOWNS.
C     NDF IS THE SAME AT EACH NODE.
C
C  SAME AS SPRSBLDW1 EXCEPT FOR THE ORDER OF THE ARGUMENTS FOR THE IN 
C     ARRAY
C 
C  NDIM IS THE FIRST PHYSICAL DIMENSION OF IN; NDIM.GE.NNE
C  NNE INDICATES THAT IN(I,L) FOR I=1 THRU NNE ARE NODE NUMBERS; 
C     THE BALANCE (IF ANY) OF IN WILL BE IGNORED.  
C  IQ(I) IS THE ADDRESS IN QV OF THE LAST NONZERO ENTRY IN ROW I OF Q
C  JQ(K) IS THE COLUMN INDEX IN Q FOR QV(K)
C  NN IS THE ACTUAL NUMBER OF NODES
C  NE IS THE ACTUAL NUMBER OF ELEMENTS
C  MN EQUALS OR EXCEEDS THE MAXIMUM NUMBER OF NONZREO ENTRYS PER ROW, 
C     FOR A SCALAR (i.e. NDF=1) PROBLEM.
C  NDF IS THE DEGREES OF FREEDOM, ASSUMED THE SAME AT EACH NODE.
C
C  IT IS ASSUMED THAT IN MAIN PROGRAM, DIMENSIONS ARE AT LEAST 
C     IQ(NN*NDF), JQ((NN*NDF)*(MN*NDF)).  NOTE THAT JQ IS OVERSIZED; 
C     BUT THE ALGORITHM MAY USE THIS EXTRA SIZE.  
C
C  CALLS ARE MADE TO SPRSINSRTD AND SPRSBLD3DD
C
      DIMENSION IN(NDIM,*),IQ(*),JQ(*)
C
C  INITIALIZE JQ
C
       DO 5 I=1,NN*MN
       JQ(I)=0
    5 CONTINUE
C
C  INITIALIZE JQ AS DIAGONAL MATRIX
C
      DO 10 I=1,NN
      J=(I-1)*MN+1
      JQ(J)=I
   10 CONTINUE
C
C  ELEMENT LOOP FOR JQ, ASSUMING NDF=1.
C  ALL ROWS ALLOWED MN ENTRIES.
C  TRAILING ZEROS WILL BE REMOVED BELOW.
C
      DO 100 L=1,NE
      DO 90 I=1,NNE
      II=IN(I,L)
      DO 90 J=1,NNE
      JJ=IN(J,L)
      IF(I.NE.J) CALL SPRSINSRTD(II,JJ,JQ,MN)
   90 CONTINUE
  100 CONTINUE
C
C FIND IQ FROM JQ
C
       K=0
       KK=0
       DO 110 I=1,NN
       DO 120 J=1,MN
       KK=KK+1
       IF(JQ(KK).NE.0) K=K+1
  120 CONTINUE
       IQ(I)=K
  110 CONTINUE      
C
C REARRANGE (CONDENSE) THE JQ -- REMOVE THE TRAILING ZEROS
C
       DO 140 I=2,NN
        KMIN=IQ(I-1)+1
        KMAX=IQ(I)
        KJMIN=(I-1)*MN+1
        DO 130 K=KMIN,KMAX
          JQ(K)=JQ(KJMIN)
          KJMIN=KJMIN+1
  130 CONTINUE
  140 CONTINUE
C
C IQ, JQ NOW CORRECT FOR NDF=1
C EXPAND FOR NDF.GT.1
C
        IF (NDF.GT.1) THEN
	CALL SPRSBLD3DD(IQ,JQ,NN,NDF)
	ELSE
	 GO TO 200
        END IF
  200 RETURN
      END
C
C***********************************************************************
C****************************************************************SPRSPAK
C
      SUBROUTINE SPRSINSRTD(I,J,JQ,MN)
C
C  INSERTS GLOBAL ADDRESS (I,J) INTO JQ(), 
C  ASSUMING EXACTLY MN ENTRIES IN EACH ROW.
C
      DIMENSION JQ(*)
      KMIN=(I-1)*MN+1
      KMAX=I*MN
C
C  SEARCH FOR INSERTION POINT KK
C
      DO 10 K=KMIN,KMAX
      KK=K
      IF(J.LT.JQ(K)) GO TO 20
      IF(J.EQ.JQ(K)) GO TO 50
      IF(JQ(K).EQ.0) GO TO 20
   10 CONTINUE
C
C  INSERT J AT JQ(KK); BUMP JQ
C
C CHECK INSERTED BY K PAULSEN 11/29...NOT VALIDATED AT TIME OF
C INSERTION !!!!
C
C Check commented out for FUNDY5, since the main routine does this.
C
c   20 IF(JQ(KMAX).NE.0) THEN
c       WRITE(*,*) 'MAX CONNECTIVITY MN=',MN,
c     *  'VIOLATED AT I,J=',I,J,'IN SPRSINSRTD'
c       STOP
c      END IF
   20 DO 30 K=KMAX-1,KK,-1
   30 JQ(K+1)=JQ(K)
      JQ(KK)=J   
C
C  DONE
C
   50 CONTINUE
      RETURN     
      END
C
C***********************************************************************
C****************************************************************SPRSPAK
C
       SUBROUTINE SPRSBLD3DD(IQ,JQ,NN,NDF)
C
C BUILD 3D IDEX IQQ AND JQQ
C
       DIMENSION IQ(*),JQ(*)
	DO 60 I=NN,1,-1
	 IF (I.NE.1) THEN
	  II=IQ(I)-IQ(I-1)
	 ELSE
	  II=IQ(I)
	 END IF
	 DO 10 J=NDF,1,-1
	  IQ(NDF*(I-1)+J)=IQ(I)*NDF*NDF-(NDF-J)*NDF*II
   10	 CONTINUE
	 DO 30 K=1,II
	   JK=IQ(NDF*I)+1-(K-1)*NDF
	   IF (I.NE.1) THEN
	    JJ=JQ(IQ(I)-K+1)*NDF+1
	   ELSE
	    JJ=JQ(IQ(I)/NDF-K+1)*NDF+1
	   END IF
	  DO 20 J=1,NDF
	   JK=JK-1
	   JJ=JJ-1
C
C THE FOLLOWING 3 LINES FROM WEIPING'S CODE
C	   JQ(JK)=JJ
C	   JQ(JK-NDF*II)=JJ
C	   IF(NDF.EQ.3)	   JQ(JK-6*II)=JJ
C WORKED FOR NDF=2 OR 3.  REPLACED BY INDUCTION  
C WITH THE DO 15 LOOP BELOW (DRL 1 DECEMBER 1990).   
C IT IS NOT OBVIOUS THAT JQ REMAINS UNCORRUPTED IN THE 
C OVERWRITING.  
C
           DO 15 L=1,NDF
           JQ(JK-(L-1)*NDF*II)=JJ
   15      CONTINUE
   20	  CONTINUE
   30	 CONTINUE
   60	CONTINUE
        RETURN
        END
C
C
C***********************************************************************
C****************************************************************SPRSPAK
C
      SUBROUTINE SPRSMLTRC(QV,IQ,JQ,X,B,NN)
C
C  MULTIPLIES QV*X; ANSWER RETURNS IN B
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C  VERSION RC: REAL QV, COMPLEX X
C
      DIMENSION QV(1),IQ(1),JQ(1)
      COMPLEX X(1),B(1)
      COMPLEX SUM
      KMAX=0
      DO 10 I=1,NN
      SUM=(0.,0.)
      KMIN=KMAX+1
      KMAX=IQ(I)
      DO 20 K=KMIN,KMAX
   20 SUM=SUM+QV(K)*X(JQ(K))
   10 B(I)=SUM
      RETURN
      END
C
C***********************************************************************
C****************************************************************SPRSPAK
C
      SUBROUTINE SPRSINVMLT(QV,IQ,D,NN)
C
C  REPLACES QV WITH INV(D)*QV
C  D IS A DIAGONAL MATRIX, LENGTH NN
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C
      DIMENSION QV(1),D(1),IQ(1)
      KMAX=0
      DO 10 I=1,NN
      KMIN=KMAX+1
      KMAX=IQ(I)
      FTR=1./D(I)
      DO 20 K=KMIN,KMAX
   20 QV(K)=FTR*QV(K)
   10 CONTINUE
      RETURN
      END
C
C***********************************************************************
C****************************************************************SPRSPAK
C
      FUNCTION KAY(I,J,IQ,JQ)
C
C  FINDS ADDRESS K IN QV FOR Q(I,J)
C
      DIMENSION IQ(1),JQ(1)
      KAY=IQ(I)+1
   10 KAY=KAY-1
      IF(JQ(KAY).NE.J)GO TO 10
      RETURN
      END
C
C***********************************************************************
C****************************************************************SPRSPAK
C
      SUBROUTINE SPRSCONVRC(I,PPX,PPY,IQ,JQ,U,V,Z,
     *  DUDX,DUDY,DVDX,DVDY,DZDX,DZDY,
     *  NNV,NNDIM)
C
C  VERSION RC: REAL PPX,PPY; COMPLEX U,V; REAL Z.
C  COMPUTES THE HORIZONTAL DERIVATIVES AT DEPTH, UNDER HORIZONTAL NODE I
C  PPX, PPY ARE THE SPRSPAK MATRICES FOR D/DX, D/DY
C  IQ,JQ ARE THE STANDARD POINTERS FOR PPX, PPY
C  U(I,J),V(I,J) ARE THE X- AND Y-VELOCITIES
C      AT HORIZONTAL NODE I, VERTICAL NODE J
C  Z(I,J) IS THE VERTICAL POSITION
C      OF HORIZONTAL NODE I, VERTICAL NODE J
C  THE HORIZONTAL DERIVATIVE TERMS DUDX(J), ETC. ARE RETURNED
C      FOR ALL NODES J UNDER THE GIVEN NODE I
C  NNV IS THE ACTUAL # OF NODES IN THE VERTICAL
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF U AND V
C  THIS ROUTINE IS NONDESTRUCTIVE
C
      DIMENSION PPX(1),PPY(1),IQ(1),JQ(1),Z(NNDIM,1),DZDX(1),DZDY(1)
      COMPLEX U(NNDIM,1),V(NNDIM,1),DUDX(1),DUDY(1),DVDX(1),DVDY(1)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO 20 J=1,NNV
      DUDX(J)=(0.,0.)
      DUDY(J)=(0.,0.)
      DVDX(J)=(0.,0.)
      DVDY(J)=(0.,0.)
      DZDX(J)=0.
      DZDY(J)=0.
      DO 10 K=KMIN,KMAX
      DUDX(J)=DUDX(J)+PPX(K)*U(JQ(K),J)
      DUDY(J)=DUDY(J)+PPY(K)*U(JQ(K),J)
      DVDX(J)=DVDX(J)+PPX(K)*V(JQ(K),J)
      DVDY(J)=DVDY(J)+PPY(K)*V(JQ(K),J)
      DZDX(J)=DZDX(J)+PPX(K)*Z(JQ(K),J)
      DZDY(J)=DZDY(J)+PPY(K)*Z(JQ(K),J)
   10 CONTINUE
   20 CONTINUE
      RETURN
      END
C
C***********************************************************************
C****************************************************************SPRSPAK
C
      SUBROUTINE SPRSGRADRC(I,PPX,PPY,IQ,JQ,F,Z,
     *  DFDX,DFDY,DZDX,DZDY,
     *  NNV,NNDIM)
C 
C  A SIMPLE VARIANT OF SPRSCONVRC
C  VERSION RC: REAL PPX,PPY; COMPLEX F; REAL Z.
C  COMPUTES THE HORIZONTAL GRADIENT OF F AT DEPTH, UNDER HORIZONTAL NODE #I.
C  PPX, PPY ARE THE SPRSPAK MATRICES FOR D/DX, D/DY
C  IQ,JQ ARE THE STANDARD POINTERS FOR PPX, PPY
C  F(I,J) IS A COMPLEX SCALAR FUNCTION DEFINED 
C      AT HORIZONTAL NODE I, VERTICAL NODE J
C  Z(I,J) IS THE VERTICAL POSITION
C      OF HORIZONTAL NODE I, VERTICAL NODE J
C  THE HORIZONTAL DERIVATIVE TERMS DFDX(J), DFDY(J) ARE RETURNED
C      FOR ALL NODES J UNDER THE GIVEN NODE I
C  NNV IS THE ACTUAL # OF NODES IN THE VERTICAL
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF F
C  THIS ROUTINE IS NONDESTRUCTIVE
C
      DIMENSION PPX(1),PPY(1),IQ(1),JQ(1),Z(NNDIM,1),DZDX(1),DZDY(1)
      COMPLEX F(NNDIM,1),DFDX(1),DFDY(1)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO 20 J=1,NNV
      DFDX(J)=(0.,0.)
      DFDY(J)=(0.,0.)
      DZDX(J)=0.
      DZDY(J)=0.
      DO 10 K=KMIN,KMAX
      DFDX(J)=DFDX(J)+PPX(K)*F(JQ(K),J)
      DFDY(J)=DFDY(J)+PPY(K)*F(JQ(K),J)
      DZDX(J)=DZDX(J)+PPX(K)*Z(JQ(K),J)
      DZDY(J)=DZDY(J)+PPY(K)*Z(JQ(K),J)
   10 CONTINUE
   20 CONTINUE
      RETURN
      END
C
C***********************************************************************
C****************************************************************SPRSPAK
C
      SUBROUTINE SPRSNORM(QV,IQ,NN)
C
C  CORRECTS QV SO THAT ROW SUMS ARE ZERO
C  THE CORRECTION IS SPREAD EVENLY OVER THE WHOLE ROW
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C  IQ(I) IS THE ADDRESS IN QV OF THE LAST NONZERO ENTRY IN ROW I OF Q
C
      DIMENSION QV(1),IQ(1)
      KMAX=0
      DO 10 I=1,NN
      KMIN=KMAX+1
      KMAX=IQ(I)
      SUM=0.0
      DO 20 K=KMIN,KMAX
   20 SUM=SUM+QV(K)
      CORR=SUM/(1.*(KMAX-KMIN+1))
      DO 25 K=KMIN,KMAX
   25 QV(K)=QV(K)-CORR
   10 CONTINUE
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
C f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5xC
C f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5xC
C
C Fixed subroutines added for f5x follow
C Fixed subroutines added for f5x follow
C
C f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5xC
C f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5x f5xC
C
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE BASIS2D(nndim,nedim,in,x,y,nn,ne,xval,yval,
     &nclosest,nelem,basis)
C-----------------------------------------------------------------------
C purpose: This subroutine determines the element on which the point
C            (xval,yval) resides and determines the value of the basis
C            functions at the point.  In the first pass, only elements 
C            containing the closest node to (xval,yval) are tested to 
C            see if (xval,yval) lies on them.  If this fails, all
C            elements are tested to see if (xval,yval) lies on them.
C
C
C restrictions: Applicable only for  2-D triangular linear elements
C
C calling sequence: CALL CLOSESTNODE(nndim,x,y,nn,xval,yval,nclosest)
C                   CALL POINTONEL(xlocal,ylocal,xval,yval,nflag)
C                   CALL PHI2D1(xlocal,ylocal,xval,yval,basis)
C
C inputs:  nn is the number of nodes in the mesh
C          x(nn),y(nn) 1-D arrays containing the are the x and y
C            coordinates of each nodes
C          ne is the number of elements in the mesh
C          in(3,ne) is the LABLIB Standard incidence list
C          xval,yval is the location of the point 
C          nedim is the size of the arrays in the calling routine
C
C outputs: basis is a 1-D array containing the value of the 3 basis
C            functions at the point xval,yval.
C          nelem is the element on which (xval,yval),resides
C          nclosest is the closest node to (xval,yval)
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           1992 APRIL 22 - Original Version
C           1997 MARCH 01 - Revised to call PHI2D1 instead of PHI2D to
C                           determine basis functions.  PHI2D1 is an
C                           improvement over PHI2D in that it less
C                           subject to round-off errors.
C-----------------------------------------------------------------------
C
C ARGUMENTS
      integer in(3,nedim)
      real x(nndim),y(nndim),basis(3)
C
C LOCAL VARIABLES
      dimension xlocal(3),ylocal(3)
C
C START OF EXECUTABLE CODE
C
C  Find the closest node to (xval,yval)
      CALL CLOSESTNODE(nndim,x,y,nn,xval,yval,nclosest)
C
C  Find if an element containing the closest node also contains the
C    point (xval,yval)
C
      do 10 i=1,ne
         n1=in(1,i)
         n2=in(2,i)
         n3=in(3,i)
         if(nclosest.EQ.n1.OR.nclosest.EQ.n2.OR.nclosest.EQ.n3)then
            xlocal(1)=x(n1)
            ylocal(1)=y(n1)
            xlocal(2)=x(n2)
            ylocal(2)=y(n2)
            xlocal(3)=x(n3)
            ylocal(3)=y(n3)
            CALL POINTONEL(xlocal,ylocal,xval,yval,nflag)
            if(nflag.EQ.1) then
               nelem=i
               CALL PHI2D1(xlocal,ylocal,xval,yval,basis)
               return
            endif
         endif
 10   continue
C
C  Find if any element contains the point (xval,yval)
C
      do 20 i=1,ne
         xlocal(1)=x(in(1,i))
         ylocal(1)=y(in(1,i))
         xlocal(2)=x(in(2,i))
         ylocal(2)=y(in(2,i))
         xlocal(3)=x(in(3,i))
         ylocal(3)=y(in(3,i))
         CALL POINTONEL(xlocal,ylocal,xval,yval,nflag)
         if(nflag.EQ.1) then
            nelem=i
            CALL PHI2D1(xlocal,ylocal,xval,yval,basis)
            return
         endif
 20   continue
      nelem=0
C
C END OF ROUTINE
C
      return
      end
C****67***************************************************************72
C***********************************************************************
      SUBROUTINE CLOSESTNODE(nndim,x,y,nn,xval,yval,nclosest)
C-----------------------------------------------------------------------
C purpose: This subroutine finds the node closest to a point specified
C            by the coordinates (xval,yval)
C  
C restrictions:  all nodes must lie in a flat plane
C
C inputs:  x(nn) and y(nn) are 1-D arrays containing the x and y
C            coordinates of each node
C          nn is the number of nodes
C          xval,yval are the point coordinates
C
C outputs: nclosest is the node who's (x,y) coordinates are closest to
C            (xval,yval)
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension x(nndim),y(nndim)
C
C START OF EXECUTABLE CODE
C
C  Find node closest to (xval,yval)
      dclosest=SQRT((x(1)-xval)**2.0+(y(1)-yval)**2.0)
      nclosest=1
      do 25 j=2,nn
         dval=SQRT((x(j)-xval)**2.0+(y(j)-yval)**2.0)
         if(dval.LT.dclosest)then
            dclosest=dval
            nclosest=j
         endif
 25   continue
C
C END OF ROUTINE
      return
      end
C****67***************************************************************72
C***********************************************************************
      SUBROUTINE POINTONEL(xlocal,ylocal,xval,yval,nflag)
C-----------------------------------------------------------------------
C purpose: This Subroutine determines if a point (xval,yval) is on a 
C            local 2-D triangular linear element
C    
C restrictions: Applicable only for  2-D triangular linear element
C
C inputs:  xlocal(3) and y(3) are 1-D arrays containing the x and y 
C            coordinates of each local node on this element.
C            NOTE:  The node locations must be entered in a
C                     counterclockwise order (LABLIB standard)
C          xval,yval is the location of the point 
C
C outputs: nflag is set equal to 1 if the point (xval,yval) is on the 
C            local element.
C          nflag is set equal to 0 if the point (xval,yval) is not on 
C            the local element.
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension xlocal(3),ylocal(3)
C
C LOCAL VARIABLES
      dimension vx(3),vy(3)
C
C  Calculate the x and y components of vectors pointing from (xval,yval)
C    to each node on the element
      do 10 i=1,3
         delx=xlocal(i)-xval
         dely=ylocal(i)-yval
         d=SQRT(delx**2.0+dely**2.0)
         theta=ATAN2(dely,delx)
         vx(i)=d*COS(theta)
         vy(i)=d*SIN(theta)
 10   continue
C
C  Determine if the point is on the element by calculating the
C    crossproducts of neighboring vectors in a direction
C    which yeild all positive numbers if x,y is on the element.
C    (ie: v1Xv2>0, v2Xv3>0, and v3Xv1>0 => x,y is on element)
      do 20 i=1,3
         j=i+1
         if(i.EQ.3)j=1
         crossprod=vx(i)*vy(j)-vy(i)*vx(j)
         if(crossprod.LT.0.0)then
            nflag=0
            go to 21
         endif
 20   continue
      nflag=1
 21   continue
C
C END OF ROUTINE
      return
      end
C****67***************************************************************72
C***********************************************************************
      SUBROUTINE PHI2D1(xl,yl,xvalu,yvalu,phi)
C-----------------------------------------------------------------------
C purpose: This subroutine evaluates the value of the linear basis
C            on a 2-D triangular linear element at a 
C            point (xvalu,yvalu) on a local element.
C
C restrictions: Applicable only for  2-D triangular linear element
C               Assumes that the fact that (xval,yval) lies on this 
C                 element has been determined a priori.
C               Attempts to eliminate round-off problems.
C               Minimum nodal area value is 0.0 => this subroutine
C                 can not be used to determine if the point in question
C                 is on an element
C
C inputs:  xl(3) and yl(3) are 1-D arrays containing the x and y
C            coordinates of each local node on this element.
C            NOTE:  The node locations must be entered in a
C                     counterclockwise order (LABLIB standard)
C          xvalu,yvalu is the location of the point 
C
C outputs: phi is a 1-D array containing the value of the 3 basis
C            functions at the point xval,yval.
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           1992 APRIL 22 - Original Version
C           1997 MARCH 01 - PHI2D1 - Revised version of PHI2D.
C                           Revised to reduce potential round-off errors
C                           which resulted in basis functions which
C                           were less than 0.0 and greater than 1.0 on
C                           some machines (this potential problem never
C                           revealed itself until early 1997).
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension xl(3),yl(3),phi(3)
C
C START OF EXECUTABLE CODE
C
C Calculate twice the nodal areas
      A12=xvalu*(yl(2)-yl(3))+xl(2)*(yl(3)-yvalu)+xl(3)*(yvalu-yl(2))
      A22=xl(1)*(yvalu-yl(3))+xvalu*(yl(3)-yl(1))+xl(3)*(yl(1)-yvalu)
      A32=xl(1)*(yl(2)-yvalu)+xl(2)*(yvalu-yl(1))+xvalu*(yl(1)-yl(2))
C
C Do not allow for negative areas (which only occur due to round-off
C errors)
      A12=AMAX1(A12,0.0)
      A22=AMAX1(A22,0.0)
      A32=AMAX1(A32,0.0)
C
C Sum twice the nodal areas to get twice the total elemental area and 
C then evaluate basis functions
      AELE2=A12+A22+A32
      phi(1)=A12/AELE2      
      phi(2)=A22/AELE2     
      phi(3)=A32/AELE2      
C
C END OF ROUTINE
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE BASIS1D1(nn,zref,zdes,b2,n2,b1,n1)
C-----------------------------------------------------------------------
C purpose: This subroutine evaluates the value of the basis functions 
C            alive on a 1-D linear element at a point
C
C restrictions: Applicable only for  1-D linear elements
C
C inputs:  nn is the number of entries in the reference z array
C          zref(nn) is the 1-D array containing the reference z array
C               NOTE: zref must increase from zref(1) to zref(nn)
C          zdes is the z value at which the basis functions are desired
C
C outputs: n1 is the reference node immediately below zdes
C          b1 is the value of the n1's basis function at zdes
C          n2 is the reference node immediately above zdes
C          b2 is the value of the n2's basis function at zdes
C
C notes:   zdes<zref(1)  => n1=nn,b1=0.0,n2=1,b2=1.0 
C          zdes>zref(nn) => n1=nn,b1=1.0,n2=1,b2=0.0 
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           26 AUGUST 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension zref(nn)
C
C START OF EXECUTABLE CODE
      do 10 k=1,nn 
         if(zdes.LT.zref(k))then
C
C if zdes<zref(1), zdes is below the reference data 
C                    => n1=nn,b1=0.0,n2=1,b2=1.0 
            n2=k
            n1=k-1
            if(n2.EQ.1)then
                n1=nn
                b1=0.0
                b2=1.0
                go to 20
             endif
             zdiff=zref(n2)-zref(n1)
             b1=(zref(n2)-zdes)/zdiff
             b2=(zdes-zref(n1))/zdiff  
             go to 20
         endif
 10   continue
C
C if zdes>zref(nn), zdes is above the reference data
C                     => n1=nn,b1=1.0,n2=1,b2=0.0
      b1=1.0
      n1=nn
      b2=0.0
      n2=1
 20   continue
C
C END OF ROUTINE
      return
      end
C
C***********************************************************************
C*****************************************************************FUND
C
      SUBROUTINE GetFreqI(NF, NF1, NF2, Freq, Infer, Amprat, NFDIM)
      INTEGER NF, NF1, NF2, Infer(*), NFDIM
      REAL    FREQ(NFDIM)
      COMPLEX Amprat(*)
      INTEGER IFREQ
c
      write(*,*)'How many frequencies (max = ',NFDIM,' )?'
      read(*,*)NF
      if(NF.gt.NFDIM) then 
	 write(*,*)'Too many frequencies.  NF, NFDIM ',NF,NFDIM
	 write(*,*)'I quit'
	 write(2,*)'Too many frequencies.  NF, NFDIM ',NF,NFDIM
	 write(2,*)'I quit'
	 close(2)
	 stop 
      endif
c      
      write(*,*)'How many primary frequencies (max = ',NF,' )?'
      read(*,*) NF1
      if(NF1.gt.NF) then 
	 write(*,*)'Too many frequencies.  NF1, NF ',NF1, NF
	 write(*,*)'I quit'
	 write(2,*)'Too many frequencies.  NF1, NF ',NF1, NF
	 write(2,*)'I quit'
	 close(2)
	 stop 
      endif
c      
      do 2000 ifreq = 1, NF1 
	   write(*,*)'Enter frequency ',ifreq,' (radians/sec)' 
	   read(*,*)freq(ifreq)
	   Infer(ifreq)=ifreq
	   Amprat(ifreq)=cmplx(1.0, 0.0)
 2000 continue
c
      nf2=nf-nf1 
      if(nf2.gt.0) then
        do ifreq=nf1+1, nf
         write(*,*)'Enter frequency ',ifreq,' (radians/sec), '  
         write(*,*)'inference frequency # from above (integer), '
         write(*,*)'and amplitude ratio(complex).'
         read(*,*) freq(ifreq), Infer(ifreq), Amprat(ifreq)
        enddo
      endif
c
      return
      end
C
C***********************************************************************
C*****************************************************************FUNDY4
C
      Subroutine BuildW(Wls,nobs,nbc,nfreq,w0,w1,
     & LRVS,x,y)
C
C  nobs2,  nbc2, nfreq
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Arrays passed:
C
      Real Wls(nbc2nfdim, nbc2nfdim)
      Real X(*), Y(*)
      Integer LRVS(*)
C
C Local Arrays:
C
      Real S(nbc2dim, nbc2dim), StS(nbc2dim, nbc2dim)
c      Real INBE(NBCdim,4) (this arrar was real when used with readbel.)
      Integer INBE(NBCdim,4)
      Real IBC(nndim)
C
      nobs2=2*nobs
      nbc2=2*nbc
C
C Read the .bel file
C
c      call ReadBel(NBE, INBE)
c      if(NBE.NE.NBC) then
c      write(2,*)'NBE read in ReadBel is different than NBC; ',
c     $          ' quitting'
c      close(2)
c      stop
c      endif
c
c June 10, 1997 drl:
c Instead of using redundant .bel file, use the LRVS list assuming
c nodes are ordered sequentially, counterclockwise.  The subroutine
c BNodeOrder is intended to guarantee this; as of 6/10/97 that subr. 
c is blank; relying on default ordering from .inp file.
c
      write(2,*)'Subr. BuildW: making boundary elements from the BC',
     &' node list LRVS.'
      write(2,*)'Sequential, CCW order is assumed.'
      nbe=nbc
      do L=1, nbc-1
      inbe(L,1)=lrvs(L)
      inbe(L,2)=lrvs(L+1)
      enddo
      inbe(nbc,1)=lrvs(nbc)
      inbe(nbc,2)=lrvs(1)
C
C Compute IBC(I), i.e. the inverse of LRVS
C IBC(i) is the boundary node number of global node number i. 
C IBC(i) = 0 if global node i is an interior node
C IBC(i) is between 1 and NBC if global node i is on the boundary
C
      do i=1,nndim
      IBC(i)=0
      enddo
      do i=1,nbc
      IBC(LRVS(i))=i
      enddo
C
C Check for flaws in IBC, INBE
C
c      write(2,*)'Subr. BuildW: Opening BuildW.out for .bel diagnosis'
c      open(unit=19,file='BuildW.out',status='unknown')
c      write(19,*) ' INBE:'
c      Do L=1,NBE
c      I1=INBE(L,1)
c      I2=INBE(L,2)
c      write(19,*) I1, I2
c      enddo
c      write(19,*) ' IBC:'
c      do I=1,nn
c      if(IBC(i).ne.0) then
c      write(19,*) IBC(i)
c      endif
c      enddo
c      write(19,*) ' LRVS:'
c      do I=1,nbc
c      write(19,*) LRVS(i)
c      enddo   
C      
C Initialize S and StS
C
      do i=1, nbc2
      do j=1, nbc2
      S(i,j)=0.0      
      StS(i,j)=0.0
      enddo
      enddo
C
C Make S, the slope matrix
C
      do L=1, nbe
      I1=IBC( INBE(L,1) )
      I2=IBC( INBE(L,2) )
      I11= INBE(L,1) 
      I22=  INBE(L,2) 
C
      if(i1.le.0.or.i2.le.0) then
c      write(19,*) 'error: element #', L
c      close(19)
      write(2,*) 'error: element #', L
      close(2)
      stop
      endif
      if(i1.gt.nbc.or.i2.gt.nbc) then     
c      write(19,*) 'error: element #', L
c      close(19)
      write(2,*) 'error: element #', L
      close(2)
      stop
      endif
C
cgh   ds= sqrt( (x(I2)-x(I1))**2 + (y(I2)-y(I1))**2 )
      ds= sqrt( (x(I22)-x(I11))**2 + (y(I22)-y(I11))**2 )
      SE=1./(sqrt(ds))
      s(2*L-1,2*I2-1) =  se
      s(2*L-1,2*I1-1) = -se
      s(2*L  ,2*I2  ) =  se
      s(2*L  ,2*I1  ) = -se
      enddo
C
C Compute StS (Wls for a single frequency)
C
      do i=1,nbc2
      do j=1,nbc2
         sum=0.
         do k=1,nbc2
           sum=sum+s(k,i)*s(k,j)
         enddo
         StS(i,j)=sum
      enddo
      enddo
C
C Assemble the multi-frequency version into Wls
C 
      do k=1, nfreq
         i0=(k-1)*nbc2
         j0=(k-1)*nbc2  
         do i=1,nbc2
         do j=1,nbc2
           Wls(i0+i,j0+j)=w1*StS(i,j)
         enddo
         enddo
      enddo 
      do i=1,nbc*2*nfreq
      Wls(i,i)=Wls(i,i)+w0
      enddo
C
C Done
C
c      close(19)        
      return
      end             
      
C
C***********************************************************************
C*****************************************************************FUNDY4
C
      Subroutine ReadBel(NBE, INBE)
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
      real INBE(NBCdim,*)
      character*120 filename
      write(*,*)'Enter name of .bel file'
      read(*,'(a)')filename
      open(unit=4, file=filename, status='old')
c      open(unit=4, file='bank150.bel', status='old') 
      i=1  
      read(4,*)
      read(4,*)     
  100 continue
      read(4,*,end=101) L,INBE(L,1), INBE(L,2),
     &                    INBE(L,3), INBE(L,4)             
      i=i+1
      go to 100
  101 continue
      NBE=i-1
      close(4)
      return
      end
C
C***********************************************************************
C*****************************************************************
C
      subroutine BuildFP1(FP, Nbc, Nmode, Lspec)
c
c Each column of FP is a Fourier mode sampled in j=space.
C The number of modes will be Lspec, beginning with k=-Nmode and working
c upward.  The intention is to build the spectrum
c    -Nmode <= k <= Nmode
c with length (width) Lspec = 2*Nmode +1. 
c
c The FP columns will be in the natural order of ascending wavenumber, 
c beginning with k=-Nmode and working up.
c
      INCLUDE 'TRUXTON6ARRAYS.DIM'
      complex FP(nbcdim, nbcdim), eye, f1, f2
      pi=acos(-1.)
      eye=cmplx(0.,1.)
      f1=eye*2.*pi/(1.*nbc)
c
c Check for Lspec out of bounds
c
      if(Lspec.gt.Nbc) then
      write(2,*) 'Lspec.gt.Nbc in subroutine BuildFP1; stopping'
      close(2)
      stop
      endif
c
      do j=1,Lspec
        k=j-(Nmode+1)
        do i=1, Nbc
          f2=k*i*f1
          FP(i,j)=cexp(f2)
        enddo
      enddo
      return
      end
C
C***********************************************************************
C*****************************************************************
C
      subroutine svd_se(se,Ccov,nbc,nf,zls,freq,
     &rcond,irank,ncol,RMSV)
      INCLUDE 'TRUXTON6ARRAYS.DIM'
      real se(*), zls(*), freq(*) 
      complex Ccov(NBCdim*NFdim,NBCdim*NFdim)
      write(2,*) ' '
      write(2,*)' SVD INVERTER: STANDARD ERROR REPORT'
      write(2,*) ' '
      write(2,68) RMSV
   68 format(/'NOTE:',/
     & 'All standard errors have been normalized ',/
     & 'by the RMS Velocity Error =', 
     & e12.4,//
     & 'The CVs are not normalized i.e. they are SE/Solution',/)
c
c      open(9,file='se_out1',status='unknown')
c      open(10,file='se_out2',status='unknown')
c      
      do IFREQ=1, nf
c      
        write(2,*) ' '
        write(2,13)
        write(2,*) 'Tidal Constituent #', IFREQ,
     &                     'Frequency', freq(ifreq)
        write(2,*) 'Node #,  Magnitude,  std. error,  CV'
        ioffset=(IFREQ-1)*2*NBC
        Noffset=(IFREQ-1)*NBC
        zmagsum=0.
        varsum=0.
        stderrsum=0.
        do ibc=1,nbc
          i1=ioffset+2*ibc-1
          i2=ioffset+2*ibc
          zmag  = sqrt( zls(i1)**2 + zls(i2)**2 )
c          stderr= sqrt(  se(i1)**2 +  se(i2)**2 )
          variance=cabs(Ccov(Noffset+ibc, Noffset+ibc))
          stderr=sqrt(variance)
          cv    = RMSV*stderr/zmag
          write(2,10) ibc, zmag, stderr, cv
c          write(9,10) ibc, zmag, stderr, cv
          zmagsum=zmagsum+zmag
          stderrsum=stderrsum+stderr
          varsum=varsum+variance
        enddo
        zmean=zmagsum/(1.*nbc)
        semean= stderrsum /(1.*nbc)
        varmean=varsum/(1.*nbc)          
        rootmeanvar=sqrt(varsum)
        cv=RMSV*semean/zmean
        write(2,*)' '
        write(2,*)'Mean Values:'
        write(2,*)'     Magnitude, std error, Root(meanvar), CV'
        write(2,12) zmean, semean, rootmeanvar, cv
c        write(10,'(2i5,6e12.4)')
c     &  irank,ncol,rcond,e2,freq(ifreq),zmean,cvmean,cvmean1
c
      enddo  
c
      write(2,13)
c      close(8)
c      close(9)
c      close(10)
   10 format(1x,i4,3(1x,e12.4))
   11 format(6(1x,e12.4))
   12 format(2x,4e12.4)
   13 format(60('-'))
      return
      end
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      subroutine fourier_se(se,Ccov,nbc,nf,zls,lspec,fp,freq,
     &rcond,irank,ncol,RMSV)
      INCLUDE 'TRUXTON6ARRAYS.DIM'
      real se(*), zls(*), freq(*)
      complex Ccov(NBCdim*NFdim,NBCdim*NFdim)      
      complex FP(nbcdim, nbcdim)
      real meansize, meanse, meanvar
      complex z1, z2, z3, zzz, zvar
c         
c      open(8,file='se.out',status='unknown')
      write(2,*) ' '
      write(2,*)'FOURIER INVERTER: STANDARD ERROR REPORT'
      write(2,*) ' '
      write(2,68) RMSV
   68 format(/'NOTE:',/
     & 'All standard errors have been normalized ',/
     & 'by the RMS Velocity Error =', 
     & e12.4,//
     & 'The CVs are not normalized i.e. they are SE/Solution',/)
c
c      open(9,file='se_out1',status='unknown')
c      open(10,file='se_out2',status='unknown')
c      
      do IFREQ=1, nf
c      
        write(2,*) ' '
        write(2,13)
        write(2,*) 'Tidal Constituent #', IFREQ, 
     &                     'Frequency', freq(ifreq)
        write(2,*) ' '
        write(2,*) 'Fourier Space:'
        write(2,*)'Mode #,  Magnitude,  std. error,  CV'
        koffset=(IFREQ-1)*2*lspec
        Noffset=(IFREQ-1)*lspec
        do k=1,lspec
          k1=koffset+2*k-1
          k2=koffset+2*k
          zmag  = sqrt( zls(k1)**2 + zls(k2)**2 )
c          stderr= sqrt(  se(k1)**2 +  se(k2)**2 )
          variance=cabs(Ccov(Noffset+k, Noffset+k))
          stderr=sqrt(variance)
          cv    = RMSV*stderr/zmag
          write(2,10)  k, zmag, stderr, cv
c          write(9,10)  k, zmag, stderr, cv
        enddo
c
c Physical Space
c
        write(2,*) ' '
        write(2,*) 'Physical Space:' 
        write(2,*)'Node #,  Magnitude,  std. error,  CV'
        meansize=0.
        meanse=0.
        meanvar=0.
c Node loop        
        do i=1,nbc
          zzz =cmplx(0.,0.)
          zvar=cmplx(0.,0.)
          koffset=(IFREQ-1)*2*lspec
          koff=(IFREQ-1)*lspec
          joff=(IFREQ-1)*lspec
c Basic fourier synthesis loop          
          do k=1,lspec
            k1=koffset+2*k-1
            k2=koffset+2*k
            z1=cmplx(zls(k1), -zls(k2))
            zzz = zzz + FP(i,k)*z1
c Extra loop for covariance
              do j=1,lspec
              zvar=zvar+fp(i,k)*conjg(fp(i,j))*Ccov(k+koff,j+joff)
              enddo
           enddo
c End synthesis and covariance loops              
          size=cabs(zzz)         
          variance=cabs(zvar)
          sigma=sqrt(variance)
          meansize = meansize + size
          meanvar  = meanvar  + variance 
          meanse   = meanse   + sigma
          cv=RMSV*sigma/size
          write(2,10) i, size, sigma, cv 
c          write(9,10) i, size, sigma, cv 
        enddo 
c end Node loop        
        meansize =meansize/(1.*nbc)
        meanvar  =meanvar/(1.*nbc)
        meanse   =meanse/(1.*nbc)
        rootmeanvar=sqrt(meanvar)
        cv=RMSV*meanse/meansize
        write(2,*)' '
        write(2,*)'Mean Values:'
        write(2,*)'     Magnitude, std error, Root(meanvar), CV'
        write(2,12) meansize, meanse, rootmeanvar, cv
c        write(10,'(2i5,6e12.4)')
c     &  irank,ncol,rcond,e2,freq(ifreq),meansize,sigma,cv
c
      enddo  
c
      write(2,13)
c      close(8)
c      close(9)
c      close(10)
   10 format(1x,i4,3(1x,e12.4))
   11 format(6(1x,e12.4))
   12 format(2x,4e12.4)
   13 format(60('-'))
      return
      end
C
C***********************************************************************
C*****************************************************************
C
      Subroutine Invert3(nrow, ncol,
     &                  als, zls, bls, els,
     &                       se, Cov, Rcond,irank,e2)            
c
c may 23 '97 drl: passes Cov back to main. 
c   Cov is normalized by mean sq. velocity;
c   se is normalized by rms vel.
c May 12 '97 drl: passes std error se() back to main
c Feb. 8 '97, drl; from Invert.f: adds Rcond to arguments passed.
c May 7 '97, drl: added calculation of standard error SE of LS solution
C
c Least Squares Als*Zls=Bls by Lapack SVD
c Nondestructive: Als, Bls untouched
c Answers:
c    Zls contains best-fit answer with minimum L2 norm
c    Els is the least-squares error = Als*zls-bls
c From lsq.f, jan 4 '97
c
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
cgh   parameter( nrowdim=nobs2dim, ncoldim=nbc2dim, nrhsdim=1)
      parameter( nrowdim=nobs2dim, ncoldim=nbc2nfdim, nrhsdim=1)
      parameter( nbdim=max(nrowdim, ncoldim),
     &           nsdim=min(nrowdim, ncoldim) )
C      parameter( lwork=10000 )
      parameter(lwork=3*nsdim +max(2*nsdim, nbdim, nrhsdim) )
C
C From main program:
C
      Real als(nrowdim,*), zls(*), bls(*), els(*)
      Real se(*), Cov(NBC2NFdim,NBC2NFdim) 
C
C Local:
C
      real a(nrowdim,ncoldim),b(nbdim,nrhsdim),
     &     s(nsdim),work(lwork)
c     &     ,se(ncoldim)
C
c      open(unit=3, file='invert.out', status='unknown')  
c
      nrhs=1
c
      write(2,*) ' '  
      write(2,*) 'Output from subroutine Invert3:'
      write(2,*) ' '      
      write(2,*) '    Number of Rows:               ', nrow 
      write(2,*) '    Number of Columns:            ', ncol 
      nb=max(nrow, ncol)
      ns=min(nrow, ncol)      
      if    (nrow.gt.nrowdim 
     &  .or. ncol.gt.ncoldim 
     &  .or. nrhs.gt.nrhsdim) then 
      write(2,*) 'dimensions nrow, ncol, nrhs:'
      write(2,*) nrowdim, ncoldim, nrhsdim      
      write(2,*) 'insufficient dimensions; terminate'
      close(2)
      stop
      endif
c      
      do i=1,nrow
      do j=1,ncol
      a(i,j)=als(i,j)
      enddo
      b(i,1)=bls(i)
      enddo
c      
      call SGELSS(nrow, ncol, nrhs, a, nrowdim, b, nbdim, s, rcond, 
     &           irank, work, lwork, info)
c
      if(info.eq.0) then
      write(2,*) ' '      
      write(2,*) '    Inversion successful with rank',irank 
      else
      write(2,*) ' '      
      write(2,*) 'INVERSION NOT SUCCESSFUL with rank =',irank, 
     &                                          'info =', info
      endif
      write(2,*) '    Condition Number:            ', rcond
c
c Solution and Error vectors
c
      do i=1, ncol
      zls(i)=b(i,1)
      se(i)=0.0
      enddo
      do i=1,nrow
      sum=0.
         do j=1,ncol
         sum=sum+als(i,j)*zls(j)
         enddo
      els(i)=sum-bls(i)
      enddo
c
c Size of error vector
c
      e2=0.
      einf=0.      
      do i=1,nrow
      e2=e2 + els(i)**2
      einf=amax1(einf, abs(els(i)))
      enddo
      write(2,*) ' '

      write(2,*) 'Errors:'
      write(2,*) 'Sum of Squares     ', e2
      e2=e2/(0.5*nrow)
      e2=sqrt(e2)
      write(2,*) 'RMS Velocity Error ', e2
      write(2,*) 'MAX |Error|        ', einf
      write(2,50)
   50 format('    Note: RMS Error combines x-and y-velocity errors.',
     & /'          MAX |Error| is for x- or y-component separately.')
c
c Size of solution
c
      s2=0.
      sinf=0.
      do i=1,ncol
      s2=s2 + zls(i)**2
      sinf=amax1(sinf, abs(zls(i)))
      enddo
      write(2,*) ' '      
      write(2,*) 'Size of Solution:' 
      write(2,*) 'Sum of Squares ', s2
      s2=s2/(0.5*ncol)
      s2=sqrt(s2)
      write(2,*) 'RMS size       ', s2
      write(2,*) 'MAX |size|     ', sinf
      write(2,51)
   51 format('    Note: RMS size combines Re and Im parts.',
     & /'          MAX |size| is for Re or Im part separately.')
c
c Alternate size of error (should be the same as above)
c
      if(nrow.gt.ncol.and.irank.eq.ncol)then
      e2alt=0.
      do i=ncol+1, nrow
      e2alt=e2alt+b(i,1)**2
      enddo      
      write(2,*) ' '      
      write(2,*) 'Sum of Squared Errors: alternate measure', e2alt
      endif      
c
c standard error SE(i) of LS solution:
c       SE(j)=sqrt[ sum{ A(i,j)^2 / S(i)^2 } ]    , j = (1,ncol)
c       se(j) normalized by rms velocity (els) here
c       Summation limited to i=1,irank
c       A(i,j) is ith right singular vector
c       S(i) is ith singular value
c       Assuming that nrow.ge.ncol
c Covariance Cov(i,j) by same formula; Cov normalized by ms vel. 
c
      SigSquare=e2*e2/2.
      if(nrow.lt.ncol) then
         write(2,*) 'nrow.lt.ncol; skipping calculation of std errors'
      else
         do j=1,ncol
           do k=1,ncol
           sum=0.0
             do i=1,irank
             sum = sum + ( A(i,j)*A(i,k) )/( s(i)**2 )
             enddo
           Cov(j,k)=SigSquare*sum           
           Cov(j,k)=Cov(j,k)/(e2*e2)
           enddo
         se(j)=sqrt(Cov(j,j))     
         enddo
c
c confirm orthogonality of right singular vectors
c 
c      write(3,*) '      '
c      write(3,*) 'orthogonality check'
c      write(3,*) ' i,    j,     v_i*v_j'
c      do i1=1,irank
c      do i2=1,irank     
c      sum=0.
c      do j=1,ncol
c      sum=sum+a(i1,j)*a(i2,j)
c      enddo
c      write(3,12) i1, i2, sum
c      enddo
c      enddo
c
c
      endif
c
      write(2,*) ' '           
      write(2,*) ' '      
      write(2,*) 'Answer: Value, standard error'
      write(2,*) ' '           
      do i=1,ncol          
        write(2,11) i, zls(i), se(i)
      enddo
      write(2,24) e2
   24 format(/'NOTE:',/
     & 'The standard errors have been normalized ',/
     & 'by the RMS Velocity Error =', 
     & e12.4/)
c      write(2,*) ' '           
c      write(2,*) ' '           
c      write(2,*) ' errors = A*zls - b:'
c      write(2,*) ' '           
c      write(2,10) (els(i),i=1,nrow)
c      write(2,*) ' '      
      write(2,*) ' '           
      write(2,*) ' Singular Values:'      
      write(2,*) ' '                      
      write(2,10) (s(i),i=1,ns)
   10 format(4(1x,e12.4))
   11 format(1x,i5, 4(1x,e12.4))
   12 format(1x, 2i5, 1x, e12.4)
c                   
c      close(3)  
      return
      end
C
C***********************************************************************
C*****************************************************************FUNDY5
C
      Subroutine InvertWLS3(nrow, ncol,
     &                  als, zls, bls, els, 
     &                   wls, se, Cov, rcond,irank,e2)            
C
C Subroutine InvertWLS3
c
c may 23 '97 drl: passes Cov back to main.
c   Cov is normalized by mean sq. velocity;
c   se is normalized by rms vel.
C May 21 '97, drl: from InvertWLS1.f
C     Adds calculation of std error SE() of estimated bc's and sends
c       same back to main
c     Std error computed from usual SVD relation as in Invert2.f
c     The WLS formulation is normalized so that
c       each RHS has the same variance as BLS:
c            Divide ea. row of  [W + Als* Als] {Zls} = [Als*] {Bls}
c            by sqrt[row_sum((Als*)^2)]  
c
C Feb. 8 '97, drl; from InvertWLS.f
C
C Same as WLS but 
C    a) [Wls] is a general symmetric positive definite (ncol x ncol) matrix 
C    b) [Wls] is passed from main rather than being formed locally
C    c) rcond is also passed, eliminating the need to edit this subroutine.
C
C Solves the Weighted Least Squares problem [Als]{Zls}={Bls}:
C   Min (  {Zls*}[Wls]{Zls} +  || [Als]{Zls}-{Bls} ||  )
C
c Nondestructive: Als, Bls, Wls untouched.
C
C Creates the enhanced normal equations:
C        [W + Als* Als] {Zls} = [Als*] {Bls}
C   then finds the inverse solution by Lapack SVD
C
C Otherwise identical to Subroutine Invert which solves 
C        [Als]{Zls}={Bls} by Lapack SVD
C
C Argument list identical to Invert except for the two
C   extra arguments wls and rcond. 
C 
c Answers:
c    Zls contains best-fit answer
c    Els is the weighted least-squares error = Als*zls-bls
c
C
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
      parameter( nrowdim=nobs2dim, ncoldim=nbc2nfdim, nrhsdim=1)
      parameter( nbdim=max(nrowdim, ncoldim),
     &           nsdim=min(nrowdim, ncoldim) )
C      parameter( lwork=10000 )
      parameter(lwork=3*nsdim +max(2*nsdim, nbdim, nrhsdim) )
C
C From main program:
C
      Real als(nrowdim,*), zls(*), bls(*), els(*)
      Real wls(ncoldim, *)       
      Real se(*), Cov(NBC2NFdim,NBC2NFdim) 
C
C Local:
C
      real a(ncoldim,ncoldim),b(nbdim,nrhsdim),
     &     s(nsdim),work(lwork)
C
c      open(unit=3, file='invertwls.out', status='unknown')  
C
      nrhs=1
      write(2,*) ' '  
      write(2,*) 'Output from subroutine Invert3:'
      write(2,*) ' '      
      write(2,*) '    Number of Rows:               ', nrow 
      write(2,*) '    Number of Columns:            ', ncol 
c      write(2,*) 'Output from subroutine InvertWLS3'
c      write(2,*) '      data nrow, ncol, nrhs:' 
c      write(2,*) nrow,    ncol,    nrhs
      nb=max(nrow, ncol)
      ns=min(nrow, ncol)      
      if    (nrow.gt.nrowdim 
     &  .or. ncol.gt.ncoldim 
     &  .or. nrhs.gt.nrhsdim) then 
      write(2,*) 'dimensions nrow, ncol, nrhs:'
      write(2,*) nrowdim, ncoldim, nrhsdim      
      write(2,*) 'insufficient dimensions; terminate'
      close(2)
      stop
      endif
c     
      do i=1,ncol
      do j=1,ncol
         a(i,j)=wls(i,j)
         do k=1,nrow
         a(i,j)=a(i,j)+als(k,i)*als(k,j)
         enddo
      enddo
      enddo
C      
      do i=1,ncol
         b(i,1)=0.0
         do j=1,nrow
         b(i,1)=b(i,1)+als(j,i)*bls(j)
         enddo
      enddo
c 
c InvertWLS2: Normalize each equation so variance of RHS is same:
c  note: want rowsum of transpose of als
c
      do i=1, ncol
        Rowsum=0. 
        do j=1, nrow
           Rowsum=Rowsum+(als(j,i))**2
        enddo
        Rowsum=sqrt(Rowsum)
        b(i,1)=b(i,1)/Rowsum        
        do j=1,ncol
        a(i,j)=a(i,j)/Rowsum
        enddo 
      enddo
c 
c Done with normalization; proceed to solve
c 
      call SGELSS(ncol, ncol, nrhs, a, ncoldim, b, nbdim, s, rcond, 
     &           irank, work, lwork, info)
c
      if(info.eq.0) then
      write(2,*) ' '      
      write(2,*) '    Inversion successful with rank',irank
      else
      write(2,*) ' '      
      write(2,*) 'INVERSION NOT SUCCESSFUL with rank =',irank, 
     &                                          'info =', info
      endif
      write(2,*) '    Condition Number:            ', rcond
c      write(2,*) 'rcond=', rcond
c
c Solution and Error vectors
c
      do i=1, ncol
      zls(i)=b(i,1)      
      se(i)=0.0
      enddo
      do i=1,nrow
      sum=0.
         do j=1,ncol
         sum=sum+als(i,j)*zls(j)
         enddo
      els(i)=sum-bls(i)
      enddo
c
c Size of error vector
c 
      e2=0.
      einf=0.      
      do i=1,nrow
      e2=e2 + els(i)**2
      einf=amax1(einf, abs(els(i)))
      enddo
      write(2,*) ' '
      write(2,*) 'Errors:'
      write(2,*) 'Sum of Squares      ', e2
      e2=e2/(0.5*nrow)
      e2=sqrt(e2)
      write(2,*) 'RMS Velocity Error  ', e2
      write(2,*) 'MAX |Error|         ', einf
      write(2,50)
   50 format('    Note: RMS Error combines x-and y-velocity errors.',
     & /'          MAX |Error| is for x- or y-component separately.')
c
c Size of solution
c
      s2=0.
      sinf=0.
      do i=1,ncol
      s2=s2 + zls(i)**2
      sinf=amax1(sinf, abs(zls(i)))
      enddo
      write(2,*) ' '      
      write(2,*) 'Size of Solution:' 
      write(2,*) 'Sum of Squares ', s2
      s2=s2/(0.5*ncol)
      s2=sqrt(s2)
      write(2,*) 'RMS size       ', s2      
      write(2,*) 'MAX |size|     ', sinf
      write(2,51)
   51 format('    Note: RMS size combines Re and Im parts.',
     & /'          MAX |size| is for Re or Im part separately.')
c
c standard SE error of LS solution:
c       SE(j)=sqrt[ sum{ A(i,j)^2 / S(i)^2 } ]    , j = (1,ncol)
c       se(j) normalized by rms velocity (sqrt(2)*els) here
c       Summation limited to i=1,irank
c       A(i,j) is ith right singular vector
c       S(i) is ith singular value
c       Assuming that nrow.ge.ncol
c Covariance Cov(i,j) by same formula; Cov normalized by ms vel.
c
      SigSquare=e2*e2/2.
      if(nrow.lt.ncol) then
         write(2,*) 'nrow.lt.ncol; skipping calculation of std errors'
      else
         do j=1,ncol
           do k=1,ncol
           sum=0.0
             do i=1,irank
             sum = sum + ( A(i,j)*A(i,k) )/( s(i)**2 )
             enddo
           Cov(j,k)=SigSquare*sum           
           Cov(j,k)=Cov(j,k)/(e2*e2)
           enddo
         se(j)=sqrt(Cov(j,j))     
         enddo
c
c confirm orthogonality of right singular vectors
c 
c      write(3,*) '      '
c      write(3,*) 'orthogonality check'
c      write(3,*) ' i,    j,     v_i*v_j'
c      do i1=1,irank
c      do i2=1,irank     
c      sum=0.
c      do j=1,ncol
c      sum=sum+a(i1,j)*a(i2,j)
c      enddo
c      write(3,12) i1, i2, sum
c      enddo
c      enddo
c
c
      endif
c
      write(2,*) ' '           
      write(2,*) ' '      
      write(2,*) ' Answer: Value, standard error'
      write(2,*) ' '           
      do i=1,ncol          
        write(2,11) i, zls(i), se(i)
      enddo
      write(2,24) e2
   24 format(/'NOTE:',/
     & 'The standard errors have been normalized ',/
     & 'by the RMS Velocity Error =', 
     & e12.4/)
c      write(2,*) ' '           
c      write(2,*) ' '           
c      write(2,*) ' errors = A*zls - b:'
c      write(2,*) ' '           
c      write(2,10) (els(i),i=1,nrow)
c      write(2,*) ' '      
      write(2,*) ' '           
      write(2,*) ' Singular Values:'
      write(2,*) ' '           
      write(2,10) (s(i),i=1,ns)
   10 format(4(1x,e12.4))
   11 format(1x,i5, 4(1x,e12.4))
   12 format(1x, 2i5, 1x, e12.4)
c
C      close(3)  
      return
      end
C
C***********************************************************************
C*****************************************************************
C
      Subroutine CmplxCov(CCov, Cov, ncol)
      INCLUDE 'TRUXTON6ARRAYS.DIM'  
c
c 23 May '97 drl: converts real cov to complex CCov         
c      
      real Cov(NBC2NFdim,NBC2NFdim)
      complex CCov(NBCdim*NFdim,NBCdim*NFdim)
c
      ncolh=ncol/2
      do i=1,ncolh
        Ireal=2*i-1
        Iimag=2*i
        do j=1,ncolh
          Jreal=2*j-1
          Jimag=2*j
          c1= cov(Ireal, Jreal) +cov(Iimag, Jimag)
c          c2= cov(Iimag, Jreal) -cov(Ireal, Jimag)
            c2=0.
          Ccov(i,j)=cmplx( c1,-c2)
c the -c2 is because the answers from the solver 
c are the complex conjugates of the true values         
        enddo
      enddo
      return
      end
C
C***********************************************************************
C*****************************************************************
C
      Subroutine BnodeCheck(NBC, NBCdim, NBE, NVN, 
     &  KST, NNDIM, NPCN, LST)
c            
      If(NBC.gt.NBCdim) then       
      WRITE(2,*) NBC, NBCdim
      WRITE(2,*)'EXECUTION TERMINATED - NBC.gt.NBCdim'
      close(2)
      STOP
      ENDIF

C      
      If(NBC.le.0) then 
      WRITE(2,*) NBC
      WRITE(2,*)'EXECUTION TERMINATED - NBC.le.0'
      close(2)
      stop
      endif
C      
      If(NBE.ne.0) then 
      WRITE(2,*) NBE
      WRITE(2,*)'EXECUTION TERMINATED - NBE.ne.0'
      close(2)
      stop
      endif
C      
      If(NVN.ne.0) then 
      WRITE(2,*) NVN
      WRITE(2,*)'EXECUTION TERMINATED - NVN.ne.0'
      close(2)
      stop
      endif
C      
      If (KST.le.NNDIM) then 
      WRITE(2,*) KST, NNDIM
      WRITE(2,*)'EXECUTION TERMINATED - KST.le NNDIM'
      close(2)
      stop
      endif
C      
      If(NPCN.ne.0) then
      WRITE(2,*) NPCN
      WRITE(2,*)'EXECUTION TERMINATED - NPCN.ne.0'
      close(2)
      stop
      endif
C      
      If(LST.le.NNDIM) then
      WRITE(2,*) LST, NNDIM
      WRITE(2,*)'EXECUTION TERMINATED - LST.le NNDIM'
      close(2)
      stop
      endif
c
      return 
      end
C
C***********************************************************************
C*****************************************************************
C
C getobs.f  REvised 8Feb97 
C Charles Hannah
C get observations from file format *.m3d 
C  t,x,y,u,v
C  t   in gmt0  (decimal days on IYEAR)
C  x,y in meters and inside mesh being used
C  u,v in m/s 
C

cDJM
C     Subroutine GETOBS(NOBS, Xobs, Yobs, Tobs,Tobsgmt0,IYEAR,
C    & Ubarobs,Vbarobs,Zmaxobs,Zminobs,Zbotobs,nobsdim)

      Subroutine GETOBS(NOBS, Xobs, Yobs, Tobs,Tobsgmt0,IYEAR,
     & Ubarobs,Vbarobs,Zmaxobs,Zminobs,Zbotobs,nobsdim,
     & nndim,nedim,in,x,y,nn,ne)

      real x(nndim),y(nndim)
      integer in(3,nedim)
      real basis(3)
C
c june 9 '97: from getobs_8.f
c
      integer NOBS
      real Xobs(nobsdim), Yobs(nobsdim)
      real Tobs(nobsdim), Tobsgmt0(nobsdim)
      real Ubarobs(nobsdim), Vbarobs(nobsdim) 
      real Zmaxobs(nobsdim), Zminobs(nobsdim), Zbotobs(nobsdim)

      character*120 filename
      character*120 header 

      write(*,*)'Enter observations *.m3d file'
      write(*,*)'format: t x y u v zmax zmin zbot *'
      read(*,'(a)')filename
      open(unit=4, file=filename, status='old') 


C   read and discard header
      HEADER=' '
      DO WHILE(HEADER(1:4).NE.'XXXX')
         READ(4,'(A)')HEADER
      ENDDO
C   read comment line, year, and number of columns
      READ(4,'(A)')HEADER
      READ(4,'(A)')HEADER
      READ(4,*)IYEAR
      READ(4,*)NCOL
        if(NCOL.lt.8) then 
           write(2,*)'GET_OBS: Not enough data columns; Stopping'
           write(2,*)'NCOL = ',NCOL
           close(2)
           stop
        endif
c        
      if(NCOL.gt.8) then 
      write(2,22) NCOL, filename
   22 Format('GetObs: There are',I3, ' data columns in file ',(A)/
     & 'Only the first 8 data columns will be used.',/)
      endif

C   read the forcing data

      i=1       
  100 continue
         read(4,*,end=101)tobsgmt0(i), xobs(i), yobs(i), 
     &   ubarobs(i), vbarobs(i), zmaxobs(i), zminobs(i), zbotobs(i)
cDJM
         CALL BASIS2D(NNDIM,NEDIM,IN,X,Y,NN,NE,xobs(i),yobs(i),
     &                NCLOSEST,NELEM,BASIS)

         if (nelem .eq. 0) then
            write(2,*)'GET_OBS: Observation point not within mesh'
            goto 100
          endif
cDJM
	 if(i.gt.nobsdim) then 
	    write(2,*)'GET_OBS: Too many observations; Stopping '
	    write(2,*)'NOBSDIM = ',nobsdim 
	    close(2)
	    stop
	 endif
         i=i+1
      go to 100
  101 continue
      nobs=i-1
      close(4)
      do i = 1, nobs
	 tobs(i) = (tobsgmt0(i)-tobsgmt0(1))*24.*3600.
      enddo
      write(2,21) nobs,filename 
   21 format('GetObs read',i5,' observations from file ',(a))
      write(2,*) ' '
c
      return
      end
C
C***********************************************************************
C*****************************************************************
C
      SUBROUTINE SAMPLER(ICALL,XOBS,YOBS,NOBS,ZMAXOBS,ZMINOBS,
     &USAMP,VSAMP,
     &X,Y,Z,IN,UBAR,VBAR,U,V,NN,NE,NNV)
C
C from sampler_3d_frac_v.f, 9 June 97
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Observational location arrays
      REAL XOBS(NOBSDIM),YOBS(NOBSDIM),ZOBS(NOBSDIM,NNVDIM),
     &ZOBSI(NNVDIM)
      REAL ZMINOBS(NOBSDIM),ZMAXOBS(NOBSDIM)
      REAL BZMIN(NOBSDIM,2),BZMAX(NOBSDIM,2)
      INTEGER NZMIN(NOBSDIM,2),NZMAX(NOBSDIM,2)
      COMPLEX USAMP(NOBSDIM),VSAMP(NOBSDIM)
      COMPLEX UOBSI(NNVDIM),VOBSI(NNVDIM)
      COMPLEX UINT,VINT,UOBSMIN,VOBSMIN,UOBSMAX,VOBSMAX
C
C FEM Arrays
      REAL X(NNDIM),Y(NNDIM),Z(NNDIM,NNVDIM)
      INTEGER IN(3,NEDIM)
      COMPLEX UBAR(NNDIM),VBAR(NNDIM)
      COMPLEX U(NNDIM,NNVDIM),V(NNDIM,NNVDIM)
C
C Local Arrays
      REAL BASIS(3),BOBS(3,NOBSDIM)
      INTEGER NEOBS(NOBSDIM)
C
C Save Appropriate Local Arguments
      SAVE BOBS,NEOBS,ZOBS,BZMIN,NZMIN,BZMAX,NZMAX
C
C START OF EXECUTABLE PROGRAM
C
C Check dimensioning
      IF(NOBS.GT.NOBSDIM)THEN
         WRITE(2,*)'NOBS exceeds NOBSDIM in SUBROUTINE SAMPLER'
         WRITE(2,*)'=> EXECUTION TERMINATED'
         STOP
      ENDIF
C
C If this is the first call, load the stationary basis function arrays
      IF(ICALL.NE.1)GO TO 10
      DO I=1,NOBS
         XVAL=XOBS(I)
         YVAL=YOBS(I)
         CALL BASIS2D(NNDIM,NEDIM,IN,X,Y,NN,NE,XVAL,YVAL,
     &                NCLOSEST,NELEM,BASIS)
         IF(NELEM.NE.0)THEN
            NEOBS(I)=NELEM
            BOBS(1,I)=BASIS(1)
            BOBS(2,I)=BASIS(2)
            BOBS(3,I)=BASIS(3)
            DO J=1,NNV
               ZOBSI(J)=Z(IN(1,NELEM),J)*BASIS(1)+
     &         Z(IN(2,NELEM),J)*BASIS(2)+Z(IN(3,NELEM),J)*BASIS(3)
               ZOBS(I,J)=ZOBSI(J)
            ENDDO
            CALL BASIS1D1(NNV,ZOBSI,ZMINOBS(I),
     &      BZMIN(I,2),NZMIN(I,2),BZMIN(I,1),NZMIN(I,1))
            IF(NZMIN(I,1).EQ.NNV)THEN
               NZMIN(I,1)=1
               BZMIN(I,1)=1.0
               NZMIN(I,2)=2
               BZMIN(I,2)=0.0
               ZMINOBS(I)=ZOBSI(1)
            ENDIF
            CALL BASIS1D1(NNV,ZOBSI,ZMAXOBS(I),
     &      BZMAX(I,2),NZMAX(I,2),BZMAX(I,1),NZMAX(I,1))     
            IF(NZMAX(I,2).EQ.1)THEN
               NZMAX(I,1)=NNV-1
               BZMAX(I,1)=0.0
               NZMAX(I,2)=NNV
               BZMAX(I,2)=1.0
               ZMAXOBS(I)=ZOBSI(NNV)
            ENDIF
         ELSE
            WRITE(2,*)'Your observational location does not lie on',
     &                ' your horizontal mesh'
            WRITE(2,*)'Executation terminated in SUBROUTINE SAMPLER'
            STOP
         ENDIF
      ENDDO
C
C Interpolate complex velocities to the observation locations
 10   DO I=1,NOBS
         NELEM=NEOBS(I)
         N1=IN(1,NELEM)
         N2=IN(2,NELEM)
         N3=IN(3,NELEM)
         BN1=BOBS(1,I)
         BN2=BOBS(2,I)
         BN3=BOBS(3,I)
         DO J=1,NNV
            UOBSI(J)=U(N1,J)*BN1+U(N2,J)*BN2+U(N3,J)*BN3
            VOBSI(J)=V(N1,J)*BN1+V(N2,J)*BN2+V(N3,J)*BN3
         ENDDO
C
C   Vertically integrate over obversational extent of water column
         NV1=NZMIN(I,1)
         BV1=BZMIN(I,1)
         NV2=NZMIN(I,2)
         BV2=BZMIN(I,2)
         UOBSMIN=UOBSI(NV1)*BV1+UOBSI(NV2)*BV2
         VOBSMIN=VOBSI(NV1)*BV1+VOBSI(NV2)*BV2
         DZ=ZOBS(I,NV2)-ZMINOBS(I)
         UINT=0.5*(UOBSMIN+UOBSI(NV2))*DZ
         VINT=0.5*(VOBSMIN+VOBSI(NV2))*DZ
         DO J=NZMIN(I,2)+1,NZMAX(I,1)
            DZ=ZOBS(I,J)-ZOBS(I,J-1)
            UINT=UINT+0.5*(UOBSI(J)+UOBSI(J-1))*DZ
            VINT=VINT+0.5*(VOBSI(J)+VOBSI(J-1))*DZ
         ENDDO
         NV1=NZMAX(I,1)
         BV1=BZMAX(I,1)
         NV2=NZMAX(I,2)
         BV2=BZMAX(I,2)
         UOBSMAX=UOBSI(NV1)*BV1+UOBSI(NV2)*BV2
         VOBSMAX=VOBSI(NV1)*BV1+VOBSI(NV2)*BV2
         DZ=ZMAXOBS(I)-ZOBS(I,NV1)
         UINT=UINT+0.5*(UOBSI(NV1)+UOBSMAX)*DZ
         VINT=VINT+0.5*(VOBSI(NV1)+VOBSMAX)*DZ
C
C   Compute vertical average
         DZTOT=ZMAXOBS(I)-ZMINOBS(I)
         USAMP(I)=UINT/DZTOT
         VSAMP(I)=VINT/DZTOT
      ENDDO
      RETURN
      END  
C      
C***********************************************************************
C***********************************************************************
C
C***********************************************************************
C***********************************************************************
C User-Specified FUNDY5 Subroutines which are restricted in Truxton
C User-Specified FUNDY5 Subroutines which are restricted in Truxton
C***********************************************************************
C***********************************************************************
c
      SUBROUTINE BC(WW,U,V,NNDIM,III)
C Fixed global variables
      COMPLEX U(NNDIM),V(NNDIM)
C Fixed local assignment
      COMPLEX VNORM
C-----------------------------------------------------------------------
C This subroutine must be modified by the user to assign proper
C   Boundary Conditions.  There are four sections which require
C   modification, beginning at statement labels 3,4,5, and 6.  
C Each section must begin with the appropriately labelled continue 
C   statement, and end with a return statement.
C================Beginning-of-user-specified-instructions===============
C
C User-specified local assignments
C
C
C Beginning of executable code
C
      IF(III.EQ.0) GO TO 3
      IF(III-2)4,5,6
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C This section assigns the frequency of the forcing in SEC**-1 to WW.
C A return statement must follow the WW=___  statement.
C
    3 CONTINUE
      WW=.1E-30
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Surface Elevation BC's are assigned in this section.  One statement of
C   the form U(I)=CVALUE must appear for each specified surface
C   elevation node, where I is the horizontal node number and CVALUE is
C   the complex elevation for node I.  
C A return statement must follow the last elevation BC assignment.
C
    4 CONTINUE
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Nonzero normal velocity BC's are assigned in this section.  One 
C   statement of the form U(I)=CVALUE must appear for each set nonzero
C   normal velocity node, where I is the horizontal node number and 
C   CVALUE is the complex value of the normal velocity for node I.  
C   Positive normal velocity indicates flow out of the system.  
C A return statement must follow the last normal velocity BC assignment.
C
    5 CONTINUE
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Nonzero velocity BC's are assigned in this section.  Two statements
C   of the following form must appear for each nonzero velocity node:
C      U(I)=CXVEL
C      V(I)=CYVEL
C   I is the horizontal node number, (CXVEL,CYVEL) are the complex
C   values of the (X,Y) velocities for node I.  
C A return statement must follow the last velocity BC assignment.
C
    6 CONTINUE
      RETURN
C
C End of routine
C
      END
c
C***********************************************************************
C***********************************************************************
c
      SUBROUTINE ATMOS(ATMX,ATMY,NNDIM)
C Fixed global variables
      COMPLEX ATMX(NNDIM),ATMY(NNDIM)
C-----------------------------------------------------------------------
C This subroutine must be modified by the user to assign proper values
C   to arrays ATMX and ATMY, the right-sides of the conservative X- and
C   Y-momentum equations, respectively.  (i.e. H*PSI)  Two statements of
C   the following form must appear for each node where the atmospheric 
C   forcing is nonzero:
C      ATMX(I)=XVALUE
C      ATMY(I)=YVALUE
C   I is the node number and XVALUE,YVALUE are the complex values of the
C   forcing.  ATMX,ATMY are initialized with the value of zero by the
C   main routine before this routine is called.  Therefore, if no value
C   is assigned for a node in this routine, the default value of zero 
C   is used.
C A return statement must follow the last assignment.
C================Beginning-of-user-specified-instructions===============
C
C User-specified local assignments
C
C
C End of routine
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C End of Restricted FUNDY5 User-Specified Subroutines (BC and ATMOS)
C End of Restricted FUNDY5 User-Specified Subroutines (BC and ATMOS)
C***********************************************************************
C***********************************************************************
C
C***********************************************************************
C***********************************************************************
C END OF FIXED SUBROUTINES.  USER-SPECIFIED SUBROUTINES FOLLOW.
C END OF FIXED SUBROUTINES.  USER-SPECIFIED SUBROUTINES FOLLOW.
C***********************************************************************
C***********************************************************************


