subroutine moleq(id,tt,an,aein,ane,energ,entt,wm,ipri) c ====================================================== c c calculation of the equilibrium state of atoms and molecules c c Input: id - depth point c tt - temperature [K] c an - total particle number density c aein - initial estimate of the electron density c c Output: 1) parameters transported through formal parameters" C ane - electron density c entt - entropy c energ - internal energy c 2) parameters transported through COMMON blocks c ELEC(ID) - electron denity c DENS(ID) - mass density c WMM(ID) - mean molecular weight c QADD(ID) - total charge of non-explicit species c ABUND(I,ID) - abundances of explicit atoms, c counting only atoms and ions c (not atoms sequestered in molecules) c c Input data for molecules given in file tsuji.molec c INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ATOMIC.FOR' c character*128 MOLEC COMMON/COMFH1/C(600,5),PPMOL(600),APMLOG(600), * XIP(100),XIP2(100),CCOMP(100),UIIDUI(100), * P(100),FP(100),XKP(100),XK2(100),EPS,SWITER, * NELEM(5,600),NATO(5,600),MMAX(600), * NELEMX(100),NMETAL,NMOLEC,NIMAX common/eospar/anmol(600,mdepth), * anato(100,mdepth), * anion(100,mdepth) common/ioniz2/anion2(30,mdepth) common/entrop/entato(100),ention(100),entmol(600) common/hmolab/anh2(mdepth),anhm(mdepth) common/terden/rhoter,anta,entrp common/adchar/qadd(mdepth) common/moldat/moltab,irwtab DIMENSION NATOMM(5),NELEMM(5), * emass(100),uelem(100),ull(100),anden(800), * aelem(100),ammol(600),cmol(600) dimension anat0(100),anio0(100),anmo0(600),pfmol(600) dimension denso(mdepth),eleco(mdepth),wmmo(mdepth) c data nmetal/92/ c data iread/1/ c if(ifmol.eq.0) return MOLEC ='data/tsuji.molec_bc2' if(moltab.eq.0) MOLEC='data/tsuji.molec_orig' c ECONST=4.342945E-1 tk=1./(tt*bolk) pgas=an/tk sahcon=1.87840e20*tt*sqrt(tt) ev2erg=1.6018e-12 nimax=3000 eps=0.001 switer=1 C if(iread.eq.1) then c C---- data for atoms ---------------- C do i=1,nmetal ia=i nelemx(i)=ia ccomp(ia)=abndd(ia,id) xip(ia)=enev(ia,1) xip2(ia)=enev(ia,2) emass(ia)=amas(ia) end do c c---- read molecular data from tsuji.molec ---------------------- c J=0 OPEN(UNIT=26,FILE=MOLEC,STATUS='OLD') 10 J=J+1 READ(26,510,END=20) CMOL(J),(C(J,K),K=1,5),MMAX(J), * (NELEMM(M),NATOMM(M),M=1,4) 510 FORMAT(A8,5E13.5,9I3) c c for now, exclude all molecules with 4 or more C atoms c do m=1,4 if(nelemm(m).eq.6.and.natomm(m).ge.5) then j=j-1 go to 10 end if end do MMAXJ=MMAX(J) IF(MMAXJ.GT.0) THEN DO M=1,MMAXJ NELEM(M,J)=NELEMM(M) NATO(M,J)=NATOMM(M) END DO GO TO 10 END IF 20 CONTINUE NMOLEC=J-1 close(26) c DO I=1,NMETAL NELEMI=NELEMX(I) P(NELEMI)=1.D-70 END DO iread=0 end if c c---- end of reading atomic and molecular data ---------------------- c p(99)= aein/tk pesave=p(99) p(99)=pesave c THET=5040./tt TEM=tt PGLOG=log10(Pgas) PG=Pgas tkln25=-log(tk)*2.5 tkln15=log(bolk*tt)*1.5 tkev=5040./tt c CALL RUSSEL(TEM,PG) c PE=P(99) ane=pe*tk PELOG=log10(PE) emass(99)=5.486e-4 uelem(99)=2. aelem(99)=pe*tk/(2.*sahcon*emass(nelemi)**1.5) ull(99)=log10(aelem(99)) idstd=nd*2/3 entcon=103.973 ann=an-ane tkk=bolk*tt*tt c c----atoms----------------------------------------------------------------- c entt=0 antt=0. energ=0. tmass=0. DO I=1,NMETAL NELEMI=NELEMX(I) FPLOG=log10(FP(NELEMI)) anden(i)=(p(nelemi)+1.d-70)*tk tmass=tmass+anden(i)*emass(nelemi) call mpartf(nelemi,1,0,tt,u0,dulog) uelem(nelemi)=u0 aelem(nelemi)=anden(i)/(u0*sahcon*emass(nelemi)**1.5) ull(nelemi)=log10(aelem(nelemi)) anat0(nelemi)=anden(i) if(dulog.lt.0.) dulog=0. entato(nelemi)=tkln15-log(anden(i))+log(u0)+ * 1.5*log(emass(nelemi))+entcon+tkk*dulog anx=anden(i)/ann antt=antt+anx entt=entt+entato(nelemi)*anden(i) energ=energ+dulog/tk*anden(i) END DO an1=anden(1) c c---- positive ions --------------------------------------------------------- c DO I=1,NMETAL NELEMI=NELEMX(I) PLOG= log10(P(NELEMI)+1.0D-70) XKPLOG=log10(XKP(NELEMI)+1.0D-70) PIONL=PLOG+XKPLOG-PELOG anden(i+nmetal)=exp(pionl/econst)*tk tmass=tmass+anden(i+nmetal)*emass(nelemi) call mpartf(nelemi,2,0,tt,u1,dulog) if(dulog.lt.0.) dulog=0. anio0(nelemi)=anden(i+nmetal) ention(nelemi)=tkln15-log(anden(i+nmetal))+log(u1)+ * 1.5*log(emass(nelemi))+entcon+tkk*dulog anx=anio0(nelemi)/ann antt=antt+anx entt=entt+ention(nelemi)*anio0(nelemi) energ=energ+(xip(nelemi)*1.6018e-12+ * dulog/tk)*anio0(nelemi) if(nelemi.ge.2.and.nelemi.le.30) then x2log=log10(XK2(NELEMI)+1.0D-70) pion2=pionl+x2log-pelog anion2(nelemi,id)=exp(pion2/econst)*tk end if END DO c c H- (incorrect in the original Tsuji's table) c jbeg=1 if(moltab.eq.0) then j=1 anmo0(1)=1.0353e-16/tt/sqrt(tt)*exp(8762.9/tt)* * anat0(1)*ane ammol(1)=emass(1) pfmol(1)=1. entmol(1)=tkln15-log(anmo0(1))+1.5*log(emass(1))+entcon anx=anmo0(1)/ann antt=antt+anx entt=entt+entmol(j)*anmo0(1) tmass=tmass+emass(1)*anmo0(1) jbeg=2 end if c c---- molecules------------------------------------------------------------- c DO J=jbeg,NMOLEC jm=j+2*nmetal PMOLL=log10(PPMOL(J)+1.0D-70) anden(jm)=exp(pmoll/econst)*tk if(pmoll.lt.-20.) go to 100 umoll=log10(anden(jm))+c(j,2)*thet amasm=0. do jjj=1,mmax(j) i=nelem(jjj,j) amasm=amasm+NATO(jjj,j)*emass(i) umoll=umoll-NATO(jjj,j)*ull(i) end do ammol(j)=amasm tmass=tmass+anden(jm)*amasm umoll=exp(umoll/econst)/(sahcon*amasm**1.5) c c replace with Irwin data whenever available c call mpartf(0,0,j,tt,um,dulog) if(um.gt.0.) umoll=um if(umoll.lt.1.) umoll=1. c anmo0(j)=anden(jm) pfmol(j)=umoll if(dulog.lt.0.) dulog=0. entmol(j)=tkln15-log(anden(jm))+log(umoll)+ * 1.5*log(amasm)+entcon+tkk*dulog anx=anden(jm)/ann antt=antt+anx entt=entt+entmol(j)*anden(jm) energ=energ+dulog/tk*anden(jm) if(j.eq.2) energ=energ-4.476*ev2erg*anden(jm) c if(ipri.gt.0.and.ipmole.gt.0.and.mod(id,10).eq.1.and.j.le.10) c * write(80,680) id,j,anmo0(j),anx,antt,pfmol(j), c * entmol(j),entt,dulog,dulog/tk*anden(jm),energ 100 continue END DO c c electrons c entel=tkln15-log(ane)+1.5*log(emass(99))+entcon entt=entt+entel*ane antt=antt+ane/ann c c final entropy, density, and mean mnolecular weight c entt=entt*bolk rhoter=tmass*hmass ahtot=anat0(1)+anio0(1)+anmo0(1)+2.*anmo0(2) wm=tmass/antt/ann c if(ipri.eq.0) return c do i=1,nmetal j=nelemx(i) anato(j,id)=anat0(j) anion(j,id)=anio0(j) end do c do j=1,nmolec anmol(j,id)=anmo0(j) end do c anhm(id)=anmol(1,id) anh2(id)=anmol(2,id) C C save new density, molecular weight, and abundances of c atomic species c denso(id)=dens(id) eleco(id)=elec(id) wmmo(id)=wmm(id) dens(id)=tmass*hmass elec(id)=ane wmm(id)=dens(id)/(an-elec(id)) qadd(id)=0. do i=1,nmetal NELEMI=NELEMX(I) ia=iatex(nelemi) if(ia.gt.0) then abund(ia,id)=(anato(nelemi,id)+anion(nelemi,id))* * wmm(id)/dens(id)*ytot(id) else qadd(id)=qadd(id)+anion(nelemi,id) end if end do if(ielhm.eq.0) qadd(id)=qadd(id)-anhm(id) c if(ipri.eq.0) return c c don't change structure if particle conservation in not solved c c IF(INPC.eq.0.and.ifryb.eq.0) THEN IF(INPC.eq.0) THEN dens(id)=denso(id) elec(id)=eleco(id) wmm(id)=wmmo(id) ENDIF c RETURN END