program atmos include 'implno.dek' include 'const.dek' include 'network.dek' include 'vector_eos.dek' include 'atmos.dek' c..this program computes grey, lte, hydrostatic stellar atmospheres. c..declare external pzero character*80 file_root,outfile logical success integer i,j,k,n,niter,ltot,ng,lop,jrem,ilop,kb,ke, 1 lenstr,lenfile double precision teff,zbrent,taubeg,tauend,lo,hi,start,stopp, 1 step,bet,natural,sum,qtau,toler parameter (toler = 1.0d-12) c..for the ionization names character*6 roman(30) data roman /'I ', 'II ', 'III ', 'IV ', 'V ', 1 'VI ', 'VII ', 'VIII ', 'IX ', 'X ', 2 'XI ', 'XII ', 'XIII ', 'XIV ', 'XV ', 3 'XVI ', 'XVII ', 'XVIII ', 'XIX ', 'XX ', 4 'XXI ', 'XXII ', 'XXIII ', 'XXIV ', 'XXV ', 5 'XXVI ', 'XXVII ', 'XXVIII', 'XXIX ', 'XXX '/ c..for the hopf function spline fit c..hopdim = number of points in the spline c..tauray = tau points c..hopray = integration values at tau c..coef = spline coefficients integer hopdim parameter (hopdim=45) double precision tauray(hopdim),hopray(hopdim),coef(hopdim) c..data structures for the optical depth and hopf function c..tau array data tauray/ 0.0d0, 1.0d-4, 1.0d-3, 5.0d-3, 0.01d0, 0.02d0, 1 0.03d0, 0.05d0, 0.075d0, 0.1d0, 0.15d0, 0.2d0, 2 0.25d0, 0.3d0, 0.35d0, 0.4d0, 0.45d0, 0.5d0, 3 0.55d0, 0.6d0, 0.65d0, 0.7d0, 0.75d0, 0.8d0, 4 0.85d0, 0.9d0, 0.95d0, 1.0d0, 1.25d0, 1.5d0, 5 2.0d0, 2.5d0, 3.0d0, 3.5d0, 4.0d0, 4.5d0, 6 5.0d0, 6.0d0, 7.0d0, 8.0d0, 9.0d0, 10.0d0, 7 15.0d0, 20.0d0, 1.0d4/ c..the hopf function integration values data (hopray(i), i=1,30) / 1 5.7734470605805144E-01, 5.7749372259025067E-01, 2 5.7878311637859092E-01, 5.8361367328205960E-01, 3 5.8824571876827791E-01, 5.9548623664237532E-01, 4 6.0127056863486483E-01, 6.1073478587583896E-01, 5 6.2017475118816212E-01, 6.2790803944660067E-01, 6 6.4012847468710965E-01, 6.4954570640122333E-01, 7 6.5711469305603321E-01, 6.6336114367761934E-01, 8 6.6861194193240192E-01, 6.7308668731845589E-01, 9 6.7694097390571328E-01, 6.8028943386138840E-01, & 6.8321909226503308E-01, 6.8579760926267874E-01, 1 6.8807862048050283E-01, 6.9010533503061888E-01, 2 6.9191303889069611E-01, 6.9353088446143718E-01, 3 6.9498319972686284E-01, 6.9629046525143556E-01, 4 6.9747005605677348E-01, 6.9853681361570130E-01, 5 7.0256943722183141E-01, 7.0512862018715849E-01/ data (hopray(i), i=31,45) / 6 7.0791569376311414E-01, 7.0919255312074059E-01, 7 7.0980740940795184E-01, 7.1011381168281651E-01, 8 7.1027039279619053E-01, 7.1035197164159514E-01, 9 7.1039513093225670E-01, 7.1043076552803253E-01, & 7.1044135817064913E-01, 7.1044459903991219E-01, 1 7.1044561325364786E-01, 7.1044593531658484E-01, 2 7.1044608923199515E-01, 7.1044608986148083E-01, 3 7.1044608986444246E-01/ c..format statements 01 format('* model parameters :',/, 1 2('* ',t10,a,' ',1pe10.3,/), 2 '* ') 02 format(t2,a,t8,a,t16,a,t26,a,t36,a,t46,a, 2 t56,a,t66,a,t76,a) 03 format(i3,1p20e10.3) c..spline interpolation coefficients of the hopf function natural = 1.0e32 call spline(tauray,hopray,hopdim,natural,natural,coef) c..initialize names and numbers call init_h_to_zinc do i=1,ionmax xmass(i) = 1.0d-20 enddo c..set the composition of the atmosphere xmass(ih1) = 0.70 xmass(ihe4) = 0.28 xmass(ic12) = 0.004 xmass(in14) = 0.003 xmass(io16) = 0.004 xmass(ina22) = 0.002 xmass(img24) = 0.002 xmass(ial27) = 0.002 xmass(ife56) = 0.003 c..normalize the input composition sum = 0.0d0 do i=1,ionmax sum = sum + xmass(i) enddo sum = 1.0d0/sum do i=1,ionmax xmass(i) = min(1.0d0,max(1.0d-20,xmass(i)*sum)) enddo c..model parameters for the sun teff = 5800.0d0 surfg = g * msol/rsol**2 taubeg = 1.0e-3 tauend = 150.0d0 file_root = 'sun' c..choose nlayer number of equally log spaced optical depth points start = log10(taubeg) stopp = log10(tauend) step = (stopp - start)/float(nlayer-1) c..set the initial density limits the root finder lo = 1.0e-12 hi = 1.0e-4 write(6,02) 'n','tau','temp','den','xne','press','opac','depth' c..for each layer of the atmosphere c do layer = 1,1 do layer = 1,nlayer c..set the value of the optical depth tau; and its differential tau(layer) = 10.0d0 ** (start + float(layer-1) * step) if (layer .eq. 1) then dtau(layer) = tau(layer) else dtau(layer) = tau(layer) - tau(layer-1) endif c..get the temperature of this layer if (tau(layer) .gt. tauray(hopdim)) then qtau = hopray(hopdim) else call splint(tauray,hopray,coef,hopdim,tau(layer),qtau) end if temp(layer) = teff *(0.75d0*(tau(layer)+qtau))**(0.25d0) c..make sure the root is bracketed call zbrac(pzero,lo,hi,success) if (.not. success) stop 'could not bracket root' c..find the denity where the pressure from hydrostatic equilibrium equals the c..equation of state pressure. den(layer) = zbrent(pzero,lo,hi,toler,niter) c..reset the root finder limits so the next layer doesn't have to work as hard lo = 0.5d0 * den(layer) hi = 2.0d0 * den(layer) c..write out the physical variables as we go write(6,03) layer,tau(layer),temp(layer),den(layer), 1 xne(layer),press(layer),kappa(layer),depth(layer) c..and back for another layer enddo c..now write the results to files c..physical variables lenfile = lenstr(file_root,80) outfile = file_root(1:lenfile)//'_physical.dat' open(unit=2,file=outfile,status='unknown') write(2,01) 'teff = ',teff,'surfg = ',surfg write(2,02) 'n','tau','temp','dens','xne', 1 'press','opac','depth' do layer =1,nlayer write(2,03) layer,tau(layer),temp(layer), 1 den(layer),xne(layer),press(layer), 2 kappa(layer),depth(layer) enddo close(unit=2) c..ionization fractions in groups of 6 for each element c..lop is how many groups of 6 exist; jrem is the remainder do n = 1,ionmax if (xmass(n) .gt. 1.0e-16) then ltot = int(zion(n)) + 1 ng = 6 lop = int(ltot/ng) jrem = ltot - ng*lop do ilop = 1,lop+1 kb = 1 + ng*(ilop-1) ke = ng + ng*(ilop-1) if (ilop .eq. lop+1 .and. jrem .eq. 0) goto 33 if (ilop .eq. lop+1) ke = ltot write(outfile,111) file_root(1:lenfile),ionam(n),ilop,'.dat' 111 format(a4,'_',a5,'_',i2.2,a4) call sqeeze(outfile) open(unit=2,file=outfile,status='unknown') write(2,01) 'teff = ',teff,'surfg = ',surfg write(2,112) 'n','tau','temp',(ionam(n),roman(k),k=kb,ke) 112 format(1x,t2,a,t6,a,t16,a,t26,a,a,t36,a,a,t46,a,a, 1 t56,a,a,t66,a,a,t76,a,a) do layer = 1,nlayer write(2,03) layer,tau(layer),temp(layer), 1 (fracs(k,n,layer), k=kb,ke) enddo close(unit=2) c..end of loop over groups of six 33 continue enddo c..end of loop over elements end if enddo c..close up shop stop 'normal termination' end double precision function pzero(mden) include 'implno.dek' include 'const.dek' include 'network.dek' include 'vector_eos.dek' include 'vector_kap.dek' include 'atmos.dek' c..given the total number density of particles, this function computes the c..gas pressure by the equation of state and then by hydrostatic equilibrium. c..a root finder drives the difference between the two expressions to zero. c..declare integer i,j,k,kend,ifirst,olayer double precision mden,phydro,pold,rold,xx data ifirst/0/, olayer/0/ c..this only needs to be done once if (ifirst .eq. 0) then ifirst = 1 jlo_eos = 1 jhi_eos = 1 jlo_kap = 1 jhi_kap = 1 niso = ionmax do j=jlo_eos,jhi_eos do i=1,ionmax xmass_row(i,j) = xmass(i) aion_row(i,j) = aion(i) zion_row(i,j) = zion(i) enddo enddo end if c..this only needs to be done once for every new layer if (olayer .ne. layer) then olayer = layer pold = 0.0d0 rold = 0.0d0 do i = 1,layer-1 xx = dtau(i)/kappa(i) pold = pold + xx rold = rold + xx/den(i) enddo end if c..set the eos input temp_row(1) = temp(layer) den_row(1) = mden c..get the eos call eosion c..set the pressure of this layer press(layer) = ptot_row(1) c..set the electron number density xne(layer) = xne_row(1) c..set the ionzation fractions do i=1,ionmax kend = int(zion(i)) + 1 do k=1,kend fracs(k,i,layer) = frac_row(k,i,1) enddo enddo c..load the opacity input from the eos pin_row(1) = pele_row(1) + ppos_row(1) dpind_row(1) = dpepd_row(1) dpint_row(1) = dpept_row(1) xin_row(1) = xne_row(1) + xnp_row(1) dxind_row(1) = dxned_row(1) dxint_row(1) = dxnet_row(1) ein_row(1) = etaele_row(1) deind_row(1) = detad_row(1) deint_row(1) = detat_row(1) c..get the hydrogen, helium and metal mass fractions c..and a few other opacity composition variables call kapbar(xmass,aion,zion,ionmax, 1 xh_row(1),xhe_row(1),xz_row(1), 2 w4_row(1),w5_row(1),w6_row(1),zalf_row(1)) call kapfxt c..set the opacity of this layer kappa(layer) = opac_row(1) c..compute the gas pressure via hydrostatic equilibrium integral xx = dtau(layer)/kappa(layer) phydro = surfg * ( pold + xx ) c..compute the physical depth at this layer depth(layer) = rold + xx/mden c..for self consistency we want the pgas = pydro, i.e find the zero of pzero = press(layer) - phydro return end subroutine init_h_to_zinc include 'implno.dek' include 'network.dek' c..this routine initializes stuff for the aprox19 network c..set the size of the network and the number of rates ionmax = 30 c..set the id numbers of the elements ih1 = 1 ihe4 = 2 ili7 = 3 ibe9 = 4 ib11 = 5 ic12 = 6 in14 = 7 io16 = 8 if19 = 9 ine20 = 10 ina22 = 11 img24 = 12 ial27 = 13 isi28 = 14 ip31 = 15 is32 = 16 icl35 = 17 iar36 = 18 ik39 = 19 ica40 = 20 isc45 = 21 iti48 = 22 iv51 = 23 icr52 = 24 imn55 = 25 ife56 = 26 ico59 = 27 ini58 = 28 icu63 = 29 izn64 = 30 c..set the names of the elements ionam(ih1) = 'h1 ' ionam(ihe4) = 'he4 ' ionam(ili7) = 'li7 ' ionam(ibe9) = 'be9' ionam(ib11) = 'b11 ' ionam(ic12) = 'c12 ' ionam(in14) = 'n14 ' ionam(io16) = 'o16 ' ionam(if19) = 'f19 ' ionam(ine20) = 'ne20' ionam(ina22) = 'na22' ionam(img24) = 'mg24' ionam(ial27) = 'al27' ionam(isi28) = 'si28' ionam(ip31) = 'p31 ' ionam(is32) = 's32 ' ionam(icl35) = 'cl35' ionam(iar36) = 'ar36' ionam(ik39) = 'k39 ' ionam(ica40) = 'ca40' ionam(isc45) = 'sc45' ionam(iti48) = 'ti48' ionam(iv51) = 'v51 ' ionam(icr52) = 'cr52' ionam(imn55) = 'mn55 ' ionam(ife56) = 'fe56' ionam(ico59) = 'co59' ionam(ini58) = 'ni58' ionam(icu63) = 'cu63' ionam(izn64) = 'zn64' c..set the number of nucleons in the element aion(ih1) = 1.0d0 aion(ihe4) = 4.0d0 aion(ili7) = 7.0d0 aion(ibe9) = 10.0d0 aion(ib11) = 11.0d0 aion(ic12) = 12.0d0 aion(in14) = 14.0d0 aion(io16) = 16.0d0 aion(if19) = 19.0d0 aion(ine20) = 20.0d0 aion(ina22) = 22.0d0 aion(img24) = 24.0d0 aion(ial27) = 27.0d0 aion(isi28) = 28.0d0 aion(ip31) = 31.0d0 aion(is32) = 32.0d0 aion(icl35) = 35.0d0 aion(iar36) = 36.0d0 aion(ik39) = 39.0d0 aion(ica40) = 40.0d0 aion(isc45) = 45.0d0 aion(iti48) = 44.0d0 aion(iv51) = 51.0d0 aion(icr52) = 52.0d0 aion(imn55) = 55.0d0 aion(ife56) = 56.0d0 aion(ico59) = 59.0d0 aion(ini58) = 58.0d0 aion(icu63) = 63.0d0 aion(izn64) = 64.0d0 c..set the number of protons in the element zion(ih1) = 1.0d0 zion(ihe4) = 2.0d0 zion(ili7) = 3.0d0 zion(ibe9) = 4.0d0 zion(ib11) = 5.0d0 zion(ic12) = 6.0d0 zion(in14) = 7.0d0 zion(io16) = 8.0d0 zion(if19) = 9.0d0 zion(ine20) = 10.0d0 zion(ina22) = 11.0d0 zion(img24) = 12.0d0 zion(ial27) = 13.0d0 zion(isi28) = 14.0d0 zion(ip31) = 15.0d0 zion(is32) = 16.0d0 zion(icl35) = 17.0d0 zion(iar36) = 18.0d0 zion(ik39) = 19.0d0 zion(ica40) = 20.0d0 zion(isc45) = 21.0d0 zion(iti48) = 22.0d0 zion(iv51) = 23.0d0 zion(icr52) = 24.0d0 zion(imn55) = 25.0d0 zion(ife56) = 26.0d0 zion(ico59) = 27.0d0 zion(ini58) = 28.0d0 zion(icu63) = 29.0d0 zion(izn64) = 30.0d0 zion(ih1) = 1.0d0 zion(ihe3) = 2.0d0 zion(ihe4) = 2.0d0 zion(ic12) = 6.0d0 zion(in14) = 7.0d0 zion(io16) = 8.0d0 zion(ine20) = 10.0d0 zion(img24) = 12.0d0 zion(isi28) = 14.0d0 zion(is32) = 16.0d0 zion(iar36) = 18.0d0 zion(ica40) = 20.0d0 zion(iti48) = 22.0d0 zion(icr48) = 24.0d0 zion(ife52) = 26.0d0 zion(ife54) = 26.0d0 zion(ini56) = 28.0d0 zion(ineut) = 0.0d0 zion(iprot) = 1.0d0 return end subroutine spline(x,y,n,yp1,ypn,y2) include 'implno.dek' c..given arrays x(n) and y(n) containing a tabulated function y=f(x), c..this routine returns the array y2(n) of spline coefficients. c..x must be a monotonically increasing c.. c..if yp1 and/or ypn are set larger than 1e30, the routine sets the boundary c..condtion for a natural spline (zero second derivative) at that boundary, c..otherwise yp1 and ypn are the first derivatives to use at the endpoints c.. c..this routine is only called once for any given x and y. c..declare integer n,i,k,nmax parameter (nmax=1000) double precision x(n),y(n),y2(n),yp1,ypn,u(nmax),sig,p,qn,un c..the lower boundary condition is set to be either natural if (yp1 .gt. 1.0d30) then y2(1) = 0.0d0 u(1) = 0.0d0 c..or it has a specified first derivative else y2(1) = -0.5d0 u(1) = (3.0d0/(x(2)-x(1))) * ((y(2)-y(1))/(x(2)-x(1)) - yp1) end if c..this is the decomposition loop of the tridiagonal algorithm. y2 and u c..are used for temporary storage of the decomposition factors. do i=2,n-1 sig = (x(i) - x(i-1))/(x(i+1) - x(i-1)) p = sig*y2(i-1) + 2.0d0 y2(i) = (sig - 1.0d0)/p u(i) = (6.0d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) / 1 (x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p enddo c..the upper boundary condition is set to be either natural if (ypn .gt. 1.0d30) then qn = 0.0d0 un = 0.0d0 c..or it has a specified first derivative else qn = 0.5d0 un = (3.0d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) end if c..this is the backsubstitution loop of the tridiagonal algorithm y2(n) = (un-qn*u(n-1)) / (qn*y2(n-1) + 1.0d0) do k=n-1,1,-1 y2(k) = y2(k)*y2(k+1) + u(k) enddo return end subroutine splint(xa,ya,y2a,n,x,y) include 'implno.dek' c..given arrays xa and ya of length n, which tablulate a monotonic function, c..and given the array y2a, which is the output of spline (above), and given c..a value of x this routine returns a cubic spline interpolated value of y. c.. c..declare integer n,klo,khi,k double precision xa(n),ya(n),y2a(n),x,y,h,a,b c..find the right place in the table by bisection. this is optimal if c..the sequential calls to this routine are at random values of x. c..if the sequential calls are in order and closely spaced, one might store c..the values of klo and khi and test if they remain appropriate on next call klo = 1 khi = n 1 if (khi-klo .gt. 1) then k = (khi+klo)/2 if (xa(k) .gt. x) then khi = k else klo = k end if goto 1 end if c..klo and khi now bracket the input value of x; the xa's must be distinct h = xa(khi) - xa(klo) if (h .eq. 0.0) stop 'bad xa input in routine splint' c..evaluate the cubic spline a = (xa(khi) - x)/h b = (x - xa(klo))/h y = a*ya(klo) + b*ya(khi) + 1 ((a*a*a-a)*y2a(klo) + (b*b*b-b)*y2a(khi)) * (h*h)/6.0d0 return end subroutine zbrac(func,x1,x2,succes) include 'implno.dek' c..given a function func and an initial guessed range x1 to x2, the routine c..expands the range geometrically until a root is bracketed by the returned c..values x1 and x2 (in which case succes returns as .true.) or until the c..range becomes unacceptably large (in which succes returns as .false.). c..success guaranteed for a function which has opposite sign for sufficiently c..large and small arguments. c.. c..declare external func logical succes integer ntry,j parameter (ntry=50) double precision func,x1,x2,factor,f1,f2 parameter (factor = 1.6d0) if (x1 .eq. x2) stop ' x1 = x2 in routine zbrac' f1 = func(x1) f2 = func(x2) succes = .true. do j=1,ntry if (f1*f2 .lt. 0.0) return if (abs(f1) .lt. abs(f2)) then x1 = x1 + factor * (x1-x2) f1 = func(x1) else x2 = x2 + factor * (x2-x1) f2 = func(x2) end if enddo succes = .false. return end double precision function zbrent(func,x1,x2,tol,niter) include 'implno.dek' c.. c..using brent's method this routine finds the root of a function func c..between the limits x1 and x2. the root is when accuracy is less than tol. c.. c..note: eps the the machine floating point precision c.. c..declare external func integer niter,itmax,iter parameter (itmax = 100) double precision func,x1,x2,tol,a,b,c,d,e,fa, 1 fb,fc,xm,tol1,p,q,r,s,eps parameter (eps=3.0d-15) c..initialize niter = 0 a = x1 b = x2 fa = func(a) fb = func(b) if ( (fa .gt. 0.0 .and. fb .gt. 0.0) .or. 1 (fa .lt. 0.0 .and. fb .lt. 0.0) ) then write(6,100) x1,fa,x2,fb 100 format(1x,' x1=',1pe11.3,' f(x1)=',1pe11.3,/, 1 1x,' x2=',1pe11.3,' f(x2)=',1pe11.3) stop 'root not bracketed in routine zbrent' end if c = b fc = fb c..rename a,b,c and adjusting bound interval d do iter =1,itmax niter = niter + 1 if ( (fb .gt. 0.0 .and. fc .gt. 0.0) .or. 1 (fb .lt. 0.0 .and. fc .lt. 0.0) ) then c = a fc = fa d = b-a e = d end if if (abs(fc) .lt. abs(fb)) then a = b b = c c = a fa = fb fb = fc fc = fa end if tol1 = 2.0d0 * eps * abs(b) + 0.5d0 * tol xm = 0.5d0 * (c-b) c..convergence check if (abs(xm) .le. tol1 .or. fb .eq. 0.0) then zbrent = b return end if c..attempt quadratic interpolation if (abs(e) .ge. tol1 .and. abs(fa) .gt. abs(fb)) then s = fb/fa if (a .eq. c) then p = 2.0d0 * xm * s q = 1.0d0 - s else q = fa/fc r = fb/fc p = s * (2.0d0 * xm * q *(q-r) - (b-a)*(r - 1.0d0)) q = (q - 1.0d0) * (r - 1.0d0) * (s - 1.0d0) end if c..check if in bounds if (p .gt. 0.0) q = -q p = abs(p) c..accept interpolation if (2.0d0*p .lt. min(3.0d0*xm*q - abs(tol1*q),abs(e*q))) then e = d d = p/q c..or bisect else d = xm e = d end if c..bounds decreasing to slowly use bisection else d = xm e = d end if c..move best guess to a a = b fa = fb if (abs(d) .gt. tol1) then b = b + d else b = b + sign(tol1,xm) end if fb = func(b) enddo stop 'too many iterations in routine zbrent' end integer function lenstr(string,istrln) include 'implno.dek' c.. c..lenstr returns the non blank length length of the string. c.. c..declare integer istrln,i character*(*) string lenstr=0 do i=istrln,1,-1 if (string(i:i).ne. ' ') then if (ichar(string(i:i)).ne. 0 )then lenstr=i goto 20 end if end if enddo 20 return end subroutine sqeeze(line) include 'implno.dek' c.. c..this routine takes line and removes all blanks, such as c..those from writing to string with fortran format statements c.. c..declare character*(*) line character*1 achar integer l,n,k,lend,lsiz,lenstr c..find the end of the line lsiz = len(line) lend = lenstr(line,lsiz) n = 0 l = 0 c..do the compression in place 10 continue l = l + 1 achar = line(l:l) if (achar .eq. ' ') goto 10 n = n + 1 line(n:n) = achar if (l .lt. lend) goto 10 c..blank the rest of the line do k=n+1,lsiz line(k:k) = ' ' enddo return end c..download the kapfxt code from my conductivities page c..and include that subroutine here include 'public_kapfxt.f' c..download the eosion code from my ionization code page c..and include that subroutine here include 'public_ionization.f'