subroutine moleq(id,tt,an,aein,ane,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 - number density c aein - initial estimate of the electron density c c Output: ane - electron density c C Output through common/atomol: c rrr(id,j,i) - N/U for the atom with atomic number i and c ion j (j=1 for neutral, and j=2 for 1st ions) c rrmol(imol,id) - N/U for the molecule with index imol c (the index is given by the ordering of c in the input file tsuji.molec c c c Input data for molecules iven in the file c tsuji.molec c INCLUDE 'PARAMS.FOR' INCLUDE 'MODELP.FOR' character*128 MOLEC COMMON/COMFH1/C(600,5),PPMOL(600),APMLOG(600),P(100), * XIP(100),XI2(100),CCOMP(100),UIIDUI(100), * FP(100),XKP(100),XK2(100),EPS,SWITER, * NELEM(5,600),NATO(5,600),MMAX(600), * NELEMX(100),NMETAL,NIMAX common/moltst/pfmol(600,mdepth),anmol(600,mdepth), * pfato(100,mdepth),anato(100,mdepth), * pfion(100,mdepth),anion(100,mdepth) common/ioniz2/anion2(30,mdepth) DIMENSION NATOMM(5),NELEMM(5), * emass(100),uelem(100),ull(100),anden(800), * aelem(100) dimension denso(mdepth),eleco(mdepth),wmmo(mdepth) c data nmetal/92/ c data iread/1/ c MOLEC ='data/tsuji.molec_bc2' if(moltab.eq.0) MOLEC='data/tsuji.molec_orig' c ECONST=4.342945E-1 AVO=0.602217E+24 SPA=0.196E-01 GRA=0.275423E+05 AHE=0.100E+00 tk=1./(tt*1.38054e-16) pgas=an/tk sahcon=1.87840e20*tt*sqrt(tt) nimax=3000 eps=1.e-5 switer=0.0 C C---- data for atoms ---------------- C if(iread.eq.1) then c do i=1,nmetal ia=i nelemx(i)=ia ccomp(ia)=abndd(ia,id) xip(ia)=enev(ia,1) xi2(ia)=enev(ia,2) emass(ia)=amas(ia) end do c c---- read molecular data from a table ---------------------- c J=0 OPEN(UNIT=26,FILE=MOLEC,STATUS='OLD') 10 J=J+1 IF(MOLTAB.GE.1) * READ (26,510,end=20) CMOL(J),(C(J,K),K=1,5),MMAX(J), * (NELEMM(M),NATOMM(M),M=1,4) IF(MOLTAB.EQ.0) * READ (26,511,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) 511 FORMAT (A8,E11.5,4E12.5,I1,(I2,I3),3(I2,I2)) 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 c MMAXJ=MMAX(J) IF(MMAXJ.EQ.0) GO TO 20 DO M=1,MMAXJ NELEM(M,J)=NELEMM(M) NATO(M,J)=NATOMM(M) END DO c write(6,680) j,cmol(j) c 680 format(i5,a10) GO TO 10 20 NMOLEC=J-1 close(26) c DO I=1,NMETAL NELEMI=NELEMX(I) P(NELEMI)=1.D-70 END DO iread=0 endif c c---- end of reading atomic and molecular data ---------------------- c p(99)= aein/tk pesave=p(99) p(99)=pesave c THETA=5040./tt TEM=tt PGLOG=log10(Pgas) PG=Pgas 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)) c c----atoms----------------------------------------------------------------- c 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 irwpf(nelemi,1,0,tt,u0) uelem(nelemi)=u0 aelem(nelemi)=anden(i)/(u0*sahcon*emass(nelemi)**1.5) ull(nelemi)=log10(aelem(nelemi)) rrr(id,1,nelemi)=anden(i)/u0 anato(nelemi,id)=anden(i) pfato(nelemi,id)=u0 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 irwpf(nelemi,2,0,tt,u1) anion(nelemi,id)=anden(i+nmetal) pfion(nelemi,id)=u1 rrr(id,2,nelemi)=anden(i+nmetal)/u1 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 anion2(1,id)=0. c c---- molecules------------------------------------------------------------- c DO J=1,NMOLEC jm=j+2*nmetal PMOLL=log10(PPMOL(J)+1.0D-70) anden(jm)=exp(pmoll/econst)*tk rrmol(j,id)=0. umoll=1. if(pmoll.gt.-30.) then umoll=log10(anden(jm))+c(j,2)*theta 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 EXOMOL data whenever available c um=0. if(ipfexo.gt.0.and.tt.le.9000.) * call exopf(j,tt,um) if(um.gt.0.) then umoll=um else c c or with modified Irwin (Barklem & Collet) data whenever available c call irwpf(0,0,j,tt,um) if(um.gt.0.) umoll=um end if c H- c if(j.eq.1) umoll=1. c c set up array RRR = number density/partition function c rrmol(j,id)=anden(jm)/umoll end if c anmol(j,id)=anden(jm) pfmol(j,id)=umoll END DO jm=2*nmetal anhm(id)=anden(1+jm) anh2(id)=anden(2+jm) anch(id)=anden(5+jm) anoh(id)=anden(4+jm) C C C save new density, molecular weight, and abundances of c atomic species c ipri1=ipri denso(id)=dens(id) eleco(id)=elec(id) wmmo(id)=wmm(id) dens(id)=tmass*hmass elec(id)=pe*tk wmm(id)=dens(id)/(an-elec(id)) ane=elec(id) c do i=1,nmetal NELEMI=NELEMX(I) ia=iatex(nelemi) if(ia.gt.0) then attot(ia,id)=(anato(nelemi,id)+anion(nelemi,id)) end if end do c if(id.eq.nd) then write(86,610) do iid=1,nd write(86,611) iid,dm(iid),temp(iid),elec(iid),eleco(iid), * dens(iid),denso(iid),wmm(iid),wmmo(iid) end do end if 610 format(/' id m T ne(old) ne(new)', * ' dens(old) dens(new) wmm(old) wmm(new)'/) 611 format(i4,1p8e10.2) C RETURN END