C***********************************************************************
C***********************************************************************
      SUBROUTINE STATIONARYC3(NN,NNV,NE,NHBW,X,Y,HDOWN,IN,NR,LRVS,G,
     &  DELT, TAU0,THETA,AR,DX,DY,NFTR,IQ,JQ,SV,SH,EMAT,AMAT,BMAT,
     &  CXMAT,CYMAT,DXMAT,DYMAT,G0xmat,G0ymat,D1,D2,D3,DZ,ZMID, 
     &  COR,AK,BCFTR,ENZM,sumx,sumy,DSBDRY,INBE,nbe,ingeo,norm,nvn,
     &  kst,npcn,lst)

C-----------------------------------------------------------------------
C Objectives: Fill stationary geometric arrays.
C             Compute stationary SPRSPAK arrays.
C             Compute stationary diagonal mass martix.
C             Compute and LU decompose transpose of stationary depth
C              matrix.
C
C Complete Array Documentation: See main program QUODDY4_1.1
C
C Index Range: I=1,NN; J=1,NNV; L=1,NE; M=1,NR: NB1=1,NHBW+1;
C              NB2=1,2*NHBW+1
C               NN - number of horizontal nodes
C               NNV - number of vertical nodes
C               NE - number of horizontal elements
C               NR - number of elevation boundary conditions
C               NHBW - half bandwidth of horizontal mesh
C
C Inputs: NN,NE,NHBW,X(I),Y(I),HDOWN(I),IN(3,L),NR,LRVS(I) - geometric
C          information which never changes throughout the simulation
C         G - gravitational acceleration
C         DELT - size of time step
C         TAU0,THETA - numerical parameters
C
C Outputs:AR(L) - elemental areas
C         DX(3,L) - element based array of deltax for each local node
C         DY(3,L) - element based array of deltay for each local node
C         NFTR - sparse factor
C         IQ(I),JQ(I*NFTR) - SPRSPAK STANDARD pointer arrays
C         SV(I) - stationary diagonal mass matrix
C         SH(NB1,I) - transpose of the stationary depth matrix
C         EMAT(NB1,M) - matrix used for the assymetric terms resulting
C          from the elevation boundary conditions.
C          (note:  The current implementation renders SH symmetric and
C           EMAT assymetric, which facilitates the storage of SH in
C           NB1=NHBW+1 [instead of NB2=2*NHBW+1] columns.  To implement
C           SH in its nonsymmetric form [and hence eliminate the need
C           for EMAT], the following evolution must be conducted in
C           this subroutine, the MAIN PROGRAM, and SUBROUTINE
C           ELEVATIONS2:
C             1. Uncomment the currently commented lines between the
C                first C+++ and the C>>> for each case of 2.)
C             2. Comment Out (or remove) lines between the C>>> and C<<<
C                comment lines.
C
C  Wave equation RHS terms
C         Amat, Bmat   - elevation terms
C         Cxmat, Cymat - tau0 and f terms
C         Dxmat, Dymat - bottom stress
c                        Dxmat, Dymat are only stationary if the 
c                        linearized bottom stress (AK) is time 
c                        independent
c         G0xmat, G0ymat-gravity term in momentum equation
C         COR - Corilios parameter
C         AK  - linearized bottom stress assumed independent of time
C         ENZM - vertical eddy viscosity 
C         BCFTR - boundary nodes for turning off Coriolis term
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, NE, NHBW, NFTR, NR
      REAL    G
      REAL    DELT, TAU0, THETA
      INTEGER IN(3,NEDIM), LRVS(NNDIM)
      REAL    X(NNDIM), Y(NNDIM), HDOWN(NNDIM)
      REAL    AR(NEDIM),DX(3,NEDIM),DY(3,NEDIM)
      INTEGER IQ(NNDIM), JQ(NFTRDIM*NNDIM)
      REAL    AMAT(NFTRDIM*NNDIM),  BMAT(NFTRDIM*NNDIM)
      REAL    CXMAT(NFTRDIM*NNDIM), CYMAT(NFTRDIM*NNDIM)
      REAL    DXMAT(NFTRDIM*NNDIM), DYMAT(NFTRDIM*NNDIM)
      Real    G0xmat(NFTRDIM*NNDIM),G0ymat(NFTRDIM*NNDIM)
      REAL    SV(NNDIM)
      REAL    COR
      REAL    AK(NNDIM)
      REAL    ENZM(NNDIM,NNVDIM)
      REAL    BCFTR(NNDIM)
      REAL SH(NB1DIM,NNDIM),EMAT(NBWE,NRDIM)

      COMPLEX D1(NNVDIM*NNDIM), D2(NNVDIM*NNDIM), D3(NNVDIM*NNDIM)
      REAL DZ(NNVDIM*NNDIM)
      REAL ZMID(NNDIM,NNVDIM)

C
C Function 
      INTEGER KAY 
C Local assignments
C
      REAL      EPSN 
      PARAMETER (EPSN = 0.5) 
      INTEGER I, J,L, I1,I2,I3, II,JJ, IJ, M, N, K
      INTEGER IDIAG
      INTEGER JP, JM 
      INTEGER KBOT 
      INTEGER NB1
      REAL   ARR, SUMHD
      REAL DUMMY(NNDIM)
      REAL    D33, D11
      REAL DZT(NNVDIM),DZB(NNVDIM),DZJ(NNVDIM)
      REAL ENZMI(NNVDIM)
      real    SUMX(NNDIM),SUMY(NNDIM), DSBDRY(NNDIM)
      integer NBE, NGEO, NVN, LST,  KST, NPCN
      integer INBE(4,NBEDIM), INGEO(2,NBEDIM)
      integer NORM(NNDIM)
      real gtang,si,co
C big dummy array for call to CTHOMAS
      COMPLEX RZ(NNVDIM*nndim)

      COMPLEX EYE

      REAL ZERO
      DATA ZERO/0.0/

      EYE=CMPLX(0.0,1.0)
C
C START OF EXECUTABLE CODE
C
C>>>Initialize for symmetric storage
      NB1=NHBW+1
      DO I=1,NN
         DO J=1,NB1
            SH(J,I)=0.0
         ENDDO
         SV(I)=0.0           
      ENDDO
C<<<
C
C Initialize sparse matrix arrays
C
      DO I=1,NFTRDIM*NNDIM
        
	 AMAT(I)  = 0.0
	 BMAT(I)  = 0.0
	 CXMAT(I) = 0.0
	 CYMAT(I) = 0.0
	 DXMAT(I) = 0.0
	 DYMAT(I) = 0.0         
	 G0xmat(I)= 0.0
         G0ymat(I)= 0.0
      ENDDO
C
C Prepare for SPRSPAK storage
C
      CALL SPRSBLD1(IN,3,3,IQ,JQ,NN,NE,NFTRDIM,1)
      NFTR=IQ(NN)/NN+1
      IF(NFTR.GT.NFTRDIM)THEN
         WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE STATIONARYQ4 ',
     &   'DUE TO INSUFFICIENT SIZE OF NFTRDIM'
         WRITE(2,*)' NFTR   =',NFTR,   ',NFTRDIM   =',NFTRDIM
         STOP
      END IF
C                
C Compute elemental areas, deltax's, and deltay's
C                
      DO L=1,NE       
         I1=IN(1,L)          
         I2=IN(2,L)          
         I3=IN(3,L)          
         DX(1,L)=X(I2)-X(I3)          
         DX(2,L)=X(I3)-X(I1)          
         DX(3,L)=X(I1)-X(I2)          
         DY(1,L)=Y(I2)-Y(I3)          
         DY(2,L)=Y(I3)-Y(I1)          
         DY(3,L)=Y(I1)-Y(I2)          
         AR(L)=0.5*(X(I1)*DY(1,L)+X(I2)*DY(2,L)+X(I3)*DY(3,L))          
         IF(AR(L).LE.0.0)WRITE(2,*)' NEGATIVE AREA IN ELEMENT', L
      ENDDO
C
C Compute stationary matricies:
C
C  SH - transpose of the stationary depth matrix
C  SV - stationary diagonal mass matrix <Pi>
C
      DO L=1,NE       
         ARR=AR(L)       
         SUMHD=HDOWN(IN(1,L))+HDOWN(IN(2,L))+HDOWN(IN(3,L))
         DO I=1,3        
            II=IN(I,L)          
            JJ=1
            SH(JJ,II)=SH(JJ,II)+ARR*(1.+TAU0*DELT*0.5)/3.
            SV(II)=SV(II)+ARR/3.0
            DO J=1,3
               IJ=IN(J,L)
               M=II       
               IF(THETA.GT.ZERO)THEN
                  N=1+(IJ-II)
                  if(N.GE.1)THEN
                     SH(N,M)=SH(N,M)+DELT*DELT*G
     &                               *(DY(I,L)*DY(J,L)+DX(I,L)*DX(J,L))
     &                               *SUMHD*THETA/24.0/ARR              
                  ENDIF
               ENDIF
               K=KAY(II,IJ,IQ,JQ)
	       idiag = 0
	       if(i.eq.j) idiag = 1
C
      Amat(k)=Amat(k) + idiag*2.*arr/3.0
     & - (1.-theta)*delt*delt*g*sumhd*
     & ( dy(i,L)*dy(j,L) + dx(i,L)*dx(j,L) )/(arr*12.)
C
      Bmat(k)=Bmat(k) + idiag*(-1.+tau0*delt*0.5)*arr/3.0
     & - theta*delt*delt*0.5*g*sumhd*
     & ( dy(i,L)*dy(j,L) + dx(i,L)*dx(j,L) )/(arr*12.)
C
      Cxmat(k)=Cxmat(k) + delt*delt*( tau0*dy(i,L) + cor*dx(i,L))/6.
C
      Cymat(k)=Cymat(k) + delt*delt*(-tau0*dx(i,L) + cor*dy(i,L))/6.
C
      Dxmat(k)=Dxmat(k) - delt*delt*( ak(ij)*dy(i,L)/6.)
C
      Dymat(k)=Dymat(k) - delt*delt*(-ak(ij)*dx(i,L)/6.)
C
      G0xmat(k)=G0xmat(k) - 0.5*g*dy(j,L)/6.     
C
      G0ymat(k)=G0ymat(k) + 0.5*g*dx(j,L)/6.
C           
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc


            ENDDO
         ENDDO
      ENDDO
      
C********************************************************************
C Enforce normal velocity bcs on land boundary nodes
C********************************************************************
      if (kst.le.nndim)then
       do i=kst,nndim
        ii=norm(i)
	si=sumx(ii)
	co=sumy(ii)
	if (ii.gt.1)then
	 k1=iq(ii-1)+1
	else
	 k1=1
	endif
	k2=iq(ii)
	do k=k1,k2
	 gtang=-si*g0xmat(k)+co*g0ymat(k)
	 g0xmat(k)=-si*gtang
	 g0ymat(k)= co*gtang
	enddo
       enddo
      endif
C********************************************************************
C
C Correct SH for the type 1 BC's and construct EMAT
C  (Symmetric storage)
C
      IF(NR.EQ.0)GO TO 10
      IF(THETA.GT.ZERO)THEN
         CALL EMATBUILD(SH,EMAT,LRVS,NN,NHBW,NR,NBWE,NB1DIM)
      ELSE
         DO I=1,NR
            II=LRVS(I)
            SH(1,II)=1.0
         ENDDO
      ENDIF
C<<<
 10   CONTINUE
C  
C LU decompose SH (DUMMY is not modified for KKK=1) 
C
C+++      CALL SOLVETR(1,SH,DUMMY,NN,NHBW,NNDIM,NBWE)
C>>>Symmetric storage
      IF(THETA.GT.ZERO)CALL BANSOLTR(1,SH,DUMMY,NN,NB1,NNDIM,NB1DIM)
C<<<
      DO I = 1, NR
         ii=lrvs(i)
         call sprsrowzap(Amat, iq, ii)
         call sprsrowzap(Bmat, iq, ii)
         call sprsrowzap(Cxmat, iq, ii)
         call sprsrowzap(Cymat, iq, ii)
         call sprsrowzap(Dxmat, iq, ii) 
         call sprsrowzap(Dymat, iq, ii) 
      ENDDO


C
C Casco: epsn is hardwired at 1/2 (Crank-Nicolson).      
C
c     epsn=0.5
      
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C PART 1: Assemble and Factor the Complex 3-Diag. Matrix for 
C          2*[average of (U+iV) over the time step]
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C k is a counter for position in the global velocity vector
C
C
C Loop over horizontal nodes
C
      k=0
      DO 100 I=1,NN
C
C  Store some arrays in local 1-D arrays
C
         DO J=1,NNV
            ENZMI(J)=ENZM(I,J)
            JP=MIN0(J+1,NNV)
            JM=MAX0(J-1,1)
            DZT(J)=Zmid(I,JP)-Zmid(I,J)
            DZB(J)=Zmid(I,J) -Zmid(I,JM)
            DZJ(J)=0.5*(DZT(J)+DZB(J))
         ENDDO
C
C Assemble tridiagonal matrix for 2*[average of (U+iV) for the time step]
C
C    LHS from d/dt, vertical diffusion, and coriolis
C    RHS from d/dt, surface pressure gradient.
C    BCFTR nullifies the coriolis term along the coast:
C      BCFTR=0. for coastal nodes
C      BCFTR=1. for all others
C      This is a convenience to allow simple BC application later
c    Casco: we are not allowing velocity BC's; BCFTR is left coded but unnecessary
C    Casco: epsn is hardwired for crank-nicolson
C
         DO J=1,NNV
            k=k+1
            dz(k)=dzj(j)
            JP=MIN0(J+1,NNV)
            JM=MAX0(J-1,1)
C
            D11=0.
            D33=0.
c      
            IF(J.NE.1 )D11=-ENZMI(J-1)/DZB(J)
            IF(J.NE.NNV) D33=-ENZMI(J)/DZT(J)
            D1(k)=0.5*D11*2.*EPSN*sv(i)
            D3(k)=0.5*D33*2.*EPSN*sv(i)
            D2(k)=-(D1(k)+D3(k))
     &       + 0.5*(2./DELT+BCFTR(I)*EYE*COR)*sv(i)*dz(k)
         ENDDO
C
C   Apply bottom stress BC
C
         kbot=nnv*(i-1) + 1
         D2(kbot)=D2(kbot)+0.5*AK(I)*sv(i)
c
c end of assembly loop
c
  100 continue
C
C Decompose
C
         CALL CTHOMAS(1,D1,D2,D3,RZ,NNV*nn)
C
C
C
C END OF ROUTINE
C
      RETURN
      END
