Initial commit

1 files changed, 513 insertions, 0 deletions

diff --git a/Setmols.r b/Setmols.r
new file mode 100755
index 0000000..149ba41
--- /dev/null
+++ b/Setmols.r
@@ -0,0 +1,513 @@
+
+      subroutine inmodel 
+c******************************************************************************
+c     This subroutine reads in the model 
+c****************************************************************************** 
+
+      implicit real*8 (a-h,o-z) 
+      include 'Atmos.com'
+      include 'Linex.com'
+      include 'Mol.com'
+      include 'Quants.com'
+      include 'Factor.com'
+      include 'Dummy.com'
+      include 'Pstuff.com'
+      real*8 element(95), logepsilon(95)
+      real*8 kaprefmass(100)
+      real*8 bmol(110)
+      integer nel(7)
+      character list*80, list2*70
+      integer append
+      data (nel(i),i=1,7)/  1,  2,  6, 12, 13, 14, 26/
+
+
+c*****Read in the key word to define the model type
+      modelnum = modelnum + 1
+      rewind nfmodel
+      read (nfmodel,2001) modtype
+      write (nf1out,1010) modtype
+      if (modtype .eq. 'begn      ' .or.  modtype .eq. 'BEGN      ') 
+     .   write (nf1out,1011)
+
+
+c*****Read a comment line (usually describing the model)
+      read (nfmodel,2002) moditle
+
+
+c*****Read the number of depth points
+      read (nfmodel,2002) list
+      list2 = list(11:)
+      read (list2,*) ntau
+      if (ntau .gt. 100) then
+         write (array,1012)
+         call prinfo (10)
+         stop
+      endif
+
+
+c*****EITHER: Read in a model from the output of the experimental new
+c     MARCS code.  This modtype is called "NEWMARCS".  On each line 
+c     the numbers are:
+c     tau(5000), t, pe, pgas, rho,  model microtrubulent velocity,
+c     and mean opacity (cm^2/gm) at the reference wavelength (5000A).
+      if (modtype .eq. 'NEWMARCS  ') then
+         read (nfmodel,*) wavref    
+         do i=1,ntau
+            read (nfmodel,*) tauref(i),t(i),ne(i),pgas(i),rho(i),
+     .                         vturb(1),kaprefmass(i)
+         enddo
+c*****OR: Read in a model from the output of the on-line new
+c     MARCS code.  This modtype is called "WEBMARCS".  On each line
+c     the numbers are:
+c     layer number (not needed), log{tau(Rosseland)} (not needed),
+c     log{tau(5000)}, depth, t, pe, pgas, prad (not read in) and
+c     pturb (not read in)
+      elseif (modtype .eq. 'WEBMARCS') then
+         read (nfmodel,*) wavref
+         do i=1,ntau
+            read (nfmodel,*) k, dummy1(k), tauref(i), dummy2(k), t(i),
+     .                       ne(i), pgas(i)
+         enddo
+c*****OR: Read in a model from an alternative form of on-line new
+c     MARCS code.  This modtype is called "WEB2MARC".  On each line
+c     the numbers are:
+c     atmospheric layer number (not needed), log{tau(5000)}, t, 
+c     log(Pe), log(Pgas), rhox
+      elseif (modtype .eq. 'WEB2MARC') then
+         read (nfmodel,*) wavref
+         do i=1,ntau
+            read (nfmodel,*) k,tauref(i),t(i),ne(i),pgas(i),rhox(i)
+         enddo
+c     OR: Read in a model from the output of the ATLAS code.  This
+c     modtype is called "KURUCZ".  On each line the numbers are:
+c     rhox, t, pgas, ne, and Rosseland mean opacity (cm^2/gm), and
+c     two numbers not used by MOOG.  
+      elseif (modtype .eq. 'KURUCZ    ') then
+         do i=1,ntau
+            read (nfmodel,*) rhox(i),t(i),pgas(i),ne(i),kaprefmass(i)
+         enddo
+c     OR: Read in a model from the output of the NEXTGEN code.  This
+c     modtype is called "NEXTGEN".  These models have tau scaled at a 
+c     specific wavelength that is read in before the model. MOOG will 
+c     need to generate the opacities internally.On each line the numbers 
+c     are: tau, t, pgas, pe, density, mean molecular weight, two numbers
+c     not used by MOOG, and Rosseland mean opacity (cm^2/gm).
+      elseif (modtype .eq. 'NEXTGEN   ') then
+         read (nfmodel,*) wavref
+         do i=1,ntau
+            read (nfmodel,*) tauref(i),t(i),pgas(i),ne(i), rho(i),  
+     .                       molweight(i), x2, x3, kaprefmass(i)
+         enddo
+c     OR: Read in a model from the output of the MARCS code.  This modtype
+c     type is called "BEGN".  On each line the numbers are:
+c     tauross, t, log(pg), log(pe), mol weight, and kappaross.
+      elseif (modtype .eq. 'BEGN      ') then
+         do i=1,ntau
+            read (nfmodel,*) tauref(i),t(i),pgas(i),ne(i),
+     .                          molweight(i),  kaprefmass(i)
+         enddo
+c     OR: Read in a model generated from ATLAS, but without accompanying
+c     opacities.  MOOG will need to generate the opacities internally,
+c     using a reference wavelength that it reads in before the model.
+      elseif (modtype .eq. 'KURTYPE') then
+         read (nfmodel,*) wavref    
+         do i=1,ntau
+            read (nfmodel,*) rhox(i),t(i),pgas(i),ne(i)
+         enddo
+c     OR: Read in a model generated from ATLAS, with output generated
+c     in Padova.  The columns are in somewhat different order than normal
+      elseif (modtype .eq. 'KUR-PADOVA') then
+         read (nfmodel,*) wavref
+         do i=1,ntau
+            read (nfmodel,*) tauref(i),t(i),kaprefmass(i),
+     .                         ne(i),pgas(i),rho(i)
+         enddo
+c     OR: Read in a generic model that has a tau scale at a specific 
+c     wavelength that is read in before the model.  
+c     MOOG will need to generate the opacities internally.
+      elseif (modtype .eq. 'GENERIC   ') then
+         read (nfmodel,*) wavref    
+         do i=1,ntau
+            read (nfmodel,*) tauref(i),t(i),pgas(i),ne(i)
+         enddo
+c     OR: quit in utter confusion if those model types are not specified
+      else
+         write (*,1001)
+         stop
+      endif
+
+
+c*****Compute other convenient forms of the temperatures
+      do i=1,ntau
+          theta(i) = 5040./t(i)
+          tkev(i) = 8.6171d-5*t(i)
+          tlog(i) = dlog(t(i))
+      enddo
+
+
+c*****Convert from logarithmic Pgas scales, if needed
+      if (pgas(ntau)/pgas(1) .lt. 10.) then
+         do i=1,ntau                                                    
+            pgas(i) = 10.0**pgas(i)
+         enddo
+      endif
+
+
+c*****Convert from logarithmic Ne scales, if needed
+      if(ne(ntau)/ne(1) .lt. 20.) then
+         do i=1,ntau                                                    
+            ne(i) = 10.0**ne(i)
+         enddo
+      endif
+
+
+c*****Convert from Pe to Ne, if needed
+      if(ne(ntau) .lt. 1.0e7) then
+         do i=1,ntau                                                    
+            ne(i) = ne(i)/1.38054d-16/t(i)
+         enddo
+      endif
+
+
+c*****compute the atomic partition functions
+      do j=1,95
+         elem(j) = dble(j)
+         call partfn (elem(j),j)
+      enddo
+
+
+c*****Read the microturbulence (either a single value to apply to 
+c     all layers, or a value for each of the ntau layers). 
+c     Conversion to cm/sec from km/sec is done if needed
+      read (nfmodel,2003) (vturb(i),i=1,6)
+      if (vturb(2) .ne. 0.) then
+         read (nfmodel,2003) (vturb(i),i=7,ntau) 
+      else
+         do i=2,ntau                                                    
+            vturb(i) = vturb(1)
+         enddo
+      endif
+      if (vturb(1) .lt. 100.) then
+         write (moditle(55:62),1008) vturb(1)
+         do i=1,ntau
+            vturb(i) = 1.0e5*vturb(i)
+         enddo
+      else
+         write (moditle(55:62),1008) vturb(1)/1.0e5
+      endif
+
+
+c*****Read in the abundance data, storing the original abundances in xabu
+c*****The abundances not read in explicity are taken from the default
+c*****solar ones contained in array xsolar.
+      read (nfmodel,2002) list
+      list2 = list(11:)
+      read (list2,*) natoms,abscale
+      write (moditle(63:73),1009) abscale
+      if(natoms .ne. 0) 
+     .         read (nfmodel,*) (element(i),logepsilon(i),i=1,natoms) 
+      xhyd = 10.0**xsolar(1)
+      xabund(1) = 1.0
+      xabund(2) = 10.0**xsolar(2)/xhyd
+      do i=3,95                                                      
+         xabund(i) = 10.0**(xsolar(i)+abscale)/xhyd
+         xabu(i) = xabund(i)
+      enddo
+      if (natoms .ne. 0) then
+         do i=1,natoms                                                  
+            xabund(idint(element(i))) = 10.0**logepsilon(i)/xhyd
+            xabu(idint(element(i))) = 10.0**logepsilon(i)/xhyd
+         enddo
+      endif
+
+c*****Compute the mean molecular weight, ignoring molecule formation
+c     in this approximation (maybe make more general some day?)
+      wtnum = 0.
+      wtden = 0.
+      do i=1,95
+         wtnum = wtnum + xabund(i)*xam(i)
+         wtden = wtden + xabund(i)
+      enddo
+      wtmol = wtnum/(xam(1)*wtden)
+      nomolweight = 0
+      if (modtype .eq. 'BEGN      ' .or. modtype .eq. 'NEXTGEN   ') then
+         nomolweight = 1
+      endif
+      if (nomolweight .ne. 1) then
+         do i=1,ntau
+             molweight(i) = wtmol
+         enddo
+      endif
+
+c*****Compute the density 
+      if (modtype .ne. 'NEXTGEN   ') then
+         do i=1,ntau                                                    
+            rho(i) = pgas(i)*molweight(i)*1.6606d-24/(1.38054d-16*t(i))
+         enddo
+      endif
+
+
+c*****Calculate the fictitious number density of hydrogen
+c     Note:  ph = (-b1 + dsqrt(b1*b1 - 4.0*a1*c1))/(2.0*a1)
+      iatom = 1
+      call partfn (dble(iatom),iatom)
+      do i=1,ntau    
+         th = 5040.0/t(i)         
+         ah2 = 10.0**(-(12.7422+(-5.1137+(0.1145-0.0091*th)*th)*th))
+         a1 = (1.0+2.0*xabund(2))*ah2
+         b1 = 1.0 + xabund(2)
+         c1 = -pgas(i)         
+         ph = (-b1/2.0/a1)+dsqrt((b1**2/(4.0*a1*a1))-(c1/a1))
+         nhtot(i) = (ph+2.0*ph*ph*ah2)/(1.38054d-16*t(i))
+      enddo
+
+
+c*****Molecular equilibrium called here.
+c     First, a default list of ions and molecules is considered. Then the 
+c     user's list is read in. A check is performed to see if any of these 
+c     species need to be added. If so, then they are appended to the 
+c     default list. The molecular equilibrium routine is then called.
+c     Certain species are important for continuous opacities and damping
+c     calculations - these are read from the equilibrium solution and 
+c     saved.
+
+
+c*****Set up the default molecule list
+      if (molset .eq. 0) then
+         nmol = 30
+      else
+         nmol = 59
+      endif
+      if     (molset .eq. 0) then
+         do i=1,110
+            amol(i) = smallmollist(i)
+         enddo
+      elseif (molset .eq. 1) then
+         do i=1,110
+            amol(i) = largemollist(i)
+         enddo
+      else
+         array = 'molset = 0 or 1 only; I quit!'
+         call putasci (29,6)
+         stop
+      endif
+
+
+c*****Read in the names of additional molecules to be used in 
+c     molecular equilibrium if needed.
+      read (nfmodel,2002,end=101) list
+      list2 = list(11:)
+      read (list2,*) moremol
+      if (moremol .ne. 0) then
+         read (nfmodel,*) (bmol(i),i=1,moremol)
+         append = 1
+         do k=1,moremol
+            do l=1,nmol
+               if (nint(bmol(k)) .eq. nint(amol(l))) 
+     .         append = 0  
+            enddo
+            if (append .eq. 1) then 
+               nmol = nmol + 1
+               amol(nmol) = bmol(k)
+            endif
+            append = 1
+         enddo  
+      endif
+
+
+c*****do the general molecular equilibrium
+101   call eqlib
+      call resetmol
+
+
+c     In the number density array "numdens", the elements denoted by
+c     the first subscripts are named in at the ends of the assignment
+c     lines; at present these are the only ones needed for continuous 
+c     opacities
+c     the desired elements are H, He, C, Mg, Al, Si, Fe, along
+c     with the H_2 molecule
+c     first do the neutrals
+      do j=1,7
+         do k=1,neq
+            if (nel(j) .eq. iorder(k)) then
+               do i=1,ntau
+                  numdens(j,1,i) = xamol(k,i)
+               enddo
+               exit
+            endif
+         enddo
+      enddo
+c     then do the ions
+      do j=1,7
+         ispec10 = nint(dble(10*nel(j)+1))
+         do k=1,nmol
+            kmol10 = nint(10.*amol(k))
+            if (ispec10 .eq. kmol10) then
+               do i=1,ntau
+                  numdens(j,2,i) = xmol(k,i)
+               enddo
+               exit
+            endif
+         enddo
+      enddo
+c     finally add in H_2
+      do k=1,nmol
+         ispec = 101
+         if (ispec .eq. nint(amol(k))) then
+            do i=1,ntau
+               numdens(8,1,i) = xmol(k,i)
+            enddo
+            exit
+         endif
+      enddo
+
+
+c*****compute partitiion functions for H_2O and CO_2; 
+      do i=1,ntau
+         h2olog = 0.
+         co2log = 0.
+         if (t(i) .gt. 5000.) then
+            uh2o(i) = 1.d8
+            uco2(i) = 1.d8
+         else
+            do j=1,5
+               h2olog = h2olog + h2ocoeff(j)*t(i)**(j-1)
+               co2log = co2log + co2coeff(j)*t(i)**(j-1)
+            enddo
+            uh2o(i) = 10**h2olog
+            uco2(i) = 10**co2log
+         endif
+      enddo
+
+
+c*****transfer H_2O and CO_2 number densities from the molecular 
+c     equilibrium output
+c     note:  HITRAN partition functions are given only for T < 5000K
+      do j=1,nmol
+         if (nint(amol(j)) .eq. 10108) ih2o = j
+         if (nint(amol(j)) .eq. 60808) ico2 = j
+      enddo
+      do i=1,ntau
+         xnh2o(i) = xmol(ih2o,i)
+         xnco2(i) = xmol(ico2,i)
+      enddo
+         
+      
+c*****SPECIAL NEEDS: for NEWMARCS models, to convert kaprefs to our units
+      if (modtype .eq. 'NEWMARCS  ') then
+         do i=1,ntau
+            kapref(i) = kaprefmass(i)*rho(i)
+         enddo
+c     SPECIAL NEEDS: for KURUCZ models, to create the optical depth array,
+c     and to convert kaprefs to our units
+      elseif (modtype .eq. 'KURUCZ    ') then
+         first = rhox(1)*kaprefmass(1)
+         tottau = rinteg(rhox,kaprefmass,tauref,ntau,first) 
+         tauref(1) = first
+         do i=2,ntau
+            tauref(i) = tauref(i-1) + tauref(i)
+         enddo
+         do i=1,ntau
+            kapref(i) = kaprefmass(i)*rho(i)
+         enddo
+c     SPECIAL NEEDS: for NEXTGEN models, to convert kaprefs to our units
+      elseif (modtype .eq. 'NEXTGEN   ') then
+         do i=1,ntau                                                    
+            kapref(i) = kaprefmass(i)*rho(i)
+         enddo
+c     SPECIAL NEEDS: for BEGN models, to convert kaprefs to our units
+      elseif (modtype .eq. 'BEGN      ') then
+         do i=1,ntau                                                    
+            kapref(i) = kaprefmass(i)*rho(i)
+         enddo
+c     SPECIAL NEEDS: for KURTYPE models, to create internal kaprefs,
+c     and to compute taurefs from the kaprefs converted to mass units
+      elseif (modtype .eq. 'KURTYPE   ') then
+         call opacit (1,wavref)
+         do i=1,ntau                                                    
+            kaprefmass(i) = kapref(i)/rho(i)
+         enddo
+         first = rhox(1)*kaprefmass(1)
+         tottau = rinteg(rhox,kaprefmass,tauref,ntau,first) 
+         tauref(1) = first
+         do i=2,ntau
+            tauref(i) = tauref(i-1) + tauref(i)
+         enddo
+c     SPECIAL NEEDS: for NEWMARCS models, to convert kaprefs to our units
+      elseif (modtype .eq. 'KUR-PADOVA') then
+         do i=1,ntau
+            kapref(i) = kaprefmass(i)*rho(i)
+         enddo
+c     SPECIAL NEEDS: for generic models, to create internal kaprefs,
+      elseif (modtype .eq. 'GENERIC   ' .or.
+     .        modtype .eq. 'WEBMARCS  ' .or.
+     .        modtype .eq. 'WEB2MARC  ') then
+         call opacit (1,wavref)
+      endif
+
+
+c*****Convert from logarithmic optical depth scales, or vice versa.
+c     xref will contain the log of the tauref
+      if(tauref(1) .lt. 0.) then
+         do i=1,ntau                                                    
+            xref(i) = tauref(i)
+            tauref(i) = 10.0**xref(i)
+         enddo
+      else
+         do i=1,ntau                                                    
+            xref(i) = dlog10(tauref(i))
+         enddo
+      endif
+
+
+c*****Write information to output files
+      if (modprintopt .lt. 1) return
+      write (nf1out,1002) moditle
+      do i=1,ntau
+         dummy1(i) = dlog10(pgas(i))
+         dummy2(i) = dlog10(ne(i)*1.38054d-16*t(i))
+      enddo
+      write (nf1out,1003) wavref,(i,xref(i),tauref(i),t(i),dummy1(i),
+     .                    pgas(i),dummy2(i),ne(i),vturb(i),i=1,ntau)
+      write (nf1out,1004)
+      do i=1,95
+         dummy1(i) = dlog10(xabund(i)) + 12.0
+      enddo
+      write (nf1out,1005) (names(i),i,dummy1(i),i=1,95)
+      write (nf1out,1006) modprintopt, molopt, linprintopt, fluxintopt
+      write (nf1out,1007) (kapref(i),i=1,ntau)
+      return
+
+
+c*****format statements
+2001  format (a10)
+2002  format (a80)
+2003  format (6d13.0)
+1001  format('permitted model types are:'/'KURUCZ, BEGN, ',
+     .       'KURTYPE, KUR-PADOVA, NEWMARCS, WEBMARCS, NEXTGEN, ',
+     .       'WEB2MARC, or GENERIC'/ 'MOOG quits!')
+1002  format (/'MODEL ATMOSPHERE HEADER:'/a80/)
+1003  format ('INPUT ATMOSPHERE QUANTITIES', 10x,
+     .        '(reference wavelength =',f10.2,')'/3x, 'i', 2x, 'xref',
+     .        3x, 'tauref', 7x, 'T', 6x, 'logPg', 4x, 'Pgas',
+     .        6x, 'logPe', 5x, 'Ne', 9x, 'Vturb'/
+     .        (i4, 0pf6.2, 1pd11.4, 0pf9.1, f8.3, 1pd11.4, 0pf8.3,
+     .        1pd11.4, d11.2))
+1004  format (/'INPUT ABUNDANCES: (log10 number densities, log H=12)'/
+     .       '      Default solar abundances: Asplund et al. 2009')
+1005  format (5(3x,a2, '(',i2,')=', f5.2))
+1006  format (/'OPTIONS: atmosphere = ', i1, 5x, 'molecules  = ', i1/
+     .        '         lines      = ', i1, 5x, 'flux/int   = ', i1)
+1007  format (/'KAPREF ARRAY:'/(6(1pd12.4)))
+1008  format ('vt=', f5.2)
+1009  format (' M/H=', f5.2)
+1010  format (13('-'),'MOOG OUTPUT FILE', 10('-'),
+     .        '(MOOG version from 01/28/09)', 13('-')//
+     .        'THE MODEL TYPE: ', a10)
+1011  format ('   The Rosseland opacities and optical depths have ',
+     .        'been read in')
+1012  format ('HOUSTON, WE HAVE MORE THAN 100 DEPTH POINTS! I QUIT!')
+
+
+      end
+