C***********************************************************************
C***********************************************************************
      PROGRAM FUNDY5
C-----------------------------------------------------------------------
C FUNDY5:  version97.1
C          4/11/98 - Added SUBROUTINE LOADLST; which had been
C                    inadvertantly left out of the most recent version.
C          1/22/97
C          Modifications from FUNDY5 version93:
C            1. Standardized complex assignments to F77 Standard 
C                (i.e. CMPLX(REAL,IMAGINARY)).
C            2. Fixed bug in SUBROUTINE READINPF5, whereby an input
C                error resulted if there were no set depth nodes.
C            3. Modified reporting of input errors, such that error
C                statement written to echo file provides indication of 
C                what type of error occured.
C            4. Reverted back to the standard SPRSPAK version of 
C                SUBROUTINE SPRSINSRTD to correct bug related to 
C                incorrect dimensioning of NFTRDIM.
C          Modifications conducted by Christopher E. Naimie
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  
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(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   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 'FUNDY5ARRAYS.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 Include file for parameter statements
C
      INCLUDE 'FUNDY5ARRAYS.DIM'
C
C Arrays filled by Subroutine F5READ
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 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 Load input deck information
C
      CALL READINPF5(NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &INBE,NBE,NORM,NVN,KST,LRVS,NR,NPCN,LST,COR,OMEGA)
C
C Load vertical structure and baroclinic forcing
C
      CALL VERTGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,X,Y,IN,
     &NNV,Z,AK,ENZ,RHOX,RHOY,HRBARELX,HRBARELY,HDOWN)
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'
      STOP
      ENDIF
      WRITE(2,110)NNV
C
C Initialize atmospheric forcing then load nonzero values
C
      DO 10 I=1,NN
         ATMX(I)=CMPLX(0.0,0.0)
         ATMY(I)=CMPLX(0.0,0.0)
 10   CONTINUE
      CALL ATMOS(ATMX,ATMY,NNDIM)
C
C Elevation Solution
C
      CALL ELEVATION5(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 Compute IQ,JQ,PPX,PPY
C
      CALL BUILDPPXPPY5(NN,NE,X,Y,IN,
     &IQ,JQ,PPX,PPY)
C
C Velocity Solution
C
      CALL VELOCITY5(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)
C
C Write the results
C
      CALL OUTPUT5(X,Y,Z,RHS,UBAR,VBAR,UVERT,VVERT,WVERT,OMEGA,
     * NN,NNV,NNDIM,NNVDIM)
C
C Write successful completion message to ECHO file
C
      WRITE(2,120)
      CLOSE(2)
C
C Format statements
C
 110  FORMAT(///,'Set in Subroutine VERTGRID:',
     &/5X,'The number of nodes in the vertical = ',I6)
 120  FORMAT(//,'FUNDY5 Run successfully completed')
      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
      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 'FUNDY5ARRAYS.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)
      IF(TITLE(1).EQ.INCI) GO TO 74
      WRITE(2,*)'INPUT ERROR: Unable to find beginning of incidence list
     &'
      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,*)'INPUT ERROR: Unable to find end of incidence list'
      STOP
      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,*)'INPUT ERROR: Unable to find beginning of nodal location
     &s'
      STOP
      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,*)'INPUT ERROR: Unable to find end of node locations'
      STOP
      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)
      IF(TITLE(1).NE.SET )THEN
      WRITE(2,*)'INPUT ERROR: Unable to find beginning of set depth node
     &s'
      STOP
      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)
      IF(TITLE(1).NE. SPEC)THEN
      WRITE(2,*)'INPUT ERROR: Unable to find beginning of specified norm
     &al nodes'
      STOP
      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)
      IF(TITLE(1).NE.KNOW)THEN
      WRITE(2,*)'INPUT ERROR: Unable to find beginning of known velociti
     &es'
      STOP
      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)
      IF(TITLE(1).NE.AVER)THEN
      WRITE(2,*)'INPUT ERROR: Unable to find beginning of average normal
     &s'
      STOP
      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,*)'INPUT ERROR: Cannot find end of average normals'
         STOP
         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,*)'INPUT ERROR: Unable to find end of average normals'
      STOP
      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)
      IF(TITLE(1).NE.GEOS)THEN
      WRITE(2,*)'INPUT ERROR: Unable to find beginning of geostrophic be
     &s'
      STOP
      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,*)'INPUT ERROR: Unable to find beginning of bathymetries'
      STOP
      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,*)'INPUT ERROR: Unable to find bathymetry for node',I
      STOP
      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,*)'INPUT ERROR: Unable to find end of bathymetries'
      STOP
      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)
         STOP
      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)
         STOP
      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)
         STOP
      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)
         STOP
      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)
         STOP
      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 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,'F U N D Y 5',28X,2(1H*)/
     &2(1H*),68X,2(1H*)/2(1H*),26X,'3-D HYDRODYNAMICS',25X,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*),26X,'JANUARY 19, 1993',26X,2(1H*)/
     &72(1H*)/72(1H*))
 3010 FORMAT(//'Parameters specified in include file FUNDY5ARRAYS.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*****************************************************************FUNDY5
C
      SUBROUTINE ELEVATION5(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 'FUNDY5ARRAYS.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
C Local Arrays
C
      REAL DELX(3),DELY(3)
      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
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,0.0)
         SUMF=CMPLX(0.0,0.0)
         SUMAX=CMPLX(0.0,0.0)
         SUMAY=CMPLX(0.0,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
 101  IF(NR.EQ.0) GO TO 111
C   
      DO 105 I=1,NN
         RBARX(I)=CMPLX(0.0,DEF)
 105  CONTINUE
C   
      CALL BC(W,RBARX,BETAY,NNDIM,1)
      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
         IF(AIMAG(RBARX(N)).NE.DEF)THEN
            RHS(N)=RBARX(N)
C   Assign elevation=0 if Subroutine BC contained no value for this node
         ELSE
            RHS(N)=CMPLX(0.0,0.0)
         ENDIF
         RBARX(N)=CMPLX(0.0,DEF)
 109  CONTINUE
C
      WRITE(2,*)'    Nodes where an elevation BC was incorrectly set:'
      DO 110 I=1,NN
         IF(REAL(RBARX(I)).EQ.0.0.AND.AIMAG(RBARX(I)).EQ.DEF)GO TO 110
         WRITE(2,*)I
         NVIOL=NVIOL+1
 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)
      CALL CSOLVET(2,QG,RHS,NN,NHBW,NNDIM,NBWE)
C
C RHS now contains the elevation solution with all BC's applied
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 'FUNDY5ARRAYS.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 VELOCITY5(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)
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 'FUNDY5ARRAYS.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
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)=CMPLX(0.0,0.0)
         VBAR(N)=CMPLX(0.0,0.0)
         DO 59 J=1,NNV
            UVERT(N,J)=CMPLX(0.0,0.0)
            VVERT(N,J)=CMPLX(0.0,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
      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,0.0)
         BCONT(1)=CMPLX(1.0,0.0)
         CCONT(1)=CMPLX(0.0,0.0)
         DCONT(1)=WVERT(I,1)
         WVERT(I,NNV)=EYE*OMEGA*RHS(I)
         ACONT(NNV)=CMPLX(0.0,0.0)
         BCONT(NNV)=CMPLX(1.0,0.0)
         CCONT(NNV)=CMPLX(0.0,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
C
C The velocity solution is complete
C
      RETURN
      END
C
C***********************************************************************
C*****************************************************************FUNDY5
      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)=CMPLX(0.0,0.0)
         HRBARELY(I)=CMPLX(0.0,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)=CMPLX(0.0,0.0)
         HRBARELY(I)=CMPLX(0.0,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 'FUNDY5ARRAYS.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=CMPLX(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,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,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,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,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
      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)=CMPLX(0.0,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 'FUNDY5ARRAYS.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 'FUNDY5ARRAYS.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 FUNDY5ARRAYS.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 FUNDY5ARRAYS.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 FUNDY5ARRAYS.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 FUNDY5ARRAYS.DIM:  Currently RHOREF =',rhoref
      write(*,*)'If another value is desired, change the value of RHOREF
     & in the FUNDY5ARRAYS.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 'FUNDY5ARRAYS.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
   20 IF(JQ(KMAX).NE.0) THEN
       WRITE(*,*) 'MAX CONNECTIVITY MN=',MN,
     *  'VIOLATED AT I,J=',I,J,'IN SPRSINSRTD'
       STOP
      END IF
      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=CMPLX(0.0,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)=CMPLX(0.0,0.0)
      DUDY(J)=CMPLX(0.0,0.0)
      DVDX(J)=CMPLX(0.0,0.0)
      DVDY(J)=CMPLX(0.0,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)=CMPLX(0.0,0.0)
      DFDY(J)=CMPLX(0.0,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 END OF FIXED SUBROUTINES.  USER-SPECIFIED SUBROUTINES FOLLOW.
C END OF FIXED SUBROUTINES.  USER-SPECIFIED SUBROUTINES FOLLOW.
C***********************************************************************
C***********************************************************************
C
