C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
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***********************************************************************
C
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 LOADS3R(nndim,nnvdim,nn,nnv,filename,scval,z)
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
c      include 'oaci.dim'
      real scval(nndim,nnvdim),z(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
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 i=1,nn
        do j=1,nnv
         read(63,*)node,z(node,j),scval(node,j)
        enddo
      enddo

 11   close(63)
C
C FORMAT STATEMENTS
 101  format(a)
C
C END OF ROUTINE
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc     
      SUBROUTINE LOADmir3d(nobsdim,nobs,filename,tobs,xobs,yobs,zobs,
     &sobs,stderr,inst)
C-----------------------------------------------------------------------
C Load a 3d scalar observation file with measurment stderr and instrument
C number.   
C-----------------------------------------------------------------------
      real tobs(nobsdim),xobs(nobsdim),yobs(nobsdim),zobs(nobsdim),
     &sobs(nobsdim),stderr(nobsdim)
      integer nobs,nobsdim,inst(nobsdim)
      character*72 filename,header
C Open file
      open(63,file=filename,status='old')
      read(63,*)header
      read(63,*)year
C Read data
      do 10 i=1,nobsdim
         read(63,*,end=11) tobs(i),xobs(i),yobs(i),zobs(i),sobs(i),
     &stderr(i),inst(i)
         nobs=i
 10   continue
 11   close(63)

      return
      end
C***********************************************************************
C***********************************************************************
C
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
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  Fixed subroutines from FEMPAK
C  Fixed subroutines from FEMPAK
C  Fixed subroutines from FEMPAK
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE BASIS2D(nndim,nedim,in,x,y,nn,ne,xval,yval,
     &nclosest,nelem,basis)
C-----------------------------------------------------------------------
C purpose: This subroutine determines the element on which the point
C            (xval,yval) resides and determines the value of the basis
C            functions at the point.  In the first pass, only elements 
C            containing the closest node to (xval,yval) are tested to 
C            see if (xval,yval) lies on them.  If this fails, all
C            elements are tested to see if (xval,yval) lies on them.
C
C
C restrictions: Applicable only for  2-D triangular linear elements
C
C calling sequence: CALL CLOSESTNODE(nndim,x,y,nn,xval,yval,nclosest)
C                   CALL POINTONEL(xlocal,ylocal,xval,yval,nflag)
C                   CALL PHI2D1(xlocal,ylocal,xval,yval,basis)
C
C inputs:  nn is the number of nodes in the mesh
C          x(nn),y(nn) 1-D arrays containing the are the x and y
C            coordinates of each nodes
C          ne is the number of elements in the mesh
C          in(3,ne) is the LABLIB Standard incidence list
C          xval,yval is the location of the point 
C          nedim is the size of the arrays in the calling routine
C
C outputs: basis is a 1-D array containing the value of the 3 basis
C            functions at the point xval,yval.
C          nelem is the element on which (xval,yval),resides
C          nclosest is the closest node to (xval,yval)
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           1992 APRIL 22 - Original Version
C           1997 MARCH 01 - Revised to call PHI2D1 instead of PHI2D to
C                           determine basis functions.  PHI2D1 is an
C                           improvement over PHI2D in that it less
C                           subject to round-off errors.
C-----------------------------------------------------------------------
C
C ARGUMENTS
      integer in(3,nedim)
      real x(nndim),y(nndim),basis(3)
C
C LOCAL VARIABLES
      dimension xlocal(3),ylocal(3)
C
C START OF EXECUTABLE CODE
C
C  Find the closest node to (xval,yval)
      CALL CLOSESTNODE(nndim,x,y,nn,xval,yval,nclosest)
C
C  Find if an element containing the closest node also contains the
C    point (xval,yval)
C
      do 10 i=1,ne
         n1=in(1,i)
         n2=in(2,i)
         n3=in(3,i)
         if(nclosest.EQ.n1.OR.nclosest.EQ.n2.OR.nclosest.EQ.n3)then
            xlocal(1)=x(n1)
            ylocal(1)=y(n1)
            xlocal(2)=x(n2)
            ylocal(2)=y(n2)
            xlocal(3)=x(n3)
            ylocal(3)=y(n3)
            CALL POINTONEL(xlocal,ylocal,xval,yval,nflag)
            if(nflag.EQ.1) then
               nelem=i
               CALL PHI2D1(xlocal,ylocal,xval,yval,basis)
               return
            endif
         endif
 10   continue
C
C  Find if any element contains the point (xval,yval)
C
      do 20 i=1,ne
         xlocal(1)=x(in(1,i))
         ylocal(1)=y(in(1,i))
         xlocal(2)=x(in(2,i))
         ylocal(2)=y(in(2,i))
         xlocal(3)=x(in(3,i))
         ylocal(3)=y(in(3,i))
         CALL POINTONEL(xlocal,ylocal,xval,yval,nflag)
         if(nflag.EQ.1) then
            nelem=i
            CALL PHI2D1(xlocal,ylocal,xval,yval,basis)
            return
         endif
 20   continue
      nelem=0
C
C END OF ROUTINE
C
      return
      end
C****67***************************************************************72
C***********************************************************************
C***********************************************************************
C***********************************************************************
      SUBROUTINE CLOSESTNODE(nndim,x,y,nn,xval,yval,nclosest)
C-----------------------------------------------------------------------
C purpose: This subroutine finds the node closest to a point specified
C            by the coordinates (xval,yval)
C  
C restrictions:  all nodes must lie in a flat plane
C
C inputs:  x(nn) and y(nn) are 1-D arrays containing the x and y
C            coordinates of each node
C          nn is the number of nodes
C          xval,yval are the point coordinates
C
C outputs: nclosest is the node who's (x,y) coordinates are closest to
C            (xval,yval)
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension x(nndim),y(nndim)
C
C START OF EXECUTABLE CODE
C
C  Find node closest to (xval,yval)
      dclosest=SQRT((x(1)-xval)**2.0+(y(1)-yval)**2.0)
      nclosest=1
      do 25 j=2,nn
         dval=SQRT((x(j)-xval)**2.0+(y(j)-yval)**2.0)
         if(dval.LT.dclosest)then
            dclosest=dval
            nclosest=j
         endif
 25   continue
C
C END OF ROUTINE
      return
      end
C***********************************************************************
C***********************************************************************
C****67***************************************************************72
C***********************************************************************
      SUBROUTINE POINTONEL(xlocal,ylocal,xval,yval,nflag)
C-----------------------------------------------------------------------
C purpose: This Subroutine determines if a point (xval,yval) is on a 
C            local 2-D triangular linear element
C    
C restrictions: Applicable only for  2-D triangular linear element
C
C inputs:  xlocal(3) and y(3) are 1-D arrays containing the x and y 
C            coordinates of each local node on this element.
C            NOTE:  The node locations must be entered in a
C                     counterclockwise order (LABLIB standard)
C          xval,yval is the location of the point 
C
C outputs: nflag is set equal to 1 if the point (xval,yval) is on the 
C            local element.
C          nflag is set equal to 0 if the point (xval,yval) is not on 
C            the local element.
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension xlocal(3),ylocal(3)
C
C LOCAL VARIABLES
      dimension vx(3),vy(3)
C
C  Calculate the x and y components of vectors pointing from (xval,yval)
C    to each node on the element
      do 10 i=1,3
c         delx=xlocal(i)-xval
c         dely=ylocal(i)-yval
c         d=SQRT(delx**2.0+dely**2.0)
c         theta=ATAN2(dely,delx)
c         vx(i)=d*COS(theta)
c         vy(i)=d*SIN(theta)
         vx(i)=xlocal(i)-xval
         vy(i)=ylocal(i)-yval
 10   continue
C
C  Determine if the point is on the element by calculating the
C    crossproducts of neighboring vectors in a direction
C    which yeild all positive numbers if x,y is on the element.
C    (ie: v1Xv2>0, v2Xv3>0, and v3Xv1>0 => x,y is on element)
      do 20 i=1,3
         j=i+1
         if(i.EQ.3)j=1
         crossprod=vx(i)*vy(j)-vy(i)*vx(j)
         if(crossprod.LT.0.0)then
            nflag=0
            go to 21
         endif
 20   continue
      nflag=1
 21   continue
C
C END OF ROUTINE
      return
      end
C****67***************************************************************72
C***********************************************************************
C****67***************************************************************72
C***********************************************************************
      SUBROUTINE PHI2D1(xl,yl,xvalu,yvalu,phi)
C-----------------------------------------------------------------------
C purpose: This subroutine evaluates the value of the linear basis
C            on a 2-D triangular linear element at a 
C            point (xvalu,yvalu) on a local element.
C
C restrictions: Applicable only for  2-D triangular linear element
C               Assumes that the fact that (xval,yval) lies on this 
C                 element has been determined a priori.
C               Attempts to eliminate round-off problems.
C               Minimum nodal area value is 0.0 => this subroutine
C                 can not be used to determine if the point in question
C                 is on an element
C
C inputs:  xl(3) and yl(3) are 1-D arrays containing the x and y
C            coordinates of each local node on this element.
C            NOTE:  The node locations must be entered in a
C                     counterclockwise order (LABLIB standard)
C          xvalu,yvalu is the location of the point 
C
C outputs: phi is a 1-D array containing the value of the 3 basis
C            functions at the point xval,yval.
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           1992 APRIL 22 - Original Version
C           1997 MARCH 01 - PHI2D1 - Revised version of PHI2D.
C                           Revised to reduce potential round-off errors
C                           which resulted in basis functions which
C                           were less than 0.0 and greater than 1.0 on
C                           some machines (this potential problem never
C                           revealed itself until early 1997).
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension xl(3),yl(3),phi(3)
C
C START OF EXECUTABLE CODE
C
C Calculate twice the nodal areas
      A12=xvalu*(yl(2)-yl(3))+xl(2)*(yl(3)-yvalu)+xl(3)*(yvalu-yl(2))
      A22=xl(1)*(yvalu-yl(3))+xvalu*(yl(3)-yl(1))+xl(3)*(yl(1)-yvalu)
      A32=xl(1)*(yl(2)-yvalu)+xl(2)*(yvalu-yl(1))+xvalu*(yl(1)-yl(2))
C
C Do not allow for negative areas (which only occur due to round-off
C errors)
      A12=AMAX1(A12,0.0)
      A22=AMAX1(A22,0.0)
      A32=AMAX1(A32,0.0)
C
C Sum twice the nodal areas to get twice the total elemental area and 
C then evaluate basis functions
      AELE2=A12+A22+A32
      phi(1)=A12/AELE2      
      phi(2)=A22/AELE2     
      phi(3)=A32/AELE2      
C
C END OF ROUTINE
      return
      end
C***********************************************************************
C***********************************************************************
C***********************************************************************
C***********************************************************************
      SUBROUTINE BASIS1D(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=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-----------------------------------------------------------------------
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=n1=0
            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
      b1=1.0
      n1=nn
      b2=0.0
      n2=0.0
 20   continue
C
C END OF ROUTINE
      return
      end
C***********************************************************************
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 zdes
C          n2 is the reference node immediately above zdes
C          b2 is the value of the n2's basis function at zdes
C
C notes:   zdes<zref(1)  => n1=nn,b1=0.0,n2=1,b2=1.0 
C          zdes>zref(nn) => n1=nn,b1=1.0,n2=1,b2=0.0 
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           26 AUGUST 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension zref(nn)
C
C START OF EXECUTABLE CODE
      do 10 k=1,nn 
         if(zdes.LT.zref(k))then
C
C if zdes<zref(1), zdes is below the reference data 
C                    => n1=nn,b1=0.0,n2=1,b2=1.0 
            n2=k
            n1=k-1
            if(n2.EQ.1)then
                n1=nn
                b1=0.0
                b2=1.0
                go to 20
             endif
             zdiff=zref(n2)-zref(n1)
             b1=(zref(n2)-zdes)/zdiff
             b2=(zdes-zref(n1))/zdiff  
             go to 20
         endif
 10   continue
C
C if zdes>zref(nn), zdes is above the reference data
C                     => n1=nn,b1=1.0,n2=1,b2=0.0
      b1=1.0
      n1=nn
      b2=0.0
      n2=1
 20   continue
C
C END OF ROUTINE
      return
      end
C***********************************************************************
C***********************************************************************

C***********************************************************************
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)
      CALL CSTRING(header,lhs,lhe)
C
C Write file
      open (unit=7,file=filename(lfs:lfe),status='new')
      write(7,101)mesh(lms:lme)
      write(7,101)header(lhs:lhe)
      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***********************************************************************
C***********************************************************************
      SUBROUTINE WRITES3R(nndim,nnvdim,filename,mesh,header,z,s,nn,nnv)
C-----------------------------------------------------------------------
C purpose: This routine writes a .s3r file
C
C index range: i=1,nn, j=1,nnv
C
C inputs:  nndim - parameter used to dimension nodal arrays in the
C            calling routine
C          nnvdim - parameter used to dimension the number of 
C            vertical nodes in the calling routine
C          filename - name of file to be read
C          mesh - name of the mesh
C          header - character string describing file contents
C          z(i,j) - z coordinate at horizontal node i
C            and vertical node j
C          s(i,j) - real scalar value at horizontal node i
C            and vertical node j
C          nn - number of horizontal nodes for which data exists
C          nnv - number of vertical nodes 
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 z(nndim,nnvdim),s(nndim,nnvdim)
C
C Check for oversize of arrays
      nflag=0
      if(nn.GT.nndim)then
	 write(*,*)'Dimensioning problem in SUBROUTINE WRITES3R:'
         write(*,*)'  nn exceeds nndim => rectify and recompile'
         nflag=1
      endif
         if(nnv.GT.nnvdim)then
	 write(*,*)'Dimensioning problem in SUBROUTINE WRITES3R:'
         write(*,*)'  nnv exceeds nnvdim => 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))
      write(7,101)mesh(lms:lme)
      write(7,101)header
      write(7,*)nnv
      write(7,103)((i,z(i,j),s(i,j),j=1,nnv),i=1,nn)
      close(7)
C
C Format Statements
 101  format(a)
 103  format(i5,2e12.4)
      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      

      subroutine compute_dmm(nobs,tobs,xobs,yobs,zobs,sobs,
     &inst,r,s,z,x,y,in,nn,ne,nnv)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c
c compute the residual of observations with respect to prior estimate-s
c
c Inputs:
c	nobs- number of observations
c	tobs(k)- time of observation k
c       xobs(k)- x location of observation k
c       yobs(k)- y location of observation k
c       zobs(k)- z location of observation k
c       sobs(k)- observation k
c	inst(k)- instrument number,r,s,z,x,y,in,nn,ne,nnv
c	 s(n,j)- background field value at node n, sigma level j
c	 z(n,j)- depth of sigma level j at node n
c          x(n)- x-cordinate of mesh node n
c          y(n)- y-cordinate of mesh node n
c	     in- element incedence list
c            nn- number of horizontal nodes
c           nnv- number of sigma levels
c Outputs: 
c       r(k)- sobs(k)-s for observation k
c 
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      implicit none
      include 'oaci.dim'
      
      real tobs(nobsdim),xobs(nobsdim),yobs(nobsdim),zobs(nobsdim),
     &     sobs(nobsdim),r(nobsdim)
      
      integer nobs,inst(nobsdim)
      real s(nndim,nnvdim),z(nndim,nnvdim),x(nndim),y(nndim),s0,
     &     hdown(nndim)
      integer in(3,nedim),ne,nn,inmesh,k,nnv
      integer id

      id=0
      do k=1,nobs
       if (mod(k,100).eq.0) write(*,*)
     &       'interpolating background field to observation #',k

       call interp3d(x,y,z,s,in,nn,nnv,ne,xobs(k),yobs(k),zobs(k),
     &               s0,inmesh)
       if (inmesh.eq.1)then
        id=id+1
	tobs(id)=tobs(k)
	xobs(id)=xobs(k)
	yobs(id)=yobs(k)
	zobs(id)=zobs(k)
        r(id)=sobs(k)-s0
       endif
      enddo
      nobs=id

      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      subroutine interphdown2obs(nobs,xobs,yobs,hdobs,x,y,hdown,in,
     &                           nn,ne)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c Interpolate bathymetry to xobs,yobs
c
c Inputs:
c	nobs: number of points to interpolate to
c	xobs(k): x-coordinate of point to interpolate to
c	yobs(k): y-coordinate of point to interpolate to
c	x(n): x-coordinate of mesh node n
c	y(n): y-coordinate of mesh node n
c	hdown(n): bathymetric depth at mesh node n
c	in: element incedence list 
c	nn: number of horizontal mesh nodes
c	ne: number of elements
c
c Outputs:
c	hdobs(k) bathymetric depth at point xobs(k),yobs(k)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c*******************************************************************
c kws 6/28/2003
c*******************************************************************
     
     

      implicit none
      include 'oaci.dim'

      real x(nndim),y(nndim),hdown(nndim)
      integer in(3,nedim),nn,ne,nobs
      real xobs(nobsdim),yobs(nobsdim),hdobs(nobsdim)
      integer nelem,nclosest,k,n
      real basis(3)


      do k=1,nobs
        call BASIS2D(nndim,nedim,in,x,y,nn,ne,xobs(k),yobs(k),nclosest,
     &               nelem,basis)
        if (nelem.ne.0) then
	  hdobs(k)=basis(1)*hdown(in(nelem,1))+
     &             basis(2)*hdown(in(nelem,2))+
     &             basis(3)*hdown(in(nelem,3))
        else
	  hdobs(k)=9999999.
	  write(*,*)'Error interpolating bathymetry for obs#',k
	endif
      enddo
      
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      

      subroutine interp3d(x,y,z,s,in,nn,nnv,ne,x0,y0,z0,s0,inmesh)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c
c interpolate 3 dimensional field on fe mesh s to points x0, y0, z0
c
c Inputs: 
c	x(n): x-coordinate of fe mesh node n
c	y(n): y-coordinate of fe mesh node n
c	z(n,j): z-coordinate sigma level j at node n
c       in: element incedence list
c       nn: number of horizontal nodes
c	nnv: number of sigma levels
c	ne: number of elements
c	x0(k): x- coordinate to interpolate to
c	y0(k): y- coordinate to interpolate to
c	z0(k): z- coordinate to interpolate to
c
c Outputs:
c	s0(k): linear interpolated value at x0(k),y0(k),z0(k)
c       inmesh(k): 1 if  x0(k),y0(k) is in the mesh 0 otherwise
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c*******************************************************************
c kws 6/28/2003
c*******************************************************************

      implicit none
      include 'oaci.dim'

      real x(nndim),y(nndim),z(nndim,nnvdim),s(nndim,nnvdim)
      integer in(3,nedim),nn,ne,nnv
      real x0,y0,z0,s0
      integer inmesh,nelem,nclosest,j
      real zt(nnvdim),st(nnvdim)
      real basis(3),d
      call BASIS2D(nndim,nedim,in,x,y,nn,ne,x0,y0,nclosest,nelem,
     &             basis)
      inmesh=0
      if (nelem.gt.0)then
       inmesh=1 
       do j=1,nnv
        zt(j)=basis(1)*z(in(nelem,1),j)+basis(2)*z(in(nelem,2),j)+
     &       basis(3)*z(in(nelem,3),j)
	st(j)=basis(1)*s(in(nelem,1),j)+basis(2)*s(in(nelem,2),j)+
     &       basis(3)*s(in(nelem,3),j)
       enddo
       if (z0.lt.zt(1)) s0=st(1)
       if (z0.gt.zt(nnv)) s0=st(nnv)
       if ((z0.gt.zt(1)).and.(z0.lt.zt(nnv))) then
        do 11 j=2,nnv
         if (zt(j).ge.z0) goto 12
   11   enddo
   12   continue
        d=zt(j)-zt(j-1)
        s0=st(j)*(z0-zt(j-1))/d + st(j-1)*(zt(j)-z0)/d 
       endif
      endif

      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c oa2mesh
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      subroutine oa2mesh(nn,nnv,x,y,z,hdown,t,nobs,tobs,xobs,yobs,zobs,
     &hdobs,stdobs,inst,r,p,v,nr)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c
c Primary driver routine for objective analysis to 3d mesh 
c
c inputs:
c 	nn number of horizontal nodes in mesh
c 	nnv number of vertical sigma layers in mesh
c 	x(n) x- coridinate for node n of the fe mesh
c       y(n) y- coridinate for node n of the fe mesh
c       z(n,j) z- coridinate for node n and sigma level j
c       hdown(n) bathymetric depth at x(n),y(n)
c 	t time of estimation
c
c 	nobs number of observations
c 	tobs,xobs,yobs,zobs,hdobs time,x,y,z, and bathymetric depth for
c						observation # k
c 	stdobs(k) standart error for obeservation #k 
c 	inst(k) insturment number for observation #k
c 	r(k) residual (observation-background estimate) for observation #k
c
c output: 
c	 p(n,j) estimate at x(n),y(n),z(n,j),t
c	 v(n,j) posterior variance of estimate at x(n),y(n),z(n,j),t
c	 nr(n,j) number of data points used in estimation at x(n),y(n),z(n,j),t
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c*******************************************************************
c kws 6/28/2003
c*******************************************************************
      implicit none
      include 'oaci.dim'

      integer nn,nnv,nobs
      real x(nndim),y(nndim),z(nndim,nnvdim),t,p(nndim,nnvdim),
     &v(nndim,nnvdim),hdown(nndim)
      real tobs(nobsdim),xobs(nobsdim),yobs(nobsdim),zobs(nobsdim),
     &r(nobsdim),stdobs(nobsdim),hdobs(nobsdim)
      
      integer inst(nobsdim)
      real nr(nndim,nnvdim)
      
      real c0
      
      real tl(nlobsdim),xl(nlobsdim),yl(nlobsdim),zl(nlobsdim),
     &rl(nlobsdim),sl(nlobsdim),hl(nlobsdim)
      integer il(nlobsdim)
      integer nlobs,hindex(nobsdim),vindex(nobsdim)
      integer n,k,l,j,kk,i
      integer ipivot(nlobsdim),info
      
      real Cdd(nlobsdim,nlobsdim),W(nlobsdim,nlobsdim),Cmd(nlobsdim),
     & A(nlobsdim,nlobsdim),sdot,B(nlobsdim,nlobsdim),
     & lambda(nlobsdim)
c      real tmp,pd(nlobsdim),pdist
      real tmp,pd(nobsdim),pdist
      integer hitmp

      do n=1,nn
       if (mod(n,100).eq.0) write(*,*)
     &       'interpolating to node #',n
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c select data points for interpolation to all sigma levels at node n
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c
c compute psudo-distance between first nlobs observations and node n 
c
      nlobs=min(nlobsdim,nobs)
      do j=1,nlobs
       pd(j)=pdist(x(n)-xobs(j),y(n)-yobs(j),hdown(n)-hdobs(j),
     &                t-tobs(j))
       hindex(j)=j
      enddo
c
c sort pd
c
      do i=1,nlobs
       do j=1,nlobs-i
c       do j=1,nlobs-i+1  ! RHE 04/12/2004
        if(pd(j+1).lt.pd(j))then
	 tmp=pd(j)
	 hitmp=hindex(j)
	 pd(j)=pd(j+1)
	 hindex(j)=hindex(j+1)
	 pd(j+1)=tmp
	 hindex(j+1)=hitmp
	endif
       enddo !i=1,lobs
      enddo !j=1,nlobs-i+1

c add in other observations to sorted list of nlobs closest observations
c
      do j=nlobs+1,nobs
        tmp=pdist(x(n)-xobs(j),y(n)-yobs(j),hdown(n)-hdobs(j),t-tobs(j))
	if (tmp.lt.pd(nlobs))then
	 do k=1,nlobs-1
	  i=nlobs-k
	  if (tmp.ge.pd(i)) then
	   pd(i+1)=tmp
	   hindex(i+1)=j
	   goto 102
	  else
	   pd(i+1)=pd(i)
	   hindex(i+1)=hindex(i)
	  endif	   
	 enddo !k=1,nlobs-1
	 if(i.eq.1)then
          pd(1)=tmp
	  hindex(1)=j
	 endif
	endif
 102   continue
      enddo !j=nlobs+1,nobs
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c end select data points for interpolation to all sigma levels at node n
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

       do k=1,nlobs
         tl(k)=tobs(hindex(k))
	 xl(k)=xobs(hindex(k))
	 yl(k)=yobs(hindex(k))
	 zl(k)=zobs(hindex(k))
	 hl(k)=hdobs(hindex(k))
	 sl(k)=stdobs(hindex(k))
	 il(k)=inst(hindex(k))
	 rl(k)=r(hindex(k))
       enddo
c
c compute  inv(H Cuu H' + W )(d- H u)  
c to use at all vertical sigma layers at node n
c
      do j=1,nlobs
        do k=1,nlobs
        call modelcov(Cdd(j,k),xl(j),yl(j),zl(j),hl(j),tl(j),
     &                         xl(k),yl(k),zl(k),hl(k),tl(k))   
        call obsercov(W(j,k),xl(j),yl(j),zl(j),hl(j),tl(j),sl(j),il(j),
     &                       xl(k),yl(k),zl(k),hl(k),tl(k),sl(k),il(k))
        enddo   
      enddo
      
      do j=1,nlobs
        do k=1,nlobs
	   A(j,k)=Cdd(j,k)+W(j,k)
        enddo
      enddo

      do j=1,nlobs
        do k=1,nlobs
	   B(j,k)=0.0
        enddo
	B(j,j)=1.
      enddo
      
      call sgesv(nlobs,nlobs,A,nlobsdim,ipivot,B,nlobsdim,info)
       
c interpolate to each vertical level
       do j =1,nnv
         nr(n,j)=real(nlobs)
	  
c         write(*,*)n,j,nlobs
	 if (nlobs .eq. 0) then
	   p(n,j)=0.
	   v(n,j)=9999999.
	 else
         do k=1, nlobs
           call modelcov(Cmd(k),xl(k),yl(k),zl(k),hl(k),tl(k),
     &                          x(n), y(n),z(n,j),hdown(n),t)
         enddo 

         do k=1,nlobs
	   lambda(k)=0.
	   do kk=1,nlobs
	    lambda(k)=lambda(k)+Cmd(kk)*B(k,kk)
	   enddo
c	   lambda(k)=sdot(nlobs,Cmd,1,B(k,1),1)
         enddo

         p(n,j)=sdot(nlobs,lambda,1,rl,1)
         call modelcov(c0,x(n),y(n),z(n,j),hdown(n),t,
     &                    x(n),y(n),z(n,j),hdown(n),t)
         v(n,j)=c0-sdot(nlobs,lambda,1,cmd,1)
         endif
       enddo !j=1,nnv
      enddo! n=1,nn

      return
      end
