subroutine fractn(iatnum) c ========================= c implicit double precision (a-h,o-z) parameter (mtemp=100, * melec= 60, * mion1=30, * mdat = 17) parameter (inp=71) dimension frac0(-1:mion1),ioo(-1:mion1),idat(mion1) dimension gg(mion1,mdat),g0(mion1),z0(-1:mion1) dimension uu(mion1,mdat),u0(mion1) 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) 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)) common/fracop/frac(mtemp,melec,mion1),fracm(mtemp,melec), * itemp(mtemp),ntt 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 gg/2.,29*0.,2.,1.,28*0., * 2.,1.,2.,1.,6.,9.,24*0.,2.,1.,2.,1.,6.,9.,4.,23*0., * 2.,1.,2.,1.,6.,9.,4.,9.,22*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,20*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,19*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1.,18*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1.,6.,17*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1.,6.,9.,16*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1.,6.,9.,4.,9.,14*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1.,6.,9.,4.,9.,6.,1., * 12*0.,2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1.,6.,9.,4.,9., * 6.,1.,2.,1.,10*0.,2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1., * 6.,9.,4.,9.,6.,1.,10.,21.,28.,25.,6.,7.,6*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1.,6.,9.,4.,9., * 6.,1.,10.,21.,28.,25.,6.,7.,6.,5*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1.,6.,9.,4.,9., * 6.,1.,10.,21.,28.,25.,6.,25.,30.,25.,4*0., * 2.,1.,2.,1.,6.,9.,4.,9.,6.,1.,2.,1.,6.,9.,4.,9., * 6.,1.,10.,21.,28.,25.,6.,25.,28.,21.,10.,21.,0.,0./ data uu(1,1)/109.6787/ data uu(1,2)/198.3108/ data uu(2,2)/438.9089/ data u6/90.82,196.665,386.241,520.178,3162.395,3952.061/ data u7/117.225,238.751,382.704,624.866,789.537,4452.758,5380.089/ data u8/109.837,283.24,443.086,624.384,918.657,1114.008,5963.135, * 7028.393/ 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(idat(iatnum).eq.0) then write(6,600) iatnum 600 format(' OP data for element no. ',i3,' do not exist') iatnum=-1 return end if c g0(iatnum+1)=1. do i=1,iatnum ig0=iatnum-i+1 g0(ig0)=gg(i,idat(iatnum)) u0(i)=uu(i,idat(iatnum))*1000. enddo c if(iatnum.eq.1) open(inp,file='ioniz.dat',status='old') do 10 it=1,mtemp do 10 ie=1,melec fracm(it,ie)=0. do 10 ion=1,mion1 frac(it,ie,ion)=0. 10 continue c read(inp,*) read(inp,*) it0,it1,itstp ntt=(it1-it0)/itstp+1 c do it=1,ntt read(inp,*) itt,ie0,ie1,iestp itemp(it)=itt net=(ie1-ie0)/iestp+1 t=exp(2.3025851*0.025*itt) safac0=sqrt(t)*t/2.07d-16 tkcm=0.69496*t do ie=1,net read(inp,601) iee,ion0,ion1, * (ioo(i),frac0(i),i=ion0,min(ion1,ion0+3)) ane=exp(2.3025851*0.25*iee) safac=safac0/ane nio=ion1-ion0 if(nio.ge.3) then nlin=nio/4 do ilin=1,nlin read(inp,602) (ioo(i),frac0(i), * i=ion0+4*ilin,min(ion1,ion0+4*ilin+3)) end do end if ieind=iee/2 do ion=ion0,ion1 if(ion.lt.iatnum) then if(ion.eq.ion0) then z0(ion)=g0(iatnum-ion) else z0(ion)=frac0(ion)/frac0(ion-1)*safac*z0(ion-1) z0(ion)=z0(ion)*exp(-u0(iatnum-ion)/tkcm) endif frac(it,ieind,iatnum-ion)=frac0(ion)/z0(ion) else u0hm=6090.5 z0hm=frac0(ion)/frac0(ion-1)*safac z0hm=z0hm*exp(-u0hm/tkcm) fracm(it,ieind)=frac0(ion)/z0hm end if end do end do end do 601 format(3i4,2x,4(i4,1x,e9.3)) 602 format(14x,4(i4,1x,e9.3)) return end