c***********************************************************************
c***********************************************************************
      subroutine verticalm4(nn,nnv,ne,in,delt,sv,
     & cxmat,cymat,dxmat,dymat,iq,jq,
     & d1,d2,d3,dz,zetamid,uzmid,vzmid,uznew,vznew,
     & basisobs,vweight,uerror,verror,w_delta,ktime,obs_ptr)
c-----------------------------------------------------------------------
c drl 27 oct. 98
c cgh 19 nov. 98 start to add error forcing
c adjoint momentum equation; from verticals4.f
c
c Objectives: Determine the solution of the vertical diffusion equations
c             for the 3-d adjoimt velocities (centered at time k+1/2).
c             The tridiagonal system is solved for 
C                (unew+umid) - i (vnew+vmid). 
c             Then unew, vnew are computed at time k. 
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               nftr - sparse factor
c
c Inputs: iq,jq  - the standard SPRSPAK pointers 
C         delt   - the time step 
c         sv(i)  - horizontal mass matrix.
c         zetamid(i),umid(i),vmid(i),uzmid(i,j),vzmid(i,j) are
c           the adjoint variables at time level k+1
C         basisobs, vweight - weight functions for observations and
c                             errors.  See weights.m
C         uerror, verror - velocity errors
c         ktime - (or k) the time node corresponding to the *new level 
C                 i.e. the end of current time step.  
C         obs_ptr(2) Pointers into the error and weight functions.
C            obs_ptr(1) keeps track of stuff in time element k
C            obs_ptr(2) keeps track of stuff in time element k+1 
C         When solving for velocities at time node k we require error
C              contributions from time elements k and k+1. 
C The arrays for the weight functions and the errors contain all the 
C weights and all the errors for all time.  The array contains the 
C time element number basisobs(5,i). 
C                    
c         
c Output:  uznew(i,j),vznew(i,j) -- adjoint velocities at time level k
C
C  Adjoint steps backwards in time:  k   --> new
C                                    k+1 --> mid
C                              k+2 --> old (no old in velocity adjoint) 
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,ne
      integer in(3,nedim)
      real    delt
      real    sv(nndim)
      real    cxmat(nftrdim*nndim), cymat(nftrdim*nndim)
      real    dxmat(nftrdim*nndim), dymat(nftrdim*nndim)
      integer iq(nndim),jq(nftrdim*nndim)
      real    zetamid(nndim)
      real    uzmid(nndim,nnvdim),vzmid(nndim,nnvdim)
      real    uznew(nndim,nnvdim),vznew(nndim,nnvdim)
      complex d1(nnvdim*nndim),d2(nnvdim*nndim), d3(nnvdim*nndim)
      real    dz(nnvdim*nndim)

      real    basisobs(6,nobsdim), vweight(nobsdim,nnvdim)
      real    uerror(nobsdim), verror(nobsdim) 
      real    W_delta(nobsdim)
      integer ktime
      integer obs_ptr(2)
c
c local dimensioning and assignments
c
      real surf(nndim)
      real qxmid(nndim), qymid(nndim)
      real ubotmid(nndim), vbotmid(nndim)
c
      complex rz(nnvdim*nndim)
c
      integer i,j,k 
      integer kbot 
      integer ig,il,ie,ell,jbase
      integer itime_e
      real    xyw,tw
      logical done

      complex eye
c      parameter (eye=cmplx(0.0,1.0))
c
c     write(*,*)'verticalm4: k,obs_ptr ',ktime,obs_ptr(1),obs_ptr(2)
c     write(*,*)nn,nnv,ne,delt
c     write(*,*)uerror(1),uerror(2),uerror(3),uerror(4)
c     write(*,*)verror(1),verror(2),verror(3),verror(4)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c part 2: make rhs and solve
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      eye=cmplx(0.0,1.0)
c 
c make right-hand side parts.  
c
      do i=1,nn
         surf(i) = -zetamid(i)
      enddo
c

      call sprspremlt(cxmat, iq, jq, surf, qxmid, nn)
      call sprspremlt(cymat, iq, jq, surf, qymid, nn)
      call sprspremlt(dxmat, iq, jq, surf, ubotmid, nn)
      call sprspremlt(dymat, iq, jq, surf, vbotmid, nn)
c
c assemble rhs
c remember we have (u-iv) in the adj. system
c
      k=0
      do i=1,nn
         do j=1,nnv
          k=k+1
          rz(k)=(2./delt*(uzmid(i,j)-eye*vzmid(i,j)))*sv(i)*dz(k)
     &         -(qxmid(i)-eye*qymid(i))*dz(k)         
         enddo
         kbot=(i-1)*nnv +1
         rz(kbot)=rz(kbot)-(ubotmid(i)-eye*vbotmid(i))
      enddo
c
c add error forcing
C
c      write(2,*)'***TESTING**** setting RZ=0 before adding errors'
c      write(*,*)'***TESTING**** setting RZ=0 before adding errors'
c      k = 0 
c      do i = 1,nn
c         do j = 1,nnv
c            k = k+1
c            rz(k) = 0
c         enddo
c      enddo
C  deal with stuff in time element ktime 
      itime_e = ktime
      done = .false.
      ie = obs_ptr(1)
c     if(ie.gt.1) then 
c     write(*,*)'ie itime_e basisobs(5,ie-1) nint(basisobs(5,ie-1))'
c     write(*,*)ie, itime_e, basisobs(5,ie-1), nint(basisobs(5,ie-1))
c     else 
c        write(*,*)'ie ',ie
c     endif
100   continue
        if(ie.le.1) then
           done = .true.
        elseif(nint(basisobs(5,ie-1)).lt.itime_e) then
           done = .true.
        elseif(nint(basisobs(5,ie-1)).eq.itime_e) then
           ie = ie - 1
           ell = nint(basisobs(1,ie))
           tw = basisobs(6,ie) 
c	   write(*,*)'Do it for k ',ie,ell,tw
           do il = 1,3
              ig = in(il,ell)
              xyw = basisobs(il+1,ie)
              jbase = (ig-1)*nnv
	      do j = 1, nnv
		 rz(j+jbase) = rz(j+jbase) + 
     &            (uerror(ie) - eye*verror(ie))*xyw*vweight(ie,j)*tw
     &            *W_delta(ie) 
	      enddo
           enddo
        else
           write(*,*)'Trouble with basis_obs and pointer'
	   write(*,*)'ie,basisobs(5,ie-1)', ie,basisobs(5,ie-1)
	   stop
        endif
        if(.not.done) goto 100
101   continue
      obs_ptr(1) = ie 

C  deal with stuff in time element k+1 
      itime_e = ktime+1
      done = .false.
      ie = obs_ptr(2)
c     write(*,*)'ie 2',ie
200   continue
        if(ie.le.1) then
           done = .true.
        elseif(nint(basisobs(5,ie-1)).lt.itime_e) then
           done = .true.
        elseif(nint(basisobs(5,ie-1)).eq.itime_e) then
           ie = ie - 1
           ell = nint(basisobs(1,ie))
           tw = 1.0 - basisobs(6,ie) 
c	   write(*,*)'Do it for k+1 ',ie,ell,tw
           do il = 1,3
              ig = in(il,ell)
              xyw = basisobs(il+1,ie)
              jbase = (ig-1)*nnv
	      do j = 1, nnv
		 rz(j+jbase) = rz(j+jbase) + 
     &             (uerror(ie) - eye*verror(ie))*xyw*vweight(ie,j)*tw
     &            *W_delta(ie) 
	      enddo
           enddo
        else
           write(*,*)'Trouble with basis_obs and pointer. part 2'
	   write(*,*)'ie,basisobs(5,ie-1)', ie,basisobs(5,ie-1)
        endif
        if(.not.done) goto 200
201   continue
      obs_ptr(2) = ie 

c     write(21,*)'ktime = ', ktime
c     k = 0
c     eps = 1.0e-10
c     do i = 1, nn
c         do j = 1,nnv 
c            k = k+1   
c            if(cabs(rz(k)).gt.eps) then 
c               write(21,*)i,j,k,rz(k)
c            endif 
c         enddo
c     enddo
c
c
c   solve for 2*[average of (u-iv) over the time step]
c
         call cthomas(2,d1,d2,d3,rz,nnv*nn)
c
c   convert to values at the end of the time step (uznew and vznew)
c
      k=0
      do i=1,nn
      do j=1,nnv
         k=k+1
         uznew(i,j) =  real(rz(k)) -uzmid(i,j)
         vznew(i,j) = -aimag(rz(k))-vzmid(i,j)
      enddo
      enddo
c
c uznew,vznew are complete and at time level 'new'  
c end of routine
c
      return
      end
