      program proc

      include 'ACADIA2d.DIM'

      integer in(3,nedim)
      integer nn,ne
      real x(nndim),y(nndim),bathy(nndim)
      character*72 casnam

      casnam = 'marmap1'

      call loadmeshgom(nndim,nedim,casnam,x,y,in,bathy,nn,ne)

      open(unit=1,file='cal_stage.dat',form='formatted',status='old')

  1   read(1,*,end=2)itile,iseas,xlon,xlat,stg0,stg5,stg4,stg3,stgt,
     .               s0,s5,s4,s3
      read(1,*)stot,nobs,avgdp,tarea

      call mercgom2(1,xlon,xlat,xpos,ypos)

      call closestnode(casnam,nn,xlon,xlat,inode)

      hmarmap = amin1(bathy(inode)-5.,200.)
      c0 = (10.**stg0 - 1.) / 10. / hmarmap
      c5 = (10.**stg5 - 1.) / 10. / hmarmap
      c4 = (10.**stg4 - 1.) / 10. / hmarmap
      c3 = (10.**stg3 - 1.) / 10. / hmarmap
      ct = (10.**stgt - 1.) / 10. / hmarmap

      write(10+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               c0,c0,c0,.1
      write(20+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               c5,c5,c5,.1
      write(30+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               c4,c4,c4,.1
      write(40+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               c3,c3,c3,.1
      write(50+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               ct,ct,ct,.1

      goto 1

  2   continue

 105  format(12e13.5)

c
c output format x,y,z,year day,t,s,st,.1
c

      end

C=======================================================================
C========= ROUTINES OF MERCATOR PROJECTION =============================
C=======================================================================
C SUBROUTINE MERCGOM1
C Provided to NML by D. Greenberg
C Latest Revision by Christopher E. Naimie
C Dartmouth College - 25 NOV 1991
C
C This set of routines was sent to Dartmouth by D. Greenberg.
C   In its present form it performs a mercator projection
C   which is valid for the gom1,gb1, and bf1 meshes to within
C   0.0025 deg (5th sig fig.) for longitude and latitude measurements
C   and to within 200 meters(4th sig fig.) for x and y measurements.
C
C The projection type is set by NOPT:
C   NOPT=1  =>  convert provided (dlong,dlat) to (x,y)
C   NOPT=2  =>  convert provided (x,y) to (dlong,dlat)
C
C Units: XLON is the longitude in degrees, positive east of Grenwich.
C        XLAT is the latitude in degrees, positive north of the Equator.
C        the x and y coordinates (XPOS,YPOS) are in meters.
C-----------------------------------------------------------------------

      SUBROUTINE MERCGOM1(NOPT,XLON,XLAT,XPOS,YPOS)
C
C      scale is specified at the equator
C      R = degrees/radian   S = cm per degree of latitude
C
      DATA A,B/-1.57079632,0.0067226/
      DATA R,S/57.2957795,1.111338e7/
C Set constants
      SCALE=1353100.3
      xoffset=17.1108
      yoffset=28.5281
C
      if(kkk.NE.64)then
         SF=SCALE/S
         RLAT=42.35
         RLON=-71.0333
      kkk=64
      endif
C
      IF(NOPT.EQ.1)then
         YPOS=((VB(XLAT/R)-VB(RLAT/R))*R/SF+yoffset)*10000.0
         XPOS=((XLON-RLON)/SF+xoffset)*10000.0
         RETURN
C
      else IF(NOPT.EQ.2)then
         XLON=((XPOS/10000.0-xoffset)*SF)+RLON
         ABA=((YPOS/10000.0-yoffset)*SF/R)+VB(RLAT/R)
         ABB=A+2.0*ATAN(EXP(ABA))
         ABC=A+2.0*ATAN(EXP(ABA+B*SIN(ABB)))
         XLAT=R*ABC
         RETURN
      else
         PRINT*,'INVALID MERCATER OPTION'
      endif
      RETURN
      END
C=======================================================================
      FUNCTION VB(X)
C *
C * INCREASING LATITUDE IN A MERCATOR PROJECTION
C * X IN RADIANS, VB IN RADIANS
C
      VB=ALOG(TAN(0.78539816+.5*X))-0.0067226*SIN(X)
      return
      END
C=======================================================================
C*****67****************************************************************
C***********************************************************************
C SUBROUTINE mercgom2
C Christopher E. Naimie
C Dartmouth College
C Latest version - 3 MAY 1993
C
C This subroutine performs a mercator projection, using 
C   Boston, MA., USA as the (x,y) origin.
C 
C Units: dlong is the longitude in degrees, positive east of Grenwich.
C        dlat is the latitude in degrees, positive north of the Equator.
C        x and y are in MKS units.
C
C The constants used are:
C      longitude at Boston = -71.03 deg E
C      latitude at Boston  =  42.35 deg N
C      radius of the earth =  6.3675E+6 meters 
C
C The projection type is set by nopt:
C   nopt=1  =>  convert provided (dlong,dlat) to (x,y)
C   nopt=2  =>  convert provided (x,y) to (dlong,dlat)
C-----------------------------------------------------------------------
      SUBROUTINE mercgom2(nopt,dlong,dlat,x,y)
C Set constants and (xo,yo) coordinates of Boston
      r=6.3675E+6                            
      reflong=rad(-71.03)                    
      reflat=rad(42.35)                      
      xo=r*COS(reflat)*reflong
      yo=r*COS(reflat)*ALOG((1.0+SIN(reflat))/COS(reflat))
C If nopt=1:  Convert (long,lat) to (x,y) and subtract (xo,yo)
      if(nopt.EQ.1)then
         rlong=rad(dlong)              
         rlat=rad(dlat)
         x=r*COS(reflat)*rlong-xo
         y=r*COS(reflat)*ALOG((1.0+SIN(rlat))/COS(rlat))-yo             
         return
C If nopt=2:  Convert (x+xo,y+xo) to (long,lat)
      else if(nopt.EQ.2)then
         rlong=(x+xo)/(r*COS(reflat))
         call newtlat(rlat,y+yo,reflat,r)
         dlong=deg(rlong)
         dlat=deg(rlat)
         return
      else
         write(*,*)'SUBROUTINE mercbos requires nopt be set to 1 or 2.'
         stop
      endif
      end      
C*****67****************************************************************
C***********************************************************************
      SUBROUTINE NEWTLAT(rlat,y,reflat,r)
C-----------------------------------------------------------------------
C Christopher E. Naimie
C Dartmouth College - 25 NOV 1991
C
C This subroutine uses the Newton-Raphson Method to calculate the 
C   latitude (in radians), given the following information:
C      r=radius of the earth (m)
C      y=y coordinate of the point who's latitude is desired (m/sec)
C      reflat=reference latitude where y is zero (radians)
C The latitude at y is returned as rlat (in radians).
C-----------------------------------------------------------------------
C Set constants
      tol=0.0001
      n=1000
C  Enter iteration loop
      po=reflat
      do 10 i=1,n
         p1=po-((1.0+SIN(po))/COS(po)-EXP(y/(r*COS(reflat))))
     &         /((1.0+SIN(po))/((COS(po))**2.0))
         if(ABS(p1-po).lt.tol)then
           rlat=p1
           return
         endif
         po=p1
 10   continue
C
C  If loop has completed, the method failed.
C
      write(*,*)' A root was not found after ',n,' iterations.'
      return
      end
C****67****************************************************************
C****67****************************************************************
C Convert degrees to radians
      real function rad(degrees)
      pi=3.14159265359
      rad=degrees*pi/180.0
      return
      end
C*****67****************************************************************
C***********************************************************************
C Convert radians to degrees
      real function deg(radians)
      pi=3.14159265359
      deg=radians*180.0/pi
      return
      end
C*****67****************************************************************
C***********************************************************************
C-----------------------------------------------------------------------
      SUBROUTINE CLOSESTNODE(filemesh,nn,xval,yval,nclosest)
C-----------------------------------------------------------------------
C
C purpose: This subroutine finds the node closest to a point specified
C            by the coordinates (xval,yval)
C  
C restrictions:  all nodes must lie in a flat plane
C
C inputs:  x(nn) and y(nn) are 1-D arrays containing the x and y
C            coordinates of each node
C          nn is the number of nodes
C          xval,yval are the point coordinates
C
C outputs: nclosest is the node who s (x,y) coordinates are closest to
C            (xval,yval)
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
c
c           revision by Wendy Gentleman
c 		Dartmouth College
c 		26 October 1994
C
C-----------------------------------------------------------------------
C
C ARGUMENTS
C
	include 'ACADIA2d.DIM'
      real x(nndim),y(nndim)
	character*72 filemesh, llfile
C
C START OF EXECUTABLE CODE
C
c  read in lat/long coords
	indx=index(filemesh,' ')-1
	llfile=filemesh(1:indx)//'ll.nod'
	open (unit=2, file=llfile)
	do i=1,nn
	  read(2,*) id, x(i), y(i)
	end do
c
c
C  Find node closest to (xval,yval)
C
      dclosest=SQRT((x(1)-xval)**2.0+(y(1)-yval)**2.0)
      nclosest=1
      do 25 j=2,nn
         dval=SQRT((x(j)-xval)**2.0+(y(j)-yval)**2.0)
         if(dval.LT.dclosest)then
            dclosest=dval
            nclosest=j
         endif
 25   continue
C
C END OF ROUTINE
C
	close(2)
      return
      end


C*********************************************************************
C
C FUNDY4 FIXED SUBROUTINE WRITTEN BY PROFESSOR DANIEL R. LYNCH
C
      FUNCTION PHASELAGD(Z)
C
C COMPUTES THE NEGATIVE ARGUMENT OF A COMPLEX NUMBER Z, IN DEGREES.
C RESULT IS REPORTED IN THE RANGE 0. TO 360. 
C 
      COMPLEX Z
      PI=3.1415926535 
      FACTOR=180./PI
      PHASELAGD = 0.0
      IF(CABS(Z).GT.0.)  PHASELAGD = -FACTOR*ATAN2(AIMAG(Z), REAL(Z))
      IF(PHASELAGD.LT.0.) PHASELAGD = PHASELAGD + 360.  
      RETURN
      END
C

C*********************************************************************
C*********************************************************************
C
C FUNDY4 FIXED SUBROUTINE WRITTEN BY PROFESSOR DANIEL R. LYNCH
C
      FUNCTION PHASELAG(Z)
C
C COMPUTES THE NEGATIVE ARGUMENT OF A COMPLEX NUMBER Z, IN RADIANS. 
C RESULT IS REPORTED IN THE RANGE 0. TO 2*PI. 
C 
      COMPLEX Z
      PI=3.1415926535 
      PHASELAG = 0.0
      IF(CABS(Z).GT.0.)  PHASELAG = -ATAN2(AIMAG(Z), REAL(Z))
      IF(PHASELAG.LT.0.) PHASELAG = PHASELAG + 2.*PI
      RETURN
      END
C
C*********************************************************************
C*********************************************************************


C***********************************************************************
C***********************************************************************
      SUBROUTINE LOADVEL(nndim,nnvdim,nn,nnv,filename,zval,u,v,w)
C-----------------------------------------------------------------------
C purpose: This routine loads a .vel file
C
C index range: i=1,nn, j=1,nnv
C
C inputs:  nndim - parameter used to dimension node arrays in the
C            calling routine
C          nnvdim - parameter used to dimension the number of 
C            vertical nodes in the calling routine
C          nn - number of horizontal nodes for which data should 
C            exist
C          nnv - number of vertical nodes for which data should 
C            exist  =>  data should exist for nn*nnv nodes
C          filename - name of file to be read
C
C outputs: zval(i,j) - z coordinate at horizontal node i
C            and vertical node j           
C          u(i,j) - x-component of the vector at horizontal 
C            node i and vertical node j
C          v(i,j) - y-component of the vector at horizontal 
C            node i and vertical node j
C          w(i,j) - z-component of the vector at horizontal 
C            node i and vertical node j
C 
C file format: See CEN 3/29/93 Memorandum
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           30 JULY 1992
C           Revised: 25 APRIL 1994
C-----------------------------------------------------------------------
C
C ARGUMENTS
      real zval(nndim,nnvdim)
      complex u(nndim,nnvdim),v(nndim,nnvdim),w(nndim,nnvdim)
      character*72 filename,header
C
C START OF EXECUTABLE CODE
C
C Open file
      open(63,file=filename,status='old')
C
C Read headers
      read(63,101)header
      read(63,101)header
      read(63,*)nnv
      read(63,101)header
      read(63,*)freq
C
C Check for oversize of arrays
      nflag=0
      if(nn.GT.nndim)then
	 write(*,*)'Dimensioning problem in SUBROUTINE LOADVEL:'
         write(*,*)'  nn exceeds nndim => rectify and recompile'
         nflag=1
      endif
         if(nnv.GT.nnvdim)then
	 write(*,*)'Dimensioning problem in SUBROUTINE LOADVEL:'
         write(*,*)'  nnv exceeds nnvdim => rectify and recompile'
         nflag=1
      endif
      if(nflag.EQ.1)stop
C
C Read data
      do 20 i=1,nndim
         do 10 j=1,nnv
            read(63,*,end=21)n,zval(i,j),uamp,upha,vamp,vpha,wamp,wpha 
            u(i,j)=CMPLX(uamp*COS(-upha),uamp*SIN(-upha))
            v(i,j)=CMPLX(vamp*COS(-vpha),vamp*SIN(-vpha))
            w(i,j)=CMPLX(wamp*COS(-wpha),wamp*SIN(-wpha))
            jj=j
 10      continue
         ii=i
 20   continue
 21   close(63)
C
C Message to standard out
      if(ii.NE.nn)then
         write(*,*)'There is a mismatch between the number of horizontal
     & nodes in your .vel file and nn in SUBROUTINE LOADVEL.'
         STOP
      else
         write(*,*)ii,'x',jj,' complex 3-D vector values loaded from ',
     &filename
      endif
C
C FORMAT STATEMENTS
 101  format(a)
C
C END OF ROUTINE
      return
      end

C***********************************************************************
C***********************************************************************
C***********************************************************************
      SUBROUTINE LOADMESHGOM(nndim,nedim,pathmesh,x,y,in,bathy,nn,ne)
C-----------------------------------------------------------------------
C purpose: This routine loads LABLIB Standard node and element files for
C            a 2-D triangular finite element mesh
C
C index range: i=1,nn, k=1,ne
C
C inputs:  nndim - parameter for dimensioning the nodal coordinates
C          nedim - parameter for dimensioning the incidence list
C          pathmesh - the leading part of the path for the .nod, .ele,
C                       and .bat files.
C                  i.e.: g2s.nod, g2s.ele, and g2s.bat => pathmesh=g2s
C
C outputs: nn - number of nodes in the mesh
C          x(i),y(i) - nodal x and y coordinates for node i
C          bathy(i) - bathymetric depth for node i
C          ne - number of elements in the mesh
C          in(3,k) - LABLIB Standard incidence list for element j
C 
C file format: LABLIB Standard 2-D triangular linear element mesh
C              See CEN 3/29/93 Memorandum
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C           Revised: 25 APRIL 1994
C           Revised: 25 FEBRUARY 1995
C-----------------------------------------------------------------------
C
C ARGUMENTS
      integer in(3,nedim)     
      dimension x(nndim),y(nndim),bathy(nndim)
      character*72 pathmesh
C
C START OF EXECUTABLE CODE
C
C Use cstring to find limits of pathmesh
      CALL CSTRING(pathmesh,lms,lme)
C
C open .nod, .ele, and .bat data files
      open(12,file=pathmesh(lms:lme)//'.nod',status='old')
      open(13,file=pathmesh(lms:lme)//'.bat',status='old')
      open(14,file=pathmesh(lms:lme)//'.ele',status='old')
C
C  Read node coordinates and the scalar value at each node from files
      do 5 i=1,nndim
        read(12,*,end=6)nn,x(i),y(i)
        read(13,*)nn,bathy(i)
 5    continue
         write(*,*)' Potential error reading .nod file =>',
     &             ' increase nndim such that it exceeds the number of',
     &             ' nodes by at least 1 and recompile.'
         stop
 6    continue
      write(*,*)nn, 'nodal coordinates read from ',
     &pathmesh(lms:lme),'.nod'
      write(*,*)nn, 'nodal bathymetries read from ',
     &pathmesh(lms:lme),'.bat'
      close(12)
      close(13)
C
C Read element data from file
      do 10 i=1,nedim
        read(14,*,end=11)ne,in(1,i),in(2,i),in(3,i)
 10   continue
         write(*,*)' Potential error reading .ele file =>',
     &             ' increase nedim such that it exceeds the number of',
     &             ' elements by at least 1 and recompile.'
         stop
 11   continue
      write(*,*)ne, 'elements read from ',pathmesh(lms:lme),'.ele'
      close(14)
C
C FORMAT STATEMENTS
 101  format(a)
C
C END OF ROUTINE
C
      return
      end

C***********************************************************************
C***********************************************************************
      SUBROUTINE WRITES2R(nndim,filename,mesh,header,s,nn)
C-----------------------------------------------------------------------
C purpose: This routine writes a .s2r file
C
C index range: i=1,nn
C
C inputs:  nndim - parameter used to dimension nodal arrays in the
C            calling routine
C          filename - name of file to be read
C          mesh - name of the mesh
C          header - character string describing file contents
C          s(i) - real scalar value at node i
C          nn - number of horizontal nodes for which data exists
C 
C file format: See CEN 3/29/93 Memorandum
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           26 APRIL 1994
C-----------------------------------------------------------------------
      character*72 filename,mesh,header
      real s(nndim)
C
C Check for oversize of arrays
      nflag=0
      if(nn.GT.nndim)then
	 write(*,*)'Dimensioning problem in SUBROUTINE WRITES2R:'
         write(*,*)'  nn exceeds nndim => rectify and recompile'
         nflag=1
      endif
      if(nflag.EQ.1)stop
C
C Use SUBROUTINE CSTRING to extract nonblank part of character strings
      CALL CSTRING(filename,lfs,lfe)
      CALL CSTRING(mesh,lms,lme)
C
C Write file
      open (unit=7,file=filename(lfs:lfe),status='new')
      write(7,101)mesh(lms:lme)
      write(7,101)header
      write(7,103)(I,s(I),I=1,nn)
      close(7)
C
C Format Statements
 101  format(a)
 103  format(i5,e12.4)
      return
      end


C***********************************************************************

C***********************************************************************
      SUBROUTINE LOADS2R(nndim,nn,filename,scval)
C-----------------------------------------------------------------------
C purpose: This routine loads a .s2r file
C
C index range: i=1,nn
C
C inputs:  nndim - parameter used to dimension node arrays in the
C            calling routine
C          nn - number of nodes for which data should exist
C          filename - name of file to be read
C
C outputs: scval(i) - real scalar value at node i
C 
C file format: See CEN 3/29/93 Memorandum
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           30 JULY 1992
C           Revised: 25 APRIL 1994
C-----------------------------------------------------------------------
C
C ARGUMENTS
      real scval(nndim)
      character*72 filename,header
C
C START OF EXECUTABLE CODE
C
C Open file
      open(63,file=filename,status='old')
C
C Read headers
      read(63,101)header
      read(63,101)header
C
C Check for oversize of arrays
      nflag=0
      if(nn.GT.nndim)then
	 write(*,*)'Dimensioning problem in SUBROUTINE LOADS2R:'
         write(*,*)'  nn exceeds nndim => rectify and recompile'
         nflag=1
      endif
      if(nflag.EQ.1)stop
C
C Read data
      do 10 i=1,nndim
         read(63,*,end=11)node,scval(node)
         ii=i
 10   continue
 11   close(63)
C
C Message to standard out
      if(ii.NE.nn)then
         write(*,*)'There is a mismatch between the number of horizontal
     & nodes in your .s2r file and nn in SUBROUTINE LOADS2R.'
         STOP
      else
         write(*,*)ii,' real 2-D scalar values loaded from ',filename
      endif
C
C FORMAT STATEMENTS
 101  format(a)
C
C END OF ROUTINE
      return
      end

C***********************************************************************

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 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,1),IQ(1),JQ(1)
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*********************************************************************
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*********************************************************************
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
      SUBROUTINE SPRSMLTIS3(I,QV1,qv2,qv3,IQ,JQ,X,B1,b2,b3,
     &                      NSTAGE,NNDIM,NSTAGEDIM)
C
C  MULTIPLIES THE Ith ROW OF QV BY THE Jth COLUMN OF THE MATRIX X,
C     FOR ALL VALUES OF J FROM 1 TO NNV
C  ANSWER RETURNS IN THE VECTOR B
C     (OR: OBTAINS ONLY Bi OF Q*X=B)
C  IQ,JQ ARE THE STANDARD POINTERS FOR QV
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF X
C  THIS ROUTINE IS NONDESTRUCTIVE
C
      DIMENSION QV1(1),qv2(1),qv3(1)
      dimension IQ(1),JQ(1),X(NSTAGEDIM,NNDIM)
      dimension B1(nstagedim),B2(nstagedim),B3(nstagedim)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
		
      DO ISTAGE=1,NSTAGE
      B1(ISTAGE)=0.
	b2(istage)=0.
	b3(istage)=0.
	end do

      DO 10 K=KMIN,KMAX
	qv1k=qv1(k)
	qv2k=qv2(k)
	qv3k=qv3(k)
	jqk=jq(k)
	do istage=1,nstage
	xij=X(ISTAGE,JQK)
      B1(ISTAGE)=B1(ISTAGE)+QV1K*Xij
      B2(ISTAGE)=B2(ISTAGE)+QV2K*Xij
      B3(ISTAGE)=B3(ISTAGE)+QV3K*Xij
	end do

   10 CONTINUE
      RETURN
      END
C  







C  
C
      SUBROUTINE GDAY(IDD,IMM,IYear,KD)
C!
C!  GIVEN DAY,MONTH,(EACH 2 DIGITS) and year (four digits), GDAY RETURNS
C!  THE DAY#, KD BASED ON THE GREGORIAN CALENDAR.
C!  THE GREGORIAN CALENDAR, CURRENTLY 'UNIVERSALLY' IN USE WAS
C!  INITIATED IN EUROPE IN THE SIXTEENTH CENTURY. NOTE THAT GDAY
C!  IS VALID ONLY FOR GREGORIAN CALENDAR DATES.
C
C   KD=1 CORRESPONDS TO JANUARY 1, 0000
C	
c 	Note that the Gregorian reform of the Julian calendar 
c	omitted 10 days in 1582 in order to restore the date
c	of the vernal equinox to March 21 (the day after
c	Oct 4, 1582 became Oct 15, 1582), and revised the leap 
c	year rule so that centurial years not divisible by 400
c	were not leap years.
c
C   THIS ROUTINE WAS WRITTEN BY EUGENE NEUFELD, AT IOS, IN JUNE 1990.
C
      INTEGER NDP(13)
      INTEGER NDM(12)
	
      DATA NDP/0,31,59,90,120,151,181,212,243,273,304,334,365/
      DATA NDM/31,28,31,30,31,30,31,31,30,31,30,31/
C!
	lp=6
c make iyy and icc variables
        icc=iyear/100
        iyy=iyear-icc*100

C!  TEST FOR INVALID INPUT:
      IF(ICC.LT.0)THEN
	 WRITE(LP,5000)ICC
	 STOP
      ENDIF
      IF(IYY.LT.0.OR.IYY.GT.99)THEN
	 WRITE(LP,5010)IYY
	 STOP
      ENDIF
      IF(IMM.LE.0.OR.IMM.GT.12)THEN
	 WRITE(LP,5020)IMM
	 STOP
      ENDIF
      IF(IDD.LE.0)THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.NE.2.AND.IDD.GT.NDM(IMM))THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.EQ.2.AND.IDD.GT.29)THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.EQ.2.AND.IDD.GT.28.AND.((IYY/4)*4-IYY.NE.0.OR.(IYY.EQ.0.AND
     .    .(ICC/4)*4-ICC.NE.0)))THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
5000  FORMAT(' INPUT ERROR. ICC = ',I7)
5010  FORMAT(' INPUT ERROR. IYY = ',I7)
5020  FORMAT(' INPUT ERROR. IMM = ',I7)
5030  FORMAT(' INPUT ERROR. IDD = ',I7)
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST CENTURY:
      KD = ICC*36524 + (ICC+3)/4
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST YEAR:
      KD = KD + IYY*365 + (IYY+3)/4
C!
C!  ADJUST FOR CENTURY RULE:
C!  (VIZ. NO LEAP-YEARS ON CENTURYS EXCEPT WHEN THE 2-DIGIT
C!  CENTURY IS DIVISIBLE BY 4.)
      IF(IYY.GT.0.AND.(ICC-(ICC/4)*4).NE.0) KD=KD-1
C!  KD NOW TRULY REPRESENTS THE DAY# OF THE LAST DAY OF LAST YEAR.
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST MONTH:
      KD = KD + NDP(IMM)
C!
C!  ADJUST FOR LEAP YEARS:
      IF(IMM.GT.2.AND.((IYY/4)*4-IYY).EQ.0.AND.((IYY.NE.0).OR.
     .   (((ICC/4)*4-ICC).EQ.0)))   KD=KD+1
C!  KD NOW TRULY REPRESENTS THE DAY# OF THE LAST DAY OF THE LAST
C!  MONTH.
C!
C!  CALCULATE THE CURRENT DAY#:
      KD = KD + IDD
      RETURN
C!
C!
      ENTRY DMY(IDD,IMM,IYear,KD)
C!
C!  GIVEN THE (GREGORIAN) DAY#, KD, AS CALCULATED ABOVE IN THIS ROUTINE,
C!  ENTRY DMY RETURNS THE (GREGORIAN) DAY, MONTH, YEAR AND CENTURY.
C!
C!  TEST FOR VALID INPUT:
      IF(KD.LE.0) WRITE(LP,5040)KD
5040  FORMAT(' KD = ',I7,'  INVALID INPUT. DMY STOP.')
C!
C!  SAVE KD
      KKD=KD
C!  CALCULATE ICC AND SUBTRACT THE NUMBER OF DAYS REPRESENTED BY ICC
C!  FROM KKD
C!  JFH IS THE NUMBER OF 400 YEAR INTERVALS UP TO KKD
C!  JCC IS THE NUMBER OF ADDITIONAL CENTURIES UP TO KKD
      JFH = KKD/146097
      KKD = KKD - JFH*146097
      IF(KKD.LT.36525)THEN
	 JCC = 0
      ELSE
	 KKD = KKD - 36525
	 JCC = 1 + KKD/36524
	 KKD = KKD - (JCC-1)*36524
      END IF
      ICC = 4*JFH + JCC
      IF(KKD.EQ.0)THEN
	 ICC = ICC-1
	 IYY = 99
	 IMM = 12
	 IDD = 31
c	 RETURN
	go to 110
      ENDIF
C!
C!  CALCULATE IYY. JFY IS THE NUMBER OF FOUR YEAR INTERVALS IN THE
C!  CURRENT CENTURY. THE FIRST FOUR YEAR INTERVAL IS SHORT (1460 DAYS
C!  RATHER THAN 1461)IF THE CURRENT CENTURY IS NOT DIVISIBLE BY 4, AND
C!  IN THIS CASE JCC.NE.0 AS CALCULATED ABOVE.
C!
C!  CALCULATE JFY:
      JFY = 0
      IF(JCC.EQ.0)GOTO 10
      IF(KKD.LT.1460)GOTO 10
      JFY = 1
      KKD = KKD - 1460
10    KK = KKD/1461
      JFY = JFY + KK
      KKD = KKD - KK*1461
C!
C!  CALCULATE JYY, THE REMAINING YEARS OF THE CURRENT CENTURY UP TO THE
C!  CURRENT DAY:
      JYY = 0
C!  THE NEXT YEAR IS NOT A LEAP YEAR IF JFY=0 AND JCC.NE.0.
      IF(JFY.EQ.0.AND.JCC.NE.0)GOTO 20
      IF(KKD.LT.366)GOTO 30
      JYY = 1
      KKD = KKD - 366
20    JYYY = KKD/365
      JYY = JYY + JYYY
      KKD = KKD - JYYY*365
30    IYY = 4*JFY + JYY
      IF(KKD.EQ.0) THEN
	 IYY=IYY-1
	 IMM=12
	 IDD=31
c	 RETURN
	go to 110
      END IF
C!
C!  SET L=1 IF WE HAVE A LEAP YEAR.
      L=0
      IF(IYY-(IYY/4)*4.NE.0)GOTO 40
      IF(IYY.EQ.0.AND.(ICC-(ICC/4)*4).NE.0)GOTO 40
      L=1
C!
C!  CALCULATE IMM AND IDD
40    IF(KKD.GT.31) GOTO 50
      IMM=1
      IDD=KKD
c      RETURN
       go to 110
C!
50    IF(KKD.GT.59)GOTO 60
      IMM = 2
      IDD = KKD-31
c      RETURN
	go to 110
C!
60    IF(KKD.GT.60)GOTO 70
      IF(L.EQ.0)GOTO 70
      IMM = 2
      IDD = 29
c      RETURN
	go to 110
C!
70    IF(L.EQ.1) KKD=KKD-1
      DO 80 I=4,13
	 IF(KKD.GT.NDP(I))GOTO 80
	 IMM = I-1
	 IDD = KKD - NDP(I-1)
c	 RETURN
	go to 110
C!
80    CONTINUE

 110 	iyear=icc*100+iyy
	return

90    WRITE(LP,5050)
5050  FORMAT(' ERROR IN DMY.')
      STOP
      END

	Subroutine Up_date(kd,ssec)
C 
C this subroutine updates the current gregorian day number kd
C and resets the time, ssec, accordingly
C it inherently assumes that the variable is updated every day!

      integer kd
      real ssec
	
      kd=kd+aint(ssec/8.64E4)
      ssec=amod(ssec,8.64E4)
 
      return
      end




C***********************************************************************
C*****************************************************************FUNDY5
C
      SUBROUTINE VERTAVGC(F,FBAR,Z,NN,NNV,NNDIM,NNVDIM)
      COMPLEX F(NNDIM,NNVDIM), FBAR(NNDIM)
      DIMENSION Z(NNDIM,NNVDIM)
C
C COMPUTES THE VERTICAL AVERAGE OF F, PUTS IT IN FBAR.
C DOES THIS FOR ALL HORIZONTAL NODES
C
      DO 100 I=1,NN
      FBAR(I)=(0.,0.)
      DO 200 J=2,NNV
      FBAR(I)=FBAR(I)+0.5*(F(I,J) + F(I,J-1))*(Z(I,J) - Z(I,J-1))
  200 CONTINUE
      FBAR(I)=FBAR(I)/(Z(I,NNV)-Z(I,1))
  100 CONTINUE
      RETURN
      END
C
C********************************************************************
C*********************************************************************
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 100 I=1,NN
      FINT=0.0
C
      DO 200 J=2,NNV
      FINT=FINT +0.5*(F(I,J) + F(I,J-1))*(Z(I,J) - Z(I,J-1))
  200 CONTINUE
      FAVG(I)=FINT/(Z(I,NNV) - Z(I,1))
C
  100 CONTINUE
      RETURN
      END
C

	subroutine loadvels(NN,IFLD,NFREQ,
     &               VELFILES,AMPU,AMPV,PHASEU,PHASEV)

C
C This subroutine loads in the requested velfiles and vertically 
C averages them
C
	include 'ACADIA2d.DIM'

C Global Variables
	integer nn,ifld,nfreq(nflddim)
	real ampu(nndim,nfreqdim),phaseu(nndim,nfreqdim)
	real ampv(nndim,nfreqdim),phasev(nndim,nfreqdim)
	character*72 velfiles(nflddim,nfreqdim)
C
C local variables
      real z(nndim,nnvdim)
	integer nnold, nnv
	complex ubar(nndim),vbar(nndim),u(nndim,nnvdim),v(nndim,nnvdim)
	complex w(nndim,nnvdim)
	nnold=nn

C
C loop over all velocity files for this field

	do k=1,nfreq(ifld)
	
c read in the vel file
	call loadvel(nndim,nnvdim,nn,nnv,velfiles(ifld,k),z,u,v,w)
	if (nnold.ne.nn) then
	write(*,*) 'problem with nodes in vel file'
	stop
	end if
	call vertavgc(u,ubar,z,nn,nnv,nndim,nnvdim)
	call vertavgc(v,vbar,z,nn,nnv,nndim,nnvdim)

	do i=1,nn
	ampu(i,k)=cabs(ubar(i))
	ampv(i,k)=cabs(vbar(i))
	phaseu(i,k)=phaselagd(ubar(i))
	phasev(i,k)=phaselagd(vbar(i))
	end do
	
	end do
	
 100  format (1x,a25,1x,i2,1x,i2,1x,i4,f12.4)

	return
	end


C --------------------------------------------------------------------
	SUBROUTINE READ_INPUT(CASNAM,filinp,nstage,
     &              DEG,
     &              DELT,NFLD,IDD,IMM,IYEAR,
     &              SSEC,NFREQ,KON,VELFILES)

	include 'ACADIA2d.DIM'
C
C This subroutine reads in a POP2D.inp file
C
C Global Variables
c	Integer nflddim,nfreqdim
	Integer nfld, nfreq(nflddim)
	Character*72 casnam, filinp
	character*72 velfiles(nflddim,nfreqdim)
	character*5 kon(nflddim,nfreqdim)
	REAL deg, delt, ssec(nflddim)
	integer idd(nflddim), imm(nflddim), iyear(nflddim), nstage

C Local Variables
      CHARACTER*72 CHACT1, COMINQ, LLABEL
	integer indx

      OPEN (UNIT=1, FILE=FILINP, STATUS='OLD')

      READ (1,600) LLABEL
      READ (1,600) COMINP
      READ (1,600) LLABEL
      READ (1,600) CHACT1
      ICAS = INDEX(CHACT1,' ') - 1
      CASNAM = CHACT1(1:ICAS)
      READ (1,600) LLABEL
      READ (1,*) nstage
	if (nstage.gt.nstagedim) then
	  write(*,*) 'error in nstagedim, plesae increase and recompile'
	  stop
	end if

      READ (1,*) DEG
      READ (1,*) DELT
	READ (1,*) NFLD
	if (nfld.gt.nflddim-1) then
	  write(*,*) 'error in nflddim, please increase and recompile'
	  stop
	end if

	DO IFLD=1,NFLD
        read(1,*) idd(IFLD),imm(IFLD),iyear(IFLD),
     &            ssec(IFLD),NFREQ(IFLD)
	  DO K=1,NFREQ(ifld)
	  READ(1,600) CHACT1  
	  INDX = INDEX(CHACT1,' ') - 1
        KON(IFLD,K) = CHACT1(1:INDX)
	  indx=indx+2
	  VELFILES(ifld,k)=chact1(indx:)
        END DO
      END DO
      read(1,*) idd(NFLD+1),imm(NFLD+1),iyear(NFLD+1),ssec(NFLD+1)
	NFREQ(NFLD+1)=0

      CLOSE (1)
	
	write(*,*) 'POP2D.inp read in'

 600  format (a)

	RETURN
	END


	SUBROUTINE VELCALC(nn,ifld,nfreq,
     &         ampu,ampv,phaseu,phasev,kon,kd0,ssec0,deg,u,v)

	include 'ACADIA2d.DIM' 

C Global variables
	integer nn, ifld,nfreq(nflddim),kd0
	real ampu(nndim,nfreqdim),ampv(nndim,nfreqdim)
      real phaseu(nndim,nfreqdim),phasev(nndim,nfreqdim)
      real ssec0,u(nndim),v(nndim),deg
	character*5 kon(nflddim,nfreqdim)


C local variables
      integer m, i,ifldold,k
      real sig(nfreqdim)
      real f(nfreqdim),vu(nfreqdim),cose
      character*5 konk, ddum	
      INTEGER kd0old, idum
      integer dayflag,kh
      REAL FAC, PI, twopi, ang, tm, dum
	
      SAVE kd0old, f, vu, sig, dayflag, m,ifldold, kh,twopi

	if (ifld.eq.1) then
	  pi=acos(-1.)
	  twopi=2.*pi
	  dayflag=0 
	  kd0old=kd0
	  ifldold=0
	  kh=kd0*24.
        call opnvuf(kh,ddum,deg,dum,dum,dum) 
	  ifld=1
	end if

      if (kd0.ne.kd0old) then
	  dayflag=dayflag+1
	  kd0old=kd0
	  if (dayflag.eq.30) then
          dayflag=0
	    kh=kd0*24.
        call setvuf(kh,ddum,deg,dum,dum,dum)
	  end if
	end if

	if (ifldold.ne.ifld) then
	  ifldold=ifld
          m=nfreq(ifld)
	  do k=1,m
	    konk=kon(ifld,k)
	    call vuf(kh,konk,deg,f(k),vu(k),sig(k))
	  end do
	end if


c tm is the time (in hours) passed since the f and vu 
c variables were updated

	tm=float(dayflag)*24.+ssec0/3600.

c synthesize tides

     	do i=1,nn
	    u(i)=0.
	    v(i)=0.
 	end do
      
	DO 811 K=1,m
	  sigk=sig(k)
	  vuk=vu(k)
	  fk=f(k)
	  do 810 i=1,nn
          ANG=sigK*tm + VUK - phaseu(i,K)/360d0
          ANG=TWOPI*(ANG-AINT(ANG))
          COSE=COS(ANG)
          u(i)=u(i)+fk*ampu(i,k)*cose
          ANG=sigK*tm + VUK - phasev(i,K)/360d0
          ANG=TWOPI*(ANG-AINT(ANG))
          COSE=COS(ANG)
          v(i)=v(i)+fk*ampv(i,k)*cose
 810	  continue

 811  continue

      return
	end


C***********************************************************************
c      SUBROUTINE VUF(kd0,KOnk,xlat,fk,vuk,freqk)
c taken from Mike Foreman's IOS subroutine Tide1.f
C alterations by Wendy Gentleman,   July 1995
C  1) ONE OR MORE RECORDS FOR EACH MAIN CONSTITUENT CONTAINING THE
C     FOLLOWING INFORMATION IN FORMAT (6X,A5,1X,6I3,F5.2,I4)
C      KON    = CONSTITUENT NAME
C   II,JJ,KK,LL,MM,NN = THE SIX DOODSON NUMBERS
C      SEMI   = PHASE CORRECTION
C      NJ     = THE NUMBER OF SATELLITES FOR THIS CONSTITUENT.
C   IF NJ>0, INFORMATION ON THE SATELLITE CONSTITUENTS IS READ , THREE
C   SATELLITES PER CARD, IN THE FORMAT (11X,3(3I3,F4.2,F7.4,1X,I1,1X)).
C   FOR EACH SATELLITE THE VALUES READ ARE
C      LDEL,MDEL,NDEL = THE CHANGES IN THE LAST THREE DOODSON NUMBERS
C                       FROM THOSE OF THE MAIN CONSTITUENT.
C      PH     = THE PHASE CORRECTION
C      EE     = THE AMPLITUDE RATIO OF THE SATELLITE TIDAL POTENTIAL TO
C               THAT OF THE MAIN CONSTITUENT.
C      IR     = 1 IF THE AMPLITUDE RATIO HAS TO BE MULTIPLIED BY THE
C                 LATITUDE CORRECTION FACTOR FOR DIURNAL CONSTITUENTS
C               2 IF THE AMPLITUDE RATIO HAS TO BE MULTIPLIED BY THE
C                 LATITUDE CORRECTION FACTOR FOR SEMI-DIURNAL CONSTI-
C                 TUENTS.
C               OTHERWISE IF NO CORRECTION IS REQUIRED TO THE AMPLITUDE
C                 RATIO.
C
C  2) ONE BLANK RECORD
C
C  3) ONE RECORD FOR EACH SHALLOW WATER CONSTITUENT SPECIFYING IN FORMAT
C     (11X,3(3I3,F4.2,F7.4,1X,I1,1X)), THE FOLLOWING INFORMATION
C      KON    = NAME  OF THE SHALLOW WATER CONSTITUENT
C      NJ     = NUMBER OF MAIN CONSTITUENTS FROM WHICH IT IS DERIVED.
C      COEF,KONCO = COMBINATION NUMBER AND NAME OF THESE MAIN
C                   CONSTITUENTS.
C
C  4) ONE BLANK RECORD
C
C=======================================================================
C
C*  MTOT           IS THE NUMBER OF CONSTITUENTS INCLUDED IN THE DATA
C*                 FILE. AT PRESENT THIS IS 146.
C*  M              IS THE NUMBER OF CONSTITUENTS INCLUDED IN THE
C*                 ANALYSIS. AT PRESENT IT WILL HAVE A MAXIMUM OF 69 +
C                  THE NUMBER OF SHALLOW WATER CONSTITUENTS SPECIALLY
C*                 DESIGNATED FOR INCLUSION (NOT IN THE STANDARD 69
C*                 CONSTITUENT PACKAGE) IN THE ANALYSIS.
C*  KONTAB(I)      IS THE ARRAY CONTAINING ALL THE CONSTITUENT NAMES AS
C*                 THEY ARE READ IN FROM THE DATA FILE. IT SHOULD HAVE
C*                 THE MINIMUM DIMENSION MTOT.
C*  KON(I)         IS THE ARRAY CONTAINING THE NAMES OF ALL THE
C*                 CONSTITUENTS THAT ARE TO BE INCLUDED IN THE LEAST
C*                 SQUARES ANALYSIS TO THE TIDAL RECORD. THEY ARE CHOSEN
C*                 FROM THE KONTAB LIST VIA THE RAYLEIGH CRITERION. THE
C*                 ARRAY SHOULD HAVE DIMENSION AT LEAST M.
C*  FREQ(I)        IS THE ARRAY OF FREQUENCIES (CYCLES/ HR.) CORRESPOND-
C*                 ING TO THE CONSTITUENTS CONTAINED IN KONTAB. IT
C*                 SHOULD HAVE DIMENSION AT LEAST MTOT.
C*  SIG(I)         IS THE ARRAY OF FREQUENCIES CORRESPONDING TO THE
C*                 CONSTITUENTS CONTAINED IN KON. IT SHOULD HAVE MINIMUM
C*                 DIMENSION M.
C*  V              ASTRONOMICAL PHASE ARGUMENT
C*  U & F          NODAL MODULATION &PHASE AND AMPLITUDE CORRECTIONS
C
C Right now these (v, u & f) are calculated mtot constituents
C
C MINIMAL DIMENSION REQUIREMENTS
C ==============================
C THE DIMENSION OF
C   KON AND NJ SHOULD BE AT LEAST EQUAL TO THE TOTAL NUMBER OF POSSIBLE
C   CONSTITUENTS PLUS ONE (PRESENTLY 146+1), THE DIMENSION OF II, JJ,
C   KK, LL, MM, NN, AND SEMI SHOULD BE AT LEAST EQUAL TO THE NUMBER OF
C   MAIN CONSTITUENTS PLUS ONE (PRESENTLY 45+1), THE DIMENSION OF EE,
C   LDEL, MDEL, NDEL, IR, AND PH SHOULD BE AT LEAST EQUAL TO THE TOTAL
C   NUMBER OF SATELLITES TO ALL THE MAIN CONSTITUENTS PLUS THE NUMBER
C   OF CONSTITUENTS WITH NO SATELLITES (PRESENTLY 162+8 IF THERE ARE NO
C   THIRD ORDER SATELLITES FOR BOTH N2 AND L2, AND 167+8 IF ALL THIRD
C   ORDER SATELLITES ARE INCLUDED FOR BOTH N2 AND L2),
C   AND THE DIMENSION OF KONCO, AND COEFF SHOULD BE AT LEAST EQUAL TO
C   THE SUM FOR ALL SHALLOW WATER CONSTITUENTS OF THE NUMBER OF MAIN
C   CONSTITUENTS FROM WHICH EACH IS DERIVED,PLUS FOUR (PRESENTLY 251+4).
C=======================================================================
C
      SUBROUTINE VUF(kh,KOnk,xlat,fk,vuk,freqk)
      CHARACTER*5 KONk,KBLANK,KONCO,kontab(170)
      REAL freq(170),f(170),vu(170),fk,vuk,sigk
      real*8 d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp,hh,tau,dtau
      real EE(205),PH(205),SEMI(50)
      real xlat
      integer kh, line
      integer II(50),JJ(50),KK(50),LL(50),MM(50),NN(50),indx
      integer LDEL(205),MDEL(205),NDEL(205),IR(205),NJ(170)
      DIMENSION KONCO(320),COEF(320)
      DATA KBLANK/'     '/

      Save kontab, freq, konco, ntotal
     &              ee, ph, semi, ii, jj, kk, ll, mm, nn,
     &             ldel, mdel, ndel, ir, nj, coef

C
c in opnvuf - reads in kontab and data and falls to setvuf
c in setvuf - evaluates f and vu for all constituents in kontab
c in vuf - finds appropriate f,vu and sig for a specified constituent

c determines the index for the specified constituents
        do  209 l=1,ntotal
	  if(konk.eq.kontab(l)) then
	    indx=l
          go to 211 
         end if
 209	continue
	  write(*,*) 'specified constituent',konk,'is not in file'
	  stop
 211    continue

        fk=f(indx)
        vuk=vu(indx)
        freqk=freq(indx)

c        write(*,*) konk,vuk*360.,fk,freqk
 98   format(' ',a5,5x,2f12.5)

	
      RETURN


      ENTRY OPNVUF(Kh,KONk,XLAT,fk,vuk,freqk)	

      KR=8
      PI=3.1415926536
      TWOPI=2.*PI

      open(unit=kr,file='IOS_tidetbl',status='old',form='formatted') 

C
C***********************************************************************
C*  HERE THE MAIN CONSTITUENTS AND THEIR DOODSON NUMBERS ARE READ IN
C*  FORMAT (6X,A5,1X,6I3,F5.2,I4). THE VALUES ARE RESPECTIVELY
C*     KON    = CONSTITUENT NAME
C*  II,JJ,KK,LL,MM,NN = THE SIX DOODSON NUMBERS
C*     SEMI   = PHASE CORRECTION
C*     NJ     = THE NUMBER OF SATELLITES FOR THIS CONSTITUENT.
C*  THE END OF ALL MAIN CONSTITUENTS IS DENOTED BY A BLANK CARD.
C

      JBASE=0
      DO 90 K=1,1000
      READ(KR,60)KONTAB(K),II(K),JJ(K),KK(K),LL(K),
     & MM(K),NN(K),SEMI(K),NJ(K)
   60 FORMAT(6X,A5,1X,6I3,F5.2,I4)
      IF(KONTAB(K).eq.KBLANK) go to 100
70    J1=JBASE+1
      IF(NJ(K).GE.1) GO TO 75
      NJ(K)=1
      JL=J1
      PH(J1)=0.
      EE(J1)=0.
      LDEL(J1)=0
      MDEL(J1)=0
      NDEL(J1)=0
      IR(J1)=0
      GO TO 90
   75 JL=JBASE+NJ(K)
C
C***********************************************************************
C*  IF NJ>0, INFORMATION ON THE SATELLITE CONSTITUENTS IS READ , THREE
C*  SATELLITES PER CARD, IN THE FORMAT (11X,3(3I3,F4.2,F7.4,1X,I1,1X)).
C*  FOR EACH SATELLITE THE VALUES READ ARE
C*     LDEL,MDEL,NDEL = THE CHANGES IN THE LAST THREE DOODSON NUMBERS
C*                      FROM THOSE OF THE MAIN CONSTITUENT.
C*     PH     = THE PHASE CORRECTION
C*     EE     = THE AMPLITUDE RATIO OF THE SATELLITE TIDAL POTENTIAL TO
C*              THAT OF THE MAIN CONSTITUENT.
C*     IR     = 1 IF THE AMPLITUDE RATIO HAS TO BE MULTIPLIED BY THE
C*                LATITUDE CORRECTION FACTOR FOR DIURNAL CONSTITUENTS
C*              2 IF THE AMPLITUDE RATIO HAS TO BE MULTIPLIED BY THE
C*                LATITUDE CORRECTION FACTOR FOR SEMI-DIURNAL CONSTI-
C*                TUENTS.
C*              OTHERWISE IF NO CORRECTION IS REQUIRED TO THE AMPLITUDE
C*                RATIO.
C

      READ(KR,80)(LDEL(J),MDEL(J),NDEL(J),PH(J),EE(J),IR(J),J=J1,JL)
   80 FORMAT((11X,3(3I3,F4.2,F7.4,1X,I1,1X)))

90    JBASE=JL
100   NTIDAL=K-1		!for us this will be 45

C
C***********************************************************************
C*  THE SHALLOW WATER CONSTITUENTS AND THE MAIN CONSTITUENTS FROM WHICH
C*  THEY ARE DERIVED ARE READ IN HERE WITH THE FORMAT
C*  (6X,A5,I1,2X,4(F5.2,A5,5X)). THE VALUES ARE RESPECTIVELY
C*     KON    = NAME  OF THE SHALLOW WATER CONSTITUENT
C*     NJ     = NUMBER OF MAIN CONSTITUENTS FROM WHICH IT IS DERIVED.
C*     COEF,KONCO = COMBINATION NUMBER AND NAME OF THESE MAIN
C*                  CONSTITUENTS.
C*  THE END OF THESE CONSTITUENTS IS DENOTED BY A BLANK CARD.
C
      JBASE=0
      K1=NTIDAL+1
      DO 160 K=K1,1000
      J1=JBASE+1
      J4=J1+3

      READ(KR,130)KONTAB(K),NJ(K),(COEF(J),KONCO(J),J=J1,J4)
  130 FORMAT(6X,A5,I1,2X,4(F5.2,A5,5X))
      IF(KONTAB(K).eq.KBLANK) go to 170
160   JBASE=JBASE+NJ(K)
170   NTOTAL=K-1		! for us this will be 146
	close(kr)

C
C***********************************************************************
C*  NTIDAL IS THE NUMBER OF MAIN CONSTITUENTS
C*  NTOTAL IS THE NUMBER OF CONSTITUENTS (MAIN + SHALLOW WATER)
C*  FOR  THE GIVEN TIME KH, THE TABLE OF F AND V+U VALUES IS
C*  CALCULATED FOR ALL THE CONSTITUENTS.
C*     F IS THE NODAL MODULATION ADJUSTMENT FACTOR FOR AMPLITUDE
C*     U IS THE NODAL MODULATION ADJUSTMENT FACTOR FOR PHASE
C*     V IS THE ASTRONOMICAL ARGUMENT ADJUSTMENT FOR PHASE.
C

       ENTRY SETVUF(kh,KONk,XLAT,fk,vuk,freqk)   

      SLAT=SIN(PI*XLAT/180.)	

C
C***********************************************************************
C*  THE ASTRONOMICAL ARGUMENTS ARE CALCULATED BY LINEAR APPROXIMATION
C*  AT THE MID POINT OF THE ANALYSIS PERIOD.
C*     DS,DH,DP,DNP,DPP ARE THEIR RESPECTIVE RATES OF CHANGE OVER A 365
C*     DAY PERIOD AS OF 000 ET 1/1/1976.
       
      d1=kh/24.d0
      call gday(31,12,1899,kd)
      d1=d1-dfloat(kd)-0.5d0
      call astr(d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp)
      INT24=24
      INTDYS=KH/INT24
      HH=dfloat(KH-INTDYS*INT24)
      TAU=HH/24.D0+H-S
      dtau=365.d0+dh-ds


C
C***********************************************************************
C*  ONLY THE FRACTIONAL PART OF A SOLAR DAY NEED BE RETAINED FOR COMPU-
C*  TING THE LUNAR TIME TAU.
C

      JBASE=0
      DO 210 K=1,ntidal

        FREQ(k)=(II(K)*DTAU+JJ(K)*DS+KK(K)*DH+LL(K)*DP+MM(K)*DNP+
     1          NN(K)*DPP)/(365.*24.)

        VDBL=II(k)*TAU+JJ(k)*S+KK(k)*H+LL(k)*P+
     &          MM(k)*ENP+NN(k)*PP+SEMI(k)

        IV=VDBL
        IV=(IV/2)*2
        V=VDBL-IV
        J1=JBASE+1
        JL=JBASE+NJ(k)
        SUMC=1.
        SUMS=0.

        DO 200 J=J1,JL
C
C***********************************************************************
C*  HERE THE SATELLITE AMPLITUDE RATIO ADJUSTMENT FOR LATITUDE IS MADE
C

          RR=EE(J)
          L=IR(J)+1
          GO TO (901,902,903),L
 902      RR=EE(J)*0.36309*(1.-5.*SLAT*SLAT)/SLAT
          GO TO 901
 903      RR=EE(J)*2.59808*SLAT
 901      CONTINUE
          UUDBL=LDEL(J)*P+MDEL(J)*ENP+NDEL(J)*PP+PH(J)
          IUU=UUDBL
          UU=UUDBL-IUU
          SUMC=SUMC+RR*COS(UU*TWOPI)
          SUMS=SUMS+RR*SIN(UU*TWOPI)

  200   CONTINUE

        F(K)=SQRT(SUMC*SUMC+SUMS*SUMS)
        VU(K)=V+ATAN2(SUMS,SUMC)/TWOPI
	
210   JBASE=JL

C
C***********************************************************************
C*  HERE F AND V+U OF THE SHALLOW WATER CONSTITUENTS ARE COMPUTED FROM
C*  THE VALUES OF THE MAIN CONSTITUENT FROM WHICH THEY ARE DERIVED.
C

      JBASE=0
      K1=NTIDAL+1
      IF(k1.gt.NTOTAL) RETURN
      DO 270 K=K1,NTOTAL
        F(K)=1.0
        VU(K)=0.0
        FREQ(k)=0.
        iflag=1
	go to 268
 268  J1=JBASE+1
      JL=JBASE+NJ(K)
      DO 260 J=J1,JL
        KM1=K-1
        DO 240 L=1,KM1
          IF(KONTAB(L).eq.KONCO(J)) go to 250
 240    CONTINUE
        WRITE(LP,241)KONCO(J)
 241    FORMAT('   SETVUF STOP.',A5)
        STOP
 250  F(K)=F(K)*F(L)**ABS(COEF(J))
      VU(K)=VU(K)+COEF(J)*VU(L)
      if (iflag.eq.1) FREQ(k)=FREQ(k)+COEF(J)*FREQ(l)
 260   continue
 270  JBASE=JL

      RETURN
      END




	subroutine astr(d1,h,pp,s,p,np,dh,dpp,ds,dp,dnp)
c	this subroutine calculates the following five ephermides
c	of the sun and moon
c	h = mean longitude of the sum
c	pp = mean longitude of the solar perigee
c	s = mean longitude of the moon
c	p = mean longitude of the lunar perigee
c	np = negative of the longitude of the mean ascending node
c	and their rates of change.
c	Units for the ephermides are cycles and for their derivatives
c	are cycles/365 days
c	The formulae for calculating this ephermides were taken from
c	pages 98 and 107 of the Explanatory Supplement to the
c	Astronomical Ephermeris and the American Ephermis and
c	Nautical Almanac (1961)
c
	implicit real*8(a-h,o-z)
	real*8 np
	d2=d1*1.d-4
	f=360.d0
	f2=f/365.d0
	h=279.696678d0+.9856473354d0*d1+.00002267d0*d2*d2
	pp=281.220833d0+.0000470684d0*d1+.0000339d0*d2*d2+
     1  .00000007d0*d2**3
	s=270.434164d0+13.1763965268d0*d1-.000085d0*d2*d2+
     1  .000000039d0*d2**3
	p=334.329556d0+.1114040803d0*d1-.0007739d0*d2*d2-
     1  .00000026d0*d2**3
	np=-259.183275d0+.0529539222d0*d1-.0001557d0*d2*d2-
     1  .00000005d0*d2**3
	h=h/f
	pp=pp/f
	s=s/f
	p=p/f
	np=np/f
	h=h-dint(h)
	pp=pp-dint(pp)
	s=s-dint(s)
	p=p-dint(p)
	np=np-dint(np)
	dh=.9856473354d0+2.d-8*.00002267d0*d1
	dpp=.0000470684d0+2.d-8*.0000339d0*d1
     1  +3.d-12*.00000007d0*d1**2
	ds=13.1763965268d0-2.d-8*.000085d0*d1+
     1  3.d-12*.000000039d0*d1**2
	dp=.1114040803d0-2.d-8*.0007739d0*d1-
     1  3.d-12*.00000026d0*d1**2
	dnp=+.0529539222d0-2.d-8*.0001557d0*d1-
     1  3.d-12*.00000005d0*d1**2
	dh=dh/f2
	dpp=dpp/f2
	ds=ds/f2
	dp=dp/f2
	dnp=dnp/f2
	return
	end


C****67***************************************************************72
C***********************************************************************
      SUBROUTINE CSTRING(string,istart,iend)
C-----------------------------------------------------------------------
C purpose: This subroutine finds the first nonblank space (istart) in a 
C            arbitrary character string and the first subsequent blank
C            space (iend).  The range of string between istart and iend
C            can then be separated from the rest of string by specifying
C            string(istart:iend)
C
C
C inputs:  string - a character string of length .LE. 72
C
C outputs: istart - first nonblank character in string
C          iend   - first blank character in string after istart
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           12 JUNE 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      character*72 string
C
C START OF EXECUTABLE CODE
C
C find beginning of string
      do 10 i=1,72
         istart=i
	 if(string(i:i).NE.' ')go to 11
 10   continue
C
C find end of string
 11   do 20 i=istart,72
         iend=i-1
	 if(string(i:i).EQ.' ')go to 21
 20   continue
C
C END OF ROUTINE
 21   return
      end
C***********************************************************************
C***********************************************************************
C*********************************************************************
C*********************************************************************
C
       SUBROUTINE SPRSBLD3DD(IQ,JQ,NN,NDF)
C
C BUILD 3D IDEX IQQ AND JQQ
C
       DIMENSION IQ(1),JQ(1)
	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
      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(1)
      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***********************************************************************
C
C    Subroutine: INTTOCHR (INTG, FILPFX, IPOS)
C   Description: This subroutine converts integer INTG (maximum of I20)
C                to character string FILPFX(1:IPOS) (maximum of 20
C                characters).
C     Variables: [Argument list]
C                INTG -- Integer (maximum I20) to be converted
C                FILPFX -- Character string of maximum 20 characters
C                        converted from integer INTG
C                IPOS -- Position of the last meaningful (non-blank)
C                        character in FILPFX
C                [Local integer variables]
C                CBLANK -- The blank space character
C                I -- Do loop index
C                INDX -- Position of the first blank space at the tail
C                        of FILPFX before editting out the blanks
C                LMAX -- Maximum number of characters allowed in the
C                        converted character string
C         Input: INTG
C        Output: FILPFX, IPOS
C
C***********************************************************************
C23456789012345678901234567890123456789012345678901234567890123456789012
C________1_________2_________3_________4_________5_________6_________7_|
C
C...Argument list of the subroutine
C
      SUBROUTINE INTTOCHR (INTG, FILPFX, IPOS)
C
C...Declare no implicit variable type
C
      IMPLICIT NONE
C
C...Variable declaration of the argument list
C
      CHARACTER*20 FILPFX
      INTEGER INTG, IPOS
C
C...Declaration of the local variables
C
      CHARACTER*1 CBLANK
C###      INTEGER I, INDX
      INTEGER I, LMAX
C
C...Define the maximum length of the character string
C
      LMAX = 20
C
C...Define a blank space character
C
      CBLANK = ' '
C
C...Convert an I20 integer INTG to a CHARACTER*72 string FILPFX using
C   interal write statement
C
      WRITE (FILPFX,'(I20)') INTG
C
C...Edit out all the blank spaces from the beginning of FILPFX and
C   reload the meaningful string back as FILPFX(1:IPOS), where IPOS
C   is the position of the last meaningful (non-blank) character in
C   FILPFX
C
      DO 10 I=1,LMAX
         IF (FILPFX(I:I) .NE. CBLANK) GOTO 20
   10 CONTINUE
      WRITE (*,*) 'Invalid file prefix for output files'
      STOP 1111
   20 CONTINUE
      IPOS = LMAX - I + 1
      FILPFX = FILPFX(I:LMAX)
C
C...Delete the tailing blanks of the string
C
C###      IPOS = INDEX(FILPFX(I:72),CBLANK)
C###      INDX = IPOS + I - 1
C###      FILPFX = FILPFX(I:INDX-1)
C###      IPOS = IPOS - 1
C
C...Terminate subroutine INTTOCHR and return to the calling program
C
      RETURN
      END


	subroutine sparsestor(nn,ne,in,x,y,iter,bathy,diff,
     &                      iq,jq,sv,ppx,ppy,agpgp)
	include 'ACADIA2d.DIM'
      INTEGER IN(3,NEDIM)
	integer nftr, nn, ne, iter, dtype
	real x(nndim), y(nndim)
	real DX(3,NEDIM),DY(3,NEDIM), AR(NEDIM)
      INTEGER IQ(NNDIM), JQ(NFTRDIM*NNDIM)
      REAL PPX(NFTRDIM*NNDIM), PPY(NFTRDIM*NNDIM), sv(nndim)
      REAL AGPGP(NFTRDIM*NNDIM)
	real diff(nndim), bathy(nndim),dhav(nedim)

	save dx, dy, ar, nftr

	if (iter.ne.0) go to 100

C initalize and assemble PPX and PPY
	do i=1,nn
	sv(i)=0.
	end do

      DO 30 I=1,NFTRDIM*NNDIM
      PPX(I)  = 0.0
      PPY(I)  = 0.0
 30   CONTINUE      

C  SPRSBLD1 is used instead of SPRSBLD for efficiency (7.8.93)

      CALL SPRSBLD1 (IN, 3, 3, IQ, JQ, NN, NE, NFTRDIM, 1)
      WRITE (*,*) 'Subroutine SPRSBLD1 executed'
C
      NFTR=IQ(NN)/NN +1
      IF(NFTR.GT.NFTRDIM)THEN
      write(*,*) 'NFTR, NFTRDIM=',NFTR,NFTRDIM,';EXECUTION TERMINATED'
      STOP 
      END IF

C  COMPUTE ELEMENT AREAS, DELX, DELY
C                
      DO 40 L=1,NE       
      I1=IN(1,L)          
      I2=IN(2,L)          
      I3=IN(3,L) 
      DX(1,L)=X(I2)-X(I3)          
      DX(2,L)=X(I3)-X(I1)          
      DX(3,L)=X(I1)-X(I2)          
      DY(1,L)=Y(I2)-Y(I3)          
      DY(2,L)=Y(I3)-Y(I1)          
      DY(3,L)=Y(I1)-Y(I2)          
      AR(L) = 0.5*(X(I1)*DY(1,L)+X(I2)*DY(2,L)+X(I3)*DY(3,L))   
      IF(AR(L).LE.0.0) WRITE(2,*) ' NEGATIVE AREA IN ELEMENT', L 
  40  CONTINUE
C                
C  COMPUTE DIAGONAL MASS MATRIX SV(I);
C  BANDED DEPTH MATRIX SH(J,i) (transposed)
C  AND DERIVATIVE MATRICES IN SPRSPAK FORM:
C     PPX = <Pi*dP/dxj>
C     PPY = <Pi*dP/dyj>

      DO 50 L=1,NE       
      ARR=AR(L)       
      DO 60 I=1,3        
      II=IN(I,L)          
      SV(II)=SV(II)+ARR/3.0
      DO 70 J=1,3
      IJ=IN(J,L)
      K=KAY(II,IJ,IQ,JQ)
      PPX(K)=PPX(K)+DY(J,L)/6.   
	PPY(K)=PPY(K)-DX(J,L)/6. 

  70  CONTINUE
  60  CONTINUE
  50  CONTINUE

C
C  PREMULTIPLY PPX AND PPY BY INVERSE OF SV
C
      CALL SPRSINVMLT (PPX, IQ, SV, NN)
      CALL SPRSINVMLT (PPY, IQ, SV, NN)
	
	write(*,*) 'ppx and ppy assembled'

 100   Continue

C Compute Laplacian

C Laplacian matrix 
C     AGPGP = -<AH*grad(phi_i) . grad(phi_j)>
C
      DO 156 I=1,NFTR*NN
      AGPGP(I)= 0.0
  156 CONTINUE         

	do l=1,ne
       I1=IN(1,L)          
       I2=IN(2,L)          
       I3=IN(3,L) 
   	    dhav(l)=(diff(i1)*bathy(i1)+
     &         diff(i2)*bathy(i2)+diff(i3)*bathy(i3))/3.   
	end do

      DO L=1,NE       
      ARR=AR(L) 
	do i=1,3
      II=IN(I,L)
      DO J=1,3
      IJ=IN(J,L)
      K=KAY(II,IJ,IQ,JQ)
      AGPGP(K) = AGPGP(K) - dhav(l)*(DX(J,L)*DX(I,L)
     &                             + DY(J,L)*DY(I,L))/4./ARR

	end do
	end do
      end do

      CALL SPRSINVMLT (AGPGP, IQ, SV, NN)         

	write(*,*) 'agpgp assembled'

      RETURN
      END

C*********************************************************************
C
      SUBROUTINE SPRSMLTIS(I,QV,IQ,JQ,X,B,NSTAGE,NNDIM,NSTAGEDIM)
C
C  MULTIPLIES THE Ith ROW OF QV BY THE Jth COLUMN OF THE MATRIX X,
C     FOR ALL VALUES OF J FROM 1 TO NNV
C  ANSWER RETURNS IN THE VECTOR B
C     (OR: OBTAINS ONLY Bi OF Q*X=B)
C  IQ,JQ ARE THE STANDARD POINTERS FOR QV
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF X
C  THIS ROUTINE IS NONDESTRUCTIVE
C
      DIMENSION QV(1)
      dimension IQ(1),JQ(1),X(NSTAGEDIM,NNDIM)
      dimension B(nstagedim)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
		
      DO ISTAGE=1,NSTAGE
      B(ISTAGE)=0.
	end do

      DO 10 K=KMIN,KMAX
	qvk=qv(k)
	jqk=jq(k)
	do istage=1,nstage
	xij=X(ISTAGE,JQK)
      B(ISTAGE)=B(ISTAGE)+QVK*Xij
	end do

   10 CONTINUE
      RETURN
      END
C  
C*********************************************************************
C
      SUBROUTINE SPRSMLTI(I,QV,IQ,JQ,X,B)
C
C  MULTIPLIES THE Ith ROW OF QV BY THE VECTOR X; ANSWER RETURNS IN THE SCALAR B
C     (OR: OBTAINS ONLY Bi OF Q*X=B)
C  IQ,JQ ARE THE STANDARD POINTERS FOR QV
C  THIS ROUTINE IS NONDESTRUCTIVE
C
      DIMENSION QV(1),IQ(1),JQ(1),X(1)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      B=0.
      DO 20 K=KMIN,KMAX
   20 B=B+QV(K)*X(JQ(K))
      RETURN
      END
C  
C***********************************************************************
C***********************************************************************
      SUBROUTINE LOADS2C(nndim,nn,filename,scval,freq)
C-----------------------------------------------------------------------
C purpose: This routine loads a .s2c file
C
C index range: i=1,nn
C
C inputs:  nndim - parameter used to dimension node arrays in the
C            calling routine
C          nn - number of nodes for which data should exist
C          filename - name of file to be read
C
C outputs: scval(nnn) - complex scalar value at node nnn
C          freq - frequency in radians/sec
C 
C file format: See CEN 3/29/93 Memorandum
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           30 JULY 1992
C           Revised: 25 APRIL 1994
C-----------------------------------------------------------------------
C
C ARGUMENTS
      complex scval(nndim)
      character*72 filename,header
C
C START OF EXECUTABLE CODE
C
C Open file
      open(63,file=filename,status='old')
C
C Read headers
      read(63,101)header
      read(63,101)header
      read(63,*)freq
C
C Check for oversize of arrays
      nflag=0
      if(nn.GT.nndim)then
	 write(*,*)'Dimensioning problem in SUBROUTINE LOADS2C:'
         write(*,*)'  nn exceeds nndim => rectify and recompile'
         nflag=1
      endif
      if(nflag.EQ.1)stop
C
C Read data
      fac=2.0*ATAN2(1.0,0.0)/180.0
      do 10 i=1,nndim
         read(63,*,end=11)node,amp,phase
         scval(i)=CMPLX(amp*COS(-phase*fac),amp*SIN(-phase*fac))
         ii=i
 10   continue
 11   close(63)
C
C Message to standard out
      if(ii.NE.nn)then
         write(*,*)'There is a mismatch between the number of horizontal
     & nodes in your .s2c file and nn in SUBROUTINE LOADS2C.'
         STOP
      else
         write(*,*)ii,' complex 2-D scalar values loaded from ',filename
      endif
C
C FORMAT STATEMENTS
 101  format(a)
C
C END OF ROUTINE
      return
      end
