c
c compositionuncertainties.f
c
c This code calculates the uncertainties in neutrino fluxes and 
c radiochemical neutrino capture rates using individual uncertainties in 
c heavy element compositions (began Friday, September 17, 2004). On 
c November 11, 2004, I added the option of combining coherently the Neon, Argon, and
c Oxygen abundances, since Ne and Ar are measured relative to O(see Asplund, Grevesse,
c and Sauval 2005).
c
c Partial derivatives are from the paper astro-ph/0412096, Bahcall and Serenelli. The relevant
c equations are given in that paper. 
c
c
c ***************************
c Precepts and numerical data
c ***************************
c
c Input composition uncertainties are 1-sigma; output uncertainties are also 1-sigma. (Gallium
c  neutrino absorption cross sections are an exception.)
c

c
c ***********************************************
c Information re: order of the 8 neutrino sources
c ***********************************************
c
c The order of the neutrino sources is: pp,pep,Be7,B8, N13,O15,F17,hep.
c
c************
c Definitions
c************
c
c Stndrd(I):  Standard solar neutrino fluxes.
c 
c
c CrossCl(I): Cross sections for the chlorine detector.  Unit: 10**-46 
c             cm**2. 
c
c CrossGa(I): Cross sections for the gallium detector.  Unit: 10**-46 
c             cm**2. 
c
c CrossLi(I): Cross sections for the Li detector.  Unit: 10**-46 
c             cm**2.   
c
c
c PartCl(I):  Product of flux times cross section for Ith neutrino 
c             source.  
c
c PartGa(I):  Same as PartCl(I) except for Ga instead of Cl.
c

c
c PartLi(I):  Product of flux times cross section for Ith neutrino 
c             source.  
c
c UncertJ:    The uncertainty in the parameter J defined by 
c             [1 + (dx sub j)/(x sub j) ] as in equation 23 
c             Bahcall-Ulrich RMP '82. 
c
c Uncertsq(J): Squares of uncertainties: Uncert(J)*Uncert(J).
c
c From Bahcall and Serenelli 2004: Composition derivatives.
c
c
c AlphaC(I):     Logariathmic Partial Derivatives of neutrino fluxes with respect 
c                to Carbon/Hydrogen abundance.
c
c AlphaN(I):     Logariathmic Partial Derivatives of neutrino fluxes with respect 
c                to Nitrogen/Hydrogen abundance.
cc
c AlphaO(I):     Logariathmic Partial Derivatives of neutrino fluxes with respect 
c                to Oxygen/Hydrogen abundance.
c
c AlphaNe(I):    Logariathmic Partial Derivatives of neutrino fluxes with respect 
c                to Neon/Hydrogen abundance.
c
c
c AlphaMg(I):    Logariathmic Partial Derivatives of neutrino fluxes with respect 
c                to Magnesium/Hydrogen abundance.
c
c AlphaSi(I):    Logariathmic Partial Derivatives of neutrino fluxes with respect 
c                to Silicon/Hydrogen abundance.
c
c AlphaS(I):     Logariathmic Partial Derivatives of neutrino fluxes with respect 
c                to Sulphur/Hydrogen abundance.
cc
c AlphaAr(I):    Logariathmic Partial Derivatives of neutrino fluxes with respect 
c                to Argon/Hydrogen abundance.
c
c AlphaFe(I):    Partial Derivatives of neutrino fluxes with respect to 3-4
c                rate. 
c UncertX:       The fractional uncertainty in the abundance of X. X can be Carbon,
c                Nitrogen, Oxygen,...,Iron.
c
c Carbon(K):     The fractional uncertainty in the k-th neutrino flux due to uncertainty
c                in the Carbon/Hydrogen abundance ratio.  Same for Nitrogen(K), ...,
c                Iron(K)
c
c*********************
c Dimension Statements
c*********************
c
      Dimension Stndrd(8)
      Dimension AlphaC(8),AlphaN(8),AlphaO(8),AlphaNe(8),AlphaMg(8)
      Dimension AlphaSi(8),AlphaS(8),AlphaAr(8),AlphaFe(8)
      Dimension CrossCl(8),CrossGa(8),CrossLi(8)
      Dimension PartCl(8),PartGa(8),PartLi(8)
      Dimension Carbon(8),Oxygen(8),Silicon(8),Sulphur(8),Argon(8)
c
c******************
c Type declarations
*******************
c     
c
      Real Nitrogen(8),Neon(8),Magnesium(8),Iron(8)
      Real N13Carbon,N13Nitrogen,N13Oxygen,N13Neon,
     *N13Magnesium,N13Sulphur,N13Silicon,N13Argon,N13Iron 
      Real O15Carbon,O15Nitrogen,O15Oxygen,O15Neon,
     *O15Magnesium,O15Sulphur,O15Silicon,O15Argon,O15Iron 
      Integer PP,PEP,BE7,B8,N13,O15,F17,HEP,CL,GA,LI
      Character*4 Source(11)
      Real N13coherent,O15coherent,F17coherent,hepcoherent
c
c ********************
c Parameter statements
c ********************
c
c Here are the  parameters used for indexing.
c
      Parameter (PP=1,PEP=2,BE7=3,B8=4,N13=5,O15=6)
      Parameter (F17=7,HEP=8,CL=9,GA=10,LI=11)
c
c
c*****************************************************************
c Input neutrino absorption cross sections and their uncertainties.
c*****************************************************************
c
c Take 37Cl cross section for 8B neutrinos from Bahcall
c  et al. (1996), Phys. Rev. C,54, 411.
c
      Data CrossCl/0.0,16.,2.4,11400.,1.7,6.8,6.9,42600./
c
c Gallium cross sections are from Phys Rev C 56, 3391,1997, see 
c  Table VI. These are 1-sigma uncertainties.
c
      Data CrossGa/11.72,204.,71.7,24000.,60.4,113.7,113.9,71400./  
c
c Have included the Li cross section for 7Be from PRD,49,3923(1994).
c
      Data CrossLi/0.0,655.,19.,39000.,42.4,246.,249.,84000./
c
c*********************************
c Open the output and input files.
c*********************************
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Open(1,file='compositionuncertainties.output',status='unknown')
c
c Input neutrino fluxes. Units of 10^10.
c
      Open(3,file = 'compositionuncertainties.dat', status = 'old')
c
c
c
c********************
c Write output title.
c********************
c
      Write(6,102)
      Write(1,102)
      Write(1,100)
      Write(6,100)
      Write(1,102)
      Write(6,102)
c
c 
c*****************************
c Read in the neutrino fluxes. 
c*****************************
c
c First four lines of compositiondata.dat are not read; they are skipped.
c
      Read(3,300)Skip
      Read(3,*)(Stndrd(I), I = 1,8)
      Write(6,501)
      Write(6,*)(Stndrd(I), I = 1,8)
c 
c****************************************
c Input neutrino composition derivatives.
c****************************************
c
c
c The tabulated quantites are the logarithmic partial derivatives
c  of neutrino fluxes with respect to composition parameters, defined by
c  Eq. (7.3) of Neutrino Astrophysics. Thus we tabulate
c  AlphaComposition_{i} = \partial (ln (flux_i) )/\partial (ln composition parameter) .
c
c
c
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c  Partial derivatives w.r.t. the solar model BP04 (Grevesse&Suaval 98 composition).
c  Reference paper is Bahcall and Pinsonneault, PRL, 92 (2004)121301.
c  Partial derivatives are from the paper astro-ph/0412096, Bahcall and Serenelli.      
c234567
c
       Data AlphaC/-0.014,-0.025,-0.002,0.030,0.845,0.826,0.033,-0.015/
       Data AlphaN/-0.003,-0.006,+0.002,0.011,0.181,0.209,0.010,-0.004/
       Data AlphaO/-0.006,-0.011,+0.052,0.121,0.079,0.093,1.102,-0.023/
       Data AlphaNe/-0.005,-0.005,+0.049,0.096,0.057,0.068,0.076,-0.017/
       Data AlphaMg/-0.005,-0.005,+0.051,0.096,0.060,0.070,0.078,-0.018/
       Data AlphaSi/-0.011,-0.014,+0.104,0.194,0.128,0.150,0.164,-0.037/
       Data AlphaS/-0.008,-0.017,+0.074,0.137,0.094,0.109,0.120,-0.028/
       Data AlphaAr/-0.002,-0.006,+0.018,0.034,0.024,0.028,0.031,-0.007/
       Data AlphaFe/-0.023,-0.065,+0.209,0.515,0.342,0.401,0.444,-0.069/
c
c  Partial derivatives w.r.t. the solar model BP04+ (lower values heavy element abundances).
c
c      Data AlphaC/-0.010,-0.018,0.005,0.035,0.846,0.824,0.035,-0.012/
c      Data AlphaN/-0.003,-0.004,+0.002,0.009,0.180,0.211,0.008,-0.003/
c      Data AlphaO/-0.005,-0.008,+0.046,0.099,0.055,0.068,1.074,-0.018/
c      Data AlphaNe/-0.003,-0.002,+0.033,0.064,0.036,0.043,0.048,-0.011/
c      Data AlphaMg/-0.005,-0.003,+0.057,0.107,0.065,0.077,0.086,-0.019/
c      Data AlphaSi/-0.010,-0.015,+0.115,0.212,0.139,0.163,0.180,-0.039/
c      Data AlphaS/-0.007,-0.015,+0.080,0.150,0.102,0.119,0.130,-0.029/
c      Data AlphaAr/-0.001,-0.003,+0.012,0.023,0.015,0.018,0.020,-0.005/
c      Data AlphaFe/-0.022,-0.062,+0.230,0.553,0.355,0.423,0.469,-0.072/
c
c Partial derivatives w.r.t. solar model same as BP04+ but uses slightly different composition.
c Uses AGS 2005 composition.
c
c
c      Data AlphaC/-0.009,-0.016,0.005,0.033,0.853,0.828,0.032,-0.011/
c      Data AlphaN/-0.002,-0.003,+0.002,0.007,0.169,0.202,0.006,-0.003/
c      Data AlphaO/-0.004,-0.005,+0.048,0.104,0.056,0.071,1.078,-0.018/
c      Data AlphaNe/-0.003,-0.001,+0.035,0.066,0.034,0.044,0.048,-0.012/
c      Data AlphaMg/-0.006,-0.004,+0.057,0.105,0.059,0.074,0.081,-0.018/
c      Data AlphaSi/-0.010,-0.014,+0.114,0.210,0.133,0.161,0.176,-0.037/
c      Data AlphaS/-0.007,-0.017,+0.081,0.149,0.097,0.116,0.128,-0.028/
c      Data AlphaAr/-0.001,-0.004,+0.013,0.025,0.016,0.020,0.022,-0.005/
c      Data AlphaFe/-0.0220,-0.059,+0.226,0.536,0.335,0.406,0.452,-0.069/ 
c
c
c Partial Derivatives for debugging purposes.
c
c      Data AlphaC/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/
c      Data AlphaN/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/
c      Data AlphaO/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/
c      Data AlphaNe/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/
c      Data AlphaMg/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/
c      Data AlphaSi/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/
c      Data AlphaS/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/
c      Data AlphaAr/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/
c      Data AlphaFe/1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/
c      
c      
c*******************************************************************************
c Write partial derivatives of fluxes with respect to each composition variable.
c*******************************************************************************
c
      Write(6,110)
      Write(1,110)
      Write(6,*)(AlphaC(I),I= 1,8)
      Write(1,*)(AlphaC(I),I= 1,8)
      Write(6,*)(AlphaN(I),I= 1,8)
      Write(1,*)(AlphaN(I),I= 1,8)
      Write(6,*)(AlphaO(I),I= 1,8)
      Write(1,*)(AlphaO(I),I= 1,8)
      Write(6,*)(AlphaNe(I),I= 1,8)
      Write(1,*)(AlphaNe(I),I= 1,8)
      Write(6,*)(AlphaMg(I),I= 1,8)
      Write(1,*)(AlphaMg(I),I= 1,8)
      Write(6,*)(AlphaSi(I),I= 1,8)
      Write(1,*)(AlphaSi(I),I= 1,8)
      Write(6,*)(AlphaS(I),I= 1,8)
      Write(1,*)(AlphaS(I),I= 1,8)
      Write(6,*)(AlphaAr(I),I= 1,8)
      Write(1,*)(AlphaAr(I),I= 1,8)
      Write(6,*)(AlphaFe(I),I= 1,8)
      Write(1,*)(AlphaFe(I),I= 1,8)
c
c********************************************
c Composition uncertainties from measurements
c********************************************
c 
c
c***************
c Uncertainties
c***************
c
c The uncertainties for Fe, S, Si, and Mag are from N. Grevesse and A. J. Sauval,
c Space Science Reviews 85, 161-174, 1998; Adv.Spac. Res. 30, 3-11 (2002).
c
c The oxygen uncertainty comes from Allende Prieto, C., Lambert, D. L. and Asplund, M. 2001, ApJ 556, L63.
c The abundance is 8.69 +- 0.05 . The same uncertainty is given in M. Asplund, 
c N. Grevesse, A. J. Sauval, C. Allende Prieto, and D. Kiselman, A&A, 417, 751, who give 8.66 +- 0.05 for O.
c This later paper also gives the uncertainties in neon and argon: Ne = 7.84 +- 0.06 and Ar = 6.18 +- 0.08
c (see their Section 4).
c
c The carbon uncertainty is from C. Allende Prieto, D. L. Lambert, and M. Asplund, ApJ 573, L137, 2002.
c They obtain C = 8.39 +- 0.04.
c
c The nitrogen uncertainty is from M. Asplund, astro-ph/0302, N = 7.80 +- 0.04.
c
c The uncertainties given below are those quoted by Asplund, Grevesse, and Sauval (2005). Corrected November 25,
c 2004. 
c
      Data UncertCarbon/1.122/
      Data UncertNitrogen/1.148/
      Data UncertOxygen/1.122/
      Data UncertNeon/1.148/
      Data UncertMagnesium/1.072/
      Data UncertSilicon/1.047/
      Data UncertSulphur/1.096/
      Data UncertArgon/1.202/
      Data UncertIron/1.072/     
c234567
c The uncertainties listed below are computed by regarding the differences between the Grevesse
c and Sauval (1998) abundances  and the recent lower abundances by Asplund, Grevesse, and Sauval (2005)
c as 1 sigma differences.
c234567
c
c       Data UncertCarbon/1.297/
c       Data UncertNitrogen/1.320/
c       Data UncertOxygen/1.387/
c       Data UncertNeon/1.539/
c       Data UncertMagnesium/1.115/
c       Data UncertSilicon/1.115/
c       Data UncertSulphur/1.092/
c       Data UncertArgon/1.496/
c       Data UncertIron/1.115/   
c
c Uncertainties for debugging purposes.
c
c       Data UncertCarbon/2.0/
c       Data UncertNitrogen/2.0/
c       Data UncertOxygen/2.0/
c       Data UncertNeon/2.0/
c       Data UncertMagnesium/2.0/
c       Data UncertSilicon/2.0/
c       Data UncertSulphur/2.0/
c       Data UncertArgon/2.0/
c       Data UncertIron/2.0/   
c
c
c********************************
c Write composition uncertainties
c********************************
c      
c
      Write(6,130)
      Write(1,130)
130   Format(//,1x,'Composition uncertainties: C,N,O,Ne,Mg,Si,S,Ar,Fe',
     *//,1x,'               1 + 1 sigma Uncertainty',//)
      Write(6,*)UncertCarbon,UncertNitrogen,UncertOxygen,UncertNeon,
     *UncertMagnesium,UncertSilicon,UncertSulphur,UncertArgon,
     *UncertIron
      Write(1,*)UncertCarbon,UncertNitrogen,UncertOxygen,UncertNeon,
     *UncertMagnesium,UncertSilicon,UncertSulphur,UncertArgon,
     *UncertIron      
c
c*************************************  
c Here are the labels used for output.
c************************************* 
c   
c23456789012345678901234567890123456789012345678901234567890123456789012
      Data Source/'pp','pep','7Be','8B','13N','15O','17F','hep','Cl',
     * 'Ga','Li'/
c
c*********************************************
c Calculate Rates: Fluxes times cross sections
c*********************************************
c
      RateCl = 0.0
      RateGa = 0.0
      RateLi = 0.0
c
c Calculate each source's contribution to each experiment.
c      
      Do I = 1,8      
        PartCl(I) = Stndrd(I)*CrossCl(I)
        PartGa(I) = Stndrd(I)*CrossGa(I)
        PartLi(I) = Stndrd(I)*CrossLi(I)
        RateCl = RateCl + PartCl(I)
        RateGa = RateGa + PartGa(I)
        RateLi = RateLi + PartLi(I)
      End do
c
c
c Write out individual contributions to each target.
c
      Write(6,103)
      Write(6,104)(Source(I),Stndrd(I),
     *PartCl(I),PartGa(I),PartLi(I),I=1,8)
c     
      Write(1,103)
      Write(1,104)(Source(I),Stndrd(I),
     *PartCl(I),PartGa(I),PartLi(I),I=1,8)
c
c**********************************************************
c Write rates  for radiochemical solar neutrino experiments
c**********************************************************
c
      Write(6,120)
      Write(6,121)RateCl
      Write(6,122)RateGa
      Write(6,123)RateLi
      Write(1,120)
      Write(1,121)RateCl
      Write(1,122)RateGa
      Write(1,123)RateLi
c
c
c ******************************************************
c Calculate the uncertainties for Cl, Ga, Li experiments 
c ******************************************************
c
c The equation used here is equation (23) of RMP 82.  I calculate here 
c  for j = p-p, 3-3, 3-4, 1-14,1-7, Z/X, luminosity, opacity, age, and 
c  neutrino cross sections.  I first calculate all of the partial 
c  uncertainties (the sum over i for each parameter separately 
c  in equation 23 ) and then sum them quadratically and take the square
c  root. At this stage, the uncertainties are still 3 sigma.
c
      Del11Cl = 0.0
      Del11Ga = 0.0
      Del11Li = 0.0
c
c Carbon
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Do I = 1,8
        DelCarbonCl = DelCarbonCl + PartCl(I)*
     1( UncertCarbon**alphaC(I) -1.0 )
        DelCarbonGa = DelCarbonGa + PartGa(I)*
     1( UncertCarbon**alphaC(I) -1.0 )
        DelCarbonLi = DelCarbonLi + PartLi(I)
     1*( UncertCarbon**alphaC(I) -1.0 )
c
      End Do
c
c Nitrogen
c
c
      Do I = 1,8
        DelNitrogenCl = DelNitrogenCl + PartCl(I)*
     1( UncertNitrogen**alphaN(I) -1.0 )
        DelNitrogenGa = DelNitrogenGa + PartGa(I)
     1*( UncertNitrogen**alphaN(I) -1.0 )
        DelNitrogenLi = DelNitrogenLi + PartLi(I)
     1*( UncertNitrogen**alphaN(I) -1.0 )
c
      End Do
c
c Oxygen
c
      Do I = 1,8
c234567
c23456789012345678901234567890123456789012345678901234567890123456789012
      DelOxygenCl = DelOxygenCl + PartCl(I)
     1*( UncertOxygen**alphaO(I) -1.0 )
      DelOxygenGa = DelOxygenGa + PartGa(I)
     1*( UncertOxygen**alphaO(I) -1.0 )
      DelOxygenLi = DelOxygenLi + PartLi(I)
     1*( UncertOxygen**alphaO(I) -1.0 )
c234567
      End Do
c
c Neon
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Do I = 1,8
        DelNeonCl = DelNeonCl + PartCl(I)*( UncertNeon**alphaNe(I) -1.0)
        DelNeonGa = DelNeonGa + PartGa(I)*( UncertNeon**alphaNe(I) -1.0)
        DelNeonLi = DelNeonLi + PartLi(I)*( UncertNeon**alphaNe(I) -1.0)
      End Do
c
c Magnesium
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Do I = 1,8
c234567
        DelMagnesiumCl = DelMagnesiumCl + PartCl(I)*
     1( UncertMagnesium**alphaMg(I) -1.0 )
        DelMagnesiumGa = DelMagnesiumGa + PartGa(I)*
     1( UncertMagnesium**alphaMg(I) -1.0 )
        DelMagnesiumLi = DelMagnesiumLi + PartLi(I)*
     1( UncertMagnesium**alphaMg(I) -1.0 )
c23456789012345678901234567890123456789012345678901234567890123456789012
      End Do
c
C Silicon
c 
c23456789012345678901234567890123456789012345678901234567890123456789012
      Do I = 1,8
        DelSiliconCl = DelSiliconCl + PartCl(I)
     1*( UncertSilicon**alphaSi(I) -1.0 )
        DelSiliconGa = DelSiliconGa + PartGa(I)
     1*( UncertSilicon**alphaSi(I) -1.0 )
        DelSiliconLi = DelSiliconLi + PartLi(I)
     1*( UncertSilicon**alphaSi(I) -1.0 )
c234567
      End Do
c
c Sulphur
c
c
      Do I = 1,8
        DelSulphurCl = DelSulphurCl + PartCl(I)
c
     1*( UncertSulphur**alphaS(I) -1.0 )
        DelSulphurGa = DelSulphurGa + PartGa(I)
     1*( UncertSulphur**alphaS(I) -1.0 )
        DelSulphurLi = DelSulphurLi + PartLi(I)
     1*( UncertSulphur**alphaS(I) -1.0 )
      End Do
c
c Argon
c
c
      Do I = 1,8
        DelArgonCl = DelArgonCl + PartCl(I)*
c234567
     1( UncertArgon**alphaAr(I) -1.0 )
        DelArgonGa = DelArgonGa + PartGa(I)
     1*( UncertArgon**alphaAr(I) -1.0 )
        DelArgonLi = DelArgonLi + PartLi(I)
     1*( UncertArgon**alphaAr(I) -1.0 )
      End Do
c
c Iron
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Do I = 1,8
        DelIronCl = DelIronCl + PartCl(I)*( UncertIron**alphaFe(I) -1.0)
        DelIronGa = DelIronGa + PartGa(I)*( UncertIron**alphaFe(I) -1.0)
        DelIronLi = DelIronLi + PartLi(I)*( UncertIron**alphaFe(I) -1.0)
      End Do
c
c
c **********************************************************************
c Add quadratically the uncertainties for the radiochemical experiments.
c ***********************************************************************
c
c****************************
c First the Cl uncertainties.
c****************************
c23456789012345678901234567890123456789012345678901234567890123456789012
      DelClComposition2 = DelCarbonCl*DelCarbonCl  
     1 + DelNitrogenCl*DelNitrogenCl + DelOxygenCl*DelOxygenCl
     1 + DelNeonCl*DelNeonCl + DelMagnesiumCl*DelMagnesiumCl  
     2 + DelSiliconCl*DelSiliconCl + DelSulphurCl*DelSulphurCl 
     3 + DelArgonCl*DelArgonCl + DelIronCl*DelIronCl 
c
c**********************************************
c Coherent combination of O,Ne,Ar uncertainties
c**********************************************
c
      DelCoherentCl = DelOxygenCl + DelNeonCl + DelArgonCl
      DelCoherentCl = DelCoherentCl*DelCoherentCl
      DelCoherentCl = DelCoherentCl + DelCarbonCl*DelCarbonCl 
     1+ DelNitrogenCl*DelNitrogenCl + DelMagnesiumCl*DelMagnesiumCl 
     1+ DelSiliconCl*DelSiliconCl + DelSulphurCl*DelSulphurCl 
     1+ DelIronCl*DelIronCl      
c      
c*************************************
c Now calculate Gallium uncertainties.
c*************************************
c
      DelGaComposition2 = DelCarbonGa*DelCarbonGa  
c234567
     1 + DelNitrogenGa*DelNitrogenGa + DelOxygenGa*DelOxygenGa
     1 + DelNeonGa*DelNeonGa +DelMagnesiumGa*DelMagnesiumGa  
     2 + DelSiliconGa*DelSiliconGa + DelSulphurGa*DelSulphurGa 
     3 + DelArgonGa*DelArgonGa + DelIronGa*DelIronGa 
c
c**********************************************
c Coherent combination of O,Ne,Ar uncertainties
c**********************************************
c
      DelCoherentGa = DelOxygenGa + DelNeonGa + DelArgonGa
      DelCoherentGa = DelCoherentGa*DelCoherentGa
      DelCoherentGa = DelCoherentGa + DelCarbonGa*DelCarbonGa 
     1+ DelNitrogenGa*DelNitrogenGa + DelMagnesiumGa*DelMagnesiumGa 
     1+ DelSiliconGa*DelSiliconGa + DelSulphurGa*DelSulphurGa 
     1+ DelIronGa*DelIronGa
      
c********************************
c Now calculate Li uncertainties.
c********************************
c
      DelLiComposition2 = DelCarbonLi*DelCarbonLi  
     1 + DelNitrogenLi*DelNitrogenLi + DelOxygenLi*DelOxygenLi
     1 + DelNeonLi*DelNeonLi +DelMagnesiumLi*DelMagnesiumLi  
     2 + DelSiliconLi*DelSiliconLi + DelSulphurLi*DelSulphurLi 
     3 + DelArgonLi*DelArgonLi + DelIronLi*DelIronLi 
c
c**********************************************
c Coherent combination of O,Ne,Ar uncertainties
c**********************************************
c
      DelCoherentLi = DelOxygenLi + DelNeonLi + DelArgonLi
      DelCoherentLi = DelCoherentLi*DelCoherentLi
      DelCoherentLi = DelCoherentLi + DelCarbonLi*DelCarbonLi 
     1+ DelNitrogenLi*DelNitrogenLi + DelMagnesiumLi*DelMagnesiumLi 
     1+ DelSiliconLi*DelSiliconLi + DelSulphurLi*DelSulphurLi 
     1+ DelIronLi*DelIronLi
c
c******************************************************
c Calculate square root of the quadratic uncertainties.
c*****************************************************
c
      UnClComposition = sqrt(DelClComposition2)
      UnGaComposition = sqrt(DelGaComposition2)
      UnLiComposition = sqrt(DelLiComposition2)
      UnCoherentCl = sqrt(DelCoherentCl)
      UnCoherentGa = sqrt(DelCoherentGa)
      UnCoherentLi = sqrt(DelCoherentLi)
c
c**************************************************************
c Write composition uncertainties for radiochemical experiments
c**************************************************************
c
c23456789012345678901234567890123456789012345678901234567890123456789012
c
      Write(1,507)
      Write(6,507)
      Write(1,503)UnClComposition, UnGaComposition, UnLiComposition
c      Write(6,502)DelClComposition2,DelGaComposition2,DelLiComposition2
      Write(6,503)UnClComposition,UnGaComposition,UnLiComposition
      Write(6,506)UnCoherentCl,UnCoherentGa,UnCoherentLi
      Write(1,506)UnCoherentCl,UnCoherentGa,UnCoherentLi
      Write(1,507)
      Write(6,507)
c
c Write Uncertainties for debugging.
c
c      Write(6,510)
c510   Format(1x,//,1x,'Debugging partial uncertainties.',/)
c      Write(6,*)DelCarbonCl,DelNitrogenCl,DelOxygenCl,DelNeonCl,
c     1DelMagnesiumCl,DelSiliconCl,DelSulphurCl,DelArgonCl,DelIronCl 
c      Write(6,510)
c      Write(6,*)DelCarbonGa,DelNitrogenGa,DelOxygenGa,DelNeonGa,
c     1DelMagnesiumGa,DelSiliconGa,DelSulphurGa,DelArgonGa,DelIronGa 
c      Write(6,510)
c      Write(6,*)DelCarbonLi,DelNitrogenLi,DelOxygenLi,DelNeonLi,
c     1DelMagnesiumLi,DelSiliconLi,DelSulphurLi,DelArgonLi,DelIronLi            
c     
c
c ********************************
c Uncertainties in neutrino fluxes
c ********************************
c
c Use improved version equation (24) of RMP 82. The version used here
c  reduces to the previous version when the individual uncertainties are
c  small. It is equation (16) of the Bahcall-Ulrich 1988 RMP paper.
c
c For all except B8, the uncertainty in the 7Be flux from the 7Be electron 
c  capture rate is negligible because the 7Be electron capture is so fast, 
c  about 10^3 times faster than the proton capture rate.
c
c Compute the sum of the squares and then take square root and divide
c  by three to get the 1 sigma uncertainties.
c
c  Use a parameter to control the index for partial derivatives
c
c pp neutrino flux
c
      K = PP
c
c
c Calculate individual contributions to pp uncertainty.
c      
c23456789012345678901234567890123456789012345678901234567890123456789012
      Carbon(K) = (UncertCarbon**alphaC(K) - 1.0)
      ppCarbon = Carbon(K)  
      ppCarbon = ppCarbon*ppCarbon
c 
      Nitrogen(K) = (UncertNitrogen**alphaN(K) - 1.0)
      ppNitrogen = Nitrogen(K)
      ppNitrogen = ppNitrogen*ppNitrogen
c    
      Oxygen(K) = (UncertOxygen**alphaO(K) - 1.0)
      ppOxygen = Oxygen(K)
      ppOxygen = ppOxygen*ppOxygen
c    
      Neon(K) = (UncertNeon**alphaNe(K) - 1.0)
      ppNeon = Neon(K)
      ppNeon = ppNeon*ppNeon
c    
      Magnesium(K) = (UncertMagnesium**alphaMg(K) - 1.0)
      ppMagnesium = Magnesium(K)
      ppMagnesium = ppMagnesium*ppMagnesium
c    
      Silicon(K) = (UncertSilicon**alphaSi(K) - 1.0)
      ppSilicon = Silicon(K)
      ppSilicon = ppSilicon*ppSilicon
c    
      Sulphur(K) = (UncertSulphur**alphaS(K) - 1.0)
      ppSulphur = Sulphur(K)
      ppSulphur = ppSulphur*ppSulphur
c    
      Argon(K) = (UncertArgon**alphaAr(K) - 1.0)
      ppArgon = Argon(K)
      ppArgon = ppArgon*ppArgon
c    
      Iron(K) = (UncertIron**alphaFe(K) - 1.0)
      ppIron = Iron(K)
      ppIron = ppIron*ppIron
c       
c Sum the squares
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Uncertsqpp = ppCarbon+ppNitrogen+ppOxygen+ppNeon+ppMagnesium
     1 +ppSilicon+ppSulphur+ppArgon+ppIron
c
       Uncertpp = sqrt(Uncertsqpp)
c
c **********************************************************
c Coherent addition of oxygen, neon, and argon uncertainties.
c ***********************************************************
c
      ppcoherent = Oxygen(K) + Neon(K) + Argon(K)
      ppcoherent = ppcoherent*ppcoherent
      Uncertppcoherent = ppcoherent+ppCarbon+ppNitrogen+ppMagnesium
     1+ppSilicon+ppSulphur+ppIron
      Uncertppcoherent = sqrt(Uncertppcoherent)
c
c234567
c***************************************
c Write pp partial uncertainties
c***************************************
c
      Write(6,621)
621   Format(1x,'Partial Squared uncertainties of pp Neutrino flux.',/)
      Write(6,*)PPCarbon,PPNitrogen,PPOxygen,PPNeon, PPMagnesium
      Write(6,*)PPSilicon,PPSulphur,PPArgon,PPIron
      Write(1,621)
      Write(1,*)PPCarbon,PPNitrogen,PPOxygen,PPNeon, PPMagnesium
      Write(1,*)PPSilicon,PPSulphur,PPArgon,PPIron
      Write(1,102)
      Write(6,102)
c      
c******************
c pep neutrino flux
c******************
c
      K = PEP
c
c Calculate individual contributions to pep uncertainty.
c  
      Carbon(K) = (UncertCarbon**alphaC(K) - 1.0) 
      pepCarbon = Carbon(K)
      pepCarbon = pepCarbon*pepCarbon
c 
      Nitrogen(K) = (UncertNitrogen**alphaN(K) - 1.0)
      pepNitrogen = Nitrogen(K)
      pepNitrogen = pepNitrogen*pepNitrogen
c    
      Oxygen(K) = (UncertOxygen**alphaO(K) - 1.0)
      pepOxygen = Oxygen(K)
      pepOxygen = pepOxygen*pepOxygen
c    
      Neon(K) = (UncertNeon**alphaNe(K) - 1.0)
      pepNeon = Neon(K)
      pepNeon = pepNeon*pepNeon
c    
      Magnesium(K) = (UncertMagnesium**alphaMg(K) - 1.0)
      pepMagnesium = Magnesium(K)
      pepMagnesium = pepMagnesium*pepMagnesium
c    
      Silicon(K) = (UncertSilicon**alphaSi(K) - 1.0)
      pepSilicon = Silicon(K)
      pepSilicon = pepSilicon*pepSilicon
c    
      Sulphur(K) = (UncertSulphur**alphaS(K) - 1.0)
      pepSulphur = Sulphur(K)
      pepSulphur = pepSulphur*pepSulphur
c    
      Argon(K) = (UncertArgon**alphaAr(K) - 1.0)
      pepArgon = Argon(K)
      pepArgon = (UncertArgon**alphaAr(K) - 1.0)
      pepArgon = pepArgon*pepArgon
c    
      Iron(K) = (UncertIron**alphaFe(K) - 1.0)
      pepIron = Iron(K)
      pepIron = pepIron*pepIron
c       
c Sum the squares
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Uncertsqpep = pepCarbon+pepNitrogen+pepOxygen+pepNeon
     1 +pepMagnesium+pepSilicon+pepSulphur+pepArgon+pepIron
c
       Uncertpep = sqrt(Uncertsqpep)
c
c
c **********************************************************
c Coherent addition of oxygen, neon, and argon uncertainties.
c ***********************************************************
c
      pepcoherent = Oxygen(K) + Neon(K) + Argon(K) 
      pepcoherent = pepcoherent*pepcoherent
      Uncertpepcoherent=pepcoherent+pepCarbon+pepNitrogen+pepMagnesium
     1+pepSilicon+pepSulphur+pepIron    
      Uncertpepcoherent = sqrt(Uncertpepcoherent)
c
c
c***************************************
c Write pep partial uncertainties
c***************************************
c
      Write(6,631)
631   Format(1x,'Partial squared uncertainties of pep Neutrino flux.',/)
      Write(6,*)PEPCarbon,PEPNitrogen,PEPOxygen,PEPNeon, PEPMagnesium
      Write(6,*)PEPSilicon,PEPSulphur,PEPArgon,PEPIron
      Write(1,631)
      Write(1,*)PEPCarbon,PEPNitrogen,PEPOxygen,PEPNeon, PEPMagnesium
      Write(1,*)PEPSilicon,PEPSulphur,PEPArgon,PEPIron

      Write(6,102)
      Write(1,102)
c
c 
c******************
c Be7 neutrino flux
c******************
c
      K = BE7
c
c 
      Carbon(K) = (UncertCarbon**alphaC(K) - 1.0) 
      Be7Carbon = Carbon(K)
      be7Carbon = be7Carbon*be7Carbon
c 
      Nitrogen(K) = (UncertNitrogen**alphaN(K) - 1.0)
      Be7Nitrogen = Nitrogen(K)
      be7Nitrogen = be7Nitrogen*be7Nitrogen
c    
      Oxygen(K) = (UncertOxygen**alphaO(K) - 1.0)
      Be7Oxygen = Oxygen(K)
      be7Oxygen = be7Oxygen*be7Oxygen
c    
      Neon(K) = (UncertNeon**alphaNe(K) - 1.0)
      Be7Neon = Neon(K)
      be7Neon = be7Neon*be7Neon
c    
      Magnesium(K) = (UncertMagnesium**alphaMg(K) - 1.0)
      Be7Magnesium = Magnesium(K)
      be7Magnesium = be7Magnesium*be7Magnesium
c    
      Silicon(K) = (UncertSilicon**alphaSi(K) - 1.0)
      Be7Silicon = Silicon(K)
      be7Silicon = be7Silicon*be7Silicon
c    
      Sulphur(K) = (UncertSulphur**alphaS(K) - 1.0)
      Be7Sulphur = Sulphur(K)
      be7Sulphur = be7Sulphur*be7Sulphur
c    
      Argon(K) = (UncertArgon**alphaAr(K) - 1.0)
      Be7Argon = Argon(K)
      be7Argon = be7Argon*be7Argon
c    
      Iron(K) = (UncertIron**alphaFe(K) - 1.0)
      Be7Iron = Iron(K)
      be7Iron = be7Iron*be7Iron
c       
c Sum the squares
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Uncertsqbe7 = be7Carbon+be7Nitrogen+be7Oxygen+be7Neon
     1 +be7Magnesium+be7Silicon+be7Sulphur+be7Argon+be7Iron
c
       Uncertbe7 = sqrt(Uncertsqbe7)     
c
c
c
c **********************************************************
c Coherent addition of oxygen, neon, and argon uncertainties.
c ***********************************************************
c
      Be7coherent = Oxygen(K) + Neon(K) + Argon(K)
      Be7coherent = Be7coherent*Be7coherent
      UncertBe7coherent=Be7coherent+Be7Carbon+Be7Nitrogen+Be7Magnesium
     1+Be7Silicon+Be7Sulphur+Be7Iron
      UncertBe7coherent = sqrt(UncertBe7coherent)
c
c
c***************************************
c Write Be7 partial uncertainties
c***************************************
c
      Write(6,641)
641   Format(1x,'Partial squared uncertainties of Be7 Neutrino flux.',/)
      Write(6,*)Be7Carbon,Be7Nitrogen,Be7Oxygen,Be7Neon, Be7Magnesium
      Write(6,*)Be7Silicon,Be7Sulphur,Be7Argon,Be7Iron
      Write(1,641)
      Write(1,*)Be7Carbon,Be7Nitrogen,Be7Oxygen,Be7Neon, Be7Magnesium
      Write(1,*)Be7Silicon,Be7Sulphur,Be7Argon,Be7Iron
      Write(6,102)
      Write(1,102)
c
c*****************
c B8 neutrino flux
c*****************
c
      K = B8
c
c Calculate individual contributions to B8 uncertainty.
c      
      Carbon(K) = (UncertCarbon**alphaC(K) - 1.0) 
      B8Carbon = Carbon(K)
      b8Carbon = b8Carbon*b8Carbon
c 
      Nitrogen(K) = (UncertNitrogen**alphaN(K) - 1.0)
      B8Nitrogen = Nitrogen(K)
      b8Nitrogen = b8Nitrogen*b8Nitrogen
c    
      Oxygen(K) = (UncertOxygen**alphaO(K) - 1.0)
      B8Oxygen = Oxygen(K)
      b8Oxygen = b8Oxygen*b8Oxygen
c    
      Neon(K) = (UncertNeon**alphaNe(K) - 1.0)
      B8Neon = Neon(K)
      b8Neon = b8Neon*b8Neon
c    
      Magnesium(K) = (UncertMagnesium**alphaMg(K) - 1.0)
      B8Magnesium = Magnesium(K)
      b8Magnesium = b8Magnesium*b8Magnesium
c    
      Silicon(K) = (UncertSilicon**alphaSi(K) - 1.0)
      B8Silicon = Silicon(K)
      b8Silicon = b8Silicon*b8Silicon
c    
      Sulphur(K) = (UncertSulphur**alphaS(K) - 1.0)
      B8Sulphur = Sulphur(K)
      b8Sulphur = b8Sulphur*b8Sulphur
c    
      Argon(K) = (UncertArgon**alphaAr(K) - 1.0)
      B8Argon = Argon(K)
      b8Argon = b8Argon*b8Argon
c    
      Iron(K) = (UncertIron**alphaFe(K) - 1.0)
      B8Iron = Iron(K)
      b8Iron = b8Iron*b8Iron
c       
c Sum the squares
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Uncertsqb8 = b8Carbon+b8Nitrogen+b8Oxygen+b8Neon
     1 +b8Magnesium+b8Silicon+b8Sulphur+b8Argon+b8Iron
c
       Uncertb8 = sqrt(Uncertsqb8)
c
c
c **********************************************************
c Coherent addition of oxygen, neon, and argon uncertainties.
c ***********************************************************
c
      B8coherent = Oxygen(K) + Neon(K) + Argon(K) 
      B8coherent = B8coherent*B8coherent
      UncertB8coherent=B8coherent+B8Carbon+B8Nitrogen+B8Magnesium
     1+B8Silicon+B8Sulphur+B8Iron    
      UncertB8coherent = sqrt(UncertB8coherent)
c
c
c***************************************
c Write B8 partial uncertainties
c***************************************
c
      Write(6,651)
651   Format(1x,'Partial squared uncertainties of B8 Neutrino flux.',/)
      Write(6,*)B8Carbon,B8Nitrogen,B8Oxygen,B8Neon, B8Magnesium
      Write(6,*)B8Silicon,B8Sulphur,B8Argon,B8Iron
      Write(1,651)
      Write(1,*)B8Carbon,B8Nitrogen,B8Oxygen,B8Neon, B8Magnesium
      Write(1,*)B8Silicon,B8Sulphur,B8Argon,B8Iron
      Write(6,102)
      Write(1,102)
c            
c******************
c N13 neutrino flux
c******************
c
      K = N13
c
c Calculate individual contributions to N13 uncertainty.
c  
      Carbon(K) = (UncertCarbon**alphaC(K) - 1.0) 
      N13Carbon = Carbon(K)
      N13Carbon = N13Carbon*N13Carbon
c 
      Nitrogen(K) = (UncertNitrogen**alphaN(K) - 1.0)
      N13Nitrogen = Nitrogen(K)
      N13Nitrogen = N13Nitrogen*N13Nitrogen
c    
      Oxygen(K) = (UncertOxygen**alphaO(K) - 1.0)
      N13Oxygen = Oxygen(K)
      N13Oxygen = N13Oxygen*N13Oxygen
c
      Neon(K) = (UncertNeon**alphaNe(K) - 1.0)
      N13Neon = Neon(K)          
      N13Neon = N13Neon*N13Neon
c    
      Magnesium(K) = (UncertMagnesium**alphaMg(K) - 1.0)
      N13Magnesium = Magnesium(K)
      N13Magnesium = N13Magnesium*N13Magnesium
c    
      Silicon(K) = (UncertSilicon**alphaSi(K) - 1.0)
      N13Silicon = Silicon(K)
      N13Silicon = N13Silicon*N13Silicon
c    
      Sulphur(K) = (UncertSulphur**alphaS(K) - 1.0)
      N13Sulphur = Sulphur(K)
      N13Sulphur = N13Sulphur*N13Sulphur
c    
      Argon(K) = (UncertArgon**alphaAr(K) - 1.0)
      N13Argon = Argon(K)
      N13Argon = N13Argon*N13Argon
c    
      Iron(K) = (UncertIron**alphaFe(K) - 1.0)
      N13Iron = Iron(K)
      N13Iron = N13Iron*N13Iron
c
c***************************************
c Write N13 partial uncertainties
c***************************************
c
      Write(6,601)
601   Format(1x,'Partial squared uncertainties of N13 neutrino flux.',/)
      Write(6,*)N13Carbon,N13Nitrogen,N13Oxygen,N13Neon, N13Magnesium
      Write(6,*)N13Silicon,N13Sulphur,N13Argon,N13Iron
      Write(1,601)
      Write(1,*)N13Carbon,N13Nitrogen,N13Oxygen,N13Neon, N13Magnesium
      Write(1,*)N13Silicon,N13Sulphur,N13Argon,N13Iron
      Write(6,102)
      Write(1,102)
c      
c *********************************************************************
c       
c Sum the squares
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      UncertsqN13 = N13Carbon+N13Nitrogen+N13Oxygen+N13Neon
     1 +N13Magnesium+N13Silicon+N13Sulphur+N13Argon+N13Iron
c
       UncertN13 = sqrt(UncertsqN13)
c
c
c **********************************************************
c Coherent addition of oxygen, neon, and argon uncertainties.
c ***********************************************************
c
      N13coherent = Oxygen(K) + Neon(K) + Argon(K) 
      N13coherent = N13coherent*N13coherent
      UncertN13coherent=N13coherent+N13Carbon+N13Nitrogen+N13Magnesium
     1+N13Silicon+N13Sulphur+N13Iron    
      UncertN13coherent = sqrt(UncertN13coherent)
c
c******************
c O15 neutrino flux
c******************
c
      K = O15
c
c Calculate individual contributions to O15 uncertainty.
c 
      Carbon(K) = (UncertCarbon**alphaC(K) - 1.0) 
      O15Carbon = Carbon(K)
      O15Carbon = O15Carbon*O15Carbon
c 
      Nitrogen(K) = (UncertNitrogen**alphaN(K) - 1.0)
      O15Nitrogen = Nitrogen(K)
      O15Nitrogen = O15Nitrogen*O15Nitrogen
c    
      Oxygen(K) = (UncertOxygen**alphaO(K) - 1.0)
      O15Oxygen = Oxygen(K)
      O15Oxygen = O15Oxygen*O15Oxygen
c    
      Neon(K) = (UncertNeon**alphaNe(K) - 1.0)
      O15Neon = Neon(K)
      O15Neon = O15Neon*O15Neon
c    
      Magnesium(K) = (UncertMagnesium**alphaMg(K) - 1.0)
      O15Magnesium = Magnesium(K)
      O15Magnesium = O15Magnesium*O15Magnesium
c    
      Silicon(K) = (UncertSilicon**alphaSi(K) - 1.0)
      O15Silicon = Silicon(K)
      O15Silicon = O15Silicon*O15Silicon
c    
      Sulphur(K) = (UncertSulphur**alphaS(K) - 1.0)
      O15Sulphur = Sulphur(K)
      O15Sulphur = O15Sulphur*O15Sulphur
c    
      Argon(K) = (UncertArgon**alphaAr(K) - 1.0)
      O15Argon = Argon(K)
      O15Argon = O15Argon*O15Argon
c    
      Iron(K) = (UncertIron**alphaFe(K) - 1.0)
      O15Iron = Iron(K)
      O15Iron = O15Iron*O15Iron
c
c***************************************
c Write O15 partial uncertainties
c***************************************
c
      Write(6,611)
611   Format(1x,'Partial squared uncertainties of O15 Neutrino flux.',/)
      Write(6,*)O15Carbon,O15Nitrogen,O15Oxygen,O15Neon, O15Magnesium
      Write(6,*)O15Silicon,O15Sulphur,O15Argon,O15Iron
      Write(1,611)
      Write(1,*)O15Carbon,O15Nitrogen,O15Oxygen,O15Neon, O15Magnesium
      Write(1,*)O15Silicon,O15Sulphur,O15Argon,O15Iron
      Write(6,102)
      Write(1,102)
c  
c *********************************************************************
c       
c       
c Sum the squares
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      UncertsqO15 = O15Carbon+O15Nitrogen+O15Oxygen+O15Neon
     1 +O15Magnesium+O15Silicon+O15Sulphur+O15Argon+O15Iron
c
       UncertO15 = sqrt(UncertsqO15)
c
c
c **********************************************************
c Coherent addition of oxygen, neon, and argon uncertainties.
c ***********************************************************
c
      O15coherent = Oxygen(K) + Neon(K) + Argon(K) 
      O15coherent = O15coherent*O15coherent
      UncertO15coherent=O15coherent+O15Carbon+O15Nitrogen+O15Magnesium
     1+O15Silicon+O15Sulphur+O15Iron    
      UncertO15coherent = sqrt(UncertO15coherent)
c
c 
c******************
c F17 neutrino flux
c******************
c
      K = F17
c
c Calculate individual contributions to F17 uncertainty.
c      
      Carbon(K) = (UncertCarbon**alphaC(K) - 1.0) 
      F17Carbon = Carbon(K)
      F17Carbon = F17Carbon*F17Carbon
c 
      Nitrogen(K) = (UncertNitrogen**alphaN(K) - 1.0)
      F17Nitrogen = Nitrogen(K)
      F17Nitrogen = F17Nitrogen*F17Nitrogen
c    
      Oxygen(K) = (UncertOxygen**alphaO(K) - 1.0)
      F17Oxygen = Oxygen(K)
      F17Oxygen = F17Oxygen*F17Oxygen
c    
      Neon(K) = (UncertNeon**alphaNe(K) - 1.0)
      F17Neon = Neon(K)
      F17Neon = F17Neon*F17Neon
c    
      Magnesium(K) = (UncertMagnesium**alphaMg(K) - 1.0)
      F17Magnesium = Magnesium(K)
      F17Magnesium = F17Magnesium*F17Magnesium
c    
      Silicon(K) = (UncertSilicon**alphaSi(K) - 1.0)
      F17Silicon = Silicon(K)
      F17Silicon = F17Silicon*F17Silicon
c    
      Sulphur(K) = (UncertSulphur**alphaS(K) - 1.0)
      F17Sulphur = Sulphur(K)
      F17Sulphur = F17Sulphur*F17Sulphur
c    
      Argon(K) = (UncertArgon**alphaAr(K) - 1.0)
      F17Argon = Argon(K)
      F17Argon = F17Argon*F17Argon
c    
      Iron(K) = (UncertIron**alphaFe(K) - 1.0)
      F17Iron = Iron(K)
      F17Iron = F17Iron*F17Iron
c       
c Sum the squares
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      UncertsqF17 = F17Carbon+F17Nitrogen+F17Oxygen+F17Neon
     1 +F17Magnesium+F17Silicon+F17Sulphur+F17Argon+F17Iron
c
       UncertF17 = sqrt(UncertsqF17)
c
c
c **********************************************************
c Coherent addition of oxygen, neon, and argon uncertainties.
c ***********************************************************
c
      F17coherent = Oxygen(K) + Neon(K) + Argon(K) 
      F17coherent = F17coherent*F17coherent
      UncertF17coherent=F17coherent+F17Carbon+F17Nitrogen+F17Magnesium
     1+F17Silicon+F17Sulphur+F17Iron    
      UncertF17coherent = sqrt(UncertF17coherent)
c
c
c***************************************
c Write F17 partial uncertainties
c***************************************
c
      Write(6,661)
661   Format(1x,'Partial squared uncertainties of F17 Neutrino flux.',/)
      Write(6,*)F17Carbon,F17Nitrogen,F17Oxygen,F17Neon, F17Magnesium
      Write(6,*)F17Silicon,F17Sulphur,F17Argon,F17Iron
      Write(1,661)
      Write(1,*)F17Carbon,F17Nitrogen,F17Oxygen,F17Neon, F17Magnesium
      Write(1,*)F17Silicon,F17Sulphur,F17Argon,F17Iron
      Write(6,102)
      Write(1,102)
c      
c******************
c hep neutrino flux
c******************
c
      K = HEP
c
c Calculate individual contributions to hep uncertainty.
c      
      Carbon(K) = (UncertCarbon**alphaC(K) - 1.0) 
      hepCarbon = Carbon(K)
      hepCarbon = hepCarbon*hepCarbon
c 
      Nitrogen(K) = (UncertNitrogen**alphaN(K) - 1.0)
      hepNitrogen = Nitrogen(K)
      hepNitrogen = hepNitrogen*hepNitrogen
c    
      Oxygen(K) = (UncertOxygen**alphaO(K) - 1.0)
      hepOxygen = Oxygen(K)
      hepOxygen = hepOxygen*hepOxygen
c    
      Neon(K) = (UncertNeon**alphaNe(K) - 1.0)
      hepNeon = Neon(K)
      hepNeon = hepNeon*hepNeon
c    
      Magnesium(K) = (UncertMagnesium**alphaMg(K) - 1.0)
      hepMagnesium = Magnesium(K)
      hepMagnesium = hepMagnesium*hepMagnesium
c    
      Silicon(K) = (UncertSilicon**alphaSi(K) - 1.0)
      hepSilicon = Silicon(K)
      hepSilicon = hepSilicon*hepSilicon
c    
      Sulphur(K) = (UncertSulphur**alphaS(K) - 1.0)
      hepSulphur = Sulphur(K)
      hepSulphur = hepSulphur*hepSulphur
c    
      Argon(K) = (UncertArgon**alphaAr(K) - 1.0)
      hepArgon = Argon(K)
      hepArgon = hepArgon*hepArgon
c    
      Iron(K) = (UncertIron**alphaFe(K) - 1.0)
      hepIron = Iron(K)
      hepIron = hepIron*hepIron
c       
c Sum the squares
c
c23456789012345678901234567890123456789012345678901234567890123456789012
      Uncertsqhep = hepCarbon+hepNitrogen+hepOxygen+hepNeon
     1 +hepMagnesium+hepSilicon+hepSulphur+hepArgon+hepIron
c
       Uncerthep = sqrt(Uncertsqhep)
c
c
c **********************************************************
c Coherent addition of oxygen, neon, and argon uncertainties.
c ***********************************************************
c
      hepcoherent = Oxygen(K) + Neon(K) + Argon(K) 
      hepcoherent = hepcoherent*hepcoherent
      Uncerthepcoherent=hepcoherent+hepCarbon+hepNitrogen+hepMagnesium
     1+hepSilicon+hepSulphur+hepIron     
      Uncerthepcoherent = sqrt(Uncerthepcoherent)
c
c
c***************************************
c Write hep partial uncertainties
c***************************************
c
      Write(6,671)
671   Format(1x,'Partial squared uncertainties of hep Neutrino flux.',/)
      Write(6,*)hepCarbon,hepNitrogen,hepOxygen,hepNeon, hepMagnesium
      Write(6,*)hepSilicon,hepSulphur,hepArgon,hepIron
      Write(1,671)
      Write(1,*)hepCarbon,hepNitrogen,hepOxygen,hepNeon, hepMagnesium
      Write(1,*)hepSilicon,hepSulphur,hepArgon,hepIron
      Write(6,102)
      Write(1,102)
c
c
c**************************************
c Write out neutrino flux uncertainties 
c**************************************
c
c *******************
c Total Uncertainties
c *******************
c
      Write(1,504)Uncertpp,Uncertpep,Uncertbe7,Uncertb8,Uncertn13,
     1Uncerto15,Uncertf17,Uncerthep
      Write(6,504)Uncertpp,Uncertpep,Uncertbe7,Uncertb8,Uncertn13,
     1Uncerto15,Uncertf17,Uncerthep
      Write(1,505)Uncertppcoherent,Uncertpepcoherent,Uncertbe7coherent,
     1UncertB8coherent,UncertN13coherent,UncertO15coherent,
     1UncertF17coherent,Uncerthepcoherent 
      Write(6,505)Uncertppcoherent,Uncertpepcoherent,Uncertbe7coherent,
     1UncertB8coherent,UncertN13coherent,UncertO15coherent,
     1UncertF17coherent,Uncerthepcoherent        
c
c ****************************
c Signed partial uncertainties
c ****************************
c
      Write(1,680)(Carbon(K),K = 1,8)
      Write(6,680)(Carbon(K),K = 1,8)
c23456789012345678901234567890123456789012345678901234567890123456789012
680   Format(1x,'Carbon Signed Uncertainties, pp,pep,7Be,8B,...,hep',//,
     11P8E10.2,//)
      Write(1,681)(Nitrogen(K),K = 1,8)
      Write(6,681)(Nitrogen(K),K = 1,8)
681   Format(1x,'N Signed Uncertainties, pp,pep,7Be,8B,...,hep',//,
     11P8E10.2,//)
      Write(1,682)(Oxygen(K),K = 1,8)
      Write(6,682)(Oxygen(K),K = 1,8)
682   Format(1x,'Oxy. Signed Uncertainties, pp,pep,7Be,8B,...,hep',//,
     11P8E10.2,//)
      Write(1,683)(Neon(K),K = 1,8)
      Write(6,683)(Neon(K),K = 1,8)
683   Format(1x,'Neon signed Uncertainties, pp,pep,7Be,8B,...,hep',//,
     11P8E10.2,//) 
      Write(1,684)(Magnesium(K),K = 1,8)
      Write(6,684)(Magnesium(K),K = 1,8)
684   Format(1x,'Mg signed Uncertainties, pp,pep,7Be,8B,...,hep',//,
     11P8E10.2,//)  
      Write(1,685)(Silicon(K),K = 1,8)
      Write(6,685)(Silicon(K),K = 1,8)
685   Format(1x,'Si signed Uncertainties, pp,pep,7Be,8B,...,hep',//,
     11P8E10.2,//)  
      Write(1,686)(Sulphur(K),K = 1,8)
      Write(6,686)(Sulphur(K),K = 1,8)
686   Format(1x,'Sulphur signed Uncertainties, pp,pep,7Be,8B,...,hep',//,
     11P8E10.2,//)    
      Write(1,687)(Argon(K),K = 1,8)
      Write(6,687)(Argon(K),K = 1,8)
687   Format(1x,'Argon signed Uncertainties, pp,pep,7Be,8B,...,hep',//,
     11P8E10.2,//)    
      Write(1,688)(Iron(K),K = 1,8)
      Write(6,688)(Iron(K),K = 1,8)
688   Format(1x,'Iron signed Uncertainties, pp,pep,7Be,8B,...,hep',//,
     11P8E10.2,//)    
c
c          
c ********************************
c Format Statements
c ********************************
c
c23456789012345678901234567890123456789012345678901234567890123456789012
100   Format(1x,'Ucertainties in Fluxes and Rates due to Composition')
102   Format(1x, //)
103   Format(//,
     *'                            SNU Contributions'//
     *'   Source    Fluxes         Cl        Ga        Li'/)
104   Format(1x,A8,2x,E10.3, F10.2, F10.1, F10.1)
110   Format(1x,'Partial derivatives fluxes wrt C,N,O,Ne,Mg,Si,S,Ar,Fe',
     *//)
120   Format(//'Capture Rates for Radiochemical Neutrino Experiments',/)
121   Format(' Cl rate = ',F7.2)
122   Format(' Ga rate = ',F7.2)
123   Format(' Li rate = ',F7.2,/)    
300   Format(///A1)
501   Format(//,1x, 'The 8 neutrino fluxes are listed below.',/)
503   Format(1x,'Radiochemical Uncertainties: DelCl  DelGa  DelLi',//
     1,25x,F10.2,F7.2, F7.2,//)
504   Format(/,'        *********************************************',/,
     18x,'Total Neutrino flux Composition Uncertainites',
     1/,'        *********************************************',//,
     11x, '    pp    pep    7Be     8B    13N    15O    17F    hep',//,
     11x, 8F7.4,//)
505   Format(/,'        ************************************************
     1********',/,
     18x,'Coherent Total Neutrino flux Composition Uncertainites',
     1/,'        *******************************************************
     1*',//,
     11x, '    pp    pep    7Be     8B    13N    15O    17F    hep',//,
     11x, 8F7.4,////)
506   Format(1x,'Coherent Radiochemical Uncertainties: DelCl  DelGa  Del
     1Li',//,34x,F10.2,2F7.2,//)
507   Format(1x,'*******************************************************
     1*****',//) 
      Stop
      End