subroutine eqlib c****************************************************************************** c This routine solves the equations of molecular dissociative c equilibrium. the equations can include neutral and ionized molecules c and atoms c****************************************************************************** implicit real*8 (a-h,o-z) include 'Atmos.com' include 'Mol.com' include 'Quants.com' include 'Dummy.com' real*8 xfic(30), xcorr(30), deltax(30), c(30,30), ans(30,30) real*8 uu(2) real*8 lth integer ident(30,110) c*****clear the arrays iorder(1) = 0 neq = iorder(1) tdel = (t(ntau)-t(1))/4.0 + 1.0 do k=1,30 iorder(k) = 0 xcorr(k) = 0. do kk=1,110 ident(k,kk) = 0 enddo enddo c*****the number of species to be considered has been defined in "Inmodel"; c*****either read in the dissociation data for a molecular species do jmol=1,nmol if (amol(jmol) >= 100.) then do k=1,110 if (datmol(1,k) == amol(jmol)) go to 11 enddo write (nf1out,1001) amol(jmol) stop 11 do kk=1,6 const(kk,jmol) = datmol(kk+1,k) enddo c*****or read the ionization data for an atomic species else iatom1 = amol(jmol) + 0.0001 atom = iatom1 const(1,jmol) = xchi1(iatom1) do kk=1,5 ti = t(1) + (kk-1)*tdel it = ti do jj=1,2 att = atom+0.1*(jj-1) if (partflag(iatom1,jj) > 0) then uu(jj) = partnew(att,jj,it) else uu(jj) = ucalc(att,it) endif enddo const(kk+1,jmol) = uu(2)/uu(1) enddo endif enddo if (molopt >= 2) . write (nf1out,1002) (amol(i),(const(j,i),j=1,6),i=1,nmol) c*****set up the molecular equilibrium array. each element of array c*****'ident' contains the identifier (subscript of an element of c*****array 'amol') of an ion or a molecule. the neutral atom is always c*****understood, but is not explicitly contained in 'ident'. nmax = 1 do jmol=1,nmol atom = amol(jmol) 2 call sunder(atom,iatom1,iatom2) do k=1,30 if (iatom1 == iorder(k)) go to 4 enddo neq = neq + 1 iorder(neq) = iatom1 ident(neq,1) = jmol go to 6 4 do kk=1,nmax if (ident(k,kk)==0 .or. ident(k,kk)== jmol) go to 7 enddo nmax = nmax + 1 kk = nmax 7 ident(k,kk) = jmol 6 if (iatom2 /= 0) then atom = iatom2 go to 2 endif enddo if (molopt >= 2) then do i=1,neq dummy1(i) = iorder(i) enddo write (nf1out,1003) (dummy1(i),i=1,neq),(amol(i),i=1,nmol) endif c*****now begin the loop that goes through all the atmosphere tau layers do 21 kev=1,ntau c*****calculate *xfic* and make a first guess at *xatom* i = ntau + 1 - kev lev = i tk = 1.38065d-16*t(i) do k=1,neq korder = iorder(k) xfic(k) = xabund(korder)*nhtot(i) enddo if (i < ntau) then do k=1,neq xatom(k) = xatom(k)*nhtot(i)/nhtot(i+1) enddo else do k=1,neq xatom(k) = xfic(k) enddo endif c*****compute the number of molecules: c*****Here is some information about the equilibrium constants. c*****The equilibrium constants, Kp, are defined as follows: c c c P(A)*P(B) N(A)*N(B) c Kp(AB) = --------- = kT--------- = c P(AB) N(AB) c c 2*pi*kT 3/2 M(A)*M(B) 3/2 Q(A)*Q(B) c = kT*(-------) * (---------) * --------- * exp(-D(AB)/kT) c h^2 M(AB) Q(AB) c c c*****MOOG uses: n(AB)/n(A)n(B) = kt/Kp c c Kp - dissociation constant, Q - partition functions, M - masses c P - partial pressures, N - number densities, T - temperature, c D - dissociation energy, h - plank constant. Remember to use c masses in grams (1 amu = 1.660540E-24 g) and energy in ergs c (1 eV = 1.60219E-12 ergs). Also, k = 1.38065E-16 erg/K, c h = 6.626076E-27 erg s, and pi = 3.1415926536. 27 do jmol=1,nmol atom = amol(jmol) if (atom >= 100.) then if (t(i) > 12000.) then xmol(jmol,i) = 1.0d-20 else xmol(jmol,i) = 1. count = 0. 37 call sunder(atom,iatom1,iatom2) count = count + 1. do k=1,neq if (iorder(k) == iatom1) . xmol(jmol,i) = xmol(jmol,i)*xatom(k) enddo if (iatom2 /= 0) then atom = iatom2 go to 37 endif hion = 10.0*(amol(jmol) - dint(amol(jmol))) th = 5040./t(i) lth = log10(th) xmol(jmol,i) = xmol(jmol,i)*(((1.38065d-16 * t(i))** . (count-1.0))/(10.0**(const(2,jmol)+(const(3,jmol)* . lth)+(const(4,jmol)*(lth**2))+(const(5,jmol)* . (lth**3))+(const(6,jmol)*(lth**4))- . (const(1,jmol)*th))))/(ne(i)**hion) endif c*****compute the number of ions: else delt = (t(i)-t(1))/tdel m = min0(idint(delt)+2,5) delt = delt - idint(delt) u1 = const(m,jmol) + . (const(m+1,jmol)-const(m,jmol))*delt iatom1 = atom do k=1,neq if (iorder(k) == iatom1) xmol(jmol,i) = . 4.825d15*u1*t(i)**1.5/ne(i)*dexp(-1.1605d4* . const(1,jmol)/t(i))*xatom(k) enddo endif enddo c*****compute matrix *c*, which is the derivative of each equation with c*****respect to each atom. do k=1,neq deltax(k) = -xfic(k) + xatom(k) do kk=1,neq c(k,kk) = 0. enddo c(k,k) = 1. korder = iorder(k) do kk=1,neq kderiv = iorder(kk) do 28 j=1,nmax jmol = ident(k,j) if (jmol == 0) go to 28 call discov(amol(jmol),kderiv,num2) if (num2 == 0) go to 28 call discov(amol(jmol),korder,num1) c(k,kk) = c(k,kk) + xmol(jmol,i)*num1*num2/xatom(kk) if (k == kk) deltax(k) = deltax(k) + num1*xmol(jmol,i) 28 continue enddo enddo c*****calculate array 'xcorr', the change in 'xatom'. array 'xcorr' is c*****'deltax' multiplied by the inverse of 'c' call invert (neq,c,ans,30) do k=1,neq x1 = xcorr(k) xcorr(k) = 0. do kk=1,neq xcorr(k) = xcorr(k) + ans(k,kk)*deltax(kk) enddo enddo c*****decide if another iteration is needed iflag = 0 do k=1,neq c*****here oscillations are damped out if (x1*xcorr(k) < -0.5*x1**2) xcorr(k) = 0.5*xcorr(k) x1 = xatom(k) if (dabs(xcorr(k)/xatom(k)) > 0.005) iflag = 1 xatom(k) = xatom(k) - xcorr(k) c*****fix element number densities which are ridiculous if (xatom(k)<=0.0 .or. xatom(k)>=1.001*xfic(k)) then iflag = 1 xatom(k) = 1.0d-2*dabs(xatom(k)+xcorr(k)) endif enddo if (iflag /= 0) go to 27 c*****print out atomic and molecular partial pressures. *xamol* is the c*****number density for each neutral atom do k=1,neq xamol(k,i) = xatom(k) patom(k) = dlog10(xatom(k)*tk) enddo do jmol=1,nmol pmol(jmol) = dlog10(xmol(jmol,i)*tk) enddo if (molopt >= 2) then pglog = dlog10(pgas(lev)) write (nf1out,1004) lev,int(t(lev)+0.001),pglog, . (patom(i),i=1,neq), (pmol(i),i=1,nmol) endif c*****here the big loop in tau ends 21 continue c*****finally, transfer some of this number density output to other arrays call setmols return c*****format statements 1001 format ('I do not know this molecule: ',f10.0) 1002 format (/'INPUT EQUILIBRIUM DATA:'/1x, 'species', 2x, . 'D0/Chi ', 3x, 'const1', 6x, 'const2', 6x, . 'const3', 6x, 'const4', 6x, 'const5'/ . (0pf8.1, f8.3, 1p5d12.4)) 1003 format (/'MOLECULAR EQUILIBRIUM SOLUTIONS: (log partial', . ' pressures listed under names)'/ . 2x,'i',5x,'T',2x,'Pgas',8f8.1/(15x,8f8.1)) 1004 format (i3,i6,f6.2,8f8.2/(15x,8f8.2)) end