	PROGRAM D1R4
C	Driver for routine CALDAT
	CHARACTER NAME(12)*10
C	Check whether CALDAT properly undoes the operation of JULDAY
	DATA NAME/'January','February','March','April','May',
     *		'June','July','August','September','October',
     *		'November','December'/
	OPEN(5,FILE='DATES.DAT',STATUS='OLD')
	READ(5,*)
	READ(5,*) N
	WRITE(*,'(/1X,A,T40,A)') 'Original Date:','Reconstructed Date:'
	WRITE(*,'(1X,A,T12,A,T17,A,T25,A,T40,A,T50,A,T55,A/)') 
     *		'Month','Day','Year','Julian Day','Month','Day','Year'
	DO 11 I=1,N
		READ(5,'(I2,I3,I5)') IM,ID,IY
		IYCOPY=IY
		J=JULDAY(IM,ID,IYCOPY)
		CALL CALDAT(J,IMM,IDD,IYY)
		WRITE(*,'(1X,A,I3,I6,4X,I9,6X,A,I3,I6)') NAME(IM),ID,
     *			IY,J,NAME(IMM),IDD,IYY
11	CONTINUE
	END
	PROGRAM D6R29
C	Driver for routine CEL
	COMMON A,B,P,AKC
	EXTERNAL FUNC
	PARAMETER(PIO2=1.5707963)
	WRITE(*,*) 'Complete Elliptic Integral'
	WRITE(*,'(1X,T7,A,T18,A,T28,A,T38,A,T47,A,T54,A)')
     *		'kc','p','a','b','CEL','Integral'
	IDUM=-55
	AGO=0.0
	ASTOP=PIO2
	DO 11 I=1,20
		AKC=0.1+RAN3(IDUM)	
		A=10.0*RAN3(IDUM)
		B=10.0*RAN3(IDUM)
		P=0.1+RAN3(IDUM)	
		CALL QSIMP(FUNC,AGO,ASTOP,S)
		WRITE(*,'(1X,6F10.6)')
     *			AKC,P,A,B,CEL(AKC,P,A,B),S
11	CONTINUE
	END

	FUNCTION FUNC(PHI)
	COMMON A,B,P,AKC
	CS=COS(PHI)
	CSQ=CS*CS
	SSQ=1.0-CSQ
	FUNC=(A*CSQ+B*SSQ)/(CSQ+P*SSQ)/SQRT(CSQ+AKC*AKC*SSQ)
	END
	PROGRAM D5R7
C	Driver for routine CHDER
	PARAMETER(NVAL=40, PIO2=1.5707963)
	EXTERNAL FUNC,FDER
	DIMENSION C(NVAL),CDER(NVAL)
	A=-PIO2
	B=PIO2
	CALL CHEBFT(A,B,C,NVAL,FUNC)
C	Test derivative
10	WRITE(*,*) 'How many terms in Chebyshev evaluation?'
	WRITE(*,'(1X,A,I2,A)') 'Enter n between 6 and ',NVAL,
     *				'. Enter n=0 to end.'
	READ(*,*) MVAL
	IF ((MVAL.LE.0).OR.(MVAL.GT.NVAL)) GOTO 20
	CALL CHDER(A,B,C,CDER,MVAL)
	WRITE(*,'(1X,T10,A,T19,A,T28,A)') 'X','Actual','Cheby. Deriv.'
	DO 11 I=-8,8,1
		X=I*PIO2/10.0
		WRITE(*,'(1X,3F12.6)') X,FDER(X),CHEBEV(A,B,CDER,MVAL,X)
11	CONTINUE
	GOTO 10
20	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END

	FUNCTION FDER(X)
C	Derivative of FUNC
	FDER=4.0*X*((X**2)-1.0)*SIN(X)+(X**2)*(X**2-2.0)*COS(X)
	END
	PROGRAM D5R5
C	Driver for routine CHEBEV
	PARAMETER(NVAL=40, PIO2=1.5707963)
	EXTERNAL FUNC
	DIMENSION C(NVAL)
	A=-PIO2
	B=PIO2
	CALL CHEBFT(A,B,C,NVAL,FUNC)
C	Test Chebyshev evaluation routine
10	WRITE(*,*) 'How many terms in Chebyshev evaluation?'
	WRITE(*,'(1X,A,I2,A)') 'Enter n between 6 and ',NVAL,
     *				'. Enter n=0 to end.'
	READ(*,*) MVAL
	IF ((MVAL.LE.0).OR.(MVAL.GT.NVAL)) GOTO 20
	WRITE(*,'(1X,T10,A,T19,A,T28,A)') 'X','Actual','Chebyshev fit'
	DO 11 I=-8,8,1
		X=I*PIO2/10.0
		WRITE(*,'(1X,3F12.6)') X,FUNC(X),CHEBEV(A,B,C,MVAL,X)
11	CONTINUE
	GOTO 10
20	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END
	PROGRAM D5R4
C	Driver for routine CHEBFT
	PARAMETER(NVAL=40, PIO2=1.5707963, EPS=1E-6)
	EXTERNAL FUNC
	DIMENSION C(NVAL)
	A=-PIO2
	B=PIO2
	CALL CHEBFT(A,B,C,NVAL,FUNC)
C	Test result
10	WRITE(*,*) 'How many terms in Chebyshev evaluation?'
	WRITE(*,'(1X,A,I2,A)') 'Enter n between 6 and ',NVAL,
     *				'. Enter n=0 to end.'
	READ(*,*) MVAL
	IF ((MVAL.LE.0).OR.(MVAL.GT.NVAL)) GOTO 20
	WRITE(*,'(1X,T10,A,T19,A,T28,A)') 'X','Actual','Chebyshev fit'
	DO 12 I=-8,8,1
		X=I*PIO2/10.0
		Y=(X-0.5*(B+A))/(0.5*(B-A))
C	Evaluate Chebyshev polynomial without using routine CHEBEV
		T0=1.0
		T1=Y
		F=C(2)*T1+C(1)*0.5
		DO 11 J=3,MVAL
			DUM=T1
			T1=2.0*Y*T1-T0
			T0=DUM
			TERM=C(J)*T1
			F=F+TERM
11		CONTINUE
		WRITE(*,'(1X,3F12.6)') X,FUNC(X),F
12	CONTINUE
	GOTO 10
20	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END
	PROGRAM D5R8
C	Driver for routine CHEBPC
	PARAMETER(NVAL=40, PIO2=1.5707963)
	EXTERNAL FUNC
	DIMENSION C(NVAL),D(NVAL)
	A=-PIO2
	B=PIO2
	CALL CHEBFT(A,B,C,NVAL,FUNC)
10	WRITE(*,*) 'How many terms in Chebyshev evaluation?'
	WRITE(*,'(1X,A,I2,A)') 'Enter n between 6 and ',NVAL,
     *				'. Enter n=0 to end.'
	READ(*,*) MVAL
	IF ((MVAL.LE.0).OR.(MVAL.GT.NVAL)) GOTO 20
	CALL CHEBPC(C,D,MVAL)
C	Test polynomial
	WRITE(*,'(1X,T10,A,T19,A,T29,A)') 'X','Actual','Polynomial'
	DO 12 I=-8,8,1
		X=I*PIO2/10.0
		Y=(X-(0.5*(B+A)))/(0.5*(B-A))
		POLY=D(MVAL)
		DO 11 J=MVAL-1,1,-1
			POLY=POLY*Y+D(J)
11		CONTINUE
		WRITE(*,'(1X,3F12.6)') X,FUNC(X),POLY
12	CONTINUE
	GOTO 10
20	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END
	PROGRAM D5R6
C	Driver for routine CHINT
	PARAMETER(NVAL=40, PIO2=1.5707963)
	EXTERNAL FUNC,FINT
	DIMENSION C(NVAL),CINT(NVAL)
	A=-PIO2
	B=PIO2
	CALL CHEBFT(A,B,C,NVAL,FUNC)
C	Test integral
10	WRITE(*,*) 'How many terms in Chebyshev evaluation?'
	WRITE(*,'(1X,A,I2,A)') 'Enter n between 6 and ',NVAL,
     *				'. Enter n=0 to end.'
	READ(*,*) MVAL
	IF ((MVAL.LE.0).OR.(MVAL.GT.NVAL)) GOTO 20
	CALL CHINT(A,B,C,CINT,MVAL)
	WRITE(*,'(1X,T10,A,T19,A,T29,A)') 'X','Actual','Cheby. Integ.'
	DO 11 I=-8,8,1
		X=I*PIO2/10.0
		WRITE(*,'(1X,3F12.6)') X,FINT(X)-FINT(-PIO2)
     *				,CHEBEV(A,B,CINT,MVAL,X)
11	CONTINUE
	GOTO 10
20	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END

	FUNCTION FINT(X)
C	Integral of FUNC
	FINT=4.0*X*((X**2)-7.0)*SIN(X)
     *		-((X**4)-14.0*(X**2)+28.0)*COS(X)
	END
	PROGRAM D13R8
C	Driver for routine CHSONE
	PARAMETER(NBINS=10,NPTS=2000)
	DIMENSION BINS(NBINS),EBINS(NBINS)
	IDUM=-15
	DO 11 J=1,NBINS
		BINS(J)=0.0
11	CONTINUE
	DO 12 I=1,NPTS
		X=EXPDEV(IDUM)
		IBIN=X*NBINS/3.0+1
		IF(IBIN.LE.NBINS) BINS(IBIN)=BINS(IBIN)+1.0
12	CONTINUE
	DO 13 I=1,NBINS
		EBINS(I)=3.0*NPTS/NBINS*EXP(-3.0*(I-0.5)/NBINS)
13	CONTINUE
	CALL CHSONE(BINS,EBINS,NBINS,-1,DF,CHSQ,PROB)
	WRITE(*,'(1X,T10,A,T25,A)') 'Expected','Observed'
	DO 14 I=1,NBINS
		WRITE(*,'(1X,2F15.2)') EBINS(I),BINS(I)
14	CONTINUE
	WRITE(*,'(/1X,T9,A,E12.4)') 'Chi-squared:',CHSQ
	WRITE(*,'(1X,T9,A,E12.4)') 'Probability:',PROB
	END
	PROGRAM D13R9
C	Driver for routine CHSTWO
	PARAMETER(NBINS=10,NPTS=2000)
	DIMENSION BINS1(NBINS),BINS2(NBINS)
	IDUM=-17
	DO 11 J=1,NBINS
		BINS1(J)=0.0
		BINS2(J)=0.0
11	CONTINUE
	DO 12 I=1,NPTS
		X=EXPDEV(IDUM)
		IBIN=X*NBINS/3.0+1
		IF(IBIN.LE.NBINS) BINS1(IBIN)=BINS1(IBIN)+1.0
		X=EXPDEV(IDUM)
		IBIN=X*NBINS/3.0+1
		IF(IBIN.LE.NBINS) BINS2(IBIN)=BINS2(IBIN)+1.0
12	CONTINUE
	CALL CHSTWO(BINS1,BINS2,NBINS,-1,DF,CHSQ,PROB)
	WRITE(*,'(1X,T10,A,T25,A)') 'Dataset 1','Dataset 2'
	DO 13 I=1,NBINS
		WRITE(*,'(1X,2F15.2)') BINS1(I),BINS2(I)
13	CONTINUE
	WRITE(*,'(/1X,T10,A,E12.4)') 'Chi-squared:',CHSQ
	WRITE(*,'(1X,T10,A,E12.4)') 'Probability:',PROB
	END
	PROGRAM D13R13
C	Driver for routine CNTAB1
C	Contingency table in file TABLE.DAT
	PARAMETER(NDAT=9,NMON=12)
	DIMENSION NMBR(NDAT,NMON)
	CHARACTER FATE(NDAT)*15,MON(NMON)*5,TEXT*64
	OPEN(5,FILE='TABLE.DAT',STATUS='OLD')
	READ(5,*)
	READ(5,'(A)') TEXT
	READ(5,'(15X,12A5/)') (MON(I),I=1,12)
	DO 11 I=1,NDAT
		READ(5,'(A15,12I5)') FATE(I),(NMBR(I,J),J=1,12)
11	CONTINUE
	CLOSE(5)
	WRITE(*,'(/1X,A/)') TEXT
	WRITE(*,'(1X,15X,12A5)') (MON(I),I=1,12)
	DO 12 I=1,NDAT
		WRITE(*,'(1X,A,12I5)') FATE(I),(NMBR(I,J),J=1,12)
12	CONTINUE
	CALL CNTAB1(NMBR,NDAT,NMON,CHISQ,DF,PROB,CRAMRV,CCC)
	WRITE(*,'(/1X,A,T20,F20.2)') 'Chi-squared',CHISQ
	WRITE(*,'(1X,A,T20,F20.2)') 'Degrees of Freedom',DF
	WRITE(*,'(1X,A,T20,F20.4)') 'Probability',PROB
	WRITE(*,'(1X,A,T20,F20.4)') 'Cramer-V',CRAMRV
	WRITE(*,'(1X,A,T20,F20.4)') 'Contingency Coeff.',CCC
	END
	PROGRAM D13R14
C	Driver for routine CNTAB2
C	Contingency table in file TABLE.DAT
	PARAMETER(NI=9,NMON=12)
	DIMENSION NMBR(NI,NMON)
	CHARACTER FATE(NI)*15,MON(NMON)*5,TEXT*64
	OPEN(5,FILE='TABLE.DAT',STATUS='OLD')
	READ(5,*)
	READ(5,'(A)') TEXT
	READ(5,'(15X,12A5/)') (MON(I),I=1,12)
	DO 11 I=1,NI
		READ(5,'(A15,12I5)') FATE(I),(NMBR(I,J),J=1,12)
11	CONTINUE
	CLOSE(5)
	WRITE(*,'(/1X,A/)') TEXT
	WRITE(*,'(1X,15X,12A5)') (MON(I),I=1,12)
	DO 12 I=1,NI
		WRITE(*,'(1X,A,12I5)') FATE(I),(NMBR(I,J),J=1,12)
12	CONTINUE
	CALL CNTAB2(NMBR,NI,NMON,H,HX,HY,HYGX,HXGY,UYGX,UXGY,UXY)
	WRITE(*,'(/1X,A,T30,F10.4)') 'Entropy of Table',H
	WRITE(*,'(1X,A,T30,F10.4)') 'Entropy of x-distribution',HX
	WRITE(*,'(1X,A,T30,F10.4)') 'Entropy of y-distribution',HY
	WRITE(*,'(1X,A,T30,F10.4)') 'Entropy of y given x',HYGX
	WRITE(*,'(1X,A,T30,F10.4)') 'Entropy of x given y',HXGY
	WRITE(*,'(1X,A,T30,F10.4)') 'Dependency of y on x',UYGX
	WRITE(*,'(1X,A,T30,F10.4)') 'Dependency of x on y',UXGY
	WRITE(*,'(1X,A,T30,F10.4/)') 'Symmetrical dependency',UXY
	END
	PROGRAM D12R6
C	Driver for routine CONVLV
	PARAMETER(N=16,N2P2=34,M=9,PI=3.14159265)
	DIMENSION DATA(N),RESPNS(M),RESP(N),ANS(N2P2)
	DO 11 I=1,N
		DATA(I)=0.0
		IF ((I.GE.(N/2-N/8)).AND.(I.LE.(N/2+N/8))) DATA(I)=1.0
11	CONTINUE
	DO 12 I=1,M
		RESPNS(I)=0.0
		IF (I.GT.2 .AND. I.LT.7) RESPNS(I)=1.0
		RESP(I)=RESPNS(I)
12	CONTINUE
	ISIGN=1
	CALL CONVLV(DATA,N,RESP,M,ISIGN,ANS)
C	Compare with a direct convolution
	WRITE(*,'(/1X,T4,A,T13,A,T24,A)') 'I','CONVLV','Expected'
	DO 14 I=1,N
		CMP=0.0
		DO 13 J=1,M/2
			CMP=CMP+DATA(MOD(I-J-1+N,N)+1)*RESPNS(J+1)
			CMP=CMP+DATA(MOD(I+J-1,N)+1)*RESPNS(M-J+1)
13		CONTINUE
		CMP=CMP+DATA(I)*RESPNS(1)
		WRITE(*,'(1X,I3,3X,2F12.6)') I,ANS(I),CMP
14	CONTINUE
	END
	PROGRAM D12R7
C	Driver for routine CORREL
	PARAMETER(N=64,N2=128,PI=3.1415927)
	DIMENSION DATA1(N),DATA2(N),ANS(N2)
	DO 11 I=1,N
		DATA1(I)=0.0
		IF ((I.GT.(N/2-N/8)).AND.(I.LT.(N/2+N/8))) DATA1(I)=1.0
		DATA2(I)=DATA1(I)
11	CONTINUE
	CALL CORREL(DATA1,DATA2,N,ANS)
C	Calculate directly
	WRITE(*,'(/1X,T4,A,T13,A,T25,A/)') 'n','CORREL','Direct Calc.'
	DO 13 I=0,16
		CMP=0.0
		DO 12 J=1,N
			CMP=CMP+DATA1(MOD(I+J-1,N)+1)*DATA2(J)
12		CONTINUE
		WRITE(*,'(1X,I3,3X,F12.6,F15.6)') I,ANS(I+1),CMP
13	CONTINUE
	END
	PROGRAM D12R5
C	Driver for routine COSFT
	PARAMETER(EPS=1.0E-3,NP=16,WIDTH=30.0,PI=3.14159)
	DIMENSION DATA(NP),SIZE(NP)
1	WRITE(*,'(1X,A,I2,A)') 'Period of cosine in channels (2-',NP,')'
	READ(*,*) PER
	IF (PER.LE.0.) STOP
	DO 11 I=1,NP
		DATA(I)=COS(2.0*PI*(I-1)/PER)
11	CONTINUE
	CALL COSFT(DATA,NP,+1)
	BIG=-1.0E10
	SMALL=1.0E10
	DO 12 I=1,NP
		IF (DATA(I).LT.SMALL) SMALL=DATA(I)
		IF (DATA(I).GT.BIG) BIG=DATA(I)
12	CONTINUE
	SCAL=WIDTH/(BIG-SMALL)
	DO 13 I=1,NP
		NLIM=SCAL*(DATA(I)-SMALL)+EPS
		WRITE(*,'(1X,I2,F6.2,1X,60A1)') I,DATA(I),('*',J=1,NLIM+1)
13	CONTINUE
	WRITE(*,*) 'press continue ...'
	READ(*,*)
	CALL COSFT(DATA,NP,-1)
	BIG=-1.0E10
	SMALL=1.0E10
	DO 14 I=1,NP
		IF(DATA(I).LT.SMALL) SMALL=DATA(I)
		IF(DATA(I).GT.BIG) BIG=DATA(I)
14	CONTINUE
	SCAL=WIDTH/(BIG-SMALL)
	DO 15 I=1,NP
		NLIM=SCAL*(DATA(I)-SMALL)+EPS
		WRITE(*,'(1X,I4,1X,60A1)') I,('*',J=1,NLIM+1)
15	CONTINUE
	GOTO 1
	END
	PROGRAM D14R3
C	Driver for routine COVSRT
	PARAMETER(MA=10,MFIT=5)
	DIMENSION COVAR(MA,MA),LISTA(MFIT)
	DO 12 I=1,MA
		DO 11 J=1,MA
			COVAR(I,J)=0.0
			IF (I.LE.5 .AND. J.LE.5) THEN
				COVAR(I,J)=I+J-1
			ENDIF
11		CONTINUE
12	CONTINUE
	WRITE(*,'(//2X,A)') 'Original matrix'
	DO 13 I=1,MA
		WRITE(*,'(1X,10F4.1)') (COVAR(I,J),J=1,MA)
13	CONTINUE
	WRITE(*,*) ' press RETURN to continue...'
	READ(*,*)
C	Test 1 - spread by 2
	WRITE(*,'(/2X,A)') 'Test #1 - Spread by two'
	DO 14 I=1,MFIT
		LISTA(I)=2*I
14	CONTINUE
	CALL COVSRT(COVAR,MA,MA,LISTA,MFIT)
	DO 15 I=1,MA
		WRITE(*,'(1X,10F4.1)') (COVAR(I,J),J=1,MA)
15	CONTINUE
	WRITE(*,*) ' press RETURN to continue...'
	READ(*,*)
C	Test 2 - reverse
	WRITE(*,'(/2X,A)') 'Test #2 - Reverse'
	DO 17 I=1,MA
		DO 16 J=1,MA
			COVAR(I,J)=0.0
			IF (I.LE.5 .AND. J.LE.5) THEN
				COVAR(I,J)=I+J-1
			ENDIF
16		CONTINUE
17	CONTINUE
	DO 18 I=1,MFIT
		LISTA(I)=MFIT+1-I		
18	CONTINUE
	CALL COVSRT(COVAR,MA,MA,LISTA,MFIT)
	DO 19 I=1,MA
		WRITE(*,'(1X,10F4.1)') (COVAR(I,J),J=1,MA)
19	CONTINUE
	WRITE(*,*) ' press RETURN to continue...'
	READ(*,*)
C	Test 3 - spread and reverse
	WRITE(*,'(/2X,A)') 'Test #3 - Spread and reverse'
	DO 22 I=1,MA
		DO 21 J=1,MA
			COVAR(I,J)=0.0
			IF (I.LE.5 .AND. J.LE.5) THEN
				COVAR(I,J)=I+J-1
			ENDIF
21		CONTINUE
22	CONTINUE
	DO 23 I=1,MFIT
		LISTA(I)=MA+2-2*I
23	CONTINUE
	CALL COVSRT(COVAR,MA,MA,LISTA,MFIT)
	DO 24 I=1,MA
		WRITE(*,'(1X,10F4.1)') (COVAR(I,J),J=1,MA)
24	CONTINUE
	END
	PROGRAM D13R17
C	Driver for routine CRANK
	PARAMETER(NDAT=20,NMON=12)
	DIMENSION DATA(NDAT),RAYS(NDAT,NMON)
	DIMENSION ORDER(NDAT),AVE(NDAT),ZLAT(NDAT)
	CHARACTER CITY(NDAT)*15,MON(NMON)*4,TEXT*64
	OPEN(5,FILE='TABLE2.DAT',STATUS='OLD')
	READ(5,*)
	READ(5,'(A)') TEXT
	READ(5,'(15X,12A4/)') (MON(I),I=1,12)
	DO 11 I=1,NDAT
		READ(5,'(A15,12F4.0,F6.0,F6.1)')
     *			CITY(I),(RAYS(I,J),J=1,12),AVE(I),ZLAT(I)
11	CONTINUE
	CLOSE(5)
	WRITE(*,'(1X,A)') TEXT
	WRITE(*,'(1X,15X,12A4)') (MON(I),I=1,12)
	DO 12 I=1,NDAT
		WRITE(*,'(1X,A,12I4)')
     *			CITY(I),(NINT(RAYS(I,J)),J=1,12)
12	CONTINUE
C	Replace solar flux in each column by rank order
1	WRITE(*,'(1X,A)') 'Number of month (1-12):'
	READ(*,*) MONTH
	IF(MONTH.EQ.0)STOP
	MONTH=MOD(MONTH+5,12)+1
	DO 13 I=1,NDAT
		DATA(I)=RAYS(I,MONTH)
		ORDER(I)=I
13	CONTINUE
	CALL SORT2(NDAT,DATA,ORDER)
	CALL CRANK(NDAT,DATA,S)
	WRITE(*,*) 'Month of',MON(MONTH)
	WRITE(*,'(1X,A,T27,A,T36,A,T49,A)') 'City','Rank',
     *			'Solar Flux','Latitude'
	DO 14 I=1,NDAT
		NN=NINT(ORDER(I))
		WRITE(*,'(1X,A,T19,3F12.1)') CITY(NN),DATA(I),
     *			RAYS(NN,MONTH),ZLAT(NN)
14	CONTINUE
	GOTO 1
	END
