subroutine eospri c ================= c c Outprint of Equation of State parameters c INCLUDE 'PARAMS.FOR' INCLUDE 'MODELP.FOR' common/moltst/pfmol(600,mdepth),anmol(600,mdepth), * pfato(100,mdepth),anato(100,mdepth), * pfion(100,mdepth),anion(100,mdepth) common/hydmol/anhmi,ahmol common/hydato/ah,anh,anp common/ioniz2/anion2(30,mdepth) dimension nelemx(38) dimension amh2(5),xml(20),insm(20) data nelemx/ 1, 2, 3, 4, 5, 6, 7, 8, 9, * 11,12,13,14,15,16,17,19,20, * 21,22,23,24,25,26,28,29,32, * 35,37,38,39,40,41,53,56,57,58,60/ data amh2/1.13390E+01,-2.97499E+00,4.10842E-02,-3.58550E-03, * 1.31844E-04/ data insm/2,3,4,5,6,7,8,12,17,25,29,30,32,34,122,126,134, * 179,198,214/ data init/1/ c c id=idstd istp=1 if(ifeos.lt.0) istp=-ifeos c do id=1,nd,istp t=temp(id) ane=elec(id) rho=dens(id) ann = dens(id)/wmm(id)+elec(id) c if(ifmol.eq.0.or.t.gt.tmolim) then it=0 10 continue ann0=ann it=it+1 call eldens(id,t,ann,ane) anmol(1,id)=anhmi anmol(2,id)=ahmol anato(1,id)=anh anion(1,id)=anp hpop=dens(id)/wmy(id)/hmass do i=1,nmetal j=nelemx(i) anato(j,id)=anato(j,id)*hpop anion(j,id)=anion(j,id)*hpop if(j.ge.2.and.j.le.30) anion2(j,id)=anion2(j,id)*hpop end do anato(1,id)=anh anion(1,id)=anp c wmm(id)=(wmy(id)+2.*anmol(2,id)/hpop)/ytot(id)*hmass wmm(id)=wmy(id)/(ytot(id)-anmol(2,id)/hpop)*hmass ann=dens(id)/wmm(id)+ane if((ann-ann0)/ann0.gt.1.e-5) go to 10 end if c nmetal=38 write(*,*) '' write(*,*) 'atomic number densities and partition functions' write(*,*) '' atot=0. do i=1,nmetal j=nelemx(i) if(j.le.28) * write(6,621) j,typat(j),anato(j,id),pfato(j,id) atot=atot+anato(j,id) end do write(*,*) '' write(*,*) 'ionic number densities and partition functions' write(*,*) '' ctot=0. do i=1,nmetal j=nelemx(i) if(j.le.28) * write(6,622) j,typat(j),anion(j,id),pfion(j,id) atot=atot+anion(j,id) ctot=ctot+anion(j,id) end do 621 format(i4,a3,3x,1p2e12.4) 622 format(i4,a3,'+',2x,1p2e12.4) c if(ifmol.gt.0.and.t.le.tmolim) then write(6,600) do i=1,nmolec if(anmol(i,id).gt.ann*1.e-15) * write(6,601) i, cmol(i), anmol(i,id), pfmol(i,id) atot=atot+anmol(i,id) end do end if 600 format(/ 'Molecular number densities and partition functions'/) 601 format(i4,1x,A8,1x,1pe12.4,1x,e12.4) c ahmi=1.0353e-16/t/sqrt(t)*exp(8762.9/t)* * anato(1,id)*ane c c original B&C H2+ c APLOGJ=amh2(5) te=5040./t DO K=1,4 KM5=5-K APLOGJ=APLOGJ*TE + amh2(KM5) END DO tk=1.38054e-16*t ph2=-aplogj+log10(anato(1,id)*anion(1,id))+2.*log10(tk) anh2b=(10.**ph2)/tk htot=anato(1,id)+anion(1,id)+anmol(1,id)+ * 2.*(anmol(2,id)+anmol(3,id))+anmol(4,id)+anmol(5,id)+ * anmol(12,id)+2.*anmol(13,id)+anmol(14,id)+ * anmol(15,id)+ * anmol(16,id)+anmol(17,id)+anmol(32,id)+anmol(34,id)+ * 4.*anmol(37,id)+2.*anmol(38,id)+3.*anmol(39,id)+ * 2.*anmol(40,id)+3.*anmol(41,id)+2.*anmol(57,id)+ * anmol(118,id)+anmol(133,id)+ * 2.*anmol(140,id)+3.*anmol(141,id)+4.*anmol(142,id)+ * anmol(148,id)+2.*anmol(149,id)+anmol(222,id) ahe= (anato(2,id)+anion(2,id)+anion2(2,id))/htot aca= (anato(6,id)+anion(6,id)+anion2(6,id))/htot acm= (anmol(5,id)+anmol(6,id)+ * anmol(7,id)+2.*(anmol(8,id)+2.*anmol(13,id))+ * anmol(14,id)+2.*anmol(15,id)+anmol(20,id)+ * anmol(37,id)+anmol(38,id)+anmol(39,id)+ * anmol(44,id)+anmol(118,id)+anmol(119,id)+ * anmol(437,id)+anmol(453,id) * )/htot ana= (anato(7,id)+anion(7,id)+anion2(7,id))/htot anm= (anmol(7,id)+2.*anmol(9,id)+anmol(11,id)+ * anmol(12,id)+anmol(14,id)+anmol(23,id)+ * anmol(24,id)+anmol(40,id)+anmol(41,id)+ * anmol(109,id)+anmol(152,id)+anmol(347,id)+ * anmol(438,id)+anmol(452,id)+anmol(454,id) * )/htot aoa= (anato(8,id)+anion(8,id)+anion2(8,id))/htot aom= (anmol(3,id)+anmol(4,id)+ * anmol(6,id)+2.*anmol(10,id)+anmol(11,id)+anmol(25,id)+ * anmol(26,id)+anmol(29,id)+anmol(30,id)+anmol(31,id)+ * anmol(35,id)+2.*anmol(44,id)+anmol(49,id)+anmol(51,id)+ * anmol(54,id)+2.*anmol(56,id)+anmol(65,id)+ * 2.*anmol(66,id)+anmol(84,id)+anmol(109,id)+ * anmol(113,id)+anmol(115,id)+anmol(118,id)+ * anmol(119,id)+anmol(126,id)+anmol(134,id)+ * anmol(153,id)+anmol(179,id)+anmol(184,id)+ * 2.*anmol(185,id)+anmol(200,id)+anmol(216,id)+ * anmol(221,id)+2.*anmol(247,id)+anmol(292,id)+ * anmol(439,id)+anmol(453,id)+anmol(454,id) * )/htot ac=aca+acm an=ana+anm ao=aoa+aom write(6,623) t,dens(id),ann,atot+ane,ane,ctot-anmol(1,id), * anato(1,id),anion(1,id), * anmol(1,id),anmol(2,id), * anmol(312,id),anmol(426,id),anh2b, * htot, * anmol(1,id),ahmi,anmol(1,id)/ahmi, * anato(6,id),anion(6,id),anmol(6,id),anmol(37,id), * anato(7,id),anion(7,id),anmol(9,id),anmol(41,id), * anato(8,id),anion(8,id),anmol(3,id),anmol(6,id), * ahe,ahe/abndd(2,id), * ac,ac/abndd(6,id), * an,an/abndd(7,id), * ao,ao/abndd(8,id) act=ac*htot ant=an*htot aot=ao*htot 623 format(/'EOS useful quantities - summary'// * 'T,rho ',f13.2,1pe13.5/ * 'N ',1p2e13.5/ * 'n_e ',1p2e13.5/ * 'H,H+,H-,H2 ',1p4e13.5/ * 'H2-,H2+,H2+b',1p3e13.5/ * 'Htot ',1pe13.5/ * 'H- ',1p3e13.5/ * 'C,C+,CO,CH4 ',1p4e13.5/ * 'N,N+,N2,NH3 ',1p4e13.5/ * 'O,O+,H2O,CO ',1p4e13.5/ * 'He/H ',1p2e13.5/ * 'C/H ',1p2e13.5/ * 'N/H ',1p2e13.5/ * 'O/H ',1p2e13.5/) c if(init.eq.1) then write(52,625) write(51,626) write(53,653) (cmol(insm(i)),i=1,20) write(54,654) (cmol(insm(i)),i=1,20) c 625 format(' T rho w_mol Ne/Ntot N(Htot) ' * 'n(H) n(H2)',6x, * 'a(He) a(C) a(N) a(O) molfr(C) molfr(N) molfr(O)'/) c * 'a(He) a(C) a(N) a(O) n(C) n(CO) n(CH4)',5x, c * 'n(N) n(N2) n(NH3) n(O) n(H2O) n(CO)'/) init=0 end if c c write(51,624) t,dens(id),wmm(id)/hmass,ane/ann, c * htot,anato(1,id)/htot,2.*anmol(2,id)/htot, c * ahe/abndd(2,id),ac/abndd(6,id),an/abndd(7,id),ao/abndd(8,id), c * anato(6,id)/act,anmol(6,id)/act,anmol(37,id)/act, c * anato(7,id)/ant,2.*anmol(9,id)/ant,anmol(41,id)/ant, c * anato(8,id)/aot,anmol(3,id)/aot,anmol(6,id)/aot write(52,624) t,dens(id),wmm(id)/hmass,ane/ann, * htot,anato(1,id),2.*anmol(2,id), * ahe/abndd(2,id),ac/abndd(6,id),an/abndd(7,id),ao/abndd(8,id), * acm/ac,anm/an,aom/ao c * anato(6,id),anmol(6,id),anmol(37,id), c * anato(7,id),anmol(9,id),anmol(41,id), c * anato(8,id),anmol(3,id),anmol(6,id) 624 format(f8.1,1pe9.2,0pf8.5,1x,1p4e9.2,1x,0p4f8.5,1x,1p3e9.2,1x, * 3e9.2,1x,3e9.2) c write(51,627) t,dens(id),wmm(id)/hmass,ann,ane,htot, * anato(1,id),anion(1,id),anmol(1,id),anmol(2,id),anmol(312,id), * anmol(426,id) c * anmol(426,id),anh2b 626 format(' T rho w_mol N Ne N(Htot) ', * 'N(H) N(H+) N(H-) N(H2) N(H2-) N(H2+)'/) c * 'N(H) N(H+) N(H-) N(H2) N(H2-) N(H2+) N(H2+b)'/) 627 format(f8.1,1pe9.2,0pf8.5,1x,1p10e9.2) c if(ifmol.gt.0.and.t.le.tmolim) then do i=1,20 im=insm(i) xml(i)=log10(anmol(im,id)/pfmol(im,id)) end do write(53,655) t,log10(dens(id)),(xml(i),i=1,20) do i=1,20 im=insm(i) xml(i)=log10(anmol(im,id)/htot) c xml(i)=log10(anmol(im,id)) end do write(54,655) t,log10(dens(id)),(xml(i),i=1,20) end if c 653 format(' log10(N/U)'/' T rho ',20a6/) 654 format(' log10[N/n(H)]'/' T rho ',20a6/) 655 format(2f6.1,1x,20f6.1) c end do return end