program readopa
* Sample programme to call subroutine readmo which reads a MARCS opacity model file
* Version 1.2, 10-Oct-2006 /BE
*
* Parameters are described in subroutine readmo
*
* compile by e.g.: f90 -o readopa readopa.f"
*
* Run by either   "cat opa_file | readopa" if unzipped
*            or  "zcat opa_file | readopa" if zipped
*    or perhaps "gzcat opa_file | readopa" if gzipped
*
      implicit none
      integer mdp,mwav,j,k
      parameter (mdp=56)
      parameter (mwav=1071)
      real f,a,s,w
**************
      character mcode*4
      integer ndp,nwav
      real abund(92),wav(mwav),swave
      common /global/ mcode,ndp,swave,nwav,wav,abund
**************
      real rad(mdp),tau(mdp),t(mdp),pe(mdp),pg(mdp),rho(mdp),xi(mdp)
      real ops(mdp)
      common /structure/ rad,tau,t,pe,pg,rho,xi,ops
**************
      real abs(mwav,mdp),sca(mwav,mdp)
      common /opacities/ abs,sca
**************

      call readmo

* sample printout
* list tau, t, pg, abs(5000 A), sca(5000 A), and ops(5000 A)
      w=5000.0
      print *,'mcode=',mcode,', test printout for wavelength=',w,' AA'
      if(w.lt.wav(1) .or. w.gt.wav(nwav)) then
        stop 'no wavelength extrapol'
      endif
      print 1000
 1000 format(' k tau(5000)    T       Pg       ops       abs       sca')
* find the wavelength-interpolation points
      do j=1,nwav
        if(w.lt.wav(j)) then
          f=(w-wav(j-1))/(wav(j)-wav(j-1))
          do k=1,ndp
* linear interpolation
            a=abs(j-1,k)+(abs(j,k)-abs(j-1,k))*f
            s=sca(j-1,k)+(sca(j,k)-sca(j-1,k))*f
            print 1010,k,tau(k),t(k),pg(k),ops(k),a,s
 1010       format(i2,1p,e10.3,0p,f8.1,1p,4e10.3,0p)
          enddo
          exit
        endif
      enddo
      end
*
*
*...v....1....v....2....v....3....v....4....v....5....v....6....v....7..
*
*
      subroutine readmo
* reads a MARCS opacity data file and stores the data into arrays in 3
* common blocks: 'global', 'structure', and 'opacities'
*
************************************************************************
* Global properties of the opacity model file: common "global"
* mcode  = the 4-byte standard model code 'MRXF'
* ndp    = number of depth points (=56)
* k      = depth index
* swave  = standard wavelength for the continuous optical depth (tau)
*          scale and for the total standard opacity (ops)
* nwav   = number of wavelengths for which continuous absorption and
*          scattering opacities are given. These are chosen so that
*          linear interpolation should suffice for any wavelength
*          (=1071)
* wav(j) = wavelengths for which opacities are given
* j      = wavelength index
* abund  = logarithmic number abundances of the 92 first chemical
*          elements on a scale where the hydrogen abundance=12.00
************************************************************************
* Model structure: common "structure"
* rad(k) = radius, normalized on the outermost point, k=1. For use with
*          spherical radiative transfer.
*          For plane-parallel models rad == 1.0.
* tau(k) = continuumm optical depth at the standard wavelength swave
* t(k)   = temperature (K)
* pe(k)  = electron pressure (dyn/cm2)
* pg(k)  = total gas pressure (dyn/cm2)
* rho(k) = densigy (g/cm3)
* xi(k)  = microturbulence parameter (km/s)
* ops(k) = continuumm opacity at the standard wavelength (cm2/g)
************************************************************************
* Wavelength dependent opacities: common "opacities"
* abs(j,k) = specific continuous absorption opacity (cm2/g)
* sca(j,k) = specific continuous scattering opacity (cm2/g)
************************************************************************
*
      implicit none
      integer mdp,mwav,i,k
      parameter (mdp=56)
      parameter (mwav=1071)
**************
      character mcode*4
      integer ndp,nwav
      real abund(92),wav(mwav),swave
      common /global/ mcode,ndp,swave,nwav,wav,abund
**************
      real rad(mdp),tau(mdp),t(mdp),pe(mdp),pg(mdp),rho(mdp),xi(mdp)
      real ops(mdp)
      common /structure/ rad,tau,t,pe,pg,rho,xi,ops
**************
      real abs(mwav,mdp),sca(mwav,mdp)
      common /opacities/ abs,sca
**************
*
      read(*,1010) mcode,ndp,swave
 1010 format(1x,a4,i5,f10.0)
      if(ndp.gt.mdp) then
        stop 'unexpected (ly large) number of depth points'
      endif
      read(*,*) nwav
      if(nwav.gt.mwav) then
        stop 'unexpected (ly large) number of wavelength points'
      endif
      read(*,*) (wav(i),i=1,nwav)  
      do k=1,ndp
        read(*,*) rad(k),tau(k),t(k),pe(k),pg(k),rho(k),xi(k),ops(k)
        read(*,*) (abs(i,k),sca(i,k),i=1,nwav)
        do i=1,nwav
          abs(i,k)=abs(i,k)*ops(k)
          sca(i,k)=sca(i,k)*ops(k)
        enddo
      enddo
      read(*,*) abund
*
      return
      end