C***********************************************************************
C***********************************************************************
      SUBROUTINE ELEVATIONm4(NN,NNV,in,NHBW,THETA,
     &   ZETAOLD,ZETAMID,UZMID,VZMID,UZNEW,VZNEW,
     &   SH,IQ,JQ,G0xmat,G0ymat,Amat,Bmat,DZ, 
     &   ZETANEW,ebasisobs,eerror,w_eps,ktime,obs_ptr)
C-----------------------------------------------------------------------
C Objective: Solve the wave equation (centered at time level K) for free
C            surface elevation at the end of this time step (i.e. at
C            KDNEW,SECNEW, time level K+1).
C            Note solves for ZETANEW directly, rather than ZETANEW-ZETAOLD
C
C Complete Array Documentation: See main program QUODDY4_1.1
C
C Index Range:  I=1,NN; J=1,NNV; 
C               NN - number of horizontal nodes
C               NNV - number of nodes in the vertical
C
C Inputs: NN,NNV,NHBW
C         SH(M,I) - fixed geometry matrix.
C         IQ,JQ   - standard SPRSPAK pointers
C         G0xmat,G0ymat,Amat,Bmat - stationary matrices

C         ZETAOLD(I) - adjoint elevation at time level K+2
C         ZETAMID(I),UZMID(I,J),VZMID(I,J) - adjoint elevation and 
C                                           velocity at time level K+1
C         UZNEW(I,j),VZNEW(I,j) -  adjoint velocity at time K 
C
C Output: ZETANEW(I) - elevation at time level K+1
C 
C-----------------------------------------------------------------------
C
C Declare all the integer program parameters via INCLUDE file
C
      INCLUDE 'CASCO.DIM'
C
C Global dimensioning and assignments
C
      INTEGER NN, NNV, NHBW
      integer in(3,nedim)
      REAL    THETA
      REAL ZETAOLD(NNDIM)
      REAL ZETAMID(NNDIM)
      REAL UZMID(NNDIM,NNVDIM),VZMID(NNDIM,NNVDIM)
      REAL UZNEW(NNDIM,NNVDIM),VZNEW(NNDIM,NNVDIM)
      REAL SH(NB1DIM,NNDIM)

      REAL    AMAT(NFTRDIM*NNDIM), BMAT(NFTRDIM*NNDIM)
      REAL    G0xmat(NFTRDIM*NNDIM), G0ymat(NFTRDIM*NNDIM)
      INTEGER IQ(NNDIM),JQ(NFTRDIM*NNDIM)
      REAL    DZ(NNVDIM*NNDIM)

      REAL ZETANEW(NNDIM)
C
C Elevation observations
C
      real    ebasisobs(6,nelevobsdim)
      real    eerror(nelevobsdim)
      real    w_eps(nelevobsdim)
      integer ktime
      integer obs_ptr(2)
C
C Local dimensioning and assignments
C
      REAL Zu(NNDIM), Zv(NNDIM)

C local variables 
      integer i,j,k,ie
      INTEGER NB1,ell,ig
      real Zui, Zvi
      logical done
      
       

      REAL ZERO
      DATA ZERO/0.0/

      NB1=NHBW+1
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C moody elevation solver
C unaltered from saco solver elevations3.f above this point
C START OF EXECUTABLE CODE for MOODY
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Initialize elevation array
C
      DO I=1,NN
         ZETANEW(I)=0.0
      ENDDO
C
C First, build the RHS for ZETANEW
C
      k=0
      do i=1,nn      
         Zui=0.
         Zvi=0.
         do j=1,nnv
            k=k+1
            Zui=Zui+dz(k)*(uznew(i,j)+uzmid(i,j))
            Zvi=Zvi+dz(k)*(vznew(i,j)+vzmid(i,j))
         enddo
         Zu(i)=Zui
         Zv(i)=Zvi      
      enddo   
C     
      call sprspremltadd(Amat,IQ,JQ,zetamid,zetanew,NN)
      call sprspremltadd(Bmat,IQ,JQ,zetaold,zetanew,NN)      
      call sprspremltadd(G0xmat,IQ,JQ,Zu,zetanew,NN)      
      call sprspremltadd(G0ymat,IQ,JQ,Zv,zetanew,NN) 
c      open(59,FILE='z0.dat')
c      do i=1,nn
c      write(59,*)zetanew(i)
c      enddo
c      close(59)
c-----------------------------------------------------------------------
c zetanew=zetanew+[e_(k)*][W_e][T_(k)]
c-----------------------------------------------------------------------
      itime_e = ktime
      done = .false.
      ie = obs_ptr(1)
100   continue
        if(ie.le.1) then
           done = .true.
        elseif(nint(ebasisobs(5,ie-1)).lt.itime_e) then
           done = .true.
        elseif(nint(ebasisobs(5,ie-1)).eq.itime_e) then
           ie = ie - 1
           ell = nint(ebasisobs(1,ie))
           tw = ebasisobs(6,ie) 
           do il = 1,3
              ig = in(il,ell)
              xyw = ebasisobs(il+1,ie)
c              write(*,*)'elevationm4:',
c     &                 ell,ig,tw,xyw,eerror(ie),w_eps(ie)
	      zetanew(ig)=zetanew(ig)+eerror(ie)*xyw*tw*w_eps(ie)
           enddo
        else
           write(*,*)'Trouble with basis_obs and pointer'
	   write(*,*)'ie,ebasisobs(5,ie-1)', ie,ebasisobs(5,ie-1)
	   stop
        endif
        if(.not.done) goto 100
101   continue
      obs_ptr(1) = ie 
c-----------------------------------------------------------------------
c zetanew=zeta_new+[e_(k+1)*][W_e][U_(k+1)]
c-----------------------------------------------------------------------
      itime_e = ktime+1
      done = .false.
      ie = obs_ptr(2)
200   continue
        if(ie.le.1) then
           done = .true.
        elseif(nint(ebasisobs(5,ie-1)).lt.itime_e) then
           done = .true.
        elseif(nint(ebasisobs(5,ie-1)).eq.itime_e) then
           ie = ie - 1
           ell = nint(ebasisobs(1,ie))
           tw = 1.0 - ebasisobs(6,ie) 
           do il = 1,3
              ig = in(il,ell)
              xyw = ebasisobs(il+1,ie)
	      zetanew(ig)=zetanew(ig)+eerror(ie)*xyw*tw*w_eps(ie)
           enddo
        else
           write(*,*)'Trouble with basis_obs and pointer. part 2'
	   write(*,*)'ie,ebasisobs(5,ie-1)', ie,ebasisobs(5,ie-1)
        endif
        if(.not.done) goto 200
201   continue
      obs_ptr(2) = ie 

c      open(59,FILE='z1.dat')
c      do i=1,nn
c      write(59,*)zetanew(i)
c      enddo
c      close(59)

C
C Solve for ZETANEW
C
      IF(THETA.GT.ZERO)THEN
         CALL BANSOLTR(2,SH,ZETANEW,NN,NB1,NNDIM,NB1DIM)
      ELSE
         DO I=1,NN
            ZETANEW(I)=ZETANEW(I)/SH(1,I)
         ENDDO
      ENDIF
C
C END OF ROUTINE
C
      RETURN
      END
