      SUBROUTINE SCRSHO(FX)
      PARAMETER (ISCR=60,JSCR=21)
      CHARACTER*1 SCR(ISCR,JSCR),BLANK,ZERO,YY,XX,FF
      DIMENSION Y(ISCR)
      DATA BLANK,ZERO,YY,XX,FF/' ','-','l','-','x'/
1     CONTINUE
      WRITE (*,*) ' Enter X1,X2 (= to stop)'
      READ (*,*) X1,X2
      IF(X1.EQ.X2) RETURN
      DO 11 J=1,JSCR
        SCR(1,J)=YY
        SCR(ISCR,J)=YY
11    CONTINUE
      DO 13 I=2,ISCR-1
        SCR(I,1)=XX
        SCR(I,JSCR)=XX
        DO 12 J=2,JSCR-1
          SCR(I,J)=BLANK
12      CONTINUE
13    CONTINUE
      DX=(X2-X1)/(ISCR-1)
      X=X1
      YBIG=0.
      YSML=YBIG
      DO 14 I=1,ISCR
        Y(I)=FX(X)
        IF(Y(I).LT.YSML) YSML=Y(I)
        IF(Y(I).GT.YBIG) YBIG=Y(I)
        X=X+DX
14    CONTINUE
      IF(YBIG.EQ.YSML) YBIG=YSML+1.
      DYJ=(JSCR-1)/(YBIG-YSML)
      JZ=1-YSML*DYJ
      DO 15 I=1,ISCR
        SCR(I,JZ)=ZERO
        J=1+(Y(I)-YSML)*DYJ
        SCR(I,J)=FF
15    CONTINUE
      WRITE (*,'(1X,1PE10.3,1X,80A1)') YBIG,(SCR(I,JSCR),I=1,ISCR)
      DO 16 J=JSCR-1,2,-1
        WRITE (*,'(12X,80A1)') (SCR(I,J),I=1,ISCR)
16    CONTINUE
      WRITE (*,'(1X,1PE10.3,1X,80A1)') YSML,(SCR(I,1),I=1,ISCR)
      WRITE (*,'(12X,1PE10.3,40X,E10.3)') X1,X2
      GOTO 1
      END
      PROGRAM SFROID
      PARAMETER(NE=3,M=41,NB=1,NCI=NE,NCJ=NE-NB+1,NCK=M+1,NSI=NE,
     *    NSJ=2*NE+1,NYJ=NE,NYK=M)
      COMMON X(M),H,MM,N,C2,ANORM
      DIMENSION SCALV(NE),INDEXV(NE),Y(NE,M),C(NCI,NCJ,NCK),S(NSI,NSJ)
      ITMAX=100
      CONV=5.E-6
      SLOWC=1.
      H=1./(M-1)
      C2=0.
      WRITE(*,*)'ENTER M,N'
      READ(*,*)MM,N
      IF(MOD(N+MM,2).EQ.1)THEN
        INDEXV(1)=1
        INDEXV(2)=2
        INDEXV(3)=3
      ELSE
        INDEXV(1)=2
        INDEXV(2)=1
        INDEXV(3)=3
      ENDIF
      ANORM=1.
      IF(MM.NE.0)THEN
        Q1=N
        DO 11 I=1,MM
          ANORM=-.5*ANORM*(N+I)*(Q1/I)
          Q1=Q1-1.
11      CONTINUE
      ENDIF
      DO 12 K=1,M-1
        X(K)=(K-1)*H
        FAC1=1.-X(K)**2
        FAC2=FAC1**(-MM/2.)
        Y(1,K)=PLGNDR(N,MM,X(K))*FAC2
        DERIV=-((N-MM+1)*PLGNDR(N+1,MM,X(K))-(N+1)*
     *      X(K)*PLGNDR(N,MM,X(K)))/FAC1
        Y(2,K)=MM*X(K)*Y(1,K)/FAC1+DERIV*FAC2
        Y(3,K)=N*(N+1)-MM*(MM+1)
12    CONTINUE
      X(M)=1.
      Y(1,M)=ANORM
      Y(3,M)=N*(N+1)-MM*(MM+1)
      Y(2,M)=(Y(3,M)-C2)*Y(1,M)/(2.*(MM+1.))
      SCALV(1)=ABS(ANORM)
      SCALV(2)=MAX(ABS(ANORM),Y(2,M))
      SCALV(3)=MAX(1.,Y(3,M))
1     CONTINUE
      WRITE (*,*) 'ENTER C**2 OR 999 TO END'
      READ (*,*) C2
      IF(C2.EQ.999.)STOP
      CALL SOLVDE(ITMAX,CONV,SLOWC,SCALV,INDEXV,NE,NB,M,Y,NYJ,NYK,
     *    C,NCI,NCJ,NCK,S,NSI,NSJ)
      WRITE (*,*) ' M = ',MM,'  N = ',N,
     *    '  C**2 = ',C2,'  LAMBDA = ',Y(3,1)+MM*(MM+1)
      GO TO 1
      END
      SUBROUTINE SHELL(N,ARR)
      PARAMETER (ALN2I=1.4426950, TINY=1.E-5)
      DIMENSION ARR(N)
      LOGNB2=INT(ALOG(FLOAT(N))*ALN2I+TINY)
      M=N
      DO 12 NN=1,LOGNB2
        M=M/2
        K=N-M
        DO 11 J=1,K
          I=J
3         CONTINUE
          L=I+M
          IF(ARR(L).LT.ARR(I)) THEN
            T=ARR(I)
            ARR(I)=ARR(L)
            ARR(L)=T
            I=I-M
            IF(I.GE.1)GO TO 3
          ENDIF
11      CONTINUE
12    CONTINUE
      RETURN
      END
      SUBROUTINE SHOOT(NVAR,V,DELV,N2,X1,X2,EPS,H1,HMIN,F,DV)
      PARAMETER (NP=20)
      DIMENSION V(N2),DELV(N2),F(N2),DV(N2),Y(NP),DFDV(NP,NP),INDX(NP)
      EXTERNAL DERIVS,RKQC
      CALL LOAD(X1,V,Y)
      CALL ODEINT(Y,NVAR,X1,X2,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
      CALL SCORE(X2,Y,F)
      DO 12 IV=1,N2
        SAV=V(IV)
        V(IV)=V(IV)+DELV(IV)
        CALL LOAD(X1,V,Y)
        CALL ODEINT(Y,NVAR,X1,X2,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
        CALL SCORE(X2,Y,DV)
        DO 11 I=1,N2
          DFDV(I,IV)=(DV(I)-F(I))/DELV(IV)
11      CONTINUE
        V(IV)=SAV
12    CONTINUE
      DO 13 IV=1,N2
        DV(IV)=-F(IV)
13    CONTINUE
      CALL LUDCMP(DFDV,N2,NP,INDX,DET)
      CALL LUBKSB(DFDV,N2,NP,INDX,DV)
      DO 14 IV=1,N2
        V(IV)=V(IV)+DV(IV)
14    CONTINUE
      RETURN
      END
      SUBROUTINE SHOOTF(NVAR,V1,V2,DELV1,DELV2,N1,N2,X1,X2,XF,EPS,H1,HMI
     *N,F,DV1,DV2)
      PARAMETER (NP=20)
      DIMENSION V1(N2),DELV1(N2),V2(N1),DELV2(N1),F(NVAR),DV1(N2),
     *    DV2(N1),Y(NP),F1(NP),F2(NP),DFDV(NP,NP),INDX(NP)
      EXTERNAL DERIVS,RKQC
      CALL LOAD1(X1,V1,Y)
      CALL ODEINT(Y,NVAR,X1,XF,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
      CALL SCORE(XF,Y,F1)
      CALL LOAD2(X2,V2,Y)
      CALL ODEINT(Y,NVAR,X2,XF,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
      CALL SCORE(XF,Y,F2)
      J=0
      DO 12 IV=1,N2
        J=J+1
        SAV=V1(IV)
        V1(IV)=V1(IV)+DELV1(IV)
        CALL LOAD1(X1,V1,Y)
        CALL ODEINT(Y,NVAR,X1,XF,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
        CALL SCORE(XF,Y,F)
        DO 11 I=1,NVAR
          DFDV(I,J)=(F(I)-F1(I))/DELV1(IV)
11      CONTINUE
        V1(IV)=SAV
12    CONTINUE
      DO 14 IV=1,N1
        J=J+1
        SAV=V2(IV)
        V2(IV)=V2(IV)+DELV2(IV)
        CALL LOAD2(X2,V2,Y)
        CALL ODEINT(Y,NVAR,X2,XF,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
        CALL SCORE(XF,Y,F)
        DO 13 I=1,NVAR
          DFDV(I,J)=(F2(I)-F(I))/DELV2(IV)
13      CONTINUE
        V2(IV)=SAV
14    CONTINUE
      DO 15 I=1,NVAR
        F(I)=F1(I)-F2(I)
        F1(I)=-F(I)
15    CONTINUE
      CALL LUDCMP(DFDV,NVAR,NP,INDX,DET)
      CALL LUBKSB(DFDV,NVAR,NP,INDX,F1)
      J=0
      DO 16 IV=1,N2
        J=J+1
        V1(IV)=V1(IV)+F1(J)
        DV1(IV)=F1(J)
16    CONTINUE
      DO 17 IV=1,N1
        J=J+1
        V2(IV)=V2(IV)+F1(J)
        DV2(IV)=F1(J)
17    CONTINUE
      RETURN
      END
      SUBROUTINE SIMP1(A,MP,NP,MM,LL,NLL,IABF,KP,BMAX)
      DIMENSION A(MP,NP),LL(NP)
      KP=LL(1)
      BMAX=A(MM+1,KP+1)
      IF(NLL.LT.2)RETURN
      DO 11 K=2,NLL
        IF(IABF.EQ.0)THEN
          TEST=A(MM+1,LL(K)+1)-BMAX
        ELSE
          TEST=ABS(A(MM+1,LL(K)+1))-ABS(BMAX)
        ENDIF
        IF(TEST.GT.0.)THEN
          BMAX=A(MM+1,LL(K)+1)
          KP=LL(K)
        ENDIF
11    CONTINUE
      RETURN
      END
      SUBROUTINE SIMP2(A,M,N,MP,NP,L2,NL2,IP,KP,Q1)
      PARAMETER (EPS=1.E-6)
      DIMENSION A(MP,NP),L2(MP)
      IP=0
      IF(NL2.LT.1)RETURN
      DO 11 I=1,NL2
        IF(A(L2(I)+1,KP+1).LT.-EPS)GO TO 2
11    CONTINUE
      RETURN
2     Q1=-A(L2(I)+1,1)/A(L2(I)+1,KP+1)
      IP=L2(I)
      IF(I+1.GT.NL2)RETURN
      DO 13 I=I+1,NL2
        II=L2(I)
        IF(A(II+1,KP+1).LT.-EPS)THEN
          Q=-A(II+1,1)/A(II+1,KP+1)
          IF(Q.LT.Q1)THEN
            IP=II
            Q1=Q
          ELSE IF (Q.EQ.Q1) THEN
            DO 12 K=1,N
              QP=-A(IP+1,K+1)/A(IP+1,KP+1)
              Q0=-A(II+1,K+1)/A(II+1,KP+1)
              IF(Q0.NE.QP)GO TO 6
12          CONTINUE
6           IF(Q0.LT.QP)IP=II
          ENDIF
        ENDIF
13    CONTINUE
      RETURN
      END
      SUBROUTINE SIMP3(A,MP,NP,I1,K1,IP,KP)
      DIMENSION A(MP,NP)
      PIV=1./A(IP+1,KP+1)
      IF(I1.GE.0)THEN
        DO 12 II=1,I1+1
          IF(II-1.NE.IP)THEN
            A(II,KP+1)=A(II,KP+1)*PIV
            DO 11 KK=1,K1+1
              IF(KK-1.NE.KP)THEN
                A(II,KK)=A(II,KK)-A(IP+1,KK)*A(II,KP+1)
              ENDIF
11          CONTINUE
          ENDIF
12      CONTINUE
      ENDIF
      DO 13 KK=1,K1+1
        IF(KK-1.NE.KP)A(IP+1,KK)=-A(IP+1,KK)*PIV
13    CONTINUE
      A(IP+1,KP+1)=PIV
      RETURN
      END
      SUBROUTINE SIMPLX(A,M,N,MP,NP,M1,M2,M3,ICASE,IZROV,IPOSV)
      PARAMETER(MMAX=100,EPS=1.E-6)
      DIMENSION A(MP,NP),IZROV(N),IPOSV(M),L1(MMAX),L2(MMAX),L3(MMAX)
      IF(M.NE.M1+M2+M3)PAUSE 'Bad input constraint counts.'
      NL1=N
      DO 11 K=1,N
        L1(K)=K
        IZROV(K)=K
11    CONTINUE
      NL2=M
      DO 12 I=1,M
        IF(A(I+1,1).LT.0.)PAUSE 'Bad input tableau.'
        L2(I)=I
        IPOSV(I)=N+I
12    CONTINUE
      DO 13 I=1,M2
        L3(I)=1
13    CONTINUE
      IR=0
      IF(M2+M3.EQ.0)GO TO 30
      IR=1
      DO 15 K=1,N+1
        Q1=0.
        DO 14 I=M1+1,M
          Q1=Q1+A(I+1,K)
14      CONTINUE
        A(M+2,K)=-Q1
15    CONTINUE
10    CALL SIMP1(A,MP,NP,M+1,L1,NL1,0,KP,BMAX)
      IF(BMAX.LE.EPS.AND.A(M+2,1).LT.-EPS)THEN
        ICASE=-1
        RETURN
      ELSE IF(BMAX.LE.EPS.AND.A(M+2,1).LE.EPS)THEN
        M12=M1+M2+1
        IF(M12.LE.M)THEN
          DO 16 IP=M12,M
            IF(IPOSV(IP).EQ.IP+N)THEN
              CALL SIMP1(A,MP,NP,IP,L1,NL1,1,KP,BMAX)
              IF(BMAX.GT.0.)GO TO 1
            ENDIF
16        CONTINUE
        ENDIF
        IR=0
        M12=M12-1
        IF(M1+1.GT.M12)GO TO 30
        DO 18 I=M1+1,M12
          IF(L3(I-M1).EQ.1)THEN
            DO 17 K=1,N+1
              A(I+1,K)=-A(I+1,K)
17          CONTINUE
          ENDIF
18      CONTINUE
        GO TO 30
      ENDIF
      CALL SIMP2(A,M,N,MP,NP,L2,NL2,IP,KP,Q1)
      IF(IP.EQ.0)THEN
        ICASE=-1
        RETURN
      ENDIF
1     CALL SIMP3(A,MP,NP,M+1,N,IP,KP)
      IF(IPOSV(IP).GE.N+M1+M2+1)THEN
        DO 19 K=1,NL1
          IF(L1(K).EQ.KP)GO TO 2
19      CONTINUE
2       NL1=NL1-1
        DO 21 IS=K,NL1
          L1(IS)=L1(IS+1)
21      CONTINUE
      ELSE
        IF(IPOSV(IP).LT.N+M1+1)GO TO 20
        KH=IPOSV(IP)-M1-N
        IF(L3(KH).EQ.0)GO TO 20
        L3(KH)=0
      ENDIF
      A(M+2,KP+1)=A(M+2,KP+1)+1.
      DO 22 I=1,M+2
        A(I,KP+1)=-A(I,KP+1)
22    CONTINUE
20    IS=IZROV(KP)
      IZROV(KP)=IPOSV(IP)
      IPOSV(IP)=IS
      IF(IR.NE.0)GO TO 10
30    CALL SIMP1(A,MP,NP,0,L1,NL1,0,KP,BMAX)
      IF(BMAX.LE.0.)THEN
        ICASE=0
        RETURN
      ENDIF
      CALL SIMP2(A,M,N,MP,NP,L2,NL2,IP,KP,Q1)
      IF(IP.EQ.0)THEN
        ICASE=1
        RETURN
      ENDIF
      CALL SIMP3(A,MP,NP,M,N,IP,KP)
      GO TO 20
      END
      SUBROUTINE SINFT(Y,N)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION Y(N)
      THETA=3.14159265358979D0/DBLE(N)
      WR=1.0D0
      WI=0.0D0
      WPR=-2.0D0*DSIN(0.5D0*THETA)**2
      WPI=DSIN(THETA)
      Y(1)=0.0
      M=N/2
      DO 11 J=1,M
        WTEMP=WR
        WR=WR*WPR-WI*WPI+WR
        WI=WI*WPR+WTEMP*WPI+WI
        Y1=WI*(Y(J+1)+Y(N-J+1))
        Y2=0.5*(Y(J+1)-Y(N-J+1))
        Y(J+1)=Y1+Y2
        Y(N-J+1)=Y1-Y2
11    CONTINUE
      CALL REALFT(Y,M,+1)
      SUM=0.0
      Y(1)=0.5*Y(1)
      Y(2)=0.0
      DO 12 J=1,N-1,2
        SUM=SUM+Y(J)
        Y(J)=Y(J+1)
        Y(J+1)=SUM
12    CONTINUE
      RETURN
      END
      SUBROUTINE SMOOFT(Y,N,PTS)
      PARAMETER(MMAX=1024)
      DIMENSION Y(MMAX)
      M=2
      NMIN=N+2.*PTS
1     IF(M.LT.NMIN)THEN
        M=2*M
      GO TO 1
      ENDIF
      IF(M.GT.MMAX) PAUSE 'MMAX too small'
      CONST=(PTS/M)**2
      Y1=Y(1)
      YN=Y(N)
      RN1=1./(N-1.)
      DO 11 J=1,N
        Y(J)=Y(J)-RN1*(Y1*(N-J)+YN*(J-1))
11    CONTINUE
      IF(N+1.LE.M)THEN
        DO 12 J=N+1,M
          Y(J)=0.
12      CONTINUE
      ENDIF
      MO2=M/2
      CALL REALFT(Y,MO2,1)
      Y(1)=Y(1)/MO2
      FAC=1.
      DO 13 J=1,MO2-1
        K=2*J+1
        IF(FAC.NE.0.)THEN
          FAC=AMAX1(0.,(1.-CONST*J**2)/MO2)
          Y(K)=FAC*Y(K)
          Y(K+1)=FAC*Y(K+1)
        ELSE
          Y(K)=0.
          Y(K+1)=0.
        ENDIF
13    CONTINUE
      FAC=AMAX1(0.,(1.-0.25*PTS**2)/MO2)
      Y(2)=FAC*Y(2)
      CALL REALFT(Y,MO2,-1)
      DO 14 J=1,N
        Y(J)=RN1*(Y1*(N-J)+YN*(J-1))+Y(J)
14    CONTINUE
      RETURN
      END
      SUBROUTINE SNCNDN(UU,EMMC,SN,CN,DN)
      PARAMETER (CA=.0003)
      LOGICAL BO
      DIMENSION EM(13),EN(13)
      EMC=EMMC
      U=UU
      IF(EMC.NE.0.)THEN
        BO=(EMC.LT.0.)
        IF(BO)THEN
          D=1.-EMC
          EMC=-EMC/D
          D=SQRT(D)
          U=D*U
        ENDIF
        A=1.
        DN=1.
        DO 11 I=1,13
          L=I
          EM(I)=A
          EMC=SQRT(EMC)
          EN(I)=EMC
          C=0.5*(A+EMC)
          IF(ABS(A-EMC).LE.CA*A)GO TO 1
          EMC=A*EMC
          A=C
11      CONTINUE
1       U=C*U
        SN=SIN(U)
        CN=COS(U)
        IF(SN.EQ.0.)GO TO 2
        A=CN/SN
        C=A*C
        DO 12 II=L,1,-1
          B=EM(II)
          A=C*A
          C=DN*C
          DN=(EN(II)+A)/(B+A)
          A=C/B
12      CONTINUE
        A=1./SQRT(C**2+1.)
        IF(SN.LT.0.)THEN
          SN=-A
        ELSE
          SN=A
        ENDIF
        CN=C*SN
2       IF(BO)THEN
          A=DN
          DN=CN
          CN=A
          SN=SN/D
        ENDIF
      ELSE
        CN=1./COSH(U)
        DN=CN
        SN=TANH(U)
      ENDIF
      RETURN
      END
      SUBROUTINE SOLVDE(ITMAX,CONV,SLOWC,SCALV,INDEXV,NE,NB,M,
     *    Y,NYJ,NYK,C,NCI,NCJ,NCK,S,NSI,NSJ)
      PARAMETER (NMAX=10)
      DIMENSION Y(NYJ,NYK),C(NCI,NCJ,NCK),S(NSI,NSJ),SCALV(NYJ),INDEXV(N
     *YJ)
      DIMENSION ERMAX(NMAX),KMAX(NMAX)
      K1=1
      K2=M
      NVARS=NE*M
      J1=1
      J2=NB
      J3=NB+1
      J4=NE
      J5=J4+J1
      J6=J4+J2
      J7=J4+J3
      J8=J4+J4
      J9=J8+J1
      IC1=1
      IC2=NE-NB
      IC3=IC2+1
      IC4=NE
      JC1=1
      JCF=IC3
      DO 16 IT=1,ITMAX
        K=K1
        CALL DIFEQ(K,K1,K2,J9,IC3,IC4,INDEXV,NE,S,NSI,NSJ,Y,NYJ,NYK)
        CALL PINVS(IC3,IC4,J5,J9,JC1,K1,C,NCI,NCJ,NCK,S,NSI,NSJ)
        DO 11 K=K1+1,K2
          KP=K-1
          CALL DIFEQ(K,K1,K2,J9,IC1,IC4,INDEXV,NE,S,NSI,NSJ,Y,NYJ,NYK)
          CALL RED(IC1,IC4,J1,J2,J3,J4,J9,IC3,JC1,JCF,KP,
     *        C,NCI,NCJ,NCK,S,NSI,NSJ)
          CALL PINVS(IC1,IC4,J3,J9,JC1,K,C,NCI,NCJ,NCK,S,NSI,NSJ)
11      CONTINUE
        K=K2+1
        CALL DIFEQ(K,K1,K2,J9,IC1,IC2,INDEXV,NE,S,NSI,NSJ,Y,NYJ,NYK)
        CALL RED(IC1,IC2,J5,J6,J7,J8,J9,IC3,JC1,JCF,K2,
     *      C,NCI,NCJ,NCK,S,NSI,NSJ)
        CALL PINVS(IC1,IC2,J7,J9,JCF,K2+1,C,NCI,NCJ,NCK,S,NSI,NSJ)
        CALL BKSUB(NE,NB,JCF,K1,K2,C,NCI,NCJ,NCK)
        ERR=0.
        DO 13 J=1,NE
          JV=INDEXV(J)
          ERMAX(J)=0.
          ERRJ=0.
          KMAX(J)=0
          VMAX=0.
          DO 12 K=K1,K2
            VZ=ABS(C(J,1,K))
            IF(VZ.GT.VMAX) THEN
               VMAX=VZ
               KM=K
            ENDIF
            ERRJ=ERRJ+VZ
12        CONTINUE
          ERR=ERR+ERRJ/SCALV(JV)
          ERMAX(J)=C(J,1,KM)/SCALV(JV)
          KMAX(J)=KM
13      CONTINUE
        ERR=ERR/NVARS
        FAC=SLOWC/MAX(SLOWC,ERR)
        DO 15 JV=1,NE
          J=INDEXV(JV)
          DO 14 K=K1,K2
            Y(J,K)=Y(J,K)-FAC*C(JV,1,K)
14        CONTINUE
15      CONTINUE
        WRITE(*,100) IT,ERR,FAC,(KMAX(J),ERMAX(J),J=1,NE)
        IF(ERR.LT.CONV) RETURN
16    CONTINUE
      PAUSE 'ITMAX exceeded'
100   FORMAT(1X,I4,2F12.6,(/5X,I5,F12.6))
      RETURN
      END
      SUBROUTINE SOR(A,B,C,D,E,F,U,JMAX,RJAC)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(JMAX,JMAX),B(JMAX,JMAX),C(JMAX,JMAX),
     *    D(JMAX,JMAX),E(JMAX,JMAX),F(JMAX,JMAX),U(JMAX,JMAX)
      PARAMETER(MAXITS=1000,EPS=1.D-5,ZERO=0.D0,HALF=.5D0,QTR=.25D0,ONE=
     *1.D0)
      ANORMF=ZERO
      DO 12 J=2,JMAX-1
        DO 11 L=2,JMAX-1
          ANORMF=ANORMF+ABS(F(J,L))
11      CONTINUE
12    CONTINUE
      OMEGA=ONE
      DO 15 N=1,MAXITS
        ANORM=ZERO
        DO 14 J=2,JMAX-1
          DO 13 L=2,JMAX-1
            IF(MOD(J+L,2).EQ.MOD(N,2))THEN
              RESID=A(J,L)*U(J+1,L)+B(J,L)*U(J-1,L)+
     *            C(J,L)*U(J,L+1)+D(J,L)*U(J,L-1)+
     *            E(J,L)*U(J,L)-F(J,L)
              ANORM=ANORM+ABS(RESID)
              U(J,L)=U(J,L)-OMEGA*RESID/E(J,L)
            ENDIF
13        CONTINUE
14      CONTINUE
        IF(N.EQ.1) THEN
          OMEGA=ONE/(ONE-HALF*RJAC**2)
        ELSE
          OMEGA=ONE/(ONE-QTR*RJAC**2*OMEGA)
        ENDIF
        IF((N.GT.1).AND.(ANORM.LT.EPS*ANORMF))RETURN
15    CONTINUE
      PAUSE 'MAXITS exceeded'
      END
      SUBROUTINE SORT(N,RA)
      DIMENSION RA(N)
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          RRA=RA(L)
        ELSE
          RRA=RA(IR)
          RA(IR)=RA(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            RA(1)=RRA
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(RA(J).LT.RA(J+1))J=J+1
          ENDIF
          IF(RRA.LT.RA(J))THEN
            RA(I)=RA(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        RA(I)=RRA
      GO TO 10
      END
      SUBROUTINE SORT2(N,RA,RB)
      DIMENSION RA(N),RB(N)
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          RRA=RA(L)
          RRB=RB(L)
        ELSE
          RRA=RA(IR)
          RRB=RB(IR)
          RA(IR)=RA(1)
          RB(IR)=RB(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            RA(1)=RRA
            RB(1)=RRB
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(RA(J).LT.RA(J+1))J=J+1
          ENDIF
          IF(RRA.LT.RA(J))THEN
            RA(I)=RA(J)
            RB(I)=RB(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        RA(I)=RRA
        RB(I)=RRB
      GO TO 10
      END
      SUBROUTINE SORT3(N,RA,RB,RC,WKSP,IWKSP)
      DIMENSION RA(N),RB(N),RC(N),WKSP(N),IWKSP(N)
      CALL INDEXX(N,RA,IWKSP)
      DO 11 J=1,N
        WKSP(J)=RA(J)
11    CONTINUE
      DO 12 J=1,N
        RA(J)=WKSP(IWKSP(J))
12    CONTINUE
      DO 13 J=1,N
        WKSP(J)=RB(J)
13    CONTINUE
      DO 14 J=1,N
        RB(J)=WKSP(IWKSP(J))
14    CONTINUE
      DO 15 J=1,N
        WKSP(J)=RC(J)
15    CONTINUE
      DO 16 J=1,N
        RC(J)=WKSP(IWKSP(J))
16    CONTINUE
      RETURN
      END
      SUBROUTINE SPARSE(B,N,ASUB,ATSUB,X,RSQ)
      PARAMETER (NMAX=500,EPS=1.E-6)
      DIMENSION B(N),X(N),G(NMAX),H(NMAX),XI(NMAX),XJ(NMAX)
      EPS2=N*EPS**2
      IRST=0
1     IRST=IRST+1
      CALL ASUB(X,XI)
      RP=0.
      BSQ=0.
      DO 11 J=1,N
        BSQ=BSQ+B(J)**2
        XI(J)=XI(J)-B(J)
        RP=RP+XI(J)**2
11    CONTINUE
      CALL ATSUB(XI,G)
      DO 12 J=1,N
        G(J)=-G(J)
        H(J)=G(J)
12    CONTINUE
      DO 19 ITER=1,10*N
        CALL ASUB(H,XI)
        ANUM=0.
        ADEN=0.
        DO 13 J=1,N
          ANUM=ANUM+G(J)*H(J)
          ADEN=ADEN+XI(J)**2
13      CONTINUE
        IF(ADEN.EQ.0.)PAUSE 'very singular matrix'
        ANUM=ANUM/ADEN
        DO 14 J=1,N
          XI(J)=X(J)
          X(J)=X(J)+ANUM*H(J)
14      CONTINUE
        CALL ASUB(X,XJ)
        RSQ=0.
        DO 15 J=1,N
          XJ(J)=XJ(J)-B(J)
          RSQ=RSQ+XJ(J)**2
15      CONTINUE
        IF(RSQ.EQ.RP.OR.RSQ.LE.BSQ*EPS2)RETURN
        IF(RSQ.GT.RP)THEN
          DO 16 J=1,N
            X(J)=XI(J)
16        CONTINUE
          IF(IRST.GE.3)RETURN
          GO TO 1
        ENDIF
        RP=RSQ
        CALL ATSUB(XJ,XI)
        GG=0.
        DGG=0.
        DO 17 J=1,N
          GG=GG+G(J)**2
          DGG=DGG+(XI(J)+G(J))*XI(J)
17      CONTINUE
        IF(GG.EQ.0.)RETURN
        GAM=DGG/GG
        DO 18 J=1,N
          G(J)=-XI(J)
          H(J)=G(J)+GAM*H(J)
18      CONTINUE
19    CONTINUE
      PAUSE 'too many iterations'
      RETURN
      END
      SUBROUTINE SPCTRM(P,M,K,OVRLAP,W1,W2)
      LOGICAL OVRLAP
      DIMENSION P(M),W1(*),W2(M)
      WINDOW(J)=(1.-ABS(((J-1)-FACM)*FACP))
C     WINDOW(J)=1.
C     WINDOW(J)=(1.-(((J-1)-FACM)*FACP)**2)
      MM=M+M
      M4=MM+MM
      M44=M4+4
      M43=M4+3
      DEN=0.
      FACM=M-0.5
      FACP=1./(M+0.5)
      SUMW=0.
      DO 11 J=1,MM
        SUMW=SUMW+WINDOW(J)**2
11    CONTINUE
      DO 12 J=1,M
        P(J)=0.
12    CONTINUE
      IF(OVRLAP)THEN
        READ (9,*) (W2(J),J=1,M)
      ENDIF
      DO 18 KK=1,K
        DO 15 JOFF=-1,0,1
          IF (OVRLAP) THEN
            DO 13 J=1,M
              W1(JOFF+J+J)=W2(J)
13          CONTINUE
            READ (9,*) (W2(J),J=1,M)
            JOFFN=JOFF+MM
            DO 14 J=1,M
              W1(JOFFN+J+J)=W2(J)
14          CONTINUE
          ELSE
            READ (9,*) (W1(J),J=JOFF+2,M4,2)
          ENDIF
15      CONTINUE
        DO 16 J=1,MM
          J2=J+J
          W=WINDOW(J)
          W1(J2)=W1(J2)*W
          W1(J2-1)=W1(J2-1)*W
16      CONTINUE
        CALL FOUR1(W1,MM,1)
        P(1)=P(1)+W1(1)**2+W1(2)**2
        DO 17 J=2,M
          J2=J+J
          P(J)=P(J)+W1(J2)**2+W1(J2-1)**2
     *        +W1(M44-J2)**2+W1(M43-J2)**2
17      CONTINUE
        DEN=DEN+SUMW
18    CONTINUE
      DEN=M4*DEN
      DO 19 J=1,M
        P(J)=P(J)/DEN
19    CONTINUE
      RETURN
      END
      SUBROUTINE SPEAR(DATA1,DATA2,N,WKSP1,WKSP2,D,ZD,PROBD,RS,PROBRS)
      DIMENSION DATA1(N),DATA2(N),WKSP1(N),WKSP2(N)
      DO 11 J=1,N
        WKSP1(J)=DATA1(J)
        WKSP2(J)=DATA2(J)
11    CONTINUE
      CALL SORT2(N,WKSP1,WKSP2)
      CALL CRANK(N,WKSP1,SF)
      CALL SORT2(N,WKSP2,WKSP1)
      CALL CRANK(N,WKSP2,SG)
      D=0.
      DO 12 J=1,N
        D=D+(WKSP1(J)-WKSP2(J))**2
12    CONTINUE
      EN=N
      EN3N=EN**3-EN
      AVED=EN3N/6.-(SF+SG)/12.
      FAC=(1.-SF/EN3N)*(1.-SG/EN3N)
      VARD=((EN-1.)*EN**2*(EN+1.)**2/36.)*FAC
      ZD=(D-AVED)/SQRT(VARD)
      PROBD=ERFCC(ABS(ZD)/1.4142136)
      RS=(1.-(6./EN3N)*(D+0.5*(SF+SG)))/FAC
      T=RS*SQRT((EN-2.)/((1.+RS)*(1.-RS)))
      DF=EN-2.
      PROBRS=BETAI(0.5*DF,0.5,DF/(DF+T**2))
      RETURN
      END
      SUBROUTINE SPLIE2(X1A,X2A,YA,M,N,Y2A)
      PARAMETER (NN=100)
      DIMENSION X1A(M),X2A(N),YA(M,N),Y2A(M,N),YTMP(NN),Y2TMP(NN)
      DO 13 J=1,M
        DO 11 K=1,N
          YTMP(K)=YA(J,K)
11      CONTINUE
        CALL SPLINE(X2A,YTMP,N,1.E30,1.E30,Y2TMP)
        DO 12 K=1,N
          Y2A(J,K)=Y2TMP(K)
12      CONTINUE
13    CONTINUE
      RETURN
      END
      SUBROUTINE SPLIN2(X1A,X2A,YA,Y2A,M,N,X1,X2,Y)
      PARAMETER (NN=100)
      DIMENSION X1A(M),X2A(N),YA(M,N),Y2A(M,N),YTMP(NN),Y2TMP(NN),YYTMP(
     *NN)
      DO 12 J=1,M
        DO 11 K=1,N
          YTMP(K)=YA(J,K)
          Y2TMP(K)=Y2A(J,K)
11      CONTINUE
        CALL SPLINT(X2A,YTMP,Y2TMP,N,X2,YYTMP(J))
12    CONTINUE
      CALL SPLINE(X1A,YYTMP,M,1.E30,1.E30,Y2TMP)
      CALL SPLINT(X1A,YYTMP,Y2TMP,M,X1,Y)
      RETURN
      END
      SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2)
      PARAMETER (NMAX=100)
      DIMENSION X(N),Y(N),Y2(N),U(NMAX)
      IF (YP1.GT..99E30) THEN
        Y2(1)=0.
        U(1)=0.
      ELSE
        Y2(1)=-0.5
        U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
      ENDIF
      DO 11 I=2,N-1
        SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
        P=SIG*Y2(I-1)+2.
        Y2(I)=(SIG-1.)/P
        U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
     *      /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P
11    CONTINUE
      IF (YPN.GT..99E30) THEN
        QN=0.
        UN=0.
      ELSE
        QN=0.5
        UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
      ENDIF
      Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.)
      DO 12 K=N-1,1,-1
        Y2(K)=Y2(K)*Y2(K+1)+U(K)
12    CONTINUE
      RETURN
      END
      SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)
      DIMENSION XA(N),YA(N),Y2A(N)
      KLO=1
      KHI=N
1     IF (KHI-KLO.GT.1) THEN
        K=(KHI+KLO)/2
        IF(XA(K).GT.X)THEN
          KHI=K
        ELSE
          KLO=K
        ENDIF
      GOTO 1
      ENDIF
      H=XA(KHI)-XA(KLO)
      IF (H.EQ.0.) PAUSE 'Bad XA input.'
      A=(XA(KHI)-X)/H
      B=(X-XA(KLO))/H
      Y=A*YA(KLO)+B*YA(KHI)+
     *      ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.
      RETURN
      END
      SUBROUTINE SVBKSB(U,W,V,M,N,MP,NP,B,X)
      PARAMETER (NMAX=100)
      DIMENSION U(MP,NP),W(NP),V(NP,NP),B(MP),X(NP),TMP(NMAX)
      DO 12 J=1,N
        S=0.
        IF(W(J).NE.0.)THEN
          DO 11 I=1,M
            S=S+U(I,J)*B(I)
11        CONTINUE
          S=S/W(J)
        ENDIF
        TMP(J)=S
12    CONTINUE
      DO 14 J=1,N
        S=0.
        DO 13 JJ=1,N
          S=S+V(J,JJ)*TMP(JJ)
13      CONTINUE
        X(J)=S
14    CONTINUE
      RETURN
      END
      SUBROUTINE SVDCMP(A,M,N,MP,NP,W,V)
      PARAMETER (NMAX=100)
      DIMENSION A(MP,NP),W(NP),V(NP,NP),RV1(NMAX)
      G=0.0
      SCALE=0.0
      ANORM=0.0
      DO 25 I=1,N
        L=I+1
        RV1(I)=SCALE*G
        G=0.0
        S=0.0
        SCALE=0.0
        IF (I.LE.M) THEN
          DO 11 K=I,M
            SCALE=SCALE+ABS(A(K,I))
11        CONTINUE
          IF (SCALE.NE.0.0) THEN
            DO 12 K=I,M
              A(K,I)=A(K,I)/SCALE
              S=S+A(K,I)*A(K,I)
12          CONTINUE
            F=A(I,I)
            G=-SIGN(SQRT(S),F)
            H=F*G-S
            A(I,I)=F-G
            IF (I.NE.N) THEN
              DO 15 J=L,N
                S=0.0
                DO 13 K=I,M
                  S=S+A(K,I)*A(K,J)
13              CONTINUE
                F=S/H
                DO 14 K=I,M
                  A(K,J)=A(K,J)+F*A(K,I)
14              CONTINUE
15            CONTINUE
            ENDIF
            DO 16 K= I,M
              A(K,I)=SCALE*A(K,I)
16          CONTINUE
          ENDIF
        ENDIF
        W(I)=SCALE *G
        G=0.0
        S=0.0
        SCALE=0.0
        IF ((I.LE.M).AND.(I.NE.N)) THEN
          DO 17 K=L,N
            SCALE=SCALE+ABS(A(I,K))
17        CONTINUE
          IF (SCALE.NE.0.0) THEN
            DO 18 K=L,N
              A(I,K)=A(I,K)/SCALE
              S=S+A(I,K)*A(I,K)
18          CONTINUE
            F=A(I,L)
            G=-SIGN(SQRT(S),F)
            H=F*G-S
            A(I,L)=F-G
            DO 19 K=L,N
              RV1(K)=A(I,K)/H
19          CONTINUE
            IF (I.NE.M) THEN
              DO 23 J=L,M
                S=0.0
                DO 21 K=L,N
                  S=S+A(J,K)*A(I,K)
21              CONTINUE
                DO 22 K=L,N
                  A(J,K)=A(J,K)+S*RV1(K)
22              CONTINUE
23            CONTINUE
            ENDIF
            DO 24 K=L,N
              A(I,K)=SCALE*A(I,K)
24          CONTINUE
          ENDIF
        ENDIF
        ANORM=MAX(ANORM,(ABS(W(I))+ABS(RV1(I))))
25    CONTINUE
      DO 32 I=N,1,-1
        IF (I.LT.N) THEN
          IF (G.NE.0.0) THEN
            DO 26 J=L,N
              V(J,I)=(A(I,J)/A(I,L))/G
26          CONTINUE
            DO 29 J=L,N
              S=0.0
              DO 27 K=L,N
                S=S+A(I,K)*V(K,J)
27            CONTINUE
              DO 28 K=L,N
                V(K,J)=V(K,J)+S*V(K,I)
28            CONTINUE
29          CONTINUE
          ENDIF
          DO 31 J=L,N
            V(I,J)=0.0
            V(J,I)=0.0
31        CONTINUE
        ENDIF
        V(I,I)=1.0
        G=RV1(I)
        L=I
32    CONTINUE
      DO 39 I=N,1,-1
        L=I+1
        G=W(I)
        IF (I.LT.N) THEN
          DO 33 J=L,N
            A(I,J)=0.0
33        CONTINUE
        ENDIF
        IF (G.NE.0.0) THEN
          G=1.0/G
          IF (I.NE.N) THEN
            DO 36 J=L,N
              S=0.0
              DO 34 K=L,M
                S=S+A(K,I)*A(K,J)
34            CONTINUE
              F=(S/A(I,I))*G
              DO 35 K=I,M
                A(K,J)=A(K,J)+F*A(K,I)
35            CONTINUE
36          CONTINUE
          ENDIF
          DO 37 J=I,M
            A(J,I)=A(J,I)*G
37        CONTINUE
        ELSE
          DO 38 J= I,M
            A(J,I)=0.0
38        CONTINUE
        ENDIF
        A(I,I)=A(I,I)+1.0
39    CONTINUE
      DO 49 K=N,1,-1
        DO 48 ITS=1,30
          DO 41 L=K,1,-1
            NM=L-1
            IF ((ABS(RV1(L))+ANORM).EQ.ANORM)  GO TO 2
            IF ((ABS(W(NM))+ANORM).EQ.ANORM)  GO TO 1
41        CONTINUE
1         C=0.0
          S=1.0
          DO 43 I=L,K
            F=S*RV1(I)
            IF ((ABS(F)+ANORM).NE.ANORM) THEN
              G=W(I)
              H=SQRT(F*F+G*G)
              W(I)=H
              H=1.0/H
              C= (G*H)
              S=-(F*H)
              DO 42 J=1,M
                Y=A(J,NM)
                Z=A(J,I)
                A(J,NM)=(Y*C)+(Z*S)
                A(J,I)=-(Y*S)+(Z*C)
42            CONTINUE
            ENDIF
43        CONTINUE
2         Z=W(K)
          IF (L.EQ.K) THEN
            IF (Z.LT.0.0) THEN
              W(K)=-Z
              DO 44 J=1,N
                V(J,K)=-V(J,K)
44            CONTINUE
            ENDIF
            GO TO 3
          ENDIF
          IF (ITS.EQ.30) PAUSE 'No convergence in 30 iterations'
          X=W(L)
          NM=K-1
          Y=W(NM)
          G=RV1(NM)
          H=RV1(K)
          F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
          G=SQRT(F*F+1.0)
          F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X
          C=1.0
          S=1.0
          DO 47 J=L,NM
            I=J+1
            G=RV1(I)
            Y=W(I)
            H=S*G
            G=C*G
            Z=SQRT(F*F+H*H)
            RV1(J)=Z
            C=F/Z
            S=H/Z
            F= (X*C)+(G*S)
            G=-(X*S)+(G*C)
            H=Y*S
            Y=Y*C
            DO 45 NM=1,N
              X=V(NM,J)
              Z=V(NM,I)
              V(NM,J)= (X*C)+(Z*S)
              V(NM,I)=-(X*S)+(Z*C)
45          CONTINUE
            Z=SQRT(F*F+H*H)
            W(J)=Z
            IF (Z.NE.0.0) THEN
              Z=1.0/Z
              C=F*Z
              S=H*Z
            ENDIF
            F= (C*G)+(S*Y)
            X=-(S*G)+(C*Y)
            DO 46 NM=1,M
              Y=A(NM,J)
              Z=A(NM,I)
              A(NM,J)= (Y*C)+(Z*S)
              A(NM,I)=-(Y*S)+(Z*C)
46          CONTINUE
47        CONTINUE
          RV1(L)=0.0
          RV1(K)=F
          W(K)=X
48      CONTINUE
3       CONTINUE
49    CONTINUE
      RETURN
      END
      SUBROUTINE SVDFIT(X,Y,SIG,NDATA,A,MA,U,V,W,MP,NP,CHISQ,FUNCS)
      PARAMETER(NMAX=1000,MMAX=50,TOL=1.E-5)
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),V(NP,NP),
     *    U(MP,NP),W(NP),B(NMAX),AFUNC(MMAX)
      DO 12 I=1,NDATA
        CALL FUNCS(X(I),AFUNC,MA)
        TMP=1./SIG(I)
        DO 11 J=1,MA
          U(I,J)=AFUNC(J)*TMP
11      CONTINUE
        B(I)=Y(I)*TMP
12    CONTINUE
      CALL SVDCMP(U,NDATA,MA,MP,NP,W,V)
      WMAX=0.
      DO 13 J=1,MA
        IF(W(J).GT.WMAX)WMAX=W(J)
13    CONTINUE
      THRESH=TOL*WMAX
      DO 14 J=1,MA
        IF(W(J).LT.THRESH)W(J)=0.
14    CONTINUE
      CALL SVBKSB(U,W,V,NDATA,MA,MP,NP,B,A)
      CHISQ=0.
      DO 16 I=1,NDATA
        CALL FUNCS(X(I),AFUNC,MA)
        SUM=0.
        DO 15 J=1,MA
          SUM=SUM+A(J)*AFUNC(J)
15      CONTINUE
        CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2
16    CONTINUE
      RETURN
      END
      SUBROUTINE SVDVAR(V,MA,NP,W,CVM,NCVM)
      PARAMETER (MMAX=20)
      DIMENSION V(NP,NP),W(NP),CVM(NCVM,NCVM),WTI(MMAX)
      DO 11 I=1,MA
        WTI(I)=0.
        IF(W(I).NE.0.) WTI(I)=1./(W(I)*W(I))
11    CONTINUE
      DO 14 I=1,MA
        DO 13 J=1,I
          SUM=0.
          DO 12 K=1,MA
            SUM=SUM+V(I,K)*V(J,K)*WTI(K)
12        CONTINUE
          CVM(I,J)=SUM
          CVM(J,I)=SUM
13      CONTINUE
14    CONTINUE
      RETURN
      END
