PROGRAM RateErrors
!	Error Propagation in Resonant Nuclear Reaction Rates
!	Version 2.1 DC (includes Direct Capture)
!	January 1999
!	William J. Thompson, University of North Carolina
!	wmjthompson@msn.com	
!
! - changed I/O so that two output files are produced, one with the "detail"
!   option (error contribution from each resonance andeach term), the other
!   one lists T9 vs. reaction rate on a standard grid. The input file has
!   also been changed.
! 
! - the masses of projectile and target have been changed from integer to
!   real numbers so that precise masses can be used 
! 
! - the projectile mass and charge can now be chosen more generally [the
!   original version of the code was specific to (p,gamma)]
! 
!   CI, October 2003

! Subroutine names: ErrorCalc,ErrorIn,ErrorOut,
!	NonResRates,SaveForPlot,SortOutput
!	ErrorCalc does the computation
!	ErrorIn and ErrorOut are for test Input and Output
!	NonResRate calculates nonresonant reaction rate
!	SaveForPlot saves output for plotting
!	SortOutput sorts output by descending order

double precision T9,T9val(51),kT,sdEn(51),En(51),sdOG(51), &
	SumSquares,OG(51),Boltz(51),Qvalue,sd(51),var(51),Atarget,Aproj, &
	GP(51),GG(51),PartialRate(51),TotalRate,SDPartialRate,fracSD(51), &
	sdOGGp(51),sdOGGg(51),sdOGBoltz(51),sdBoltz(51),JJ(51),jp,jT, &
	rateFac,OGdata(51),alphaBB,alphaTB,fgBB,fgTB,sdOGSpec(51)
integer resonances,numT9,Ztarget,Zproj,Lvalue(51),Detail
character*80 FileHeading,SaveFileName
common /Reals/T9,T9val,kT,OG,sdOG,En,sdEn,Qvalue, &
	Boltz,sd,var,sumSquares,GP,GG,PartialRate,Atarget,Aproj, &
	TotalRate,fracSD,sdOGGp,sdOGGg,sdOGBoltz,sdBoltz, &
	JJ,jp,jT,rateFac,OGdata,SDPartialRate,fracErrSpec,sdOGSpec
common /Integers/Ztarget,Zproj,Lvalue,Detail
common /Characters/FileHeading,SaveFileName

open(4,FILE="RateErrors_rates.out",STATUS='unknown')

write(4,*) " T9    ","  Lower L  ","  Rec. Rate   "," Upper L  ","  C.L. (68%,>16%!)"

numT9 = 51

T9val(1)=0.010
T9val(2)=0.011
T9val(3)=0.012
T9val(4)=0.013
T9val(5)=0.014
T9val(6)=0.015
T9val(7)=0.016
T9val(8)=0.018
T9val(9)=0.020
T9val(10)=0.025
T9val(11)=0.030
T9val(12)=0.040
T9val(13)=0.050
T9val(14)=0.060
T9val(15)=0.070
T9val(16)=0.080
T9val(17)=0.09
T9val(18)=0.10
T9val(19)=0.11
T9val(20)=0.12
T9val(21)=0.13
T9val(22)=0.14
T9val(23)=0.15
T9val(24)=0.16
T9val(25)=0.18
T9val(26)=0.20
T9val(27)=0.25
T9val(28)=0.30
T9val(29)=0.35
T9val(30)=0.40
T9val(31)=0.45
T9val(32)=0.50
T9val(33)=0.60
T9val(34)=0.70
T9val(35)=0.80
T9val(36)=0.90
T9val(37)=1.0
T9val(38)=1.25
T9val(39)=1.5
T9val(40)=1.75
T9val(41)=2.0
T9val(42)=2.5
T9val(43)=3.0
T9val(44)=3.5
T9val(45)=4.0
T9val(46)=5.0
T9val(47)=6.0
T9val(48)=7.0
T9val(49)=8.0
T9val(50)=9.0
T9val(51)=10.0

write (*,*) &
"RateErrors: Error Propagation in Resonant Nuclear Reaction Rates (v. 2.1)"
write (*,*) " (Includes Direct Capture rates with errors)" 
!resonances = 1
!do while (resonances > 0 )
!		Input quantities and standard deviations
		call ErrorIn(resonances)
!		Do calculations for numT9 temperatures
		do nT9 = 1, numT9
!			Error calculation
			call ErrorCalc(resonances,nT9)
!			Error output
			call ErrorOut(resonances)
!			Saving file choice
			call SaveFile(resonances,nT9)
		end do
!end do
END PROGRAM RateErrors
!*********************************************************************
SUBROUTINE ErrorCalc(resonances,nT9)
!	Calculate the errors in reaction rates
double precision T9,T9val(51),kT,sdEn(51),En(51),sdOG(51), &
	SumSquares,OG(51),Boltz(51),Qvalue,sd(51),var(51), &
	GP(51),GG(51),PartialRate(51),TotalRate,SDPartialRate,fracSD(51), &
	sdOGGp(51),sdOGGg(51),sdOGBoltz(51),sdBoltz(51),JJ(51),jp,jT, &
	rateFac,OGdata(51),alphaBB,alphaTB,fgBB,fgTB,sdOGSpec(51), &
	gammaFac,BB,fracGPerr,fracGGerr, &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate
double precision mu,fZtarget,fZproj,derOG,omegaFac,Atarget,Aproj, &
	BoltzShift,fracErrGmw,fracErrGG,fracErrBoltz,sdDkT,fracSerr
integer resonances,nres,nT9,Ztarget,Zproj,Lvalue(51),Detail		
common /Reals/T9,T9val,kT,OG,sdOG,En,sdEn,Qvalue, &
	Boltz,sd,var,sumSquares,GP,GG,PartialRate,Atarget,Aproj, &
	TotalRate,fracSD,sdOGGp,sdOGGg,sdOGBoltz,sdBoltz, &
	JJ,jp,jT,rateFac,OGdata,SDPartialRate,fracErrSpec,sdOGSpec, &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate
common /Integers/Ztarget,Zproj,Lvalue,Detail

mu = Atarget*Aproj/(Atarget+Aproj)
T9 = T9val(nT9)
kT = 86.170d00*T9
!	Overall rate factor if energies are in keV
rateFac = 1.540d05/(mu*T9)**1.5
!	Compute standard errors in all resonances
do nres = 1,resonances
!	Boltzmann factor; eq.(2.2)
	Boltz(nres) = 	exp(-En(nres)/kT)
	!	Boltzmann energy shift in eq.(3.5)
	BoltzShift = sdEn(nres)**2/(2.0d00*kT)
!	Standard deviation in Boltzmann factor; use eq.(3.8)
!	First check whether lower bound on integral is O.K.
	if ( En(nres) < sdEn(nres)*(4.0d00+sdEn(nres)/kT) ) then
		write(*,"(/a,i2,a,e9.3)") "!! ErrorCalc: For resonance ",nres," at T9 = ",T9
		write(*,"(a)") &
			"   error/energy is too large for reliable Boltzmann error."
	end if
!	Fractional error from eq.(3.8)
	sdDkT = (sdEn(nres)/kT)**2
	fracErrBoltz = exp(0.5d00*sdDkT)*sqrt(exp(sdDkT)-1.0d00)
	sdBoltz(nres) = rateFac*Boltz(nres)*fracErrBoltz
!	Error from omega-gamma and energy uncertainties
	if (OGdata(nres) /= 0.0d00 ) then
		OG(nres) = OGdata(nres)
!		Errors from omega-gamma and energy are independent; use eq.(4.2)
		sdOGBoltz(nres) = rateFac*sdOG(nres)*Boltz(nres)
		sdBoltz(nres) = OG(nres)*sdBoltz(nres)
		sd(nres) = sqrt(sdOGBoltz(nres)**2+sdBoltz(nres)**2)
!		Reaction rate from eq. (2.2)
		PartialRate(nres) = rateFac*OG(nres)*Boltz(nres)
		fracSD(nres) = sd(nres)/PartialRate(nres)
	else
!		Errors are estimated from energy errors
!		Omega-gamma; eq.(2.3)
		omegaFac = (2.0d00*JJ(nres)+1.0d00)/((2.0d00*jp+1.0d00)*(2.0d00*jT+1.0d00))
		gammaFac = GP(nres)*GG(nres)/(GP(nres)+GG(nres))
		OG(nres) = omegaFac*gammaFac
!		Gamow factor standard deviation
		fZtarget = Ztarget
		fZproj=Zproj
		BB = 31.29d00*sqrt(mu)*fZtarget*fZproj
!		Energy width from eq.(3.13)
		alphaBB = sdEn(nres)/ &
	 		sqrt(1.0d00+0.75d00*(sdEn(nres))**2*BB/((En(nres))**2.5))
!		then repeat at 2b for use in eq.(3.15)
		alphaTB = sdEn(nres)/ &
	 		sqrt(1.0d00+1.50d00*(sdEn(nres))**2*BB/((En(nres))**2.5))
!		Fractional error in Gamow factor; eq.(3.12) then eq.(3.15) 
		fgBB = alphaBB*exp((alphaBB*BB)**2/(8.0d00*En(nres)**3))/sdEn(nres)
		fgTB = alphaTB*exp((alphaTB*2.0*BB)**2/(8.0d00*En(nres)**3))/sdEn(nres)
		fracErrGmw = sqrt(fgTB - fgBB**2)
!		Gamma-g errors from energy error and Q-value
!		Fractional S.D. in gamma-g; eq.(4.8)
		fracErrGG = (2.0d00*float(Lvalue(nres))+1.0d00)*sdEn(nres)/(En(nres)+Qvalue)
!		Fractional S.D. in reaction rate; use ingredients in eq.(4.6)
		fracGPerr = gammaFac*fracErrGmw/GP(nres)
		fracGGerr = gammaFac*fracErrGG/GG(nres)
		fracSerr = gammaFac*fracErrSpec/GP(nres)
		fracSD(nres) = sqrt((fracGPerr-fracErrBoltz+fracGGerr)**2+fracSerr**2)
!		Reaction rate from eq. (2.2)
		PartialRate(nres) = rateFac*OG(nres)*Boltz(nres)	
!		Finally, S.D. in resonance reaction rate
		sd(nres) = fracSD(nres)*PartialRate(nres)
!		with contributions from 4 parts; first gamma-p
		sdOGGp(nres) = PartialRate(nres)*fracGPerr
!		then Boltzmann factor
		sdOGBoltz(nres) = PartialRate(nres)*fracErrBoltz	
!		then gamma-g		
		sdOGGg(nres) = PartialRate(nres)*fracGGerr
!		then spectroscopic factor
		sdOGSpec(nres) = PartialRate(nres)*fracSerr
	end if
end do
!	Direct-capture rate and its error
call NonResRate
!	Sum reaction rates and variances
!	in all resonances; eqs.(2.1) and (4.1)
TotalRate = 0.0d00
SumSquares = 0.0d00
do nres = 1,resonances
!	Total reaction rate is TotalRate from eq.(2.1)
	TotalRate = TotalRate+PartialRate(nres)
!	Total variance in this resonance
	var(nres) = sd(nres)**2
!	Sum of squared errors over the resonances; eq.(4.1)
	SumSquares = SumSquares+var(nres)
end do
!	Include direct-capture rate and its variance
TotalRate = TotalRate+DCRate
SumSquares = SumSquares+(DCRate*fracErrDC)**2
END SUBROUTINE ErrorCalc
!*********************************************************************
SUBROUTINE ErrorIn(resonances)
!	Input quantities, including standard deviations
!	Input is from console or from a named file
double precision T9,T9val(51),kT,sdEn(51),En(51),sdOG(51), &
	SumSquares,OG(51),Boltz(51),Qvalue,sd(51),var(51),Atarget,Aproj, &
	GP(51),GG(51),PartialRate(51),TotalRate,SDPartialRate,fracSD(51), &
	sdOGGp(51),sdOGGg(51),sdOGBoltz(51),sdBoltz(51),JJ(51),jp,jT, &
	rateFac,OGdata(51),alphaBB,alphaTB,fgBB,fgTB,sdOGSpec(51), &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate
integer resonances,numT9,nT9,Detail, &
	Ztarget,Zproj,Lvalue(51)
character*80 FileName,FileHeading,Unused,SaveFileName
common /Reals/T9,T9val,kT,OG,sdOG,En,sdEn,Qvalue, &
	Boltz,sd,var,sumSquares,GP,GG,PartialRate,Atarget,Aproj, &
	TotalRate,fracSD,sdOGGp,sdOGGg,sdOGBoltz,sdBoltz, &
	JJ,jp,jT,rateFac,OGdata,SDPartialRate,fracErrSpec,sdOGSpec, &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate
common /Integers/Ztarget,Zproj,Lvalue,Detail
common /Characters/FileHeading,SaveFileName
 
Detail=1
write(*,*)
write(*,"(a)") "Type input file name:"
read(*,"(a)") FileName
!	Input is from file whose name is in FileName
open(2,FILE=FileName,STATUS='unknown',RECL=80)
!	Heading of input file
read(2,"(a)") FileHeading
read(2,*) Zproj
read(2,*) Ztarget
read(2,*) Aproj
read(2,*) Atarget
read(2,*) jp
read(2,*) jt
read(2,*) Qvalue
read(2,*) fracErrSpec
read(2,*) S,fracErrS
read(2,*) DS,fracErrDS
read(2,*) DDS,fracErrDDS
read(2,*) Ecutoff
read(2,*) resonances
read(2,*) Unused
do nres = 1,resonances
!	Input quantities and their S.D.s for each resonance
	read(2,*) En(nres),sdEn(nres),OGdata(nres),sdOG(nres), &
		JJ(nres),GP(nres),GG(nres),Lvalue(nres)
	if (OGdata(nres) /= 0.0d00 ) then
		OG(nres) = OGdata(nres)
	end if
end do
close(2)
END SUBROUTINE ErrorIn
!*********************************************************************
SUBROUTINE ErrorOut(resonances)
!	Output of errors in reaction rates
double precision T9,T9val(51),kT,sdEn(51),En(51),sdOG(51), &
	SumSquares,OG(51),Boltz(51),Qvalue,sd(51),var(51),Atarget,Aproj, &
	GP(51),GG(51),PartialRate(51),TotalRate,SDPartialRate,fracSD(51), &
	sdOGGp(51),sdOGGg(51),sdOGBoltz(51),sdBoltz(51),JJ(51),jp,jT, &
	rateFac,OGdata(51),alphaBB,alphaTB,fgBB,fgTB,sdOGSpec(51), &
	SDplus,SDminus,Lower,Upper,FracErrSq,Confidence, &
	sdSorted(51),sdInput(51),sdOut(4),sdTotalRate, &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate
integer resonances,Ztarget,Zproj,Lvalue(51),Detail, &
	KeepRes(51),Keep(4)
character*2 big(4)
character*5 name(4),bigKeep(51)
character*80 FileHeading,SaveFileName
common /Reals/T9,T9val,kT,OG,sdOG,En,sdEn,Qvalue, &
	Boltz,sd,var,sumSquares,GP,GG,PartialRate,Atarget,Aproj, &
	TotalRate,fracSD,sdOGGp,sdOGGg,sdOGBoltz,sdBoltz, &
	JJ,jp,jT,rateFac,OGdata,SDPartialRate,fracErrSpec,sdOGSpec, &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate
common /Integers/Ztarget,Zproj,Lvalue,Detail
common /Characters/FileHeading,SaveFileName
data name(1)/"Gp   "/,name(2)/"Boltz"/,name(3)/"Gg   "/,name(4)/"Spec "/

!	Save Direct Capture rate after last resonance
sd(resonances+1) = DCRate*fracErrDC
PartialRate(resonances+1) = DCRate
bigKeep(resonances+1) = "     "
do nres = 1,resonances
!	Find and mark with # largest S.D. contribution in resonance
	big(1) = '  '
	big(2) = '  '
	big(3) = '  '
	big(4) = '  '
	if (OGdata(nres) /= 0.0d00 ) then
!		Omega-gamma was input
		if ( sdOGBoltz(nres) > sdBoltz(nres) ) then
			big(1) = ' #'
			bigKeep(nres) = "Omega"
		else
			big(2) = ' #'
			bigKeep(nres) = "Boltz"
		end if	
	else
!		Omega-gamma is estimated
		sdInput(1) = sdOGGp(nres)
		sdInput(2) = sdOGBoltz(nres)
		sdInput(3) = sdOGGg(nres)
		sdInput(4) = sdOGSpec(nres)
		call SortOutput(sdInput,sdOut,Keep,4)
		big(Keep(1)) = ' #'
		bigKeep(nres) = name(Keep(1))	
	end if
	if ( Detail == 1 ) then
		if (OGdata(nres) /= 0.0d00 ) then
		else
!			Omega-gamma is estimated
		end if
		if ( nres == KeepRes(1) ) then
			big(1) = '##'
		else
			big(1) = '  '
		end if
!		Symmetric standard deviations as in eq.(4.6)
!		Asymmetric standard deviations use eq.(4.13)
		SDplus = PartialRate(nres)*(exp(fracSD(nres)) - 1.0d00)
		SDminus = PartialRate(nres)*(1.0d00 - exp(-fracSD(nres)))
	end if
end do
!	Order resonances and DC by total S.D. in reaction rate
call SortOutput(sd,sdSorted,KeepRes,resonances+1)
!	Totals over all resonances plus direct capture
SDTotalRate = sqrt(SumSquares)
!	Symmetric errors
!	Asymmetric total rates by using eq.(4.13)
Lower = TotalRate*exp(-SDTotalRate/TotalRate)
Upper = TotalRate*exp(SDTotalRate/TotalRate)
SDminus = TotalRate - Lower
SDPlus = Upper - TotalRate
!	Confidence level by eq.(4.17)
FracErrSq = (SDTotalRate/TotalRate)**2
Confidence = 68.3d00 - FracErrSq*(23.1d00 - 2.48d00*FracErrSq)
END SUBROUTINE ErrorOut
!*********************************************************************
SUBROUTINE NonResRate
!		Calculate nonresonant reaction rate
double precision T9,T9val(51),kT,sdEn(51),En(51),sdOG(51), &
	SumSquares,OG(51),Boltz(51),Qvalue,sd(51),var(51),Atarget,Aproj, &
	GP(51),GG(51),PartialRate(51),TotalRate,SDPartialRate,fracSD(51), &
	sdOGGp(51),sdOGGg(51),sdOGBoltz(51),sdBoltz(51),JJ(51),jp,jT, &
	rateFac,OGdata(51),alphaBB,alphaTB,fgBB,fgTB,sdOGSpec(51), &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate,Energy
integer resonances,numT9,Ztarget,Zproj,Lvalue(51),Detail
common /Reals/T9,T9val,kT,OG,sdOG,En,sdEn,Qvalue, &
	Boltz,sd,var,sumSquares,GP,GG,PartialRate,Atarget,Aproj, &
	TotalRate,fracSD,sdOGGp,sdOGGg,sdOGBoltz,sdBoltz, &
	JJ,jp,jT,rateFac,OGdata,SDPartialRate,fracErrSpec,sdOGSpec, &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate
common /Integers/Ztarget,Zproj,Lvalue,Detail
double precision C1,C2,C3,C4,C5,C6,C7,Help,TCUT,R1,R2,R3,R4,mu

if ( S == 0.0 ) then
	DCRate = 0.0
	fracErrDC = 0.0
else
	mu = Atarget*Aproj/(Atarget+Aproj)
	Help = mu*(Ztarget**2)*(Zproj**2)
	CubRtHelp  =  Help**(1.0/3.0)
	TCUT = 19.92*(Ecutoff**1.5)/Sqrt(HELP)
	C1 = 7.8324E09/Sqrt(mu/CubRtHelp)
	C2 = 4.2475*CubRtHelp
	C3 = 9.810E-2/CubRtHelp
	C4 = 0.122*DS*CubRtHelp
	C5 = 8.377E-2*DS
	C6 = 7.442E-3*DDS*CubRtHelp**2
	C7 = 1.299E-2*DDS*CubRtHelp
	R1 = C1/(T9**(2.0/3.0))
	R2 = EXP((-C2/(T9**(1.0/3.0)))-(T9/TCUT)**2.0)
	R3 = S+C3*S*(T9**(1.0/3.0))+C4*(T9**(2.0/3.0))+C5*T9
	R4 = C6*(T9**(4.0/3.0))+C7*(T9**(5.0/3.0))
	DCRate = R1*R2*(R3+R4)
	!	Fractional error in direct-capture reaction rate
	! (kT is in keV, S derivatives are in /MeV)
	Energy = 0.001*kT
	fracErrDC = sqrt(fracErrS**2+(fracErrDS*Energy)**2+ &
		(0.5*fracErrDDS*Energy**2)**2)/ &
		(1.0+DS*Energy/S+0.5*DDS*Energy**2/S)
end if
END SUBROUTINE NonResRate
!*********************************************************************
SUBROUTINE SaveFile(resonances,nT9)
!	Saves output to a file
double precision T9,T9val(51),kT,sdEn(51),En(51),sdOG(51), &
	SumSquares,OG(51),Boltz(51),Qvalue,sd(51),var(51),Atarget,Aproj, &
	GP(51),GG(51),PartialRate(51),TotalRate,SDPartialRate,fracSD(51), &
	sdOGGp(51),sdOGGg(51),sdOGBoltz(51),sdBoltz(51),JJ(51),jp,jT, &
	rateFac,OGdata(51),alphaBB,alphaTB,fgBB,fgTB,sdOGSpec(51), &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate, &
	SDplus,SDminus,Lower,Upper,FracErrSq,Confidence, &
	sdSorted(51),sdInput(51),sdOut(4)
integer resonances,Ztarget,Zproj,Lvalue(51),Detail, &
	KeepRes(51),Keep(4)
character*2 big(4)
character*5 name(4),bigKeep(51)
character*80 FileHeading,SaveFileName
common /Reals/T9,T9val,kT,OG,sdOG,En,sdEn,Qvalue, &
	Boltz,sd,var,sumSquares,GP,GG,PartialRate,Atarget,Aproj, &
	TotalRate,fracSD,sdOGGp,sdOGGg,sdOGBoltz,sdBoltz, &
	JJ,jp,jT,rateFac,OGdata,SDPartialRate,fracErrSpec,sdOGSpec, &
	S,DS,DDS,Ecutoff,fracErrS,fracErrDS,fracErrDDS,fracErrDC,DCRate
common /Integers/Ztarget,Zproj,Lvalue,Detail
common /Characters/FileHeading,SaveFileName
data name(1)/"Gp   "/,name(2)/"Boltz"/,name(3)/"Gg   "/,name(4)/"Spec "/

open(3,FILE="RateErrors_detail.out",STATUS='OLD')
			
write (3,*) &
"RateErrors: Error Propagation in Resonant Nuclear Reaction Rates. Version 2.1 DC"
write(3,"(a)") "Zproj, Ztarget, Aproj, Atarget,  jp,  jT, Q (keV),  frac. S.D. in C2S"
write(3,"(2i8,2f8.4,f5.1,f4.1,e11.4,e9.2)") Zproj,Ztarget,Aproj,Atarget,jp,jT,Qvalue,fracErrSpec
write(3,"(a)") "Direct-capture data:"
write(3,"(a)") &
	" S     &     err,     DS     &     err,      DDS     &     err,   Ecutoff"
write(3,"(3(2e10.3,a),e10.3)") &
	S,fracErrS,"  ",DS,fracErrDS,"  ",DDS,fracErrDDS,"  ",Ecutoff
!write(3,"(/a)") &
write(3,"(a,a,e9.3,a,e9.3,a)") "Reaction rate error estimates for ", &
	"T9(GK) = ",T9," [kT(keV) = ",kT,"]" 
!	Save Direct Capture rate after last resonance
sd(resonances+1) = DCRate*fracErrDC
PartialRate(resonances+1) = DCRate
bigKeep(resonances+1) = "     "
if ( Detail == 1 ) then
	write(3,"(a)") "(For each resonance largest contribution to S.D. has  # )"
	write(3,"(a)") "(Over all resonances largest S.D. has  ## )"
end if
do nres = 1,resonances
!	Find and mark with # largest S.D. contribution in resonance
	big(1) = '  '
	big(2) = '  '
	big(3) = '  '
	big(4) = '  '
	if (OGdata(nres) /= 0.0d00 ) then
!		Omega-gamma was input
		if ( sdOGBoltz(nres) > sdBoltz(nres) ) then
			big(1) = ' #'
			bigKeep(nres) = "Omega"
		else
			big(2) = ' #'
			bigKeep(nres) = "Boltz"
		end if	
	else
!		Omega-gamma is estimated
		sdInput(1) = sdOGGp(nres)
		sdInput(2) = sdOGBoltz(nres)
		sdInput(3) = sdOGGg(nres)
		sdInput(4) = sdOGSpec(nres)
		call SortOutput(sdInput,sdOut,Keep,4)
		big(Keep(1)) = ' #'
		bigKeep(nres) = name(Keep(1))	
	end if
	if ( Detail == 1 ) then
		write(3,"(/a,i2)") "For resonance ",nres
		if (OGdata(nres) /= 0.0d00 ) then
			write(3,"(a,e10.4,a,e8.2,2(a,e9.3))") &
				"En +- S.D. ",En(nres)," +- ",sdEn(nres), &
				"  Omega-gamma +- S.D. ",OG(nres)," +- ",sdOG(nres)
			write(3,"(a,e10.3,a,a,e10.3,a)") &
				"S.D.s: omega-gamma ",sdOGBoltz(nres),big(1), &
				"  Boltzmann ",sdBoltz(nres),big(2)
		else
!			Omega-gamma is estimated
			write(3,"(a,e10.4,a,e8.2,a,f3.1,i2,2(a,e9.3))") &
				"En +- S.D. ",En(nres)," +- ",sdEn(nres), &
				"  J and L  ",JJ(nres),Lvalue(nres), &
				"  Gamma-p ",GP(nres),"  Gamma-g ",GG(nres)
			write(3,"(a,e11.3,a,a,e11.3,a)") &
				"S.D.s: Gp ",sdOGGp(nres),big(1), &
				"    Boltzmann ",sdOGBoltz(nres),big(2)
			write(3,"(a,e11.3,a,a,e11.3,a)") &
				"       Gg ",sdOGGg(nres),big(3), &
				"    Spec. fac.",sdOGSpec(nres),big(4)
		end if
		if ( nres == KeepRes(1) ) then
			big(1) = '##'
		else
			big(1) = '  '
		end if
		write(3,"(a,e10.3)") "Partial reaction rate ",PartialRate(nres)
!		Symmetric standard deviations as in eq.(4.6)
		write(3,"(a,e9.3,a,a,a,e9.3)") "Symmetric S.D.s (actual and fraction): ", &
			sd(nres)," ",big(1),"   ",fracSD(nres)
!		Asymmetric standard deviations use eq.(4.13)
		SDplus = PartialRate(nres)*(exp(fracSD(nres)) - 1.0d00)
		SDminus = PartialRate(nres)*(1.0d00 - exp(-fracSD(nres)))
		write(3,"(2(a,e9.3))") "Asymmetric S.D.s: + ",SDplus,"  - ",SDminus
	end if
end do
!	Order resonances and DC by total S.D. in reaction rate
call SortOutput(sd,sdSorted,KeepRes,resonances+1)
!	Totals over all resonances plus direct capture
SDTotalRate = sqrt(SumSquares)
write(3,"(/a)") &
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"
write(3,"(a,e11.3,a,e11.3)") &
	"  DCRate: ",DCRate," with fractional error: ",fracErrDC
!	Symmetric errors

write(3,"(2(a,e11.3))") &
	"  Total reaction rate +- S.D: ", TotalRate," +- ",SDTotalRate
!	Asymmetric total rates by using eq.(4.13)
Lower = TotalRate*exp(-SDTotalRate/TotalRate)
Upper = TotalRate*exp(SDTotalRate/TotalRate)
SDminus = TotalRate - Lower
SDPlus = Upper - TotalRate
!	Confidence level by eq.(4.17)
FracErrSq = (SDTotalRate/TotalRate)**2
Confidence = 68.3d00 - FracErrSq*(23.1d00 - 2.48d00*FracErrSq)

write(3,"(2(a,e9.3))") "  Asymmetric S.D.s: + ",SDplus,"  - ",SDminus
write(3,"(2(a,e11.3))") "  Asymmetric rates: Upper ",Upper,";   Lower ",Lower
write(3,"(a,f6.2,a)") "  Confidence level: ",Confidence,"%. (Unreliable if < 16%)"
write(3,"(/a)") "S.D. ordered by resonance"
write(3,*) "Resonance       S.D.     Partial Rate    Biggest"

write(4,"(f6.3,3(1x,1pe11.3),2x,1pe8.1)") T9,Lower,TotalRate, &
   Upper,Confidence

do nres = 1,resonances+1
	if ( KeepRes(nres) == resonances+1 ) then
		write(3,"(2(a,e11.3))") &
	  "  DC       ",sdSorted(nres),"    ",PartialRate(KeepRes(nres))
	else
	  write(3,"(a,i2,2(a,e11.3),2a)") &
	  "  ",KeepRes(nres),"       ",sdSorted(nres),"    ", &
		PartialRate(KeepRes(nres)),"      ",bigKeep(KeepRes(nres))
	end if
end do
write(3,*) &
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"
END SUBROUTINE SaveFile
!*********************************************************************
SUBROUTINE SortOutput(Input,Sorted,KeepRes,resonances)
!	Sorts aray  Input  into array  Sorted  by descending order
!	Original array is maintained
double precision Input(51),Sorted(51),biggest,temp
integer KeepRes(51),resonances,keep,nres,nsub

!	Copy original array
do nres = 1,resonances
	Sorted(nres) = Input(nres)
	KeepRes(nres) = nres
end do
if ( resonances == 1) then
 return
end if
!	Bubble sort; biggest to top
do nres = 1,resonances-1
	biggest = Sorted(nres)
	do nsub = nres+1,resonances
		if ( biggest < Sorted(nsub) ) then
!			swap values
			temp = biggest
			biggest = Sorted(nsub)
			Sorted(nsub) = temp
!			swap resonances
			keep = KeepRes(nres)
			KeepRes(nres) = KeepRes(nsub)
			KeepRes(nsub) = keep
		end if
	end do
	Sorted(nres) = biggest
end do
END SUBROUTINE SortOutput