C
C Charles,
C
C The attached subroutines have passed all my tests.
C Please call me when you want to discuss their implementation.
C
C Cheers,
C Chris
C

      SUBROUTINE SACO_BASIS(KD0,SECKD0,ITER0,DELT,IN,X,Y,NN,NNV,NE,Z)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'CASCO.DIM'
      INCLUDE 'observes.h'
      INCLUDE 'weights.h'
C
C Argument Declaration
C
C   FEM Arrays
      INTEGER IN(3,NEDIM)
      REAL X(NNDIM),Y(NNDIM),Z(NNDIM,NNVDIM)
C
C   Horizontal, temporal, and vertical weights at observation locations
      REAL ZMODOBSII(NNVDIM)
      REAL BXY(3)
C
C Ensure routine is called with ITER=0
      IF(ITER0.NE.0)THEN
         WRITE(*,'(3A)')'SUBROUTINE SACO_BASIS MUST BE CALLED WITH ',
     &   'TIMING PARAMETERS WHICH CORRESPOND TO THE BEGINNING OF THE ',
     &   'FORWARD MODEL'
         STOP
      ENDIF
C
C Observation Loop
      DO II=1,NOBS
C
C   Determine horizontal weights
         CALL BASIS2D(NNDIM,NEDIM,IN,X,Y,NN,NE,XOBS(II),YOBS(II),
     &   NCLOSEST,NELEM,BXY)
         IF(NELEM.EQ.0)THEN
            WRITE(*,*)'OBSERVATIONS LIE OUTSIDE HORIZONTAL DOMAIN'
            STOP
         ENDIF
         BASISOBS(1,II)=FLOAT(NELEM)
         BASISOBS(2,II)=BXY(1)
         BASISOBS(3,II)=BXY(2)
         BASISOBS(4,II)=BXY(3)
C
C   Determine vertical integration weights
         DO J=1,NNV
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,Z,
     &      ZMODOBSII(J))
         ENDDO
         CALL VERTWEIGHTS(II,NNV,ZMODOBSII)
      ENDDO
C
C Initialize for temporal loop
      KD=KD0
      SECKD=SECKD0
      ITER=ITER0
      IOBS=1
C
C Temporal loop to determine temporal weights
      DO WHILE(IOBS.LE.NOBS) 
         DTSEC=(KDOBS(IOBS)-KD)*24.0*3600.0+SECOBS(IOBS)-SECKD
         IF(DTSEC.LT.0.0)THEN
         ELSEIF(DTSEC.GE.DELT)THEN
            ITER=ITER+1
            SECKD=SECKD+DELT
            IF(SECKD.GE.86400.)CALL UP_DATE(KD,SECKD)
         ELSE
            BASISOBS(5,IOBS)=ITER+1
            BASISOBS(6,IOBS)=DTSEC/DELT
            IOBS=IOBS+1
         ENDIF
      ENDDO
c      open(59,FILE='sacobasis.vel.dat')
c      do i=1,nobs
c       write(59,'(i6,f10.2,f10.2,f10.2,i6,f10.2)')int(basisobs(1,i)),
c     &basisobs(2,i),basisobs(3,i),
c     &basisobs(4,i),int(basisobs(5,i)),basisobs(6,i)
c      enddo
c      close(59)
           
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTWEIGHTS(II,NNV,ZMODOBSII)
C-----------------------------------------------------------------------
      INCLUDE 'CASCO.DIM'
      INCLUDE 'observes.h'
      INCLUDE 'weights.h'
C
C Dependent variables for observation locations
      REAL ZMODOBSII(NNVDIM)
C
C   Initialize weights
         DO J=1,NNV
            VWEIGHT(II,J)=0.0
         ENDDO
C
C   Determine observational range over which averaging should occur.
C      ZMAXOBS, ZMINOBS are the bottom and top of the range in the
C            observed vertical coordinate system (positive downwards).
C      ZBOTTOMOBS(II) is the observed water depth. 
C
C      ZZMIN, ZZMAX are the bottom and top of the range in the 
C            model vertical coordinate system (postive upwards).
C      ZMODOBSII(I) are the models vertical nodes interpolated to 
C            the observation location. 
C
C The model vertical coordinate system is stretched/contracted so that
C the same fraction of the water column is sampled in the model and the 
C observations. 
C      ZZMIN <---> ZMAXOBS
C      ZZMAX <---> ZMINOBS 
C   
C   (Note: ZZMIN and ZZMAX are the bottom and top of the range,
C   respectively)
         ZZMIN=ZMODOBSII(1)+(ZBOTTOMOBS(II)-ZMAXOBS(II))*
     &   (ZMODOBSII(NNV)-ZMODOBSII(1))/ZBOTTOMOBS(II)
         ZZMAX=ZMODOBSII(1)+(ZBOTTOMOBS(II)-ZMINOBS(II))*
     &   (ZMODOBSII(NNV)-ZMODOBSII(1))/ZBOTTOMOBS(II)
         IF(ZZMIN.GT.ZZMAX)THEN
            WRITE(*,*)'VERTICAL INTEGRATION INTERVAL IS NEGATIVE FOR',
     &      ' OBSERVATION NUMBER ',II
	    write(*,*)ZZMIN, ZZMAX
	    write(*,*)ZMODOBSII(1),ZMODOBSII(NNV), 
     &                ZMAXOBS(II),ZMINOBS(II), ZBOTTOMOBS(II)
            STOP
         ENDIF
C
C   Evaluate vertical weights extremes of averaging interval
         CALL BASIS1D1(NNV,ZMODOBSII,ZZMIN,B2MIN,N2MIN,B1MIN,N1MIN)
         CALL BASIS1D1(NNV,ZMODOBSII,ZZMAX,B2MAX,N2MAX,B1MAX,N1MAX)
C
C   Compute vertical weights for vertical integration
         IF(N1MIN.GT.N1MAX)THEN
            WRITE(*,*)'VERTICAL INTEGRATION INTERVAL IS NEGATIVE FOR',
     &      ' OBSERVATION NUMBER ',II
            STOP
         ELSEIF(N1MIN.EQ.N1MAX)THEN
            DZ=ZZMAX-ZZMIN
            VWEIGHT(II,N1MIN)=0.5*DZ*(B1MIN+B1MAX)
            VWEIGHT(II,N2MIN)=0.5*DZ*(B2MIN+B2MAX)
         ELSE
            DZ=ZMODOBSII(N2MIN)-ZZMIN
            VWEIGHT(II,N1MIN)=0.5*DZ*B1MIN
            VWEIGHT(II,N2MIN)=0.5*DZ*(B2MIN+1.0)
            DO J=N2MIN,N1MAX-1
               DZ=ZMODOBSII(J+1)-ZMODOBSII(J)
               VWEIGHT(II,J)=VWEIGHT(II,J)+0.5*DZ
               VWEIGHT(II,J+1)=VWEIGHT(II,J+1)+0.5*DZ
            ENDDO
            DZ=ZZMAX-ZMODOBSII(N1MAX)
            VWEIGHT(II,N1MAX)=VWEIGHT(II,N1MAX)+0.5*DZ*(B1MAX+1.0)
            VWEIGHT(II,N2MAX)=VWEIGHT(II,N2MAX)+0.5*DZ*B2MAX
         ENDIF
C
C   Normalize for vertical averaging
         DZTOTAL=ZZMAX-ZZMIN
         DO J=1,NNV
            VWEIGHT(II,J)=VWEIGHT(II,J)/DZTOTAL
         ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SACO_UVERROR(ITER,IN,NNV,U,V)
C-----------------------------------------------------------------------
C Reports the unexplained observations. 
C    Observations - Model      (I think CGH Nov 98) 
C-----------------------------------------------------------------------
      INCLUDE 'CASCO.DIM'
      INCLUDE 'observes.h'
      INCLUDE 'errors.h'
      INCLUDE 'weights.h'
C
C Argument Declaration
C
C   FEM Arrays
      INTEGER IN(3,NEDIM)
      REAL U(NNDIM,NNVDIM),V(NNDIM,NNVDIM)
C
C   Horizontal basis functions
      REAL BXY(3)
C
C Save necessary arguments
      SAVE IIOBS,IIOBSDONE
C
C Initialize IIOBS during first iteration
      IF(ITER.EQ.0)THEN
         IIOBS=1
         IIOBSDONE=0
      ENDIF
C
C Return if iteration number precedes that required for the first
C observational location or is after that required for the last
C observational location 
      IF(ITER.LT.INT(BASISOBS(5,1))-1)RETURN
      IF(IIOBSDONE.EQ.NOBS)RETURN
C
C Perform calculations
C
C  SACO time level just after observational location
      DO II=IIOBSDONE+1,IIOBS-1
         USUM=0.0
         VSUM=0.0
         NELEM=INT(BASISOBS(1,II))
         BXY(1)=BASISOBS(2,II)
         BXY(2)=BASISOBS(3,II)
         BXY(3)=BASISOBS(4,II)
         TWEIGHT=BASISOBS(6,II)
         DO J=1,NNV
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,U,UIIJK)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,V,VIIJK)
            USUM=USUM+TWEIGHT*VWEIGHT(II,J)*UIIJK
            VSUM=VSUM+TWEIGHT*VWEIGHT(II,J)*VIIJK
         ENDDO
         UERROR(II)=UOBS(II)-(UERROR(II)+USUM)
         VERROR(II)=VOBS(II)-(VERROR(II)+VSUM)
      ENDDO
      IIOBSDONE=IIOBS-1
C
C  SACO time level just prior to observational location
 2    IF(ITER.EQ.INT((BASISOBS(5,IIOBS)-1)))THEN
C
C   Compute contribution to integral
         USUM=0.0
         VSUM=0.0
         NELEM=INT(BASISOBS(1,IIOBS))
         BXY(1)=BASISOBS(2,IIOBS)
         BXY(2)=BASISOBS(3,IIOBS)
         BXY(3)=BASISOBS(4,IIOBS)
         TWEIGHT=1.0-BASISOBS(6,IIOBS)
         DO J=1,NNV
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,U,UIIJK)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,V,VIIJK)
            USUM=USUM+TWEIGHT*VWEIGHT(IIOBS,J)*UIIJK
            VSUM=VSUM+TWEIGHT*VWEIGHT(IIOBS,J)*VIIJK
         ENDDO
         UERROR(IIOBS)=USUM
         VERROR(IIOBS)=VSUM
C
C   Procede to next observation
         IIOBS=IIOBS+1
         GO TO 2
      ENDIF
      RETURN
      END
C***********************************************************************
C***********************************************************************

