C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSBLD1(IN,NDIM,NNE,IQ,JQ,NN,NE,MN,NDF)
C-----------------------------------------------------------------------
C  NOVEMBER 30, 1990
C  D.LYNCH REVISION OF K.PAULSEN REVISION OF WEIPING VERSION.
C  SIMPLE ADAPTATION OF SPRSBLDW1 BY REVERSING ORDER OF IN(,) ARGUMENTS.
C
C  BUILDS SPRSPAK ARRAYS IQ,JQ FROM ELEMENT INCIDENCE LIST IN(K,L)
C     WITH L=GLOBAL ELEMENT NUMBER, K=LOCAL NODE NUMBER
C     FOR NDF DEGREES OF FREEDOM AT EACH NODE -- i.e. 
C     NDF EQUATIONS AND NDF UNKNOWNS.
C     NDF IS THE SAME AT EACH NODE.
C
C  SAME AS SPRSBLDW1 EXCEPT FOR ORDER OF ARGUMENTS FOR THE IN ARRAY
C 
C  NDIM IS THE FIRST PHYSICAL DIMENSION OF IN; NDIM.GE.NNE
C  NNE INDICATES THAT IN(I,L) FOR I=1 THRU NNE ARE NODE NUMBERS; 
C     THE BALANCE (IF ANY) OF IN WILL BE IGNORED.  
C  IQ(I) IS THE ADDRESS IN QV OF THE LAST NONZERO ENTRY IN ROW I OF Q
C  JQ(K) IS THE COLUMN INDEX IN Q FOR QV(K)
C  NN IS THE ACTUAL NUMBER OF NODES
C  NE IS THE ACTUAL NUMBER OF ELEMENTS
C  MN EQUALS OR EXCEEDS THE MAXIMUM NUMBER OF NONZREO ENTRYS PER ROW, 
C     FOR A SCALAR (i.e. NDF=1) PROBLEM.
C  NDF IS THE DEGREES OF FREEDOM, ASSUMED THE SAME AT EACH NODE.
C
C  IT IS ASSUMED THAT IN MAIN PROGRAM, DIMENSIONS ARE AT LEAST 
C     IQ(NN*NDF), JQ((NN*NDF)*(MN*NDF)).  NOTE THAT JQ IS OVERSIZED; 
C     BUT THE ALGORITHM MAY USE THIS EXTRA SIZE.  
C
C  CALLS ARE MADE TO SPRSINSRTD AND SPRSBLD3DD
C-----------------------------------------------------------------------
      DIMENSION IN(NDIM,*),IQ(*),JQ(*)
C
C  INITIALIZE JQ
C
      DO I=1,NN*MN
         JQ(I)=0
      ENDDO
C
C  INITIALIZE JQ AS DIAGONAL MATRIX
C
      DO I=1,NN
         J=(I-1)*MN+1
         JQ(J)=I
      ENDDO
C
C  ELEMENT LOOP FOR JQ, ASSUMING NDF=1.
C  ALL ROWS ALLOWED MN ENTRIES.
C  TRAILING ZEROS WILL BE REMOVED BELOW.
C
      DO L=1,NE
         DO I=1,NNE
            II=IN(I,L)
            DO J=1,NNE
               JJ=IN(J,L)
               IF(I.NE.J)CALL SPRSINSRTD(II,JJ,JQ,MN)
            ENDDO
         ENDDO
      ENDDO
C
C FIND IQ FROM JQ
C
       K=0
       KK=0
       DO I=1,NN
          DO J=1,MN
             KK=KK+1
             IF(JQ(KK).NE.0)K=K+1
          ENDDO
          IQ(I)=K
       ENDDO      
C
C REARRANGE (CONDENSE) THE JQ -- REMOVE THE TRAILING ZEROS
C
      DO I=2,NN
         KMIN=IQ(I-1)+1
         KMAX=IQ(I)
         KJMIN=(I-1)*MN+1
         DO K=KMIN,KMAX
            JQ(K)=JQ(KJMIN)
            KJMIN=KJMIN+1
         ENDDO
      ENDDO
C
C IQ, JQ NOW CORRECT FOR NDF=1
C EXPAND FOR NDF.GT.1
C
      IF(NDF.GT.1)CALL SPRSBLD3DD(IQ,JQ,NN,NDF)
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSBLD3DD(IQ,JQ,NN,NDF)
C-----------------------------------------------------------------------
C BUILD 3D IDEX IQQ AND JQQ
C-----------------------------------------------------------------------
      DIMENSION IQ(*),JQ(*)
      DO I=NN,1,-1
         IF (I.NE.1) THEN
            II=IQ(I)-IQ(I-1)
         ELSE
            II=IQ(I)
         ENDIF
         DO J=NDF,1,-1
            IQ(NDF*(I-1)+J)=IQ(I)*NDF*NDF-(NDF-J)*NDF*II
         ENDDO
         DO K=1,II
            JK=IQ(NDF*I)+1-(K-1)*NDF
            IF (I.NE.1) THEN
               JJ=JQ(IQ(I)-K+1)*NDF+1
            ELSE
               JJ=JQ(IQ(I)/NDF-K+1)*NDF+1
            END IF
            DO J=1,NDF
               JK=JK-1
               JJ=JJ-1
C
C THE FOLLOWING 3 LINES FROM WEIPING'S CODE
C	   JQ(JK)=JJ
C	   JQ(JK-NDF*II)=JJ
C	   IF(NDF.EQ.3)	   JQ(JK-6*II)=JJ
C WORKED FOR NDF=2 OR 3.  REPLACED BY INDUCTION  
C WITH THE DO 15 LOOP BELOW (DRL 1 DECEMBER 1990).   
C IT IS NOT OBVIOUS THAT JQ REMAINS UNCORRUPTED IN THE 
C OVERWRITING.  
C
               DO L=1,NDF
                  JQ(JK-(L-1)*NDF*II)=JJ
               ENDDO
            ENDDO
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSCONV(I,PPX,PPY,IQ,JQ,U,V,Z,
     &  DUDX,DUDY,DVDX,DVDY,DZDX,DZDY,
     &  NNV,NNDIM)
C-----------------------------------------------------------------------
CCOMPUTES THE HORIZONTAL DERIVATIVES AT DEPTH, UNDER HORIZONTAL NODE #I.
C  PPX, PPY ARE THE SPRSPAK MATRICES FOR D/DX, D/DY
C  IQ,JQ ARE THE STANDARD POINTERS FOR PPX, PPY
C  U(I,J),V(I,J) ARE THE X- AND Y-VELOCITIES 
C      AT HORIZONTAL NODE I, VERTICAL NODE J
C  Z(I,J) IS THE VERTICAL POSITION
C      OF HORIZONTAL NODE I, VERTICAL NODE J
C  THE HORIZONTAL DERIVATIVE TERMS DUDX(J), ETC. ARE RETURNED 
C      FOR ALL NODES J UNDER THE GIVEN NODE I
C  NNV IS THE ACTUAL # OF NODES IN THE VERTICAL
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF U AND V 
C  THIS ROUTINE IS NONDESTRUCTIVE
C-----------------------------------------------------------------------
      DIMENSION PPX(*),PPY(*),IQ(*),JQ(*),
     &U(NNDIM,*),V(NNDIM,*),Z(NNDIM,*),
     &DUDX(*),DUDY(*),DVDX(*),DVDY(*),DZDX(*),DZDY(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO J=1,NNV
         DUDX(J)=0.
         DUDY(J)=0.
         DVDX(J)=0.
         DVDY(J)=0.
         DZDX(J)=0.
         DZDY(J)=0.
      ENDDO
      DO K=KMIN,KMAX
         PPXK=PPX(K)
         PPYK=PPY(K)
         DO J=1,NNV
            UJQKJ=U(JQ(K),J)
            VJQKJ=V(JQ(K),J)
            ZJQKJ=Z(JQ(K),J)
            DUDX(J)=DUDX(J)+PPXK*UJQKJ
            DUDY(J)=DUDY(J)+PPYK*UJQKJ
            DVDX(J)=DVDX(J)+PPXK*VJQKJ
            DVDY(J)=DVDY(J)+PPYK*VJQKJ
            DZDX(J)=DZDX(J)+PPXK*ZJQKJ
            DZDY(J)=DZDY(J)+PPYK*ZJQKJ
         ENDDO
      ENDDO
      RETURN
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE SPRSINSRTD(I,J,JQ,MN)
C-----------------------------------------------------------------------
C  INSERTS GLOBAL ADDRESS (I,J) INTO JQ(), 
C  ASSUMING EXACTLY MN ENTRIES IN EACH ROW.
C-----------------------------------------------------------------------
      DIMENSION JQ(*)
      KMIN=(I-1)*MN+1
      KMAX=I*MN
C
C  SEARCH FOR INSERTION POINT KK
C
      DO K=KMIN,KMAX
         KK=K
         IF(J.LT.JQ(K))GO TO 20
         IF(J.EQ.JQ(K))RETURN
         IF(JQ(K).EQ.0)GO TO 20
      ENDDO
C
C  INSERT J AT JQ(KK); BUMP JQ
C
C CHECK INSERTED BY K PAULSEN 11/29...NOT VALIDATED AT TIME OF
C INSERTION !!!!
C
   20 IF(JQ(KMAX).NE.0)THEN
         WRITE(*,*) 'MAX CONNECTIVITY MN=',MN,
     &   'VIOLATED AT I,J=',I,J,'IN SPRSINSRTD'
         STOP
      ENDIF
      DO K=KMAX-1,KK,-1
         JQ(K+1)=JQ(K)
      ENDDO
      JQ(KK)=J   
C
C  DONE
C
      RETURN     
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE SPRSINVMLT(QV,IQ,D,NN)
C-----------------------------------------------------------------------
C  REPLACES QV WITH INV(D)*QV
C  D IS A DIAGONAL MATRIX, LENGTH NN 
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C-----------------------------------------------------------------------
      DIMENSION QV(*),D(*),IQ(*)
      KMAX=0
      DO I=1,NN
         KMIN=KMAX+1
         KMAX=IQ(I)
         FTR=1./D(I)
         DO K=KMIN,KMAX
            QV(K)=FTR*QV(K)
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSMLT(QV,IQ,JQ,X,B,NN)
C-----------------------------------------------------------------------
C  MULTIPLIES QV*X; ANSWER RETURNS IN B
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C-----------------------------------------------------------------------
      DIMENSION QV(*),X(*),B(*),IQ(*),JQ(*)
      KMAX=0
      DO I=1,NN
         SUM=0.
         KMIN=KMAX+1
         KMAX=IQ(I)
         DO K=KMIN,KMAX
            SUM=SUM+QV(K)*X(JQ(K))
         ENDDO
         B(I)=SUM
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSMLTIN(I,QV,IQ,JQ,X,B,NNV,NNDIM)
C-----------------------------------------------------------------------
C  MULTIPLIES THE Ith ROW OF QV BY THE Jth COLUMN OF THE MATRIX X,
C     FOR ALL VALUES OF J FROM 1 TO NNV
C  ANSWER RETURNS IN THE VECTOR B
C     (OR: OBTAINS ONLY Bi OF Q*X=B)
C  IQ,JQ ARE THE STANDARD POINTERS FOR QV
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF X
C  THIS ROUTINE IS NONDESTRUCTIVE
C-----------------------------------------------------------------------
      DIMENSION QV(*),IQ(*),JQ(*),X(NNDIM,*),B(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO J=1,NNV
         B(J)=0.
      ENDDO
      DO K=KMIN,KMAX
         QVK=QV(K)
         JQK=JQ(K)
         DO J=1,NNV
            B(J)=B(J)+QVK*X(JQK,J)
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSMLTIN2(I,QV,IQ,JQ,X1,B1,X2,B2,NNV,NNDIM)
C-----------------------------------------------------------------------
C  MULTIPLIES THE Ith ROW OF QV BY THE Jth COLUMN OF THE MATRIX X,
C     FOR ALL VALUES OF J FROM 1 TO NNV
C  ANSWER RETURNS IN THE VECTOR B
C     (OR: OBTAINS ONLY Bi OF Q*X=B)
C  IQ,JQ ARE THE STANDARD POINTERS FOR QV
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF X
C  THIS ROUTINE IS NONDESTRUCTIVE
C-----------------------------------------------------------------------
      INTEGER IQ(*),JQ(*)
      REAL QV(*),X1(NNDIM,*),B1(*),X2(NNDIM,*),B2(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO J=1,NNV
         B1(J)=0.
         B2(J)=0.
      ENDDO
      DO K=KMIN,KMAX
         JQK=JQ(K)
         QVK=QV(K)
         DO J=1,NNV
            B1(J)=B1(J)+QVK*X1(JQK,J)
            B2(J)=B2(J)+QVK*X2(JQK,J)
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSMLTIN3(I,IQ,JQ,X,QV1,B1,QV2,B2,QV3,B3,NNV,NNDIM)
C-----------------------------------------------------------------------
C  MULTIPLIES THE Ith ROW OF QV BY THE Jth COLUMN OF THE MATRIX X,
C     FOR ALL VALUES OF J FROM 1 TO NNV
C  ANSWER RETURNS IN THE VECTOR B
C     (OR: OBTAINS ONLY Bi OF Q*X=B)
C  IQ,JQ ARE THE STANDARD POINTERS FOR QV
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF X
C  THIS ROUTINE IS NONDESTRUCTIVE
C-----------------------------------------------------------------------
      INTEGER IQ(*),JQ(*)
      REAL X(NNDIM,*),QV1(*),B1(*),QV2(*),B2(*),QV3(*),B3(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO J=1,NNV
         B1(J)=0.
         B2(J)=0.
         B3(J)=0.
      ENDDO
      DO K=KMIN,KMAX
         JQK=JQ(K)
         QVK1=QV1(K)
         QVK2=QV2(K)
         QVK3=QV3(K)
         DO J=1,NNV
            XJQKJ=X(JQK,J)
            B1(J)=B1(J)+QVK1*XJQKJ
            B2(J)=B2(J)+QVK2*XJQKJ
            B3(J)=B3(J)+QVK3*XJQKJ
         ENDDO
      ENDDO
      RETURN
      END
C  
C*********************************************************************
C*********************************************************************
C
      SUBROUTINE SPRSMLTadd(QV,IQ,JQ,X,B,NN)
C
C drl 24 oct 98 from sprsmlt
C  MULTIPLIES QV*X; ANSWER is **added to** B
C   (same as sprsmlt except initial values in B are not overwritten)
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C
c     DIMENSION QV(1),X(1),B(1),IQ(1),JQ(1)
      DIMENSION QV(*),X(*),B(*),IQ(*),JQ(*)
      KMAX=0
      DO 10 I=1,NN
      SUM=0.
      KMIN=KMAX+1
      KMAX=IQ(I)
      DO 20 K=KMIN,KMAX
   20 SUM=SUM+QV(K)*X(JQ(K))
   10 B(I)=b(i) + SUM
      RETURN
      END
C
C*********************************************************************
C*********************************************************************
C
      SUBROUTINE SPRSROWZAP(QV,IQ,I)
C
C  REPLACES ROW I OF QV WITH ZEROS
C
      DIMENSION QV(*),IQ(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO 20 K=KMIN,KMAX
   20 QV(K)=0.0
      RETURN
      END
C
C*********************************************************************
C*********************************************************************
C
C
C*********************************************************************
C*********************************************************************
C
      SUBROUTINE SPRSBLDW1(IN,NEDIM,NNE,IQ,JQ,NN,NE,MN,NDF)
C
C  NOVEMBER 30, 1990
C  D.LYNCH REVISION OF K.PAULSEN REVISION OF WEIPING VERSION
C
C  BUILDS SPRSPAK ARRAYS IQ,JQ FROM ELEMENT INCIDENCE LIST IN(L,K)
C     WITH L=GLOBAL ELEMENT NUMBER, K=LOCAL NODE NUMBER
C     FOR NDF DEGREES OF FREEDOM AT EACH NODE -- i.e. 
C     NDF EQUATIONS AND NDF UNKNOWNS.
C     NDF IS THE SAME AT EACH NODE.
C
C  SAME AS SPRSBLD1 EXCEPT FOR THE ORDER OF THE ARGUMENTS FOR THE IN ARRAY
C 
C  NEDIM IS THE FIRST PHYSICAL DIMENSION OF IN; NEDIM.GE.NE
C  NNE INDICATES THAT IN(L,I) FOR I=1 THRU NNE ARE NODE NUMBERS; 
C     THE BALANCE (IF ANY) OF IN WILL BE IGNORED.  
C  IQ(I) IS THE ADDRESS IN QV OF THE LAST NONZERO ENTRY IN ROW I OF Q
C  JQ(K) IS THE COLUMN INDEX IN Q FOR QV(K)
C  NN IS THE ACTUAL NUMBER OF NODES
C  NE IS THE ACTUAL NUMBER OF ELEMENTS
C  MN EQUALS OR EXCEEDS THE MAXIMUM NUMBER OF NONZREO ENTRYS PER ROW, 
C     FOR A SCALAR (i.e. NDF=1) PROBLEM.
C  NDF IS THE DEGREES OF FREEDOM, ASSUMED THE SAME AT EACH NODE.
C
C  IT IS ASSUMED THAT IN MAIN PROGRAM, DIMENSIONS ARE AT LEAST 
C     IQ(NN*NDF), JQ((NN*NDF)*(MN*NDF)).  NOTE THAT JQ IS OVERSIZED; 
C     BUT THE ALGORITHM MAY USE THIS EXTRA SIZE.  
C
C  CALLS ARE MADE TO SPRSINSRTD AND SPRSBLD3DD
C
      DIMENSION IN(NEDIM,*),IQ(*),JQ(*)
C
C  INITIALIZE JQ
C
       DO 5 I=1,NN*MN
       JQ(I)=0
    5 CONTINUE
C
C  INITIALIZE JQ AS DIAGONAL MATRIX
C
      DO 10 I=1,NN
      J=(I-1)*MN+1
      JQ(J)=I
   10 CONTINUE
C
C  ELEMENT LOOP FOR JQ, ASSUMING NDF=1.
C  ALL ROWS ALLOWED MN ENTRIES.
C  TRAILING ZEROS WILL BE REMOVED BELOW.
C
      DO 100 L=1,NE
      DO 90 I=1,NNE
      II=IN(L,I)
      DO 90 J=1,NNE
      JJ=IN(L,J)
      IF(I.NE.J) CALL SPRSINSRTD(II,JJ,JQ,MN)
   90 CONTINUE
  100 CONTINUE
C
C FIND IQ FROM JQ
C
       K=0
       KK=0
       DO 110 I=1,NN
       DO 120 J=1,MN
       KK=KK+1
       IF(JQ(KK).NE.0) K=K+1
  120 CONTINUE
       IQ(I)=K
  110 CONTINUE      
C
C REARRANGE (CONDENSE) THE JQ -- REMOVE THE TRAILING ZEROS
C
       DO 140 I=2,NN
        KMIN=IQ(I-1)+1
        KMAX=IQ(I)
        KJMIN=(I-1)*MN+1
        DO 130 K=KMIN,KMAX
          JQ(K)=JQ(KJMIN)
          KJMIN=KJMIN+1
  130 CONTINUE
  140 CONTINUE
C
C IQ, JQ NOW CORRECT FOR NDF=1
C EXPAND FOR NDF.GT.1
C
        IF (NDF.GT.1) THEN
	CALL SPRSBLD3DD(IQ,JQ,NN,NDF)
	ELSE
	 GO TO 200
        END IF
  200 RETURN
      END
C
C*********************************************************************
C*********************************************************************
C
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE SPRSpreMLT(QV,IQ,JQ,X,B,NN)
C-----------------------------------------------------------------------
C  MULTIPLIES x*QV; ANSWER RETURNS IN B.  
C  Initial data in B is overwritten
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C-----------------------------------------------------------------------
C  28 Oct 98, drl: for Casco.  From SPRSMLT.
C-----------------------------------------------------------------------
c     DIMENSION QV(*),X(*),B(*), IQ(*),JQ(*)
      REAL QV(*),X(*),B(*)
      INTEGER IQ(*),JQ(*)
      INTEGER NN
      INTEGER KMAX, KMIN, I, K

      KMAX=0
      do i=1,nn
      b(i)=0.
      enddo
      DO I=1,NN
         KMIN=KMAX+1
         KMAX=IQ(I)
         DO K=KMIN,KMAX
           b(jq(k)) = b(jq(k)) + x(i)*QV(K)
         ENDDO
      ENDDO
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE SPRSpreMLTadd(QV,IQ,JQ,X,B,NN)
C-----------------------------------------------------------------------
C  MULTIPLIES x*QV; adds answer to B
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C-----------------------------------------------------------------------
C  28 Oct 98, drl: for Casco.  From SPRpreSMLT.
C-----------------------------------------------------------------------
c     DIMENSION QV(*),X(*),B(*),IQ(*),JQ(*)
      REAL  QV(*),X(*),B(*)
      INTEGER IQ(*),JQ(*)
      INTEGER NN

      INTEGER KMAX, KMIN, I, K

      KMAX=0
      DO I=1,NN
         KMIN=KMAX+1
         KMAX=IQ(I)
         DO K=KMIN,KMAX
           b(jq(k)) = b(jq(k)) + x(i)*QV(K)
         ENDDO
      ENDDO
      RETURN
      END
