      program wwwpad

C                  histpar is the number of bins in the [Fe/H] histogram
      integer      histpar
      parameter    (histpar = 28)
c      real         metalref(histpar),ageref(histpar),shist(histpar)
c      character*80 text
      real        ages(histpar),metals(histpar),shist(histpar)
      integer     status, npop
      character dirname*80, sline*40

C     interpret user input (stored in the temporary file 'tmp_input')
C     example: "12 0. 100. 
C              "DISP
C              "Wo__vanilla
C            : this means the user wants a 12 Gyr, [Fe/H]=0, pop with
C            : a relative weight of 100. 
C            : Then there may be one or two carriage returns.
C            : DISP is "display results" SEND is "send results in a file"
C            : 
C            : Finally, the "Wo__vanilla" string is a directory name

      open(unit=33,file='/home/www/root/dial/tmp_pad_input',
     &     status='old',err=6)
      goto 7
 6    status = 5
      print*, 'Could not find crucial input file tmpinput'
      open(unit=34,file='/home/www/root/dial/tmp_pad_flag',
     z    status='unknown')
      write(34,'(i2)') status
      close(34)
      stop

C     store directory string in 'dirname'
 7    npop = 0
      do i=1,34
         read(33,'(a40)',end=8) sline
         if (sline(1:4).eq.'DISP') then
            status = 1
            goto 15
         end if
         if (sline(1:4).eq.'SEND') then
            status = 2
            goto 15
         end if
         if (sline(1:1).eq.'P') then
            read(sline,'(a11)') dirname
            goto 8
         end if
         if (sline.eq.'') goto 15
C        I think we cleared all of the non-numeric cases. Now read numbers.
         npop = npop + 1
         read(sline,*,err=16) ages(npop),metals(npop),shist(npop)
         if (ages(npop).le.0.0 .or. shist(npop).le.0.0) then
            goto 16
         end if
 15   end do
 16   status = 5
 8    close(33)
      print*, 'Status = ',status

C     if things are already sour, abort
      if (status.eq.5) then
         open(unit=34,
     z   file='/home/www/root/dial/tmp_pad_flag',
     z        status='unknown')
         write(34,'(i2)') status
         close(34)
         stop
      end if

C     if things are OK, call subroutine
      dirname = '/home/www/root/dial/exportpad2/'//dirname
      nchar = 42
      print*, dirname
      call ask (metals,ages,shist,npop,status,dirname,nchar)

c      print*, 'npop = ',npop
      print*, 'status = ',status
c      do i=1,npop
c         print*, i,ages(i),metals(i),shist(i)
c      end do

C     write status flag to alert C-shell script
      open(unit=34,file='/home/www/root/dial/tmp_pad_flag',
     z    status='unknown')
      write(34,'(i2)') status
      close(34)

      stop
      end



      subroutine ask(metalref,ageref,shist,npop,status,dirname,nchar)
C     noninteractively 
C     interpolates in a model grid and combines populations.

C                  nind is the number of indices
C                  histpar is the max number of populations to combine
C                  npop is the number of populations to combine
C                  nclrs is the number of colors
      integer      nind,histpar,status,npop,nclrs
      parameter    (nind = 71, histpar = 28, nclrs=11)
      real         metalref(histpar),ageref(histpar),shist(histpar)
C     dirname is the directory for the model grid, nchar how many characters
C     in the directory string.
      character    dirname*80
      integer      nchar

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         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)
      integer      showindex

C     other variables
      character*24   ctime,t
      integer        time
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     do not output Lick indices if there is a young age pop request
      showindex = 0
      do j=1,npop
         if (ageref(j).lt.0.4) showindex = -1
      end do
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

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)

         if (status.eq.5) then
            open(unit=34,
     &         file='/home/www/root/dial/tmp_pad_flag',status='unknown')
            write(34,'(i2)') status
            close(34)
            stop
         end if

C        store important results
         do j=1,1221
            ahlam(j,i) = hlam(j)
         end do
         abolo(i) = bolo
         amassinit(i) = massinit
         amassnow(i) = massnow
c         print*, ' bolo, mass = ',bolo,mass
         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
               aind(j,i) = -2.5*alog10(Fline(j)/Fcont(j))
            end if
         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
      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
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 output file
      open(unit=45,
     z file='/home/www/root/dial/wwwpad_.out',status='unknown')
      t = ctime(time())
      write(45,'(a,a24)') 
     & ' W94 (with Padova isos) model interpolation output ',t
      write(45,'(x)')
      write(45,'(a)') 'Composite population comprised of: '
      write(45,'(a)') ' Age (Gyr)    [Fe/H]    weight'
      do i=1,npop
         write(45,'(2f10.3,f10.5)') ageref(i), metalref(i), shist(i)
      end do
      write(45,'(x)')
      write(45,'(x)')

      if (showindex.eq.0) then
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
            index(k) = -2.5*alog10(Fline(k)/Fcont(k))
         end if
      end do
      write(45,'(a)') 'Spectral Indices'
      write(45,'(7(a10,1x))') (labels(k),k=1,7)
      write(45,'(7(2x,f8.4,1x))') (index(k),k=1,7)
      write(45,'(x)')
      write(45,'(7(a10,1x))') (labels(k),k=8,14)
      write(45,'(7(2x,f8.4,1x))') (index(k),k=8,14)
      write(45,'(x)')
      write(45,'(7(a10,1x))') (labels(k),k=15,21)
      write(45,'(7(2x,f8.4,1x))') (index(k),k=15,21)
      write(45,'(x)')
      write(45,'(7(a10,1x))') (labels(k),k=22,25)
      write(45,'(7(2x,f8.4,1x))') (index(k),k=22,25)
      write(45,'(x)')
      write(45,'(x)')
      else
      write(45,'(x)')
      write(45,'(a)') 
     & 'Due to young age component, no Lick indices are output'        
      write(45,'(x)')         
      end if
C     form Lbar and output

      do k=1,11
         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.3)')  
     z         (Lbar(k),k=1,11)
      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.3)') 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.3,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
      if (status.eq.2) then
      write(45,'(a)')'Spectral energy distribution in wavelength units.'
      write(45,'(x)')
      write(45,'(a)')'Normalized so that the integral over wavelength'
      write(45,'(a)')'in nm is equal to the bolometric luminosity.'
      write(45,'(a)')'Interpolation error can cause up to 1% mismatch.'
      write(45,'(x)')
      write(45,'(a)')'Column 1: wavelength in Angstroms '
      write(45,'(a)')'Column 2: log wavelength (in Angstroms) '
      write(45,'(a)')'Column 3: Flux in wavelength units'
      write(45,'(a)')'Column 4: log of column 3'
      write(45,'(x)')
      write(45,'(x)')
      do j=1,1221
         if ( hlam(j).gt.0.0 ) then
            x1 = alog10(hlam(j))
         else
            x1 = 0.0
         end if
         write(45,'(4e14.5)') 10.0*wave(j),alog10(wave(j))+1.0,
     &        hlam(j),x1
      end do
      end if

      close(45)

      return
      end


