SUBROUTINE STATE(MODE,ID,T,ANE) C =============================== C C For MODE=0 - initialization of the basic parameters for the C chemical species C for MODE=2 - solves the set of LTE Saha equations to determine C the total charge due to ionization of non-explicit C chemical species C for MODE=3 - similar as MODE=2, but also with derivatives wrt C temperature and electron density (called from BPOP, C ie. from the linearization step) C for MODE=1 - similar as MODE=2, but the total charge is evaluated C summing the contributions of both non-explicit and C explicit chemical species (called from LTEGR) C C Input for MODE > 0: C C T - temperature C ANE - electron density C C Output for MODE > 0 (through COMMON/STATEP) C C C Q - total charge, relative to the reference atom C QM - charge of H- (in LTE), evaluated only if H- is not C explicit ion, and if desired (if IHM=1) C DQT - derivative of Q wrt temperature C DQN - derivative of Q wrt electron density C DQM - derivative of QM wrt temperature C ENER - internal energy (of all species in all considered C ionization stages (needed only for evaluating C thermodynamic derivatoves if convection is considerd) C QREF - total charge due to the reference species C DQTR - derivative of QREF wrt temperature C DQNR - derivative of QREF wrt electron density C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' INCLUDE 'MODELQ.FOR' CHARACTER*4 TYPAT,DYP COMMON/PFSTDS/PFSTD(matom,30),MODPF(matom) common/terden/rhoter,anta,entrp DIMENSION TYPAT(matom),ABND(matom),D(3,matom), * abnref(mdepth),DYP(matom), * xio(8,matom), * abun0(matom),abun1(matom) dimension xio2(9,22), xio3(9,13) dimension ffi(matom),pfstu(matom),pfstt(matom),pfstn(matom), * entot(matom) dimension idat(30),uu(30,17) dimension u6(6),u7(7),u8(8),u10(10),u11(11) dimension u12(12),u13(13),u14(14),u16(16),u18(18),u20(20) dimension u24(24),u25(25),u26(26),u28(28) PARAMETER (TH0=5.0404D3,XMX0=2.154D4,THL0=2.3025851,FI0=3.6113D1, * TRHA=1.5D0,c1qm=1.0353d-16,c2qm=8762.9, * ev2erg=1.6018d-12) character*80 dum C equivalence (u6(1),uu(1,3)),(u7(1),uu(1,4)),(u8(1),uu(1,5)) equivalence (u10(1),uu(1,6)),(u11(1),uu(1,7)),(u12(1),uu(1,8)) equivalence (u13(1),uu(1,9)),(u14(1),uu(1,10)),(u16(1),uu(1,11)) equivalence (u18(1),uu(1,12)),(u20(1),uu(1,13)),(u24(1),uu(1,14)) equivalence (u25(1),uu(1,15)),(u26(1),uu(1,16)),(u28(1),uu(1,17)) C DATA DYP/' H ',' He ',' Li ',' Be ',' B ',' C ', * ' N ',' O ',' F ',' Ne ',' Na ',' Mg ', * ' Al ',' Si ',' P ',' S ',' Cl ',' Ar ', * ' K ',' Ca ',' Sc ',' Ti ',' V ',' Cr ', * ' Mn ',' Fe ',' Co ',' Ni ',' Cu ',' Zn ', * ' Ga ',' Ge ',' As ',' Se ',' Br ',' Kr ', * ' Rb ',' Sr ',' Y ',' Zr ',' Nb ',' Mo ', * ' Tc ',' Ru ',' Rh ',' Pd ',' Ag ',' Cd ', * ' In ',' Sn ',' Sb ',' Te ',' I ',' Xe ', * ' Cs ',' Ba ',' La ',' Ce ',' Pr ',' Nd ', * ' Pm ',' Sm ',' Eu ',' Gd ',' Tb ',' Dy ', * ' Ho ',' Er ',' Tm ',' Yb ',' Lu ',' Hf ', * ' Ta ',' W ',' Re ',' Os ',' Ir ',' Pt ', * ' Au ',' Hg ',' Tl ',' Pb ',' Bi ',' Po ', * ' At ',' Rn ',' Fr ',' Ra ',' Ac ',' Th ', * ' Pa ',' U ',' Np ',' Pu ',' Am ',' Cm ', * ' Bk ',' Cf ',' Es '/ C C Standard atomic constants for first 99 species C Abundances for the first 30 from Grevesse & Sauval, C (1998, Space Sci. Rev. 85, 161) C C Element Atomic Solar Std. C weight abundance highest C C ionization stage DATA D/ 1.008, 1.0D0, 2., * 4.003, 1.00D-1, 3., * 6.941, 1.26D-11, 3., * 9.012, 2.51D-11, 3., * 10.810, 5.0D-10, 4., * 12.011, 3.31D-4, 5., * 14.007, 8.32D-5, 5., * 16.000, 6.76D-4, 5., * 18.918, 3.16D-8, 4., * 20.179, 1.20D-4, 4., * 22.990, 2.14D-6, 4., * 24.305, 3.80D-5, 4., * 26.982, 2.95D-6, 4., * 28.086, 3.55D-5, 5., * 30.974, 2.82D-7, 5., * 32.060, 2.14D-5, 5., * 35.453, 3.16D-7, 5., * 39.948, 2.52D-6, 5., * 39.098, 1.32D-7, 5., * 40.080, 2.29D-6, 5., * 44.956, 1.48D-9, 5., * 47.900, 1.05D-7, 5., * 50.941, 1.00D-8, 5., * 51.996, 4.68D-7, 5., * 54.938, 2.45D-7, 5., * 55.847, 3.16D-5, 5., * 58.933, 8.32D-8, 5., * 58.700, 1.78D-6, 5., * 63.546, 1.62D-8, 5., * 65.380, 3.98D-8, 5., * 69.72 , 1.34896324e-09 , 3., * 72.60 , 4.26579633e-09 , 3., * 74.92 , 2.34422821e-10 , 3., * 78.96 , 2.23872066e-09 , 3., * 79.91 , 4.26579633e-10 , 3., * 83.80 , 1.69824373e-09 , 3., * 85.48 , 2.51188699e-10 , 3., * 87.63 , 8.51138173e-10 , 3., * 88.91 , 1.65958702e-10 , 3., * 91.22 , 4.07380181e-10 , 3., * 92.91 , 2.51188630e-11 , 3., * 95.95 , 9.12010923e-11 , 3., * 99.00 , 1.00000000e-24 , 3., * 101.1 , 6.60693531e-11 , 3., * 102.9 , 1.23026887e-11 , 3., * 106.4 , 5.01187291e-11 , 3., * 107.9 , 1.73780087e-11 , 3., * 112.4 , 5.75439927e-11 , 3., * 114.8 , 6.60693440e-12 , 3., * 118.7 , 1.38038460e-10 , 3., * 121.8 , 1.09647810e-11 , 3., * 127.6 , 1.73780087e-10 , 3., * 126.9 , 3.23593651e-11 , 3., * 131.3 , 1.69824373e-10 , 3., * 132.9 , 1.31825676e-11 , 3., * 137.4 , 1.62181025e-10 , 3., * 138.9 , 1.58489337e-11 , 3., * 140.1 , 4.07380293e-11 , 3., * 140.9 , 6.02559549e-12 , 3., * 144.3 , 2.95120943e-11 , 3., * 147.0 , 1.00000000e-24 , 3., * 150.4 , 9.33254366e-12 , 3., * 152.0 , 3.46736869e-12 , 3., * 157.3 , 1.17489770e-11 , 3., * 158.9 , 2.13796216e-12 , 3., * 162.5 , 1.41253747e-11 , 3., * 164.9 , 3.16227767e-12 , 3., * 167.3 , 8.91250917e-12 , 3., * 168.9 , 1.34896287e-12 , 3., * 173.0 , 8.91250917e-12 , 3., * 175.0 , 1.31825674e-12 , 3., * 178.5 , 5.37031822e-12 , 3., * 181.0 , 1.34896287e-12 , 3., * 183.9 , 4.78630102e-12 , 3., * 186.3 , 1.86208719e-12 , 3., * 190.2 , 2.39883290e-11 , 3., * 192.2 , 2.34422885e-11 , 3., * 195.1 , 4.78630036e-11 , 3., * 197.0 , 6.76082952e-12 , 3., * 200.6 , 1.23026887e-11 , 3., * 204.4 , 6.60693440e-12 , 3., * 207.2 , 1.12201834e-10 , 3., * 209.0 , 5.12861361e-12 , 3., * 210.0 , 1.00000000e-24 , 3., * 211.0 , 1.00000000e-24 , 3., * 222.0 , 1.00000000e-24 , 3., * 223.0 , 1.00000000e-24 , 3., * 226.1 , 1.00000000e-24 , 3., * 227.1 , 1.00000000e-24 , 3., * 232.0 , 1.20226443e-12 , 3., * 231.0 , 1.00000000e-24 , 3., * 238.0 , 3.23593651e-13 , 3., * 237.0 , 1.00000000e-24 , 3., * 244.0 , 1.00000000e-24 , 3., * 243.0 , 1.00000000e-24 , 3., * 247.0 , 1.00000000e-24 , 3., * 247.0 , 1.00000000e-24 , 3., * 251.0 , 1.00000000e-24 , 3., * 254.0 , 1.00000000e-24 , 3./ c data abun0 / * 12.00,10.93, 1.05, 1.38, 2.70, 8.39, 7.78, 8.66, 4.56, 7.84, * 6.17, 7.53, 6.37, 7.51, 5.36, 7.14, 5.50, 6.18, 5.08, 6.31, * 3.05, 4.90, 4.00, 5.64, 5.39, 7.45, 4.92, 6.23, 4.21, 4.60, * 2.88, 3.58, 2.29, 3.33, 2.56, 3.28, 2.60, 2.92, 2.21, 2.59, * 1.42, 1.92,-9.99, 1.84, 1.12, 1.69, 0.94, 1.77, 1.60, 2.00, * 1.00, 2.19, 1.51, 2.27, 1.07, 2.17, 1.13, 1.58, 0.71, 1.45, * -9.99, 1.01, 0.52, 1.12, 0.28, 1.14, 0.51, 0.93, 0.00, 1.08, * 0.06, 0.88,-0.17, 1.11, 0.23, 1.45, 1.38, 1.64, 1.01, 1.13, * 0.90, 2.00, 0.65,-9.99,-9.99,-9.99,-9.99,-9.99, 9.99, 0.06, * -9.99,-0.52,-9.99,-9.99,-9.99,-9.99,-9.99,-9.99,-9.99/ c data abun1 / * 12.00,10.93, 3.26, 1.38, 2.79, 8.43, 7.83, 8.69, 4.56, 7.93, * 6.24, 7.60, 6.45, 7.51, 5.41, 7.12, 5.50, 6.40, 5.08, 6.34, * 3.15, 4.95, 3.93, 5.64, 5.43, 7.50, 4.99, 6.22, 4.19, 4.56, * 3.04, 3.65, 2.30, 3.34, 2.54, 3.25, 2.36, 2.87, 2.21, 2.58, * 1.46, 1.88,-9.99, 1.75, 1.06, 1.65, 1.20, 1.71, 0.76, 2.04, * 1.01, 2.18, 1.55, 2.24, 1.08, 2.18, 1.10, 1.58, 0.72, 1.42, * -9.99, 0.96, 0.52, 1.07, 0.30, 1.10, 0.48, 0.92, 0.10, 0.92, * 0.10, 0.85,-0.12, 0.65, 0.26, 1.40, 1.38, 1.62, 0.80, 1.17, * 0.77, 2.04, 0.65,-9.99,-9.99,-9.99,-9.99,-9.99,-9.99, 0.06, * -9.99,-0.54,-9.99,-9.99,-9.99,-9.99,-9.99,-9.99,-9.99/ C C C Ionization potentials for first 99 species: DATA XIo/ C C Element Ionization potentials (eV) C I II III IV V VI VII VIII C * 13.595, 0. , 0. , 0. , 0. , 0. , 0. , 0. , * 24.580, 54.400, 0. , 0. , 0. , 0. , 0. , 0. , * 5.392, 75.619,122.451, 0. , 0. , 0. , 0. , 0. , * 9.322, 18.206,153.850,217.713, 0. , 0. , 0. , 0. , * 8.296, 25.149, 37.920,259.298,340.22, 0. , 0. , 0. , * 11.264, 24.376, 47.864, 64.476,391.99,489.98, 0. , 0. , * 14.530, 29.593, 47.426, 77.450, 97.86,551.93,667.03, 0. , * 13.614, 35.108, 54.886, 77.394,113.87,138.08,739.11,871.39, * 17.418, 34.980, 62.646, 87.140,114.21,157.12,185.14,953.6 , * 21.559, 41.070, 63.500, 97.020,126.30,157.91,207.21,239.0 , * 5.138, 47.290, 71.650, 98.880,138.37,172.09,208.44,264.16, * 7.664, 15.030, 80.120,102.290,141.23,186.49,224.9 ,265.96, * 5.984, 18.823, 28.440,119.960,153.77,190.42,241.38,284.53, * 8.151, 16.350, 33.460, 45.140,166.73,205.11,246.41,303.07, * 10.484, 19.720, 30.156, 51.354, 65.01,220.41,263.31,309.26, * 10.357, 23.400, 35.000, 47.290, 72.50, 88.03,280.99,328.8 , * 12.970, 23.800, 39.900, 53.500, 67.80, 96.7 ,114.27,348.3 , * 15.755, 27.620, 40.900, 59.790, 75.00, 91.3 ,124.0 ,143.46, * 4.339, 31.810, 46.000, 60.900, 82.6 , 99.7 ,118.0 ,155.0 , * 6.111, 11.870, 51.210, 67.700, 84.39,109.0 ,128.0 ,147.0 , * 6.560, 12.890, 24.750, 73.900, 92.0 ,111.1 ,138.0 ,158.7 , * 6.830, 13.630, 28.140, 43.240, 99.8 ,120.0 ,140.8 ,168.5 , * 6.740, 14.200, 29.700, 48.000, 65.2 ,128.9 ,151.0 ,173.7 , * 6.763, 16.490, 30.950, 49.600, 73.0 , 90.6 ,161.1 ,184.7 , * 7.432, 15.640, 33.690, 53.000, 76.0 , 97.0 ,119.24,196.46, * 7.870, 16.183, 30.652, 54.800, 75.0 , 99.1 ,125.0 ,151.06, * 7.860, 17.060, 33.490, 51.300, 79.5 ,102.0 ,129.0 ,157.0 , * 7.635, 18.168, 35.170, 54.900, 75.5 ,108.0 ,133.0 ,162.0 , * 7.726, 20.292, 36.830, 55.200, 79.9 ,103.0 ,139.0 ,166.0 , * 9.394, 17.964, 39.722, 59.400, 82.6 ,108.0 ,134.0 ,174.0 , * 6.000, 20.509, 30.700, 99.99,99.99,99.99,99.99,99.99, * 7.89944,15.93462, 34.058, 45.715,99.99,99.99,99.99,99.99, * 9.7887, 18.5892, 28.351, 99.99,99.99,99.99,99.99,99.99, * 9.750,21.500, 32.000, 99.99,99.99,99.99,99.99,99.99, * 11.839,21.600, 35.900, 99.99,99.99,99.99,99.99,99.99, * 13.995,24.559, 36.900, 99.99,99.99,99.99,99.99,99.99, * 4.175,27.500, 40.000, 99.99,99.99,99.99,99.99,99.99, * 5.692,11.026, 43.000, 99.99,99.99,99.99,99.99,99.99, * 6.2171,12.2236, 20.5244,60.607,99.99,99.99,99.99,99.99, * 6.63390,13.13,23.17,34.418,80.348,99.99,99.99,99.99, * 6.879,14.319, 25.039, 99.99,99.99,99.99,99.99,99.99, * 7.099,16.149, 27.149, 99.99,99.99,99.99,99.99,99.99, * 7.280,15.259, 30.000, 99.99,99.99,99.99,99.99,99.99, * 7.364,16.759, 28.460, 99.99,99.99,99.99,99.99,99.99, * 7.460,18.070, 31.049, 99.99,99.99,99.99,99.99,99.99, * 8.329,19.419, 32.920, 99.99,99.99,99.99,99.99,99.99, * 7.574,21.480, 34.819, 99.99,99.99,99.99,99.99,99.99, * 8.990,16.903, 37.470, 99.99,99.99,99.99,99.99,99.99, * 5.784,18.860, 28.029, 99.99,99.99,99.99,99.99,99.99, * 7.342,14.627, 30.490,72.3,99.99,99.99,99.99,99.99, * 8.639,16.500, 25.299,44.2,55.7,99.99,99.99,99.99, * 9.0096,18.600, 27.96, 37.4,58.7,99.99,99.99,99.99, * 10.454,19.090, 32.000, 99.99,99.99,99.99,99.99,99.99, * 12.12984,20.975,31.05,45.,54.14,99.99,99.99,99.99, * 3.893,25.100, 35.000, 99.99,99.99,99.99,99.99,99.99, * 5.210,10.000, 37.000, 99.99,99.99,99.99,99.99,99.99, * 5.580,11.060, 19.169, 99.99,99.99,99.99,99.99,99.99, * 5.650,10.850, 20.080, 99.99,99.99,99.99,99.99,99.99, * 5.419,10.550, 23.200, 99.99,99.99,99.99,99.99,99.99, * 5.490,10.730, 20.000, 99.99,99.99,99.99,99.99,99.99, * 5.550,10.899, 20.000, 99.99,99.99,99.99,99.99,99.99, * 5.629,11.069, 20.000, 99.99,99.99,99.99,99.99,99.99, * 5.680,11.250, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.159,12.100, 20.000, 99.99,99.99,99.99,99.99,99.99, * 5.849,11.519, 20.000, 99.99,99.99,99.99,99.99,99.99, * 5.930,11.670, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.020,11.800, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.099,11.930, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.180,12.050, 23.700, 99.99,99.99,99.99,99.99,99.99, * 6.250,12.170, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.099,13.899, 19.000, 99.99,99.99,99.99,99.99,99.99, * 7.000,14.899, 23.299, 99.99,99.99,99.99,99.99,99.99, * 7.879,16.200, 24.000, 99.99,99.99,99.99,99.99,99.99, * 7.86404,17.700, 25.000, 99.99,99.99,99.99,99.99,99.99, * 7.870,16.600, 26.000, 99.99,99.99,99.99,99.99,99.99, * 8.500,17.000, 27.000, 99.99,99.99,99.99,99.99,99.99, * 9.100,20.000, 28.000, 99.99,99.99,99.99,99.99,99.99, * 8.95868,18.563,33.227, 99.99,99.99,99.99,99.99,99.99, * 9.220,20.500, 30.000, 99.99,99.99,99.99,99.99,99.99, * 10.430,18.750, 34.200, 99.99,99.99,99.99,99.99,99.99, * 6.10829,20.4283,29.852,50.72,99.99,99.99,99.99,99.99, * 7.416684,15.0325,31.9373,42.33,69.,99.99,99.99,99.99, * 7.285519,16.679, 25.563,45.32,56.0,88.,99.99,99.99, * 8.430,19.000, 27.000, 99.99,99.99,99.99,99.99,99.99, * 9.300,20.000, 29.000, 99.99,99.99,99.99,99.99,99.99, * 10.745,20.000, 30.000, 99.99,99.99,99.99,99.99,99.99, * 4.000,22.000, 33.000, 99.99,99.99,99.99,99.99,99.99, * 5.276,10.144, 34.000, 99.99,99.99,99.99,99.99,99.99, * 6.900,12.100, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99, * 6.000,12.000, 20.000, 99.99,99.99,99.99,99.99,99.99/ C c additional ionization potentials c ix x xi xii xiii xiv xv xvi xvii c c energies added for Sc, Ti, V, Cr, Co,Ni, Cy,Zn c source G. C. Rodrigues et al. , c Systematic calculation of total atomic energies of ground state configurations, c Atomic Data and Nuclear Data Tables, Vol. 86, Issue 2, March 2004, P.117-233. data xio2/ * 1103., 0., 0., 0., 0., 0., 0., 0., 0., ! F * 1196., 1362., 0., 0., 0., 0., 0., 0., 0., ! Ne * 300., 1465.,1649., 0., 0., 0., 0., 0., 0., ! Na * 328., 367.,1762.,1963., 0., 0., 0., 0., 0., ! Mg * 330., 398., 442.,2085.,2304., 0., 0., 0., 0., ! Al * 351., 401., 476., 523.,2438.,2673., 0., 0., 0., ! Si * 372., 424., 480., 560., 612.,2816.,3069., 0., 0., ! P * 379., 447., 505., 565., 652., 707.,3223.,3494., 0., ! S * 400., 456., 529., 592., 657., 750., 809.,3658.,3946., ! Cl * 422., 479., 539., 618., 686., 756., 854., 918.,4121., ! Ar * 176., 503., 564., 629., 714., 787., 862., 968.,1034., ! K * 188., 211., 591., 656., 726., 817., 895., 974.,1087., ! Ca * 180., 225., 250., 686., 755., 830., 926.,1010.,1094., ! Sc * 193., 216., 265., 291., 787., 861., 940.,1042.,1132., ! Ti * 206., 230., 255., 308., 336., 896., 974.,1060.,1165., ! V * 209., 244., 271., 298., 355., 384.,1010.,1095.,1185., ! Cr * 222., 248., 286., 314., 344., 404., 435.,1136.,1222., ! Mn * 235., 262., 290., 331., 361., 392., 457., 490.,1266., ! Fe * 186., 276., 305., 336., 379., 411., 444., 512., 547., ! Co * 193., 224., 321., 352., 384., 430., 464., 499., 571., ! Ni * 199., 232., 266., 369., 401., 435., 484., 520., 557., ! Cu * 203., 238., 274., 311., 420., 454., 490., 542., 579./ ! Zn c c even higher ionization potentials c 18 19 20 21 22 23 24 25 26 c data xio3/ * 4426., 0., 0., 0., 0., 0., 0., 0., 0., ! Ar * 4611., 4934., 0., 0., 0., 0., 0., 0., 0., ! K * 1158., 5129.,5470., 0., 0., 0., 0., 0., 0., ! Ca * 1206., 1288.,5675.,6034., 0., 0., 0., 0., 0., ! Sc * 1222., 1346.,1425.,6249., 0., 0., 0., 0., 0., ! Ti * 1261., 1356.,1480.,1569., 0., 0., 0., 0., 0., ! V * 1294., 1397.,1497.,1627., 0., 0., 0., 0., 0., ! Cr * 1318., 1431.,1540.,1645.,1782., 0., 0., 0., 0., ! Mn * 1357., 1459.,1574.,1689.,1799.,1958.,2346.,8828.,9278., ! Fe * 1396., 1496.,1603.,1723.,1847.,1963.,2112., 0., 0., ! Co * 606., 1538.,1643.,1756.,1880.,2011.,2133.,2288., 0., ! Ni * 628., 670.,1688.,1797.,1915.,2043.,2183.,2310.,2472., ! Cu * 616., 693., 737.,1844.,1958.,2082.,2214.,2362.,2494./ ! Zn C C C data for additional ionization potentials for the Opacity C project species (IDAT sets the internal OP indexing) C data idat / 1, 2, 0, 0, 0, 3, 4, 5, 0, 6, * 7, 8, 9,10, 0,11, 0,12, 0,13, * 0, 0, 0,14,15,16, 0,17,0,0/ data u10/173.93,330.391,511.8,783.3,1018.,1273.8,1671.792, * 1928.462,9645.005,10986.876/ data u11/41.449,381.395,577.8,797.8,1116.2,1388.5,1681.5,2130.8, * 2418.7,11817.061,13297.676/ data u12/61.671,121.268,646.41,881.1,1139.4,1504.3,1814.3,2144.7, * 2645.2,2964.4,14210.261,15829.951/ data u13/48.278,151.86,229.446,967.8,1239.8,1536.3,1947.3,2295.4, * 2663.4,3214.8,3565.6,16825.022,18584.138/ data u14/65.748,131.838,270.139,364.093,1345.1,1653.9,1988.4, * 2445.3,2831.9,3237.8,3839.8,4222.4,19661.693,21560.63/ data u16/83.558,188.2,280.9,381.541,586.2,710.184,2265.9,2647.4, * 3057.7,3606.1,4071.4,4554.3,5255.9,5703.6,26002.663, * 28182.535/ data u18/127.11,222.848,328.6,482.4,605.1,734.04,1002.73,1157.08, * 3407.3,3860.9,4347.,4986.6,5533.8,6095.5,6894.2,7404.4, * 33237.173,35699.936/ data u20/49.306,95.752,410.642,542.6,681.6,877.4,1026.,1187.6, * 1520.64,1704.047,4774.,5301.,5861.,6595.,7215.,7860., * 8770.,9338.,41366.,44177.41/ data u24/54.576,132.966,249.7,396.5,560.2,731.02,1291.9,1490., * 1688.,1971.,2184.,2404.,2862.,3098.52,8151.,8850., * 9560.,10480.,11260.,12070.,13180.,13882.,60344.,63675.9/ data u25/59.959,126.145,271.55,413.,584.,771.1,961.44,1569., * 1789.,2003.,2307.,2536.,2771.,3250.,3509.82,9152., * 9872.,10620.,11590.,12410.,13260.,14420.,15162., * 65660.,69137.4/ data u26/63.737,130.563,247.22,442.,605.,799.,1008.,1218.38, * 1884.,2114.,2341.,2668.,2912.,3163.,3686.,3946.82, * 10180.,10985.,11850.,12708.,13620.,14510.,15797., * 16500.,71203.,74829.6/ data u28/61.6,146.542,283.8,443.,613.5,870.,1070.,1310.,1560., * 1812.,2589.,2840.,3100.,3470.,3740.,4020.,4606., * 4896.2,12430.,13290.,14160.,15280.,16220.,17190., * 18510.,19351.,82984.,86909.4/ C IF(MODE.NE.0) GO TO 50 C C For MODE=0, STATE serves as an auxiliary procedure for START C Input of basic parameters for individual chemical species C C An element (hydrogen through zinc) can be considered in one of C the three following options: C 1. explicitly - some of energy levels of some of its ionization C states are considered explicitly, ie. their C populations are determined by solving statistical C equilibrium C 2. implicitly - the atom is assumed not to contribute to C opacity; but is allowed to contribute to the C total number of particles and to the total charge; C the latter is evaluated assuming LTE ionization C balance, ie. by solving a set of Saha equations C 3. not considered at all C C Input: C NATOMS - the highest atomic number of an element that is C considered (explicitly or non-explicitly) C < 0 - then NATOMS=-NATOMS, and all Opacity Project C species are treated as with MODPF>0, i.e. C their partition functions are evaluated from C the Opacity Project ionization fractions C C For each element from 1 (hydrogen) to NATOMS, the following C parameters: C C MA = 0 - if the element is not considered (option 3) C = 1 - if the element is non-explicit (option 2) C = 2 - if the element is explicit (option 1) C ABN - if ABN=0, solar abundance is assumed (given above; C abundance here is assumed as relative C to hydrogen by number C if ABN>0, non-solar abundance ABN is assumed; in an C arbitrary scale C if ABN<0, non-solar abundance ABN is assumed; C (-ABN times the solar value) C if ABN>1.e6, depth-dpendent abundance (additional input) C iabset=0 read(ibuff,'(a80)') dum read(dum,*,iostat=kstat) natoms,iabset if(kstat.ne.0) READ(dum,*) NATOMS if(natoms.eq.0) then do id=1,nd ytot(id)=1.11 wmy(id)=1.41 wmm(id)=2.17e-24 end do iatref=0 return end if c do i=1,matom iadop(i)=0 if(iabset.eq.1) then d(2,i)=10.**(abun1(i)-12.) else if(iabset.ne.2) then d(2,i)=10.**(abun0(i)-12.) end if end do c WRITE(6,600) IAT=0 IREFA=0 IFOPPF=0 IF(NATOMS.LT.0) THEN NATOMS=-NATOMS IFOPPF=2 END IF C natms=natoms if(ifmol.gt.0) natms=92 DO 20 I=1,NATMS LGR(I)=.TRUE. LRM(I)=.TRUE. IATEX(I)=-1 if(i.le.natoms) then READ(IBUFF,*) MA,ABN,MODPF0 else ma=1 abn=0. modpf0=0 end if c IF(MA.EQ.0) GO TO 20 TYPAT(I)=DYP(I) AMAS(I)=D(1,I) ABND(I)=D(2,I) IF(MA.EQ.0.or.i.gt.natoms) then abnd(i)=d(2,i)*1.e-10 enev(i,1)=xio(1,i) enev(i,2)=xio(2,i) amas(i)=d(1,i) DO ID=1,ND ABNDD(I,ID)=ABND(I) END DO GO TO 20 end if IONIZ(I)=int(D(3,I)) MODPF(I)=MODPF0 C C increase the standard highest ionization for Teff larger C than 50000 K for N, O, Ne, and Fe C IF(TEFF.GT.5.D4) THEN IF(I.EQ.7) IONIZ(I)=6 IF(I.EQ.8) IONIZ(I)=7 IF(I.EQ.10) IONIZ(I)=9 IF(I.EQ.26) IONIZ(I)=9 END IF C DO J=1,8 ENEV(I,J)=xio(J,I) END DO c if(i.ge.9.and.i.le.30) then do j=9,17 enev(i,j)=xio2(j-8,i-8) end do endif if(i.ge.18.and.i.le.30) then do j=18,26 enev(i,j)=xio3(j-17,i-17) end do end if c LGR(I)=.FALSE. IF(MODPF(I).GT.0) IFOPPF=1 IF(IFOPPF.EQ.2.and.idat(i).gt.0) MODPF(I)=1 if(modpf(i).gt.0) then if(idat(i).eq.0) modpf(i)=0 if(modpf(i).gt.0) then ioniz(i)=i+1 if(i.ge.10) then do j=9,i enev(i,j)=uu(j,idat(i))*0.1239529d0 end do end if end if end if IF(ABN.GT.0) ABND(I)=ABN IF(ABN.LT.0) ABND(I)=ABS(ABN)*D(2,I) IF(ION.NE.0) IONIZ(I)=ION IF(ABN.GT.1.E6) THEN READ(IBUFF,*) (ABNDD(I,ID),ID=1,ND) ELSE DO ID=1,ND ABNDD(I,ID)=ABND(I) END DO END IF IF(MA.EQ.1) THEN LRM(I)=.FALSE. IATEX(I)=0 ELSE IAT=IAT+1 IIFIX(IAT)=0 IF(MA.LE.-2) IIFIX(IAT)=1 IATEX(I)=IAT IF(IAT.EQ.IATREF) THEN IREFA=I DO ID=1,ND ABNREF(ID)=ABNDD(I,ID) END DO END IF C C store parameters for explicit atoms C AMASS(IAT)=AMAS(I)*HMASS NUMAT(IAT)=I if(iabs(ma).eq.3) iadop(i)=1 END IF 20 CONTINUE C C renormalize abundances to have the standard element abundance C equal to unity C if(abnref(1).eq.0.or.ioptab.gt.0) then do id=1,nd abnref(id)=1. end do end if c DO 30 I=1,NATOMS IAT=IATEX(I) IF(IAT.LT.0) GO TO 30 DO ID=1,ND ABNDD(I,ID)=ABNDD(I,ID)/ABNREF(ID) YTOT(ID)=YTOT(ID)+ABNDD(I,ID) WMY(ID)=WMY(ID)+ABNDD(I,ID)*AMAS(I) END DO ABNR=ABND(I)/D(2,I) IF(IAT.EQ.0) THEN WRITE(6,601) I,TYPAT(I),ABND(I),ABNR ELSE DO ID=1,ND ABUND(IAT,ID)=ABNDD(I,ID) END DO WRITE(6,602) I,TYPAT(I),ABND(I),ABNR,IAT END IF 30 CONTINUE DO ID=1,ND WMM(ID)=WMY(ID)*HMASS/YTOT(ID) END DO c c initialization of the Opacity Project ionization fractions c (if required) c if(ifoppf.gt.0) call opfrac(0,0,t,ane,pf,fra) c 600 FORMAT(/' '/' CHEMICAL ELEMENTS INCLUDED'/ * ' --------------------------'// * ' NUMBER ELEMENT ABUNDANCE'/1H ,16X, * 'A=N(ELEM)/N(H) A/A(SOLAR)'/) 601 FORMAT(1H ,I4,3X,A5,1P2E14.2) 602 FORMAT(1H ,I4,3X,A5,1P2E14.2,3X, * 'EXPLICIT: IAT=',I3) RETURN C 50 TLN=LOG(T)*TRHA TK=BOLK*T TKLN15=TRHA*LOG(TK) ENTCON=103.973 THET=TH0/T THL=THL0*THET XMX=XMX0*SQRT(SQRT(T/ANE)) DCH=EH/(XMX*XMX*TK) Q=0. QM=0. DQT=0. DQN=0. DQM=0. ENER=0. ENTR=0. hpop=dens(id)/wmm(id)/ytot(id) DO 70 I=1,NATOMS IF(MODE.GT.1.AND.LRM(I).OR.MODE.EQ.1.AND.LGR(I)) GO TO 70 ION=IONIZ(I) DRQT=0. DRQN=0. DRST=0. DRSN=0. DFT=0. DFN=0. ENTOT(1)=0. RS=UN CALL PARTF(I,1,T,ANE,XMX,UM,DUTM,DUNM) if(i.eq.1) pfhyd=max(um,two) pfstu(1)=um pfstt(1)=dutm pfstn(1)=dunm JMAX=1 DO J=2,ION J1=J-1 DCHT=DCH*J1 TE=ENEV(I,J1)*THL ENTOT(J)=ENTOT(J1)+TE dcht=0. FI=FI0+TLN-TE+DCHT X=J XMAX=XMX*SQRT(X) CALL PARTF(I,J,T,ANE,XMAX,U,DUT,DUN) pfstu(j)=u pfstt(j)=dut pfstn(j)=dun FFI(J)=0. IF(FI.GT.-20.) FFI(J)=EXP(FI)*U/UM/ANE IF(FFI(J).GT.UN) JMAX=J UM=U END DO RQ=JMAX-1 RI=ENTOT(JMAX) RE=PFSTT(JMAX)/PFSTU(JMAX)*T if(jmax.lt.ion) then R=UN DO J=JMAX+1,ION J1=J-1 DCHT=DCH*J1 TE=ENEV(I,J1)*THL R=R*FFI(J) c RR(I,J)=R/pfstu(j) RR(I,J)=R RS=RS+R RQ=RQ+J1*R RI=RI+R*ENTOT(J) RE=RE+R*PFSTT(J)/PFSTU(J)*T DFIT=pfstt(j)/pfstu(j)-pfstt(j1)/pfstu(j1) . +(TRHA+TE-TRHA*DCHT)/T DFIN=pfstn(j)/pfstu(j)-pfstn(j1)/pfstu(j1) . +(HALF*DCHT-UN)/ANE DFT=DFT+DFIT DFN=DFN+DFIN DFIT=DFT*R DFIN=DFN*R DRST=DRST+DFIT DRSN=DRSN+DFIN DRQT=DRQT+J1*DFIT DRQN=DRQN+J1*DFIN END DO end if if(jmax.gt.1) then R=UN DFT=0. DFN=0. jmin=min(4,jmax-1) DO JJ=1,JMIN J=JMAX-JJ J1=J-1 JP1=J+1 DCHT=DCH*J TE=ENEV(I,J)*THL R=R/FFI(JP1) C RR(I,J)=R/pfstu(j) RR(I,J)=R RS=RS+R RQ=RQ+J1*R RI=RI+R*ENTOT(J) RE=RE+R*PFSTT(J)/PFSTU(J)*T DFIT=pfstt(jp1)/pfstu(jp1)-pfstt(j)/pfstu(j) * +(TRHA+TE-TRHA*DCHT)/T DFIN=pfstn(jp1)/pfstu(jp1)-pfstn(j)/pfstu(j) * +(HALF*DCHT-UN)/ANE DFT=DFT-DFIT DFN=DFN-DFIN DFIT=DFT*R DFIN=DFN*R DRST=DRST+DFIT DRSN=DRSN+DFIN DRQT=DRQT+J1*DFIT DRQN=DRQN+J1*DFIN END DO endif X=RQ/RS ABND(I)=ABNDD(I,ID) X1=ABND(I)/RS RR(I,JMAX)=X1 DO J=1,ION IF(J.NE.JMAX) RR(I,J)=RR(I,J)*X1 END DO c RR(I,JMAX)=RR(I,JMAX)/PFSTU(JMAX) c c internal energy and entropy (per 1 hydrogen atom) c chip=0 c antm=(anta-ane)/ytot(id) c do j=1,ion do j=1,2 dulog=0. c aden=rr(i,j)*antm if(aden.lt.1.e-20) aden=1.e-20 if(pfstu(j).lt.un) pfstu(j)=un if(pfstt(j).gt.0.) dulog=pfstt(j)*t/pfstu(j) ener=ener+(chip*ev2erg+tk*dulog)*aden entr=entr+(tkln15-log(aden)+log(pfstu(j))+ * trha*log(amas(i))+dulog+entcon)*aden c entr=entr+(tkln15+log(pfstu(j))+ c * trha*log(amas(i))+dulog+entcon)*aden chip=chip+enev(i,j) end do c aref=dens(id)/wmm(id)/ytot(id) c rr(i,99)=0. if(irefa.eq.0) irefa=1 IF(I.EQ.IREFA) THEN QREF=X*ABND(I) c DQTR=(DRQT-X*DRST)*X1 c DQNR=(DRQN-X*DRSN)*X1 ELSE Q=X*ABND(I)+Q c DQT=DQT+(DRQT-X*DRST)*X1 c DQN=DQN+(DRQN-X*DRSN)*X1 END IF do j=2,ion rr(i,99)=rr(i,99)+rr(i,j)*aref end do 70 CONTINUE c c c entropy of electrons c c entel=tkln15-log(ane)+1.5*log(emass(99))+entcon entel=tkln15-log(ane)-11.2622+entcon entr=entr+entel*ane C C Negative hydrogen ion C c IF(IHM.EQ.1) THEN tinv=un/t QM=C1QM*Tinv/SQRT(T)*EXP(C2QM*Tinv) DQM=-QM*Tinv*(TRHA+C2QM*Tinv) c END IF RETURN END