C
C***********************************************************************
C***********************************************************************
C
C SYSTEM FORTRAN-77
C
C     TRANSPOSED VERSION of BANSOL (DRL 8 MARCH 94):
C            B(column, row)
C            Sequence of calling arguments unchanged from BANSOL
C
C B = LHS MATRIX , DIMENSIONED (mDIM,nDIM) IN CALLING PROGRAM
C R = RIGHT-HAND-SIDE, DIMENSIONED (NDIM) IN CALLING PROGRAM
C NEQ= # OF EQUATIONS
C IBAND = DIAGONAL(COLUMN 1) + HALF-BANDWIDTH
C SOLUTION RETURNS IN R
C KKK = 1 PERFORMS LU DECOMPOSITION DESTRUCTIVELY
C KKK = 2 PERFORMS BACK SUBSTITUTION
C KKK = 3 PERFORMS OPTIONS 1 AND 2
C
      SUBROUTINE BANSOLtr(KKK,B,R,NEQ,IBAND,NDIM,MDIM)
C
C     SYMMETRIC BAND MATRIX EQUATION SOLVER,   (REF, 2)
C     DOCTORED TO IGNORE ZEROS IN LU DECOMP. STEP (DRL 8 MARCH 94)
C     TRANSPOSED VERSION: B(column, row) (DRL 8 MARCH 94)
C
      DIMENSION B(mDIM,nDIM), R(NDIM)
      NRS = NEQ - 1
      NR = NEQ
      IF(KKK .EQ. 2) GO TO 200
C
C  TRIANGULARIZE MATRIX A USING DOOLITTLE METHOD
C
      DO 120 N = 1,NRS
      M= N-1
      MR=MIN0(IBAND,NR-M)
      PIVOT = B(1,N)
      DO  120 L=2,MR
      CP= B(l,N)/PIVOT
      IF (CP.EQ.0.0) GO TO 120
      I = M+L
      J = 0
      DO 110 K=L,MR
      J = J + 1
  110 B(j,I)= B(j,I) -CP*B(k,N)
  120 B(l,N) = CP
      IF(KKK.EQ.1)GO TO 400
C
C  MODIFY LOAD VECTOR R
C
  200 DO 220 N=1,NRS
      M= N-1
      MR = MIN0(IBAND,NR-M)
      CP= R(N)
      R(N)=CP/B(1,N)
      DO 220  L=2,MR
      I = M + L
  220 R(I)=R(I)-B(l,N)*CP
C
C   BACK SOLUTION
C
      R(NR) = R(NR)/B(1,NR)
      DO 320 I = 1,NRS
      N=NR-I
      M= N-1
      MR=MIN0(IBAND,NR-M)
      DO 320 K = 2,MR
      L = M+K
C     STORE COMPUTED DISPLACEMENTS IN LOAD VECTOR R
  320 R(N)= R(N)- B(k,N)*R(L)
  400 RETURN
      END
C***********************************************************************
C***********************************************************************
C
C Miscellaneous:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE BASIS1D1(nn,zref,zdes,b2,n2,b1,n1)
C-----------------------------------------------------------------------
C purpose: This subroutine evaluates the value of the basis functions 
C            alive on a 1-D linear element at a point
C
C restrictions: Applicable only for  1-D linear elements
C
C inputs:  nn is the number of entries in the reference z array
C          zref(nn) is the 1-D array containing the reference z array
C               NOTE: zref must increase from zref(1) to zref(nn)
C          zdes is the z value at which the basis functions are desired
C
C outputs: n1 is the reference node immediately below zdes
C          b1 is the value of the n1's basis function at zres
C          n2 is the reference node immediately above zdes
C          b2 is the value of the n2's basis function at zres
C
C notes:   zdes<zref(1)  => n1=nn,b1=0.0,n2=1,b2=1.0 
C          zdes>zref(nn) => n1=nn,b1=1.0,n2=1,b2=0.0 
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           26 AUGUST 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension zref(nn)
C
C START OF EXECUTABLE CODE
      do 10 k=1,nn 
         if(zdes.LT.zref(k))then
C
C if zdes<zref(1), zdes is below the reference data => set b1=0,n1=nn
            n2=k
            n1=k-1
            if(n2.EQ.1)then
                n1=nn
                b1=0.0
                b2=1.0
                go to 20
             endif
             zdiff=zref(n2)-zref(n1)
             b1=(zref(n2)-zdes)/zdiff
             b2=(zdes-zref(n1))/zdiff  
             go to 20
         endif
 10   continue
C
C if zdes>zref(nn), zdes is above the reference data => set b2=n2=0
      b1=1.0
      n1=nn
      b2=0.0
      n2=1
 20   continue
C
C END OF ROUTINE
      return
      end
C  
C*********************************************************************
C*********************************************************************
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***********************************************************************
      SUBROUTINE ELEMCOEFS(I,NNDIM,NNV,ENZMI,ENZHI,ENZQI,
     &                          ENZM,ENZH,ENZQ)
C-----------------------------------------------------------------------
C...Store the elemental based values of the diffusion coefficients
C   into 3-D arrays
C-----------------------------------------------------------------------
      REAL ENZMI(*),ENZHI(*),ENZQI(*),
     &     ENZM(NNDIM,*),ENZH(NNDIM,*),ENZQ(NNDIM,*)
      NNV1=NNV-1
      ENZM(I,1)=ENZMI(1)
      ENZH(I,1)=ENZHI(1)
      ENZQ(I,1)=ENZQI(1)
      DO J=2,NNV1
         ENZM(I,J)=ENZMI(J)
         ENZH(I,J)=ENZHI(J)
         ENZQ(I,J)=ENZQI(J)
      ENDDO
      ENZM(I,NNV)=0.
      ENZH(I,NNV)=0.
      ENZQ(I,NNV)=0.
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE HALFBW(IN,NLOC,NE,NB)
C-----------------------------------------------------------------------
C  NB IS THE HALF-BANDWIDTH
C  NE IS THE NUMBER OF ELEMENTS
C  NOC IS THE FIRST DIMENSION OF IN
C  IN CONVENTION:IN(3,NEDIM)
C-----------------------------------------------------------------------
      DIMENSION IN(NLOC,*)      
      NB=0       
      DO L=1,NE       
         MI=IN(1,L)               
         MA=0       
         DO I=1,3        
            J=IN(I,L)           
            IF(J.LT.MI) MI=J             
            IF(J.GT.MA) MA=J             
         ENDDO 
         IDIFF=MA-MI         
         IF(IDIFF.GT.NB) NB=IDIFF     
      ENDDO 
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C Sparse matrix manipulation:
C
C***********************************************************************
C***********************************************************************
      FUNCTION KAY(I,J,IQ,JQ)
C-----------------------------------------------------------------------
C  FINDS ADDRESS K IN QV FOR Q(I,J)
C-----------------------------------------------------------------------
      DIMENSION IQ(*),JQ(*)
      KAY=IQ(I)+1
   10 KAY=KAY-1
      IF(JQ(KAY).NE.J)GO TO 10
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE NODECOEFS(I,NNDIM,NNV,ZMIDI,ENZMI,ENZHI,ENZQI,
     &                          ENZM,ENZH,ENZQ)
C-----------------------------------------------------------------------
C...Convert the elemental based values of Nz to nodal based using
C   Nz(i) = (dz(I)*Nz(I) + dz(II)*Nz(II))/(dz(I) + dz(II))
C-----------------------------------------------------------------------
      REAL ENZMI(*),ENZHI(*),ENZQI(*),ZMIDI(*),
     &    ENZM(NNDIM,*),ENZH(NNDIM,*),ENZQ(NNDIM,*)
      NNV1=NNV-1
      ENZM(I,1)=ENZMI(1)
      ENZH(I,1)=ENZHI(1)
      ENZQ(I,1)=ENZQI(1)
      DO J=2,NNV1
         DZ1=ZMIDI(J)-ZMIDI(J-1)
         DZ2=ZMIDI(J+1)-ZMIDI(J)
         DZTMP=DZ1+DZ2
         ENZM(I,J)=(DZ1*ENZMI(J-1)+DZ2*ENZMI(J))/DZTMP
         ENZH(I,J)=(DZ1*ENZHI(J-1)+DZ2*ENZHI(J))/DZTMP
         ENZQ(I,J)=(DZ1*ENZQI(J-1)+DZ2*ENZQI(J))/DZTMP
      ENDDO
      ENZM(I,NNV)=ENZMI(NNV1)
      ENZH(I,NNV)=ENZHI(NNV1)
      ENZQ(I,NNV)=ENZQI(NNV1)
      RETURN
      END
C*********************************************************************
C*********************************************************************
C
C   Symmetric banded matrix setup:
C  
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE EMATBUILD(A,AA,IDBC,NN,NB,NDBC,NB2DIM,NB1DIM)
C 
C DRL, 22 Nov, 1998 - From RHS build (see CDRL below) for CASCO 
C                     goes with EMATMULT
C DRL, 9 March 1994
C
C Condensation of a symmetric banded matrix in BANSOLtr (transposed addressing) 
C   form for symmetry-preserving Dirichlet BC enforcement 
C
C Global Columns [IDBC(I), I=1,NDBC] of A are copied into AA and replaced with zeros in A.
C Global Rows    [IDBC(I), I=1,NDBC] of A are set to zero with diagonal set to unity.
C Corresponding rows in AA are set to zero.
C 
C A is a matrix stored in transpose, column-banded, symmetric form
C     A(column, row) is A(nb+1, nn) (BANSOLtr form)
C     A(1,*) is the diagonal
C AA is a matrix in column-banded, column-condensed form:
C     A(row, column) is A((2*nb+1), nx)
C     Note this is conventional, not transposed, addressing;
C          and AA is not symmetric
C
C Column condensation of AA:
C   IDBC(k) is a list of global node numbers for Dirichlet BC's
C   IDBC(k) is the global column number of column k in AA
C       AA(*,k) multiplies the BC data for the kth Dirichlet BC
C   NDBC is the number of Dirichlet BC's (number of columns in AA)
C
C Column-banding of AA:
C   Columns stored have total band height 2NB+1
C   The column is stored beginning in row 1
C   The main diagonal is at row NB+1
C
C NB is the half-bandheight (Full BW=2*NB+1)
C NDBC is the number of Dirichlet BC's
C            (number of columns to be copied into AA)
C NN number of nodes globally
C NB1DIM is the first physical dimension of A in the main program
C NB2DIM is the first physical dimension of AA in the main program
C    --- NB1DIM.ge.(  NB+1) ---
C    --- NB2DIM.ge.(2*NB+1) ---
C
      DIMENSION A(NB1DIM,*),AA(NB2DIM,*),IDBC(*)
C
C Check dimensions passed
C
      NB2=2*NB+1
      IF (NB2DIM.LT.NB2) THEN
         WRITE(*,*) ' INSUFFICIENT DIMENSION FOR MATRIX AA IN'//
     &              ' SUBROUTINE RHSMULT'
         WRITE(*,*) ' NB2DIM =', NB2DIM, ' NB2 =', NB2
         WRITE(*,*) ' EXECUTION TERMINATED'
         STOP
      ENDIF
C
      NB1=NB+1
      IF (NB1DIM.LT.NB1) THEN
         WRITE(*,*) ' INSUFFICIENT DIMENSION FOR MATRIX A IN'//
     &              ' SUBROUTINE RHSMULT'
         WRITE(*,*) ' NB1DIM =', NB1DIM, ' NB1 =', NB1
         WRITE(*,*) ' EXECUTION TERMINATED'
         STOP
      ENDIF
C
C Zero the AA matrix
C
      DO 5 I=1,NB2
      DO 5 J=1,NDBC
    5 AA(I,J)=0.0
C
C Begin the copying/zeroing
C
      DO 30 I=1,NDBC
      II=IDBC(I)
C
C On and Above the Diagonal:
C
      DO 10 L=0,MIN0(NB,II-1)
   10 AA(NB+1-L,I)=A(L+1,II-L)
C 
C Below the Diagonal:(exploiting symmetry, use row of A)
C
      DO 20 J=2,NB1
   20 AA(NB+J,I)=A(J,II)
C
   30 CONTINUE
C
C Now the AA matrix is complete.  
C Note that AA transpose is the matrix of Galerkin equations for Dirichlet nodes.
C Begin zeroing the A matrix
C
      DO 60 I=1,NDBC
      II=IDBC(I)
C
C On and Above the Diagonal:
C  
      A(1,II)=1.0
      IF (II .GT. 1) THEN
      DO 40 L=1,MIN0(NB,II-1)
   40 A(L+1,II-L)=0.
      ENDIF
C 
C Below the Diagonal:(exploiting symmetry)
C
      DO 50 J=2,NB1
   50 A(J,II)=0.
C
   60 CONTINUE
C
C The A matrix is now ready with rows and columns zeroed and unity on diagonal.
C Now zero the Dirichlet rows of AA
C Note AA transpose will no longer hold the Dirichlet Galerkin eqns.
C
      DO 90 J=1,NDBC
      JJ=IDBC(J)
      DO 90 I=1,NDBC
      II=IDBC(I)
      IA=(NB+1)+(II-JJ)
      IF(IA.GE.1 .AND. IA.LE.NB2) AA(IA,J)=0.0
   90 CONTINUE
C
C Finished
C
CDRL  22 Nov 1998  No changes above this line
C    Convert AA to I-AA.  
C    This is EMAT for CASCO
      DO I = 1, NB2 
	 DO J = 1, NDBC
	    AA(I,J) = - AA(I,J)
	 ENDDO
      ENDDO
      DO J = 1, NDBC
	 AA(NB+1,J) = 1.0 
      ENDDO
CDRL End of changes
C
      RETURN
      END
C
C In Quoddy:
C      
C    LRVS(I) - node numbers of set depth and set velocity nodes
C                I=1,NR         are set depth nodes (NCD=5,6,9)
C                I=NR+1,NR+NPCN are V=0.0           (NCD=10)
C                I=LST,NNDIM    are V .NE. 0.0      (NCD=8,11)
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE EMATMULT(A,X,R,JX,NN,NB,NX,NB2DIM)
C 
CDRL 22 Nov 1998 : based on RHSMULT (see CDRL below) for CASCO
C  only change: r=r-A*x becomes r=r+A*x
C DRL, 9 March 1994
C
CDRL Performs the matrix algebra r=r+A*x for known r, A, and x
C A and x are undisturbed; r is overwritten as above
C
C r is a n-length vector
C A is a matrix in column-banded, column-condensed form:
C     A(row, column) is A((2*nb+1), nx)
C x is a n-length vector
C
C Column condensation:
C   JX(k) is a pointer from local to global node number
C   JX(k) is the global column number of column k in A
C       A(*,j) multiplies x(jx(j))
C
C Column-banding:
C   Columns stored have total band height 2NB+1
C   The column is stored beginning in row 1
C   The main diagonal is at row NB+1
C
C NB is the half-bandheight (Full BW=2*NB+1)
C NX is the number of columns in A (Number of Type I BC's)
C NN is the length of the r and x vectors (Number of nodes globally)
C NB2DIM is the first physical dimension of A in the main program
C    --- NB2DIM.ge.(2*NB+1) ---
C
      DIMENSION A(NB2DIM,*), R(*), X(*),JX(*)
C
C Check dimensions passed
C
      NB2=2*NB+1
      IF (NB2DIM.LT.NB2) THEN
         WRITE(*,*) ' INSUFFICIENT DIMENSION FOR MATRIX A IN'//
     &              ' SUBROUTINE RHSMULT'
         WRITE(*,*) ' NB2DIM =', NB2DIM, ' NB2 =', NB2
         WRITE(*,*) ' EXECUTION TERMINATED'
         STOP
      ENDIF
C
      DO 10 J=1,NX
C
C IR is the global row number (Position in r)
C I is the local row number in A (column-banded)
C J  is the local column number in A (columns are condensed)
C JX is the global column number (Position in x)
C Row correspondence:  I = (NB+1) + (IR-JX(J))    
C
      IR1 = MAX0(1,JX(J)-NB)
      IR2 = MIN0(NN,JX(J)+NB)
      I  = NB+1 + (IR1-JX(J))
C
CDRL r = r+A*x
      DO 20 IR=IR1,IR2
      R(IR)=R(IR)+A(I,J)*X( JX(J) )
      I=I+1
   20 CONTINUE
C
   10 CONTINUE
      RETURN
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE SINEGRID (NNDIM, NN, NNV, HDOWN, HMID, DZBL, Z)
C-----------------------------------------------------------------------
C Sinusoidal mesh and coefficients version of vertgrid
C-----------------------------------------------------------------------
C
C Global Arrays
C
      REAL Z(NNDIM,*), HDOWN(*), HMID(*)
C
C Set sinusoidal grid for the QUODDY series
C
      NEV=NNV-1
      PI=ACOS(-1.)
      DO I=1,NN
         ZETA=HMID(I)-HDOWN(I)
         A=(HDOWN(I)+ZETA-DZBL*NEV)/(NEV*SIN(2.*PI/NEV))
         IF(A .LT.0.)A=0.
         Z(I,1)=-HDOWN(I)
         DO J=2,NNV
            EPS=(J-1.0)/(NNV-1.0)
            Z(I,J)=-HDOWN(I)*(1.-EPS)-A*SIN(2.*PI*EPS)+EPS*ZETA
         ENDDO
         Z(I,NNV)=ZETA
      ENDDO
      RETURN
      END
C  
C*********************************************************************
C*********************************************************************
C
      SUBROUTINE SOLVETR(KKK,B,R,NEQ,IHALFB,NDIM,MDIM)        
C
C   16 FEB 1992 -- SOLVES FOR B TRANSPOSED -- (COLUMN, ROW)
C                
C   ASYMMETRIC BAND MATRIX EQUATION SOLVER           
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                
      DIMENSION B(MDIM,ndim),R(NDIM)        
      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             
      IF (C.EQ.0.0) GO TO 21       
      B(KC,i)=C           
      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 In Quoddy:
C      
C    LRVS(I) - node numbers of set depth and set velocity nodes
C                I=1,NR         are set depth nodes (NCD=5,6,9)
C                I=NR+1,NR+NPCN are V=0.0           (NCD=10)
C                I=LST,NNDIM    are V .NE. 0.0      (NCD=8,11)
C
C***********************************************************************
C***********************************************************************
C
C...Banded matrix solvers:
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE THOMAS(KKK,A,B,C,R,N)            
C
C THOMAS ALGORITHM FOR TRIDIAGONAL MATRICES  -  SYSTEM FORTRAN-77
C REAL 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
      DIMENSION 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
C Vertical meshing:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE UNISIGMA (NNDIM, NN, NNV, BATHY, ELEVN, Z)
C-----------------------------------------------------------------------
C ASSIGNS UNIFORM VERTICAL MESH Z(I,J) UNDER EACH HORIZONRAL NODE
C  (UNIFORM DZ AMONG NNV NODES)
C TOP:     J=NNV,  Z(I,J) =  ELEVN(I)
C BOTTOM:  J=1,    Z(I,J) = -BATHY(I)
C NNDIM IS THE FIRST PHYSICAL DIMENSION OF Z(I,J) IN MAIN.
C-----------------------------------------------------------------------
      DIMENSION Z(NNDIM,*),ELEVN(*),BATHY(*)
C
C Set uniform grid for FUNDY series
C
      DO I=1,NN
         DZ=BATHY(I)/(NNV-1.)
         ZJ=-BATHY(I)
         DO J=1,NNV
            Z(I,J)=ZJ
            ZJ=ZJ+DZ
         ENDDO
         Z(I,NNV)=0.
      ENDDO
C
C Adjust for QUODDY series
C
      DO I=1,NN
         ZETA=ELEVN(I)-BATHY(I)
         DO J=1,NNV
            ZFUNDY=Z(I,J)
            ZFACTR=(ZFUNDY+BATHY(I))/BATHY(I)
            Z(I,J)=ZFUNDY+ZETA*ZFACTR
         ENDDO
      ENDDO
      RETURN
      END

C********************************************************************
C*********************************************************************
      SUBROUTINE VERTAVG(F,FAVG,Z,NN,NNV,NNDIM)
      REAL F(NNDIM,*), FAVG(*), Z(NNDIM,*)
C-----------------------------------------------------------------------
C COMPUTES THE VERTICAL AVERAGE OF F, RETURNS IT IN FAVG.
C DOES THIS FOR ALL HORIZONTAL NODES
C DOES NOT REQUIRE UNIFORM MESH
C INTEGRATION IS PERFORMED FROM NODE 1 TO NODE NNV
C NNDIM IS THE FIRST PHYSICAL DIMENSION OF F AND Z IN MAIN.
C-----------------------------------------------------------------------
      DO I=1,NN
         FINT=0.0
         DO J=2,NNV
            FINT=FINT+0.5*(F(I,J)+F(I,J-1))*(Z(I,J)-Z(I,J-1))
         ENDDO
         FAVG(I)=FINT/(Z(I,NNV)-Z(I,1))
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTSUM(F,FSUM,NN,NEV,NNDIM)
C-----------------------------------------------------------------------
C COMPUTES THE VERTICAL SUM OF F, RETURNS IT IN FSUM.
C DOES THIS FOR ALL HORIZONTAL NODES
C SUMMATION IS PERFORMED FROM ELEMENT 1 TO ELEMENT NEV
C NNDIM IS THE FIRST PHYSICAL DIMENSION OF F IN MAIN.
C-----------------------------------------------------------------------
      REAL F(NNDIM,*), FSUM(*)
      DO I=1,NN
         FS=0.0
         DO J=1,NEV
            FS=FS+F(I,J)
         ENDDO
         FSUM(I)=FS
      ENDDO
      RETURN
      END
