C     program ask    --- driver program for interpolation in the W94 grids
      program ask
C     USER! - please edit lines 35 and 37 below before compiling.

C                  histpar is the number of bins the user can combine
      integer      histpar
      parameter    (histpar = 28)
      real         metalref(histpar),ageref(histpar),shist(histpar)
      character*80 text
      integer      ndir,idir
      parameter    (ndir=11)
      character*11 dirnames(ndir)
      character*48 explan(ndir)
      character*80 rootdir,dirname
      integer      nchar
      real         redshift

      data dirnames / 'Wo__vanilla','Wo_Salp+0.5','Wo_Salp+1.0',
     &                'Wo_Salp-0.5','Wo_Salp-1.0','Wo_MillerSc',
     &                'Wo_Yshallow','Wo_Yminimum','Wo_Ymaximum',
     &                'Wo_NoLateSt','Wo_original' /
      data explan /'Zero model: Salpeter IMF, dY/dZ=2.7             ',
     &             'IMF exponent 2.85                               ',
     &             'IMF exponent 3.35                               ',
     &             'IMF exponent 1.85                               ',
     &             'IMF exponent 1.35                               ',
     &             'Miller & Scalo (1979) piecewise IMF             ',
     &             'Y=.228+6Z below Z=.0077, Y=.274 above that value',
     &             'Y=.2 below [Fe/H]=-.5, Y=.25 above [Fe/H]=-.25  ',
     &             'Y=.3 below [Fe/H]=-.5, Y=.35 above [Fe/H]=-.25  ',
     &             'No horizontal branch (or clump) or later stages ',
     &             'Untweaked model grid exactly as in W94          '/

C     USER EDIT NEXT FEW LINES
C     Let your program know where the models are. Leave trailing slash. 
 35   rootdir = '/usr/a/www/export2/'
C     Change next line to reflect the number of characters in "rootdir"
 37   nchar = 19
C     That is all. It should compile and run now. (Cross your fingers!)

      print*, ' '
      print*, ' '
      print*, 'ASK returns a flux, colors, SBF mags, and indices'
      print*, 'for combinations of single-burst populations. '
      print*, ' '
      print*, 'First, choose the model flavor (default: vanilla)'
      do i=1,ndir
         write(*,'(i2,2x,a11,4x,a48)') i, dirnames(i), explan(i)
      end do
 5    print*, ' '
      print*, 'Enter a number: '
      read(*,'(a80)') text
      if (text.eq.' ') then
         idir = 1
         goto 7
      end if
      read(text,*,err=5) idir
      if (idir.lt.1 .or. idir.gt.ndir) goto 5
 7    print*, ' '
      print*, 'Directory: ',dirnames(idir)
      dirname = rootdir(1:nchar)//dirnames(idir)
      nchar = nchar + 11
      print*, ' '
      print*, 'You can enter up to ',histpar,
     &        ' separate populations, along '
      print*, 'with weights (which the program will normalize so that '
      print*, 'the total population mass is a million Mo). '
      print*, 'Metallicities are restricted to -2 < [Fe/H] < +0.5, and '
      print*, 'ages to between 1 and 18 Gyr. In addition, the age<8, '
      print*, '[Fe/H]<-0.225 corner is not allowed. '
      print*, ' '
      print*, 'Enter trios of numbers separated by spaces: '
      print*, '( age  [Fe/H]  weight ). '
      print*, 'Enter a blank to exit the loop '
      print*, ' ' 

      do i=1,histpar
 11      print*, 'Enter (age [Fe/H] weight) no. ',i
         read(*,'(a80)') text
         if (text.eq.' ') then
            npop =  i-1
            goto 22
         end if
         read(text,*,err=99)  ageref(i),metalref(i),shist(i)
         goto 12
 99      print*, 'I had trouble reading that. Try again.'
         goto 11
 12   end do
      npop = 28

C     normalize mass histogram to integral one
 22   xsum = 0.0
      print*, ' npop = ',npop
      do j=1,npop
         xsum = xsum + shist(j)
      end do
      do j=1,npop
         shist(j) = shist(j)/xsum
      end do

      call asksub(metalref,ageref,shist,npop,dirname,nchar,redshift)

      stop 
      end

C     SUBROUTINE ASKSUB -------------------------------------------------
      subroutine asksub(metalref,ageref,shist,npop,dirname,nchar,
     &                                                   redshift)
C                  histpar is the number of bins in the [Fe/H] histogram
C                  nind is the number of indices
C                  nclrs is the number of colors
C                  npop is the number of populations to combine
      integer      nind,histpar,nclrs,npop
      parameter    (nind = 71, histpar = 28, nclrs=11)
      real         metalref(histpar),ageref(histpar),shist(histpar)
      character    dirname*80
      integer      nchar
      real         redshift

C     primary passed variables
      real         hlam(1221),wave(1221), bolo, z, y, age
      real         mags(nclrs),Lsq(nclrs),Lavg(nclrs)
      real         widths(nind),Fcont(nind),Fline(nind)
      integer      indexunits(nind)
      character*10 labels(nind)
      character*20 imftype
      real         massinit,massnow
      integer      nsyn,nsynmax
      parameter    (nsynmax=3300)
      real         synwave(nsynmax),synflx(nsynmax)
C     other passed variables
      real         imfpar, zsolar, ysolar, yprimo, imflow, imfupp
      real         hbmorph, pagb, eta

C    
      real         zin, agein
C     variables for storing things.
      real         ahlam(1221,histpar), abolo(histpar)
      real         asynflx(nsynmax,histpar)
      real         amags(nclrs,histpar)
      real         aLsq(nclrs,histpar),aLavg(nclrs,histpar)
      real         aFcont(nind,histpar),aFline(nind,histpar)
      real         amassinit(histpar),amassnow(histpar)
      real         aind(nind,histpar)
C     variables for output
      real         index(nind), Lbar(nclrs)

C     other variables
      character*24   ctime,t
      integer        time,status
C     array of assumed solar absolute magnitudes for computation of
C     M/L ratios. Alternatively, I could have read the M/L directly from the
C     FITS files, but this way seems more transparent. UBVRcIcJHKLL'M.
      real           solmags(nclrs),mtol(nclrs)
      data      solmags  / 5.60, 5.51, 4.84, 4.48, 4.13, 3.70, 3.37,
     z            3.33, 3.27, 3.28, 3.30 /

C     get and store the interpolated population models
      do i=1,npop
         zin = 0.0169*10**(metalref(i))
         agein = ageref(i)

         call getmodel(zin,agein,dirname,nchar,
     &          hlam,wave,bolo,z,y,age,mags,Lsq,Lavg,
     &       widths,Fcont,Fline,indexunits,labels,massinit,
     &       massnow,imfpar,zsolar,
     &       ysolar,yprimo,imflow,imfupp,hbmorph,pagb,eta,imftype,
     &       synwave,synflx,nsyn,
     &       status)


C        store important results
         do j=1,1221
            ahlam(j,i) = hlam(j)
         end do
         do j=1,nsyn
            asynflx(j,i) = synflx(j)
         end do
         abolo(i) = bolo
         amassinit(i) = massinit
         amassnow(i)  = massnow

         do j=1,nclrs
            amags(j,i) = mags(j)
            aLsq(j,i) = Lsq(j)
            aLavg(j,i) = Lavg(j)
c            print*, ' mag, Lsq, Lavg = ',mags(j),Lsq(j),Lavg(j)
         end do
         do j=1,nind
            aFcont(j,i) = Fcont(j)
            aFline(j,i) = Fline(j)
            if (indexunits(j).eq. 1) then
               aind(j,i) = widths(j)*(1.0-Fline(j)/Fcont(j))
            else if ( indexunits(j).eq.0 ) then
               aind(j,i) = -2.5*alog10(Fline(j)/Fcont(j))
            else if ( indexunits(j).eq.2) then
               aind(j,i) = Fline(j)/Fcont(j)
            end if
c            print*, 'index, Fline, Fcont = ',labels(j),aind(j,i),
c     &                      Fline(j),Fcont(j)
         end do
C     end of loop over ingredient populations
      end do


C     prepare to sum by clearing arrays
      do j=1,1221
         hlam(j) = 0.0
      end do
      do j=1,nsyn
         synflx(j) = 0.0
      end do
      bolo = 0.0
      massinit = 0.0
      massnow = 0.0
      do j=1,nclrs
         mags(j) = 0
         Lsq(j) = 0.0
         Lavg(j) = 0.0
      end do
      do j=1,nind
         Fcont(j) = 0.0
         Fline(j) = 0.0
      end do

C     loop over histogram bins and sum
      do iz = 1,npop
C        if there are few stars in a bin, skip it.
         if (shist(iz).lt.0.0001) goto 55
         do j=1,1221
            hlam(j) = hlam(j) + shist(iz)*ahlam(j,iz)
         end do
         do j=1,nsyn
            synflx(j) = synflx(j) + shist(iz)*asynflx(j,iz)
         end do
C        note that we go to non-log units
c         print*, iz,bolo, shist(iz),abolo(iz)
         bolo = bolo + shist(iz)*10**(abolo(iz))
         massinit = massinit + shist(iz)*amassinit(iz)
         massnow  = massnow  + shist(iz)*amassnow(iz)
         do j=1,nclrs
            mags(j) = mags(j) +shist(iz)*10**(-0.4*amags(j,iz))
            Lsq(j) = Lsq(j) + shist(iz)*aLsq(j,iz)
            Lavg(j) = Lavg(j) + shist(iz)*aLavg(j,iz)
         end do
         do j=1,nind
            Fcont(j) = Fcont(j) + shist(iz)*aFcont(j,iz)
            Fline(j) = Fline(j) + shist(iz)*aFline(j,iz)
         end do
 55      continue
      end do



C     open some output files
      open(unit=45,file='ask_most.out',status='unknown')
      open(unit=48,file='ask_sed.out',status='unknown')
      open(unit=49,file='ask_spect.out',status='unknown')
      t = ctime(time())
      write(45,'(a,a24)') ' W94 model interpolation output ',t
      write(45,'(x)')
      write(45,'(a)') 'Composite population comprised of: '
      write(45,'(a)') ' Age (Gyr)    [Fe/H]    weight'
      write(48,'(a,a24)') ' W94 model interpolation output ',t
      write(48,'(x)')
      write(48,'(a)') 'Composite population comprised of: '
      write(48,'(a)') ' Age (Gyr)    [Fe/H]    weight'
      write(49,'(a,a24)') ' W94 model interpolation output ',t
      write(49,'(x)')
      write(49,'(a)') 'Composite population comprised of: '
      write(49,'(a)') ' Age (Gyr)    [Fe/H]    weight'
      do i=1,npop
         write(45,'(2f10.3,f10.5)') ageref(i), metalref(i), shist(i)
         write(48,'(2f10.3,f10.5)') ageref(i), metalref(i), shist(i)
         write(49,'(2f10.3,f10.5)') ageref(i), metalref(i), shist(i)
      end do
      write(45,'(x)')
      write(45,'(x)')
      write(48,'(x)')
      write(48,'(x)')
      write(49,'(x)')
      write(49,'(x)')

C     form indices and output them to a file
      do k=1,nind
         if (indexunits(k).eq. 1) then
            index(k) = widths(k)*(1.0-Fline(k)/Fcont(k))
         else if ( indexunits(k).eq.0) then
            index(k) = -2.5*alog10(Fline(k)/Fcont(k))
         else if ( indexunits(k).eq.2) then
            index(k) = Fline(k)/Fcont(k)
         else
            print*, 'Error in indexunits - yikes!'
         end if
      end do
      write(45,'(a)') 'Spectral Indices'
      write(45,'(7(a10,1x))') (labels(k),k=1,7)
      write(45,'(7(2x,f8.3,1x))') (index(k),k=1,7)
      write(45,'(x)')
      write(45,'(7(a10,1x))') (labels(k),k=8,14)
      write(45,'(7(2x,f8.3,1x))') (index(k),k=8,14)
      write(45,'(x)')
      write(45,'(7(a10,1x))') (labels(k),k=15,21)
      write(45,'(7(2x,f8.3,1x))') (index(k),k=15,21)
      write(45,'(x)')
      write(45,'(7(a10,1x))') (labels(k),k=22,25)
      write(45,'(7(2x,f8.3,1x))') (index(k),k=22,25)
      write(45,'(x)')
      write(45,'(a,a)') 'Talk with Andy Leonardi',
     z      ' before using these synthetic Rose indices:'
      write(45,'(3x,7(a10,1x))') 'Hd/4045   ','Hd/4063   ','Sr/4045   ',
     z           'Sr/4063   ','pGband    ','Hg/4325   ','4289/271  '
      write(45,'(7(2x,f8.3,1x))') index(33)/index(30),
     z      index(33)/index(31),index(32)/index(30),index(32)/index(31),
     z      index(48)/index(36),index(38)/index(37),index(35)/index(34)
      write(45,'(x)')
      write(45,'(3x,7(a10,1x))') '4384/4352 ','P[Fe/H]   ','Ca II     ',
     z           '3888/3859 ','p220/p209 ','D(Balmer) ','          '
      write(45,'(6(2x,f8.3,1x))') index(40)/index(39),
     z      (index(41)+index(43))/(2.0*index(42)),
     z      index(29)/index(28),index(27)/index(26),index(45)/index(44),
     z      index(71)/index(70)
      write(45,'(x)')
      write(45,'(x)')

C     form Lbar and output
      do k=1,nclrs
         Lbar(k) = -2.5 * alog10( Lsq(k)/Lavg(k) )
      end do
      write(45,'(a)') ' SBF magnitudes '
      write(45,'(a,a,a)') 
     &      '   Ubar   Bbar   Vbar   Rbar   Ibar',
     &      '   Jbar   Hbar   Kbar   Lbar   Lpbar',
     &       '  Mbar'
      write(45,'(11f7.2)')  
     z         (Lbar(k),k=1,nclrs)
      write(45,'(x)')
      write(45,'(x)')

C     calculate and output mass-to-light ratios
C     divide the mass by the population luminosity in solar units
      do k=1,nclrs
         mags(k) = -2.5*alog10(mags(k))
         mtol(k) = massnow / (10.0**(0.4*(solmags(k)-mags(k))))
      end do
      write(45,'(a)') ' Mass-to-light ratios '
      write(45,'(a,a,a)') 
     &      '  M/L_U  M/L_B  M/L_V  M/L_Rc M/L_Ic',
     &       ' M/L_J  M/L_H  M/L_K  M/L_L  M/L_Lp',
     &       ' M/L_M'
      write(45,'(11f7.2)') (mtol(k),k=1,nclrs)
      write(45,'(x)')
      write(45,'(x)')

C     output magnitudes and colors
      write(45,'(a)') ' Ordinary magnitudes '
      write(45,'(a,a,a)') 
     &      '      U      B      V      Rc     Ic',
     &       '     J      H      K      L      Lp',
     &       '     M'
      write(45,'(11f7.2)') (mags(k),k=1,nclrs)
      write(45,'(x)')
      write(45,'(x)')
      write(45,'(a)') 'Some colors'
      write(45,'(a,a)') 
     &      '     U-V     B-V     V-Rc    V-Ic    V-K',
     &      '     J-K     H-K'
      write(45,'(7f8.2)') mags(1)-mags(3),mags(2)-mags(3),
     &  mags(3)-mags(4),mags(3)-mags(5),mags(3)-mags(8),
     &  mags(6)-mags(8),mags(7)-mags(8)
      
      write(45,'(x)')
      write(45,'(a,f6.2,a)') 'B-V =',mags(2)-mags(3)-0.05,
     &    ' is probably more accurate. See W94.'
      write(45,'(x)')
      write(45,'(x)')
      write(45,'(a,f9.0)') 'Initial population mass   : ',massinit
      write(45,'(a,f9.0)') 'Mass in stars and remnants: ',massnow
      write(45,'(a,f7.3)') 'log Lbol/Lbol_sun         : ',alog10(bolo)
      if (imftype(2:4).eq.'SAL' .or. imftype(2:4).eq.'sal') then
         write(45,'(a,f7.3)') 'IMF power-law slope       : ',imfpar
      else if (imftype(2:4).eq.'MIL'.or.imftype(2:4).eq.'mil') then
         write(45,'(a)') 'We have a Miller-Scalo (1979) piecewise IMF'
      else
         write(45,'(a)') 'I do not recognise this IMF type!'
      end if
      write(45,'(a,f7.3)') 'IMF lower mass cut        : ',imflow
      write(45,'(a,f7.3)') 'IMF upper mass cut        : ',imfupp
      write(45,'(x)')


C     output SED
      write(48,'(a)')'Spectral energy distribution in wavelength units.'
      write(48,'(x)')
      write(48,'(a)')'Normalized so that the integral over wavelength'
      write(48,'(a)')'in nm is equal to the bolometric luminosity.'
      write(48,'(a)')'Interpolation error can cause up to 1% mismatch.'
      write(48,'(x)')
      write(48,'(a)')'Column 1: wavelength in Angstroms '
      write(48,'(a)')'Column 2: log wavelength (in Angstroms) '
      write(48,'(a)')'Column 3: Flux in wavelength units'
      write(48,'(a)')'Column 4: log of column 3'
      write(48,'(x)')
      write(48,'(x)')
      do j=1,1221
         if ( hlam(j).gt.0.0 ) then
            x1 = alog10(hlam(j))
         else
            x1 = 0.0
         end if
         write(48,'(4e14.5)') 10.0*wave(j),alog10(wave(j))+1.0,
     &        hlam(j),x1
      end do

C     write synthetic spectrum to a file
      write(49,'(a)') '3500-5500 Angstrom Spectrum'
      write(49,'(a)') ' '
      write(49,'(a)') ' Lambda     F_lambda'
      do j=1,nsyn
         write(49,'(f9.2,e14.5)') synwave(j),synflx(j)
      end do

C     end matter
      print*, ' '
      print*, 
     & 'Output: see "ask_most.out", "ask_sed.out", and "ask_spect.out"'
      print*, ' '

      close(45)
      close(48)
      close(49)

      return
      end
