C***********************************************************************
C***********************************************************************
      SUBROUTINE ELEVATIONS4(NN,NHBW,NNV,DELT,TAU0,
     &   KDNEW,SECNEW,ITER,DEG,NR,NVN,THETA,HDOWN,
     &   ZETAOLD,UOLD,VOLD,ZETAMID,UMID,VMID,ZMID,UZMID,VZMID,
     &   LRVS,NORM,SUMX,SUMY,DSBDRY,SH,EMAT,
     &   Amat,Bmat,Cxmat,Cymat,Dxmat,Dymat,IQ,JQ,
     &   ZETANEW)

c     SUBROUTINE ELEVATIONS2(NN,NE,NHBW,NNV,DELT,G,TAU0,COR,
c    & KDNEW,SECNEW,ITER,DEG,NR,NVN,THETA,IN,DX,DY,AR,HDOWN,
c    & ZETAOLD,UOLD,VOLD,ZETAMID,UMID,VMID,ZMID,UZMID,VZMID,
c    & AK, LRVS,NORM,SUMX,SUMY,DSBDRY,SV,SH,EMAT,ZETANEW)
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; M=1,NHBW+1; N=1,2*NHBW+1;
C               O=1,NR
C               NN - number of horizontal nodes
C               NNV - number of nodes in the vertical
C               NHBW - half bandwidth
C               NR - number of elevation boundary conditions
C
C Inputs: NN,NE,NHBW,NNV,NR,NVN,IN(3,L),DX(3,L),DY(3,L),AR(L),HDOWN(I),
C          LRVS(I),NORM(I),SUMX(I),SUMY(I),DSBDRY(I) - geometric
C          information which never changes throughout the simulation
C         SH(M,I) - fixed geometry matrix.
C         EMAT(N,O) - used for the assymetric terms resulting
C          from the boundary conditions.
C         DEG are physical parameters
C         TAU0,THETA are numerical parameters
C         KDNEW,SECNEW,ITER,DELT - time stepping parameters
C         ZETAOLD(I),UOLD(I),VOLD(I) - dependent variable values at time 
C          level K-1
C         ZETAMID(I),UMID(I),VMID(I),ZMID(I,J),UZMID(I,J),VZMID(I,J) -
C          dependent variable values at time level K-1
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, NHBW, NNV, NR, NVN 
      REAL    DELT, TAU0, THETA 
      INTEGER KDNEW        
      REAL    SECNEW
      INTEGER ITER
      REAL    DEG 
      REAL    CO, SI 
      REAL    VNT, VNORM 

      REAL HDOWN(NNDIM)
      REAL ZETAOLD(NNDIM),UOLD(NNDIM),VOLD(NNDIM)
      REAL ZETAMID(NNDIM),UMID(NNDIM),VMID(NNDIM)
      REAL ZMID(NNDIM,NNVDIM),UZMID(NNDIM,NNVDIM),VZMID(NNDIM,NNVDIM)
      INTEGER LRVS(NNDIM),NORM(NNDIM)
      REAL SUMX(NNDIM),SUMY(NNDIM),DSBDRY(NNDIM)
      REAL SH(NB1DIM,NNDIM),EMAT(NBWE,NRDIM)
      REAL    AMAT(NFTRDIM*NNDIM), BMAT(NFTRDIM*NNDIM)
      REAL    CXMAT(NFTRDIM*NNDIM), CYMAT(NFTRDIM*NNDIM)
      REAL    DXMAT(NFTRDIM*NNDIM), DYMAT(NFTRDIM*NNDIM)
      INTEGER IQ(NNDIM),JQ(NFTRDIM*NNDIM)

      REAL ZETANEW(NNDIM)
C
C Local dimensioning and assignments
C
      REAL DUMMY1(NNDIM),DUMMY2(NNDIM)

C local variables 
      INTEGER I, II, II1, I1
      REAL Qxmid(NNDIM), Qymid(NNDIM), ubotmid(NNDIM), vbotmid(NNDIM)
      INTEGER NB1
      REAL BINT




      REAL ZERO
      DATA ZERO/0.0/
      NB1=NHBW+1
C
C START OF EXECUTABLE CODE
C
C Initialize elevation array
C
      DO I=1,NN
         ZETANEW(I)=0.0
      ENDDO
C
C First, build the RHS for ZETANEW
C
      do i=1,nn
         Qxmid(i)=hdown(i)*umid(i)
         Qymid(i)=hdown(i)*vmid(i)
         ubotmid(i)=uzmid(i,1)
         vbotmid(i)=vzmid(i,1)  
      enddo   
C     
      call sprsmltadd(Amat, IQ,JQ,zetamid,zetanew,NN)
      call sprsmltadd(Bmat, IQ,JQ,zetaold,zetanew,NN)      
      call sprsmltadd(Cxmat,IQ,JQ,Qxmid,  zetanew,NN)      
      call sprsmltadd(Cymat,IQ,JQ,Qymid,  zetanew,NN)      
      call sprsmltadd(Dxmat,IQ,JQ,ubotmid,zetanew,NN)      
      call sprsmltadd(Dymat,IQ,JQ,vbotmid,zetanew,NN)
      
C
C
C Apply boundary conditions (call SUBROUTINE BC at time level
C K+1 to determine BC's for ZETANEW)
C  
C  Apply nonzero V.N BC's => boundary integral contributions
C  (NORM(I),I=1,NVN)
C  (WARNING:  Corners with V.NE.0 not implemented here)
C
      IF(NVN.NE.0) THEN 
         CALL BCS2(KDNEW,SECNEW,ITER,DEG,DUMMY1,DUMMY2,2)
         DO I1=1,NVN
            II1=NORM(I1)
            SI=SUMX(II1)
            CO=SUMY(II1)
            VNORM=DUMMY1(II1)-((UOLD(II1)*CO+VOLD(II1)*SI))*HDOWN(II1)
            VNT=(UOLD(II1)*CO+VOLD(II1)*SI)*HDOWN(II1)
            BINT=-DELT*DELT*(VNORM/2./DELT+TAU0*VNT)*DSBDRY(II1)
            BINT=BINT/2.
            ZETANEW(II1)=ZETANEW(II1)+BINT
         ENDDO
      ENDIF 
C
C  Apply set depth BC's (LRVS(I),I=1,NR) and solve for ZETANEW
CGH  The explicit part (theta=0) has been recoded but not tested.
C 
      IF(THETA.GT.ZERO) THEN
         IF(NR.NE.0)THEN 
            CALL BCS2(KDNEW,SECNEW,ITER,DEG,DUMMY1,DUMMY2,1)
            CALL EMATMULT(EMAT,DUMMY1,ZETANEW,LRVS,NN,NHBW,NR,NBWE)
         ENDIF 
         CALL BANSOLTR(2,SH,ZETANEW,NN,NB1,NNDIM,NB1DIM)
      ELSE 
         IF(NR.NE.0)THEN 
            CALL BCS2(KDNEW,SECNEW,ITER,DEG,DUMMY1,DUMMY2,1)
	    DO I = 1,NR
	       II = LRVS(I)
	       ZETANEW(II) = DUMMY1(II)
	    ENDDO
         ENDIF 
         DO I=1,NN
            ZETANEW(I)=ZETANEW(I)/SH(1,I)
         ENDDO
      ENDIF
C
C END OF ROUTINE
C
      RETURN
      END
