C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTICALS5(NN,NNV,DELT,SV,G0xmat,G0ymat,IQ,JQ,
     & D1,D2,D3,DZ, ZETAMID,UZMID,VZMID,ZETANEW,UZNEW,VZNEW)
      
C Remnants to do with velocity open boundaries 
C      G,COR,KST,NR,NPCN, NORM,LRVS,SUMX,SUMY,
C-----------------------------------------------------------------------
C All terms are centered at time K+1/2, except for the viscous 
C  terms which are weighted with EPSN. 
C
C
C Objectives: Determine the solution of the vertical diffusion
C             equations for the 3-d velocities (centered at time K+1/2).
C             The tridiagonal system solves for (Vnew+Vmid)
C             and then converts the answer to Vnew. 
C
C Note changes from Quoddy standard: 
C             Solving for (Vnew+Vmid) rather than (Vnew+Vmid)/2
C             EPSN (time weight for vertical viscosity) 
C               is hardwired to 0.5 herein 
C             G0xmat, G0ymat here replace PPX, PPY.  
C               They differ by the factor -g/2.
C
C Complete Array Documentation: See CASCO documentation
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: NN,NNV,KST,NR,NPCN,NORM(I),LRVS(I),HDOWN(I),SUMX(I),SUMY(I),
C          IQ(I),JQ(NFTR*I) contain geometric
C          information which never changes throughout the simulation
C         SV(I) - horizontal mass matrix.
C         DELT - time stepping parameters
C         ZETAMID(I),UMID(I),VMID(I),UZMID(I,J),VZMID(I,J) are
C          dependent variable values at time level K
C         ZETANEW(I) - elevation at time level K+1
C         
C Output:  UZNEW(I,J),VZNEW(I,J)
C 
C History:
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
      REAL    DELT
      REAL    SV(NNDIM)
      REAL    G0xmat(NFTRDIM*NNDIM),G0ymat(NFTRDIM*NNDIM)
      INTEGER IQ(NNDIM),JQ(NFTRDIM*NNDIM)
      REAL    ZETAMID(NNDIM),ZETANEW(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)
C
C Local dimensioning and assignments
C
      REAL gx(NNDIM),gy(NNDIM),SURF(NNDIM)
C
      COMPLEX RZ(NNVDIM*nndim)
      integer k
      real gxi, gyi
C
      INTEGER I,J

      COMPLEX EYE
      EYE=CMPLX(0.0,1.0)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C PART 2: Make RHS and solve
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C Determine the depth independent pressure gradient at the center of the
C   time step; save in gx, gy
C
C      
      DO I=1,NN
         SURF(I) = ZETANEW(I)+ZETAMID(I)
      ENDDO
      CALL SPRSMLT(G0xmat,IQ,JQ,SURF,gx,NN)
      CALL SPRSMLT(G0ymat,IQ,JQ,SURF,gy,NN)
C
c Assemble RHS
c
      k=0
      do i=1,nn
         GXi=gx(I)
         GYi=gy(I)
         do j=1,nnv
          k=k+1
          RZ(k)=(2./DELT*(UZMID(i,J)+EYE*VZMID(i,J)))*sv(i)*dz(k)
     &         +(GXi+EYE*GYi)*dz(k)
         enddo
      enddo
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 End of main momentum loop over horizontal nodes:
C UZNEW,VZNEW are complete and at time level K+1 except for BC's
C  
C END OF ROUTINE
C
      RETURN
      END
