      module colteff_subs

      use bessell
      use ken_hart
      use phoenix_sloan
      use phoenix_bessell

      implicit none

      ! The maximum number of points in the model isochrones.
      integer, parameter :: mmod=1000

      type isochron
        real, dimension(mmod) :: mass
        ! The colour or temperature.
        real, dimension(mmod) :: colour
        ! The absoloute magnitude or luminosity.
        real, dimension(mmod) :: magnitude
        real, dimension(mmod) :: logg
        real :: age
        integer :: npts
      end type isochron

      ! These parameters convert from colflag into the actual names for
      ! the colours in question.  It gives us a good way of naming the
      ! files.
      ! Isochrones > 20 use their own colours and bolometric corrections.
      integer, public, parameter :: n_names=50
      character(len=3), dimension(n_names), parameter, public :: 
     &col_name=(/
     &'V-I','R-I','I-K','B-V','B-V','   ','   ','   ','   ','   ',
     &'   ','   ','   ','B-V','B-V','   ','   ','   ','   ','   ',
     &'V-I','R-I','I-K','B-V','B-V','   ','   ','   ','   ','   ',
     &'   ','   ','   ','B-V','   ','   ','   ','   ','   ','   ',
     &'V-I','   ','   ','   ','   ','   ','   ','   ','   ','   '/)
      character(len=3), dimension(n_names), public, parameter :: 
     &mag_name=(/
     &'V  ','I  ','I  ','V  ','U-B','   ','   ','   ','   ','   ',
     &'   ','   ','   ','V  ','U-B','   ','   ','   ','   ','   ',
     &'V  ','I  ','I  ','V  ','U-B','   ','   ','   ','   ','   ',
     &'   ','   ','   ','V  ','   ','   ','   ','   ','   ','   ',
     &'V  ','   ','   ','   ','   ','   ','   ','   ','   ','   '/)
      character(len=40), dimension(n_names), parameter, public ::
     &pht_comm=(/
     &'Jeffries/Naylor                         ', 
     &'Jeffries/Naylor                         ', 
     &'Jeffries/Naylor                         ', 
     &'Jeffries/Flower                         ', 
     &'Jeffries/Jeffries                       ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
c
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'Bessell (solar abundance)               ', 
     &'Bessell (solar abundance)               ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
c 
     &'As supplied with isochrone.             ', 
     &'As supplied with isochrone.             ', 
     &'As supplied with isochrone.             ', 
     &'As supplied with isochrone.             ', 
     &'As supplied with isochrone.             ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ',
c 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'Bessell Tycho  (solar abundance)        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ',
c 
     &'Kenyon & Hartmann (1995).               ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        '/)
      ! Similarly this parameter allows us to convert from modflag to
      ! a name for the model.
      integer, public, parameter :: n_models=30
      character(len=13), dimension(0:n_models), parameter, public ::
     &mod_name=(/
     &'linear       ', 'dam97        ', '             ', ! 'burrows      ', 
    !&'baraffe      ', 'siessz02     ', 'newbaraffe   ', 
     &'             ', 'siessz02     ', 'newbaraffe   ', 
     &'newbaraffe_19', 'siessz01     ', 'padova       ',
     &'geneva       ', 'user         ', '             ', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             '/)
      
      interface reddening
        module procedure reddening_point
        module procedure reddening_array
      end interface

      contains
      
      character(len=30) function data_dir()
      
      ! A little function to return where the datafiles are stored.
      character(len=30) :: here
      
      call get_environment_variable('CMDDATA', here)
      data_dir=here
      
      end function data_dir
      
      subroutine colteff_calc(
     &empiso, dmod, red, age, colflag, modflag, npt, newcol, temp)
      
      ! empirical isochrone 
      character(len=*), intent(in) ::  empiso                
      ! Calibration cluster distance modulus, reddening and age.
      real, intent(in) :: dmod, red, age
      ! What colours do you want? 1 for V-I, 2 for R-I, 3 for I-K, 4 for B-V.
      integer, intent(in) ::  colflag
      ! Which model do you want?
      ! 1 for DAM, 2 for Burrows, 3 for Baraffe
      ! 4 for Siess z=0.02, 5 for new Baraffe,
      ! 6 for new Baraffe with mixing length 1.9
      ! 7 for Siess with z=0.01
      integer, intent(in) :: modflag
      
      ! no. of points in empirical isochrone
      integer, intent(out) :: npt
      ! Teff corresponding to colours for empirical isochrone.
      real, dimension(:), allocatable, intent(out) :: newcol, temp
      

      real dtemp
      
      real lum,bc               ! luminosity/bolometric corr. at a point
      ! Input colours and magnitudes.      
      real, allocatable, dimension(:) ::  col, mag
      integer :: nempiso
      real, dimension(mmod) :: lbol, teff, mass, logg    ! theoretical isochrone
      integer :: nmod                        ! no. of points in theoretical isochrone
      ! Working versions of output arrays
      real, dimension(:), allocatable :: w_newcol, w_temp
      integer i, j, k, iostat
      integer :: bcflag
      
1     format(1x,a)

      ! Read the empirical isochrone.

      open (unit=10, file=empiso, status='old')
      nempiso=0
      count_iso: do
        read(10,*,iostat=iostat)
        if (iostat < 0) exit count_iso
        nempiso=nempiso+1
      end do count_iso
      rewind(10)
      allocate(col(nempiso), mag(nempiso))
      do i=1, nempiso
        read(10,*) col(i), mag(i)
      end do
      close(10)
      
      ! Reverse the effects of distance and reddening on the CMD.
      call reddening(colflag, -dmod, -red, col, mag)

cnow read in the theoretical isochrones
c a different procedure may need to be adopted for each set of models
c so have several calls
c expects lbol,teff to be returned in log units, age to be given in Myr 
c now also returns the mass of each point in solar units

      if (modflag.eq.1) then
        call read_dam97(age,lbol,teff,mass,nmod)
        logg=4.5
      else if (modflag.eq.2) then
        call read_burrows(age,lbol,teff,mass,nmod)
        logg=4.0
      else if (modflag.eq.3) then
        call read_baraffe(age,lbol,teff,mass,nmod)
        logg=4.0
      else if (modflag.eq.4) then
        call read_siessz02(age,lbol,teff,mass,logg,nmod)
      else if (modflag.eq.5) then
        call read_newbaraffe(age,lbol,teff,mass,nmod)
        logg=4.0
      else if (modflag.eq.6) then
        call read_newbaraffe_19(age,lbol,teff,mass,nmod)
        logg=4.0
      else if (modflag.eq.7) then
        call read_siessz01(age,lbol,teff,mass,logg,nmod)
      else if (modflag.eq.9) then
        call read_padgen(modflag, age,lbol,teff,mass,logg,nmod)
      else
        print*, 'Unkown model flag ', modflag
        stop
      end if


      print*,'Lbol-Teff of calibration cluster'
      print*,'Lbol            Teff'
      print*,'--------------------'
      do i=1,nmod
        print*,10.0**lbol(i),10.0**teff(i)
      end do

c loop through the empirical isochrone
c for each point calculate the lbol and then the Teff corresponding to this.
c But need to trap instances where the empirical data lies well beyond the range 
c of the model

       
      ! The largest the output arrays can be is nempiso.
      allocate(w_temp(nempiso), w_newcol(nempiso))
      npt=0
      create: do i=1, nempiso
        if (colflag==1) bc = bcvi(col(i), logg(i), teff(i), bcflag)
        if (colflag==2) bc = bcri(col(i), bcflag)
        if (colflag==3) bc = bcik(col(i), bcflag)
        if (colflag==4) bc = bcbv(col(i), bcflag)
         
        if (bcflag /= 0) cycle create

        lum = (4.74 - bc - mag(i))/2.5

c       print*,'lum,B-V,BC=',lum,col(i),bc

c use numerical recipes subroutine - returns j, where lum is between
c lbol(j) and lbol(j+1)

        call locate(lbol,nmod,lum,j)

c       print*, 'index ',j

c this does not allow any extrapolation, but might be good to have a little bit.

        if ((j.gt.0.and.j.lt.nmod).or.
     :      (j.ge.nmod.and.(lum-lbol(nmod)).lt.0.3)) then

          npt=npt+1
c piece of black magic from numerical recipes p.113

          k = min(max(j-1,1),nmod-2)

c now interpolate (parabolic so n=3), returns temp and dtemp (error)

          call polint(lbol(k),teff(k),3,lum,w_temp(npt),dtemp)

          w_newcol(npt)=col(i)

c         print*,newcol(npt),10.0**temp(npt)

          end if

        end do create
        
        ! Now pack this into arrays which are just the right length.
        allocate(newcol(npt), temp(npt))
        newcol=w_newcol(1:npt)
        temp=w_temp(1:npt)
        deallocate(w_newcol, w_temp, col, mag)
      
c       print*,'npt=',npt

      end subroutine colteff_calc

      subroutine iso_calc(cage, colflag, modflag,
     &  npt, newcol, temp, nmod, isocol, colflg, isomag, magflg, 
     &  isomass, massflg, isoteff, teffflg)
      
        ! Given an isochrone model number and colour flag this routine
        ! interpolates the isochrone to a given age, and then pastes
        ! a model atmosphere on top.

        ! The age to be modelled.
        real, intent(in) :: cage
        ! What colours do you want? 
        integer, intent(in) ::  colflag
        ! Which model do you want?
        integer, intent(in) :: modflag

        ! If the effective temperature colour relationship is to be
        ! tweaked, you will also need the following.
        ! no. of points in empirical isochrone
        integer, intent(in) :: npt
        ! log10(teff) corresponding to colours for empirical isochrone.
        real, dimension(:), intent(in) :: newcol, temp      
        
        ! Number of points in output theoretical isochrone.
        integer, intent(out) :: nmod                        
        ! Colour, magnitude and mass of output isochrone.
        real, dimension(:), allocatable, intent(out) :: 
     &  isocol, isomag, isomass, isoteff
        ! The flag for the isochrone.
        character(len=*), allocatable, dimension(:), intent(out) :: 
     &  colflg, magflg, massflg, teffflg
        
        integer :: imod, i, j, k, iflag
        real :: bc               ! luminosity/bolometric corr. at a point
        real, dimension(mmod) :: lbol, teff, mass, logg  ! theoretical isochrone
        integer :: bcflag
        character(len=30) :: ifname

                                ! Read the appropriate models
        
        if (modflag.eq.1) then
          call read_dam97(cage,lbol,teff,mass,nmod)
          logg=4.0
        else if (modflag.eq.2) then
          call read_burrows(cage,lbol,teff,mass,nmod)
          logg=4.0
        else if (modflag.eq.3) then
          call read_baraffe(cage,lbol,teff,mass,nmod)
          logg=4.0
        else if (modflag.eq.4) then
          call read_siessz02(cage,lbol,teff,mass,logg,nmod)
        else if (modflag.eq.5) then
          call read_newbaraffe(cage,lbol,teff,mass,nmod)
          logg=4.0
        else if (modflag.eq.6) then
          call read_newbaraffe_19(cage,lbol,teff,mass,nmod)
          logg=4.0
        else if (modflag.eq.7) then
          call read_siessz01(cage,lbol,teff,mass,logg,nmod)
        else if (modflag.eq.8) then
          call read_padgen(modflag, cage,lbol,teff,mass,logg,nmod)
        else if (modflag.eq.9) then
          call read_padgen(modflag, cage,lbol,teff,mass,logg,nmod)
        else if (modflag.eq.10) then
          ifname='user_00000.00.dat'
          write(ifname(6:13), '(f8.2)') cage
          do i=6, 13
            if (ifname(i:i) == ' ') ifname(i:i)='0'
          end do
          call read_one_geneva(ifname, 0, nmod, lbol, teff, mass, logg)
        else
          print*, 'Unknown model flag ', modflag
          stop
        end if

        ! Only use isochrone points within the range of the colour-Teff
        ! relationship.  Rob's version used to allow extrapolation, but 
        ! if the extrapolation was so bad it became double valued it 
        ! confused Tim's fitting procedures.
        allocate(isocol(nmod), isomag(nmod), isomass(nmod), 
     &  isoteff(nmod))
        allocate(colflg(nmod), magflg(nmod), massflg(nmod),
     &  teffflg(nmod))
        
        ! lbol and teff derived for that isochrone, now convert to colour 
        !and mag
        imod=0
        colflg='OK'
        magflg='OK'
        massflg='OK'
        teffflg='OK'
        create: do imod=1, nmod
c          print*, 'Point ', imod, 'Lum', lbol(imod), 'teff ',
c     &    10.0**teff(imod), 'logg ', logg(imod), 'mass ', mass(imod)
          ! First deal with the mass and effective temperature.
          isomass(imod) = mass(imod)
          isoteff(imod) = teff(imod)
          ! Now the colour.
          if (colflag < 11) then
            isocol(imod)=quadint(temp, newcol, teff(imod), colflg(imod))
          else if (colflag == 11) then
            isocol(imod)=bessell_col_vi(logg(imod), teff(imod), iflag)
            if (iflag /= 0) colflg(imod)='1'
          else if (colflag==14 .or. colflag==15) then
            isocol(imod)=bessell_col_bv(logg(imod), teff(imod), iflag)
            if (iflag /= 0) colflg(imod)='1'
          else if (colflag == 34) then
            isocol(imod)=bessell_col_bv(logg(imod), teff(imod), iflag)
            call tycho_col(isocol(imod), colflg(imod))
            if (iflag /= 0) colflg(imod)='1'
         else if (colflag == 41) then
            isocol(imod)=ken_hart_col_vi(teff(imod), colflg(imod))
          else if (colflag == 51) then
            isocol(imod)=phoenix_col_gi(logg(imod), teff(imod), 
     &      colflg(imod))
          else if (colflag == 52) then
            isocol(imod)=phoenix_col_iz(logg(imod), teff(imod), 
     &      colflg(imod))
          else if (colflag == 61 .or. colflag == 66) then
            isocol(imod)=phoenix_col_vi(logg(imod), teff(imod),
     &      colflg(imod))
          else if (colflag == 64) then
            isocol(imod)=phoenix_col_bv(logg(imod), teff(imod),
     &      colflg(imod))
          else
            print*, 'The requested colour ', colflag,
     &      ' is not available' 
          end if
c          print*, 'Colour is ', isocol(imod), colflg(imod)
          if (colflg(imod) == '1') then
            colflg(imod)='Colour-Teff'
            !cycle create
          end if
          ! Now the magnitude, but note that for U-B vs B-V this is 
          ! a colour.
          if (colflag == 15) then
            isomag(imod)=bessell_col_ub(logg(imod), teff(imod), iflag)
            if (iflag /= 0) then
              magflg(imod)='Colour-Teff in U-B'
            end if
          else
            if (colflag==1) then
              bc = bcvi(isocol(imod), logg(imod), teff(imod), bcflag)
            else if (colflag==2) then
              bc = bcri(isocol(imod), bcflag)
            else if (colflag==3) then
              bc = bcik(isocol(imod), bcflag)
            else if (colflag==4) then
              bc = bcbv(isocol(imod), bcflag) 
            else if (colflag==11 .or. colflag==14) then
              bc=bessell_bc_v(logg(imod), teff(imod), bcflag)
            else if (colflag==34) then
              bc=bessell_bc_v(logg(imod), teff(imod), iflag)
              bc=bc+tycho_bc(isocol(imod), bcflag)
              if (iflag /= 0) bcflag=1
            else if (colflag==41) then
              bc=ken_hart_bc_v(teff(imod), bcflag)
            else if (colflag==51) then
              bc=phoenix_bc_g(logg(imod), teff(imod), magflg(imod))
            else if (colflag==52) then
              bc=phoenix_bc_i(logg(imod), teff(imod), magflg(imod))
            else if (colflag == 61 .or. colflag == 64) then
              bc=phoenix_blc_v(logg(imod), teff(imod), magflg(imod))
            else if (colflag == 66) then
              bc=phoenix_blc_i(logg(imod), teff(imod), magflg(imod))
            else
              print*, 'Cannot cope with colour', colflag
              stop
            end if
            !print*, 'BC is ', bc, magflg(imod)
            if (bcflag/=0 .and. colflag<51) then
              magflg(imod)='Bolometric correction'
              !cycle create
            end if
            if (magflg(imod)=='1' .and. colflag>50) then
              magflg(imod)='Bolometric correction'
              !cycle create
            end if
            isomag(imod)=4.74-bc-2.5*lbol(imod)
          end if

c          print*, isomag(imod), trim(colflg(imod)), 
c     &    isocol(imod), trim(magflg(imod))
          
        end do create

      end subroutine iso_calc



c *****************************************************************************
      subroutine reddening_array(colflag, dmod, red, col, mag)
c *****************************************************************************
      
        ! Changes an absolute magnitude sequence into a reddened one at
        ! at given distance reddening (dmod, col positve) or the reverse.
    
        ! Note that red is the reddening in the colour (e.g. B-V).
    
        integer, intent(in) :: colflag
        real, intent(in) :: dmod, red
        real, dimension(:), intent(inout) :: col, mag

        integer :: icol
        logical, save :: tycho_warn=.false.
      
        ! Change to un-reddened colour (order important for Bessell).
        if (red < 0.0) col=col+red
        
        ! Reverse the extinction appropriate for each CMD
        if (colflag==1 .or. colflag==21 .or. colflag==41 .or.
     &       colflag==61) then 
          ! Rob's code had 1.8 in here (which was a bug).
          ! Av=2.40E(V-I)
          mag=mag+2.40*red
        else if (colflag==2 .or. colflag==22) then
          ! Implication here is that E(R-I)=0.23Av
          ! Tim re-calculated this, and got the same answer by the following
          ! route. E(V-I)/B(B-V) is about 1.250 (Dean et al 1978), 
          ! Av/E(B-V)=3.1 (trad.) and E(R-I)/E(B-V) is about 0.72 (Taylor,
          ! 1986).  Magically, this gave the answer already in the code!
          ! Ai=2.57E(R-I)
          mag=mag+2.57*red
        else if (colflag==3 .or. colflag==23) then 
          ! Ai=1.2E(I-K)
          mag=mag+1.20*red
        else if (colflag==4 .or. colflag==24 .or. colflag==64) then
          ! Av=3.1E(B-V)
          mag=mag+3.10*red
        else if (colflag==5 .or. colflag==25) then
          ! E(U-B)=0.73*E(B-V) (Nathan says everyone uses this.)
          mag=mag+0.73*red
        else if (colflag == 11) then
          if (abs(red) > tiny(red)) then
            print*, 'Cant deal with reddening in Bessell V-I.'
            stop
          end if
        else if (colflag == 14) then
          ! Av=(3.26+0.22*(B_V)o)*E(B-V), for Bessell.
          mag=mag+(3.26+0.22*col)*red
          ! Special for Cep OB3b.  Works for B-V < 0.2 to about 0.01 mags.
          ! mag=mag+(3.0802126+0.20603374*col)*red
        else if (colflag == 15) then
          ! E(U-B)=(0.71+0.24*(B_V)o)*E(B-V), for Bessell.
          mag=mag+(0.71+0.24*col)*red
          ! Special for Cep OB3b.  For B-V < 0 works to about 0.01 mags.
          !mag=mag+(0.6475283-0.170563*col)*red
        else if (colflag == 34) then
          ! Colour 34 is Tycho.
          ! Vector from Mayne & Naylor.
          do icol=1, size(col,1)
            if (col(icol) < 0.065) then
              mag = mag + (3.358 + 0.237*col(icol))*red
            else if (col(icol) < 0.5) then
              mag = mag + (3.387 - 0.207*col(icol))*red
            else
              if (.not. tycho_warn) then
                print*, 'Colour B-V = ', col(icol), ' is outside'
                print*, 'the range the Tycho reddening vector is'
                print*, 'defined (B-V=0.5).'
                tycho_warn=.true.
              end if
              mag = mag + (3.515 - 0.267*0.5)*red
            end if
          end do
        else if (colflag == 51) then
           ! Ag=2.26*E(g-i), where 2.26 is average value for 3000 - 10000 K.
           mag = mag + 2.26*red
        else if (colflag == 52) then
           ! Ai=3.65*E(i-z), where 3.65 is average value for 3000 - 10000 K.
           mag = mag + 3.65*red
        else
          print*, 'Error, unknown colour flag', colflag, 
     &    ' in reddening_array.'
        end if

        ! And finally deal with the distance modulus.
        mag=mag+dmod
      
        ! Change to un-reddened colour (order important for Bessell).
        if (red > 0.0) col=col+red

      end subroutine reddening_array

c *****************************************************************************
      subroutine reddening_point(colflag, dmod, red, col, mag)
c *****************************************************************************

        integer, intent(in) :: colflag
        real, intent(in) :: dmod, red
        real, intent(inout) :: col, mag
        
        real, dimension(1) :: a_col, a_mag
        
        a_col(1)=col; a_mag(1)=mag
        
        call reddening_array(colflag, dmod, red, a_col, a_mag)

        col=a_col(1); mag=a_mag(1)
      
      end subroutine reddening_point  

c *****************************************************************************
      subroutine read_dam97(mage,lbol,teff,mass,nmod)
c *****************************************************************************

      implicit none
c import
      real mage
c export
      integer nmod
      integer nmax1
      parameter (nmax1=40)
      real lbol(nmax1),teff(nmax1),mass(nmax1)

c local
      integer nmax2,i,j,k,index
      parameter (nmax2=220)
      integer n(nmax1)         ! no. of points in each model
      
      character*32 junk

      real dlbol,dteff,age
      real logAge(nmax1,nmax2),logL(nmax1,nmax2),logT(nmax1,nmax2)
      real dummyage(nmax2),dummyl(nmax2),dummyt(nmax2)
      real dmass(nmax1)

2     format(7x,f7.4,2x,f7.4,3x,f6.4)

      data n/94,104,128,161,134,158,79,82,84,87,88,92,84,79,74,
     :      64,63,65,67,69,69,70,69,72,72,76,83,92,106,139,176,
     :     188,186,185,184,182,220,218,214,212/
      data dmass/0.017,0.018,0.02,0.03,0.04,0.045,0.05,0.055,0.06,
     :           0.065,0.07,0.075,0.08,0.085,0.09,0.10,0.12,0.14,0.16,
     :           0.18,0.2,0.25,0.3,0.35,0.4,0.5,0.6,0.7,0.8,0.9,1.0,
     :           1.1,1.2,1.3,1.4,1.5,1.7,2.0,2.5,3.0/


      age=mage
      nmod=nmax1
      age = log10(age)+6.0

c      print*,'age=',age

c tracks are all in one big file
c idea is to read in a file for each mass and interpolate to get
c log L and log T at the appropriate log(age)
      
      open (1,file=trim(data_dir())//'/dam97.dat',form='formatted',
     &access='sequential', status='old')


      do i=1,nmod
c         print*,'n(i)=',n(i)

         read(1,'(A)') junk   
         read(1,'(A)') junk         
         read(1,'(A)') junk
         read(1,2)
     :      (logAge(i,j),logL(i,j),logT(i,j), j=1,n(i))

c         print*,logAge(i,1),logL(i,1),logT(i,1)
c         print*,logAge(i,94),logL(i,94),logT(i,94)
      enddo

c now have the models in the 2d grids
c for each row in the grid, need to find the index closest to the
c required age

      index=0

      do i=1,nmod

         do j=1,n(i)
            dummyage(j)=logAge(i,j)
         if (i.eq.13) then
c            print*,j,dummyage(j)
         endif
         enddo




c use numerical recipes subroutine - returns j, where age is between
c logAge(i,j) and logAge(i,j+1)

         call locate(dummyage,n(i),age,j)

c         if (j.gt.1.and.j.lt.n(i)-1) then
         if (j.gt.0.and.j.lt.n(i)) then

c piece of black magic from numerical recipes p.113

            index=index+1

            k = min(max(j-1,1),n(i)-2)

         
c now snip out the bits of the arrays we want into dummy arrays

            do j=1,3
               dummyage(j)=logAge(i,k+j-1)
               dummyl(j)=logL(i,k+j-1)
               dummyt(j)=logT(i,k+j-1)
            enddo

c now interpolate (parabolic so n=3), returns temp and dtemp (error)

            call polint(dummyage,dummyl,3,age,lbol(index),dlbol)
            call polint(dummyage,dummyt,3,age,teff(index),dteff)

            mass(index)=dmass(i)

         endif

      enddo

      nmod=index

      close(1)

      end subroutine read_dam97



c*****************************************************************************
      subroutine read_burrows(mage,lbol,teff,mass,nmod)
c*****************************************************************************

      implicit none
c import
      real mage
c export
      integer nmod
      real lbol(100),teff(100),mass(100)


      ! Nothing is set by this routine, but lets at least initialise
      ! the outputs just in case its accessed.

      lbol=0.0
      teff=0.0
      mass=0.0
      nmod=0

      end subroutine read_burrows

c*****************************************************************************
      subroutine read_baraffe(mage,lbol,teff,mass,nmod)
c*****************************************************************************

      implicit none
c import
      real mage          ! in Myr
c export
      integer nmod
      integer nmax1
      parameter (nmax1=25)
      real lbol(nmax1),teff(nmax1),mass(nmax1)   ! need to be in log units

c local
      integer nmax2,i,j,k,index
      parameter (nmax2=53)
      integer n(nmax1)         ! no. of points in each model
      
      character*32 junk

      real dlbol,dteff,age
      real logAge(nmax1,nmax2),logL(nmax1,nmax2),logT(nmax1,nmax2)
      real dmass(nmax1)
      real dummyage(nmax2),dummyl(nmax2),dummyt(nmax2)

2     format(1x,f5.3,2x,f7.5,1x,f6.1,8x,f5.2)

      data n/49,53,53,53,53,53,53,53,53,53,53,53,53,53,53,
     :       53,53,53,53,53,53,52,50,49,49/

      age=mage
      nmod=nmax1
      age = log10(age)+6.0

c      print*,'age=',age

c tracks are all in one big file
c idea is to read in a file for each mass and interpolate to get
c log L and log T at the appropriate log(age)

      open (1,file=trim(data_dir())//'/oldbaraffe98.dat',
     &form='formatted',access='sequential',status='old')

      read(1,'(A)') junk   
      read(1,'(A)') junk         
      read(1,'(A)') junk


      do i=1,nmod
c         print*,'n(i)=',n(i)

         read(1,2)
     :      (dmass(i),logAge(i,j),logT(i,j),logL(i,j), j=1,n(i))

         do j=1,n(i)

c convert to logs

            logAge(i,j)=log10(logAge(i,j))+9.0
            logT(i,j)=log10(logT(i,j))
         
         enddo

c          print*,dmass(i),logAge(i,1),logL(i,1),logT(i,1)
c          print*,dmass(i),logAge(i,40),logL(i,40),logT(i,40)
      enddo

c now have the models in the 2d grids
c for each row in the grid, need to find the index closest to the
c required age

      index=0

      do i=1,nmod

         do j=1,n(i)
            dummyage(j)=logAge(i,j)
c         if (i.eq.13) then
c            print*,j,dummyage(j)
c         endif
         enddo




c use numerical recipes subroutine - returns j, where age is between
c logAge(i,j) and logAge(i,j+1)

         call locate(dummyage,n(i),age,j)

c         if (j.gt.1.and.j.lt.n(i)-1) then
          if (j.gt.0.and.j.lt.n(i)) then
c piece of black magic from numerical recipes p.113

            index=index+1

            k = min(max(j-1,1),n(i)-2)

         
c now snip out the bits of the arrays we want into dummy arrays

            do j=1,3
               dummyage(j)=logAge(i,k+j-1)
               dummyl(j)=logL(i,k+j-1)
               dummyt(j)=logT(i,k+j-1)
            enddo

c now interpolate (parabolic so n=3), returns temp and dtemp (error)

            call polint(dummyage,dummyl,3,age,lbol(index),dlbol)
            call polint(dummyage,dummyt,3,age,teff(index),dteff)

            mass(index)=dmass(i)
   
         endif

      enddo

      nmod=index

      close(1)

      end subroutine read_baraffe

c*****************************************************************************


c *****************************************************************************
      subroutine read_siessz02(mage,lbol,teff,mass,logg,nmod)
c *****************************************************************************

      implicit none
c import
      real mage
c export
      integer nmod
      integer nmax1
      parameter (nmax1=32)
      real lbol(nmax1),teff(nmax1),mass(nmax1),logg(nmax1)

c local
      integer nmax2,i,j,k,index
      parameter (nmax2=900)
      integer n(nmax1)         ! no. of points in each model
      
      character*32 junk

      real dlbol,dteff,age,dlogg
      real grav(nmax1,nmax2)
      real logAge(nmax1,nmax2),logL(nmax1,nmax2),logT(nmax1,nmax2)
      real dummyage(nmax2),dummyl(nmax2),dummyt(nmax2),dummyg(nmax2)
      real dmass(nmax1)
      real dummy(10)

2     format(7x,f7.4,2x,f7.4,3x,f6.4)

      data n/228,900,792,372,552,468,540,372,300,264,384,360,324,204,
     :       300,312,300,312,324,348,384,432,432,468,456,432,420,420,
     :       432,408,432,444/

      data dmass/0.1,0.13,0.16,0.2,0.25,0.3,0.4,0.5,0.6,0.7,0.8,0.9,
     :           1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.5,
     :           2.7,3.0,3.5,4.0,5.0,6.0,7.0/

      age=mage
      nmod=nmax1
      age = log10(age)+6.0

c      print*,'age=',age

c tracks are all in one big file
c idea is to read in a file for each mass and interpolate to get
c log L and log T at the appropriate log(age)

      open (1,file=trim(data_dir())//'/siessz02.dat',status='old')

      do i=1,nmod
c         print*,'n(i)=',n(i)

         read(1,'(A)') junk   
         read(1,'(A)') junk         
         read(1,'(A)') junk

         read(1,*)
     :      (dummy(1),dummy(2),logL(i,j),dummy(3),dummy(4),dummy(5),
     :      logT(i,j),dummy(6),grav(i,j),dmass(i),logAge(i,j), j=1,n(i))

         do j=1,n(i)

c convert to logs

            logAge(i,j)=log10(logAge(i,j))
            logT(i,j)=log10(logT(i,j))
            logL(i,j)=log10(logL(i,j))
         enddo

c          print*,dmass(i),logAge(i,1),logL(i,1),logT(i,1)
c          print*,dmass(i),logAge(i,40),logL(i,40),logT(i,40)
      enddo

c now have the models in the 2d grids
c for each row in the grid, need to find the index closest to the
c required age

      index=0

      do i=1,nmod

         do j=1,n(i)
            dummyage(j)=logAge(i,j)
c         if (i.eq.13) then
c            print*,j,dummyage(j)
c         endif
         enddo




c use numerical recipes subroutine - returns j, where age is between
c logAge(i,j) and logAge(i,j+1)

         call locate(dummyage,n(i),age,j)

c         if (j.gt.1.and.j.lt.n(i)-1) then
          if (j.gt.0.and.j.lt.n(i)) then
c piece of black magic from numerical recipes p.113

            index=index+1

            k = min(max(j-1,1),n(i)-2)

         
c now snip out the bits of the arrays we want into dummy arrays

            do j=1,3
               dummyage(j)=logAge(i,k+j-1)
               dummyl(j)=logL(i,k+j-1)
               dummyt(j)=logT(i,k+j-1)
               dummyg(j)=grav(i,k+j-1)
            enddo

c now interpolate (parabolic so n=3), returns temp and dtemp (error)

            call polint(dummyage,dummyl,3,age,lbol(index),dlbol)
            call polint(dummyage,dummyt,3,age,teff(index),dteff)
            call polint(dummyage,dummyg,3,age,logg(index),dlogg)

            mass(index)=dmass(i)
   
         endif

      enddo

      nmod=index

      close(1)

      end subroutine read_siessz02

c*****************************************************************************



c *****************************************************************************
      subroutine read_siessz01(mage,lbol,teff,mass,logg,nmod)
c *****************************************************************************

      implicit none
c import
      real mage
c export
      integer nmod
      integer nmax1
      parameter (nmax1=29)
      real lbol(nmax1),teff(nmax1),mass(nmax1), logg(nmax1)

c local
      integer nmax2,i,j,k,index
      parameter (nmax2=900)
      integer n(nmax1)         ! no. of points in each model
      
      character*32 junk

      real dlbol,dteff,age, dlogg
      real grav(nmax1, nmax2)
      real logAge(nmax1,nmax2),logL(nmax1,nmax2),logT(nmax1,nmax2)
      real dummyage(nmax2),dummyl(nmax2),dummyt(nmax2),dummyg(nmax2)
      real dmass(nmax1)
      real dummy(10)

2     format(7x,f7.4,2x,f7.4,3x,f6.4)

      data n/204,204,216,264,192,204,228,264,264,264,252,252,276,
     :       312,336,384,372,396,396,408,420,408,420,408,408,360,
     :       384,492,528/

      data dmass/0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,
     :           1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.5,2.7,
     :           3.0,3.5,4.0,5.0,6.0,7.0/

      age=mage
      nmod=nmax1
      age = log10(age)+6.0

c      print*,'age=',age

c tracks are all in one big file
c idea is to read in a file for each mass and interpolate to get
c log L and log T at the appropriate log(age)

      open (1,file=trim(data_dir())//'/siessz01.dat',status='old')



      do i=1,nmod
c         print*,'n(i)=',n(i)

         read(1,'(A)') junk   
         read(1,'(A)') junk         
         read(1,'(A)') junk

         read(1,*)
     :      (dummy(1),dummy(2),logL(i,j),dummy(3),dummy(4),dummy(5),
     :      logT(i,j),dummy(6),grav(i,j),dmass(i),logAge(i,j), j=1,n(i))

         do j=1,n(i)

c convert to logs

            logAge(i,j)=log10(logAge(i,j))
            logT(i,j)=log10(logT(i,j))
            logL(i,j)=log10(logL(i,j))
         enddo

c          print*,dmass(i),logAge(i,1),logL(i,1),logT(i,1)
c          print*,dmass(i),logAge(i,40),logL(i,40),logT(i,40)
      enddo

c now have the models in the 2d grids
c for each row in the grid, need to find the index closest to the
c required age

      index=0

      do i=1,nmod

         do j=1,n(i)
            dummyage(j)=logAge(i,j)
c         if (i.eq.13) then
c            print*,j,dummyage(j)
c         endif
         enddo




c use numerical recipes subroutine - returns j, where age is between
c logAge(i,j) and logAge(i,j+1)

         call locate(dummyage,n(i),age,j)

c         if (j.gt.1.and.j.lt.n(i)-1) then
          if (j.gt.0.and.j.lt.n(i)) then
c piece of black magic from numerical recipes p.113

            index=index+1

            k = min(max(j-1,1),n(i)-2)

         
c now snip out the bits of the arrays we want into dummy arrays

            do j=1,3
               dummyage(j)=logAge(i,k+j-1)
               dummyl(j)=logL(i,k+j-1)
               dummyt(j)=logT(i,k+j-1)
               dummyg(j)=grav(i,k+j-1)
            enddo

c now interpolate (parabolic so n=3), returns temp and dtemp (error)

            call polint(dummyage,dummyl,3,age,lbol(index),dlbol)
            call polint(dummyage,dummyt,3,age,teff(index),dteff)
            call polint(dummyage,dummyg,3,age,logg(index),dlogg)

            mass(index)=dmass(i)
   
         endif

      enddo

      nmod=index

      close(1)

      end subroutine read_siessz01

c*****************************************************************************







c*****************************************************************************
      subroutine read_grenoble(mage,lbol,teff,mass,nmod)
c*****************************************************************************

c grenoble data is in the form of a number of isochrones
c need to ask for the metallicity
c at the moment this is bodged so that you can only pick isochrones
c which are available - i.e. no interpolation done

      implicit none
c import
      real mage
c export
      integer nmod
      integer nmax1
      parameter (nmax1=40)
      real lbol(nmax1),teff(nmax1),mass(nmax1)

c local
      integer i,zflag,ageint
      character*128 grfile
      character*3 agechar
      character*2 zchar
      character*128 filenam
      parameter (grfile = '../files/grenobleiso')
      real rad(nmax1)
      real age
      
      ! This routine cannot set the mass, but just in case an unitialised
      ! value cases a problem...
      mass=0.0
      
2     format(1x,e9.3,2x,e9.3,2x,f6.0)

      age=mage

      ageint = nint(age)
      write(agechar,'(i3.3)')ageint
      age = log10(age)+6.0

c      print*,'age=',age

201   print*,'What metallicity? 1 is Z=0.01, 2 is Z=0.02'
      read(*,*)zflag
      if (zflag.ne.1.and.zflag.ne.2) go to 201

      write(zchar,'(i2.2)')zflag

c      print*,agechar,zchar

      filenam = trim(data_dir())//'/grenobleiso'//
     :          agechar//'z'//zchar//'.dat'

      open (1,file=filenam,form='formatted',access='sequential',
     :  status='old')


      i=1
202      read(1,2,err=203)lbol(i),rad(i),teff(i)
         print*,lbol(i),rad(i),teff(i)
         lbol(i)=log10(lbol(i))
         teff(i)=log10(teff(i))
      i=i+1
      goto 202
203   nmod=i-1

      close(1)

      end subroutine read_grenoble







c*****************************************************************************
      subroutine read_newbaraffe(mage,lbol,teff,mass,nmod)
c*****************************************************************************

c models 0.02-1.4 solar masses and alpha = 1.0

      implicit none
c import
      real mage          ! in Myr
c export
      integer nmod
      integer nmax1
      parameter (nmax1=39)
      real lbol(nmax1),teff(nmax1),mass(nmax1)   ! need to be in log units

c local
      integer nmax2,i,j,k,index
      parameter (nmax2=81)
      integer n(nmax1)         ! no. of points in each model
      
      character*32 junk

      real dlbol,dteff,age
      real logAge(nmax1,nmax2),logL(nmax1,nmax2),logT(nmax1,nmax2)
      real dmass(nmax1)
      real dummyage(nmax2),dummyl(nmax2),dummyt(nmax2)

2     format(1x,f5.3,1x,f8.5,1x,f6.1,8x,f5.2)

      data n/37,45,51,56,59,61,70,80,81,81,81,81,81,81,81,81,
     :       81,81,81,81,81,81,81,81,81,81,81,81,81,81,81,81,
     :       80,78,77,75,73,72,69/


      age=mage
      nmod=nmax1
      age = log10(age)+6.0

c       print*,'age=',age

c tracks are all in one big file
c idea is to read in a file for each mass and interpolate to get
c log L and log T at the appropriate log(age)

      open (1,file=trim(data_dir())//'/baraffe98alpha1.0',
     &form='formatted',access='sequential',status='old')

      read(1,'(A)') junk   
      read(1,'(A)') junk         
      read(1,'(A)') junk


      do i=1,nmod
c       print*,'n(i)=',n(i)

         read(1,2)
     :      (dmass(i),logAge(i,j),logT(i,j),logL(i,j), j=1,n(i))

         do j=1,n(i)

c convert to logs

            logAge(i,j)=log10(logAge(i,j))+9.0
            logT(i,j)=log10(logT(i,j))
         
         enddo

c          print*,dmass(i),logAge(i,1),logL(i,1),logT(i,1)
c          print*,dmass(i),logAge(i,40),logL(i,40),logT(i,40)
      enddo

c now have the models in the 2d grids
c for each row in the grid, need to find the index closest to the
c required age

      index=0

      do i=1,nmod

         do j=1,n(i)
            dummyage(j)=logAge(i,j)
c         if (i.eq.13) then
c            print*,j,dummyage(j)
c         endif
         enddo




c use numerical recipes subroutine - returns j, where age is between
c logAge(i,j) and logAge(i,j+1)

         call locate(dummyage,n(i),age,j)

c         if (j.gt.1.and.j.lt.n(i)-1) then
          if (j.gt.0.and.j.lt.n(i)) then
c piece of black magic from numerical recipes p.113

            index=index+1

            k = min(max(j-1,1),n(i)-2)

         
c now snip out the bits of the arrays we want into dummy arrays

            do j=1,3
               dummyage(j)=logAge(i,k+j-1)
               dummyl(j)=logL(i,k+j-1)
               dummyt(j)=logT(i,k+j-1)
            enddo

c now interpolate (parabolic so n=3), returns temp and dtemp (error)

            call polint(dummyage,dummyl,3,age,lbol(index),dlbol)
            call polint(dummyage,dummyt,3,age,teff(index),dteff)

            mass(index)=dmass(i)
c            print*,index,mass(index),lbol(index),teff(index)
            
          else
            
            print*, 'Could not find appropriate age.'
            print*, 'Required log10(age in years) is ', age
            print*, 'Available range is ', dummyage(1), 
     &            dummyage(n(i))
            stop
            
         endif

      enddo

      nmod=index

      close(1)

      end subroutine read_newbaraffe

c ******************************************************************************
      subroutine read_padgen(modflag, mage, lbol, teff, mass, logg, 
     &nmod)
c ******************************************************************************

      implicit none

      ! The required model.
      integer, intent(in) :: modflag
      ! The required age in Myr.
      real, intent(in) :: mage
      real, intent(out), dimension(:) :: lbol, teff, mass, logg
      ! The actual number of points in the isochrone. 
      integer, intent(out) :: nmod

      type(isochron), dimension(3) :: isoc
      integer :: imod, imass, jpt, iage
      real, dimension(3) :: teff_in, lbol_in, logg_in, age_in
      real :: uncer, min_mass, max_mass

      character(len=2) :: flag

      if (modflag == 8) then
        call read_padova_col(mage, 0, isoc)
      else if (modflag == 9) then
        call read_geneva_col(mage, 0, isoc)
      else
        print*, 'S/R read padgen called in model ', modflag
        stop
      end if

      ! What mass scale to use?
      min_mass=maxval(isoc%mass(1))
      max_mass=min(isoc(1)%mass(isoc(1)%npts), 
     &isoc(2)%mass(isoc(2)%npts))
      max_mass=min(max_mass, isoc(3)%mass(isoc(3)%npts))

      nmod=0
      ! At each mass point in the middle isochrone.
      each_mass: do imass=1, isoc(2)%npts
        if (isoc(2)%mass(imass) < min_mass) cycle each_mass
        if (isoc(2)%mass(imass) > max_mass) cycle each_mass
        nmod=nmod+1
        if (nmod > mmod) then
          print*, 'mmod needs enlarging.'
          stop
        end if
        mass(nmod)=isoc(2)%mass(imass)
        
        ! Set up the middle of the three values we will interpolate accross.
        teff_in(2)=isoc(2)%colour(imass)
        lbol_in(2)=isoc(2)%magnitude(imass)
        logg_in(2)=isoc(2)%logg(imass)
        age_in(2)=isoc(2)%age

        ! Now interpolate the other two points.
        each_age: do iage=1, 3, 2

          age_in(iage)=isoc(iage)%age

          teff_in(iage)=quadint(isoc(iage)%mass(1:isoc(iage)%npts), 
     &    isoc(iage)%colour(1:isoc(iage)%npts), mass(nmod), flag)
          lbol_in(iage)=quadint(isoc(iage)%mass(1:isoc(iage)%npts), 
     &    isoc(iage)%magnitude(1:isoc(iage)%npts), mass(nmod), flag)
          logg_in(iage)=quadint(isoc(iage)%mass(1:isoc(iage)%npts), 
     &    isoc(iage)%logg(1:isoc(iage)%npts), mass(nmod), flag)

        end do each_age

        ! Now interpolate over the three points to get the right age.
        teff(nmod)=quadint(age_in, teff_in, mage, flag)
        lbol(nmod)=quadint(age_in, lbol_in, mage, flag)
        logg(nmod)=quadint(age_in, logg_in, mage, flag)

      end do each_mass

      

      end subroutine read_padgen

      subroutine read_geneva_col(age, colflag, isoc)

      ! Reads the Geneva isochrones file.  Colflag follows the usual rules,
      ! except that if it is zero the magnitude is lbol, and the colour
      ! is teff.

      real, intent(in) :: age
      integer, intent(in) :: colflag
      type(isochron), intent(out), dimension(3) :: isoc

      real :: teff, lbol
      integer :: n_mass
      integer :: nages, iage, jage
      character(len=30), dimension(:), allocatable :: isofiles
      real, dimension(:), allocatable :: log10_ages

      character(len=20), parameter :: isodir='Geneva'

      ! First create a list of ages and file names.
      open(unit=2, file=trim(data_dir())//'/'//trim(isodir)// 
     &'/ages.dat',status='old')
      read(2,*) nages
      allocate(log10_ages(nages), isofiles(nages))
      read(2,'(/)')
      do iage=1, nages
        read(2,*) log10_ages(iage), isofiles(iage)
      end do
      close(2)

      ! Find the nearest age, and one on either side.
      jage=minloc(abs(log10_ages-(log10(age)+6.0)),1)
      ! If the nearest age is the first one, but its above that age
      ! then take the second age as the mid-point of the interpolation.
      if (jage==1 .and. log10_ages(1) < log10(age)+6.0) jage=2
      if (jage == 1) then
        print*, 'Youngest isochrone available is ', 10.0**log10_ages(1)
        print*, 'So cannot do an age of ', age, 'Myr'
        stop
      else if (jage > nages) then
        print*, 'Oldest isochrone available is ', 
     &  10.0**log10_ages(nages)
        print*, 'So cannot do an age of ', age, 'Myr'
        stop
      end if
      
      isoc%age=10.0**(log10_ages(jage-1:jage+1) - 6.0)

      do iage=1, 3
        call read_one_geneva(trim(data_dir())//'/'//trim(isodir)//'/'//
     &  trim(isofiles(jage-2+iage)), colflag, isoc(iage)%npts, 
     &  isoc(iage)%magnitude, isoc(iage)%colour, isoc(iage)%mass,
     &  isoc(iage)%logg)

        
      end do

      end subroutine read_geneva_col

      subroutine read_one_geneva(filnam, colflag, npts, lbol, teff, 
     &mass, logg)

      character(len=*), intent(in) :: filnam
      integer, intent(in) :: colflag

      integer, intent(out) :: npts
      real, dimension(:), intent(out) :: lbol, teff, mass, logg 


      character :: ch_junk
      real :: junk
      real :: v, u_b, b_v, v_r, v_i
      integer :: kage, i


      open(unit=2, status='old', file=filnam)
      read(2,*)
      read(2,*)
      read(2,*) ch_junk, npts
      read(2,*)
      read(2,*)
      read(2,*)
      if (npts > size(lbol)) then
         print*, 'mmod needs expanding to ', npts
         stop
      end if
      do i=1, npts
         read(2,*) junk, mass(i), junk, teff(i), 
     &        logg(i), lbol(i), v, u_b, b_v, v_r, v_i
         if (colflag == 0) then
           ! Do nothing, everything is fine.
         else if (colflag == 21) then
            lbol(i)=v
            teff(i)=v_i
         else if (colflag == 24) then
            lbol(i)=v
            teff(i)=b_v
         else if (colflag == 25) then
            lbol(i)=u_b
            teff(i)=b_v
         else
            print*,'For the Geneva isochrones only the colours V/V-I'
            print*,'V/B-V and U-B/B-V are read from the files supplied.'
            stop
         end if
      end do
      ! Sometimes (just to be a pest) there are no values for the last
      !few datapoints, which are set to -99.
      kage=npts
      do i=kage, 1, -1
         if (teff(i) < -98.0) npts=kage-1
         if (lbol(i) < -98.0) npts=kage-1
      end do
      close(2)

      end subroutine read_one_geneva 


      subroutine read_padova_col(age, colflag, isoc)

      real, intent(in) :: age
      integer, intent(in) :: colflag
      type(isochron), dimension(3), intent(out) :: isoc

      character(len=150) :: test_str
      integer :: n_ages, iostat, i_ages, i, jage, j
      real :: old_age, new_age, junk
      integer, allocatable, dimension(:) :: n_masses
      real, allocatable, dimension(:) :: log10_ages
      integer, dimension(3) :: icnt
      real :: teff, logg, mag_bol, lbol
      real :: mag_u, mag_b, mag_v, mag_r, mag_i, mag_j, mag_h, mag_k

      !character(len=50), parameter :: isofile='/Padova_MC/isoc_z004.dat'
      !character(len=50), parameter :: isofile='/Padova_MC/isoc_z001.dat'
      !character(len=50), parameter :: isofile='/Padova_MC/isoc_z030.dat'
      !character(len=50), parameter :: isofile='/Padova/isoc_z0.dat'
      !character(len=50), parameter :: isofile='/Padova/isoc_z0001.dat'
      !character(len=50), parameter :: isofile='/Padova/isoc_z0004.dat'
      !character(len=50), parameter :: isofile='/Padova/isoc_z001.dat'
      !character(len=50), parameter :: isofile='/Padova/isoc_z004.dat'
      !character(len=50), parameter :: isofile='/Padova/isoc_z008.dat'
      character(len=50), parameter :: isofile='/Padova/isoc_z019.dat'
      !character(len=50), parameter :: isofile='/Padova/isoc_z030.dat'

      real :: old_teff, new_teff, old_bol, new_bol
      logical :: first=.true.

      if (first) then
        print*, 'Using Padova datafile ', trim(isofile)
        first=.false.
      end if

      open (unit=2, file=trim(data_dir())//trim(isofile), 
     &status='old')
      ! First get a feel for the ages available.
      n_ages=0
      each_line1: do
        read(2,'(a50)',iostat=iostat) test_str
        if (iostat < 0) then
          exit each_line1
        else if (iostat > 0) then
          print*, 'Error ', iostat, ' reading file in each_line1.'
          stop
        end if
        if (test_str(1:1) == '#') cycle each_line1
        read(test_str,*) new_age
        if (n_ages == 0) then
          n_ages=1
          old_age=new_age
        else if (abs(new_age-old_age) > 2.0*tiny(old_age)) then
          n_ages=n_ages+1
          old_age=new_age
        end if
      end do each_line1
      !print*, 'Number of ages available is ', n_ages

      ! Now count the number of mass points at each age.
      rewind(2)
      allocate(n_masses(n_ages), log10_ages(n_ages))
      n_masses=0
      i_ages=0
      each_line2: do
        read(2,'(a50)',iostat=iostat) test_str
        if (iostat < 0) then
          exit each_line2
        else if (iostat > 0) then
          print*, 'Error ', iostat, ' reading file in each_line2.'
          stop
        end if
        if (test_str(1:1) == '#') cycle each_line2
        read(test_str,*) new_age, junk, junk, new_bol, new_teff
        if (i_ages == 0) then
          i_ages=1
          old_age=new_age
          n_masses(i_ages)=n_masses(i_ages)+1
          log10_ages(i_ages)=new_age
          old_teff=new_teff
          old_bol=new_bol
        else if (abs(new_age-old_age) > 2.0*tiny(old_age)) then
          i_ages=i_ages+1
          old_age=new_age
          old_teff=new_teff
          old_bol=new_bol
          log10_ages(i_ages)=new_age
        else
          if (abs(old_teff-new_teff)>0.01 .and. 
     &      abs(old_bol-new_bol)>0.01) then
            ! Some Padova masses are very close, which confuses the
            ! quadratic interpolation.
            n_masses(i_ages)=n_masses(i_ages)+1
          end if
          old_teff=new_teff
          old_bol=new_bol
        end if
      end do each_line2

      ! Find the nearest age.
      jage=minval(minloc(abs(log10_ages-(log10(age)+6.0))))
!      print*, 'Nearest age is ', 10.0**(log10_ages(jage)-6.0), 
!     &'Myr number ', jage
      ! If the nearest age is the first one, but its above that age
      ! then take the second age as the mid-point of the interpolation.
      if (jage==1 .and. log10_ages(1) < log10(age)+6.0) jage=2
      if (jage == 1) then
        print*, 'Youngest isochrone available is ', 10.0**log10_ages(1)
        print*, 'So cannot do an age of ', age, 'Myr'
        stop
      else if (jage > n_ages) then
        print*, 'Oldest isochrone available is ', 
     &  10.0**log10_ages(n_ages)
        print*, 'So cannot do an age of ', age, 'Myr'
        stop
      end if
      isoc%age=10.0**(log10_ages(jage-1:jage+1) - 6.0)
      isoc%npts=n_masses(jage-1:jage+1)

      ! Now get the three ages we want.
      rewind(2)
      icnt=0
      old_teff=0.0
      old_bol=0.0
      each_line3: do
        read(2,'(a150)',iostat=iostat) test_str
        if (iostat < 0) then
          exit each_line3
        else if (iostat > 0) then
          print*, 'Error ', iostat, ' reading file in each_line3.'
          stop
        end if
        if (test_str(1:1) == '#') cycle each_line3
        read(test_str,*) new_age, junk, junk, new_bol, new_teff
        j=minval(minloc(abs(new_age-(log10(isoc%age)+6.0))))
        if (abs(new_age-(log10(isoc(j)%age)+6.0)) 
     &    < 2.0*tiny(new_age)) then
          if (abs(old_teff-new_teff)>0.01 .and. 
     &      abs(old_bol-new_bol)>0.01) then
            icnt(j)=icnt(j)+1
            read(test_str,*) junk, isoc(j)%mass(icnt(j)), junk, lbol, 
     &      teff, isoc(j)%logg(icnt(j)), mag_bol, mag_u, mag_b, mag_v, 
     &      mag_r, mag_i, mag_j, mag_h, mag_k
            if (colflag == 0) then
              isoc(j)%magnitude(icnt(j))=lbol
              isoc(j)%colour(icnt(j))=teff
            else if (colflag == 21) then
              isoc(j)%magnitude(icnt(j))=mag_v
              isoc(j)%colour(icnt(j))=mag_v-mag_i
            else if (colflag == 22) then
              isoc(j)%magnitude(icnt(j))=mag_i
              isoc(j)%colour(icnt(j))=mag_r-mag_i
            else if (colflag == 23) then
              isoc(j)%magnitude(icnt(j))=mag_i
              isoc(j)%colour(icnt(j))=mag_i-mag_k
            else if (colflag == 24) then
              isoc(j)%magnitude(icnt(j))=mag_v
              isoc(j)%colour(icnt(j))=mag_b-mag_v
            else if (colflag == 25) then
              isoc(j)%magnitude(icnt(j))=mag_u-mag_b
              isoc(j)%colour(icnt(j))=mag_b-mag_v
            !else if (colflag == 5) then
            !  isoc(j)%magnitude(icnt(j))=mag_k
            !  isoc(j)%colour(icnt(j))=mag_j-mag_k
            else
              print*, 'For the Padova isochrones only the colours V/V-I'
              print*, 
     &        'I/R-I, I/I-K and V/B-V are read from the files supplied.'
              stop
            end if
            !print*, new_age,icnt(j), isoc(j)%mass(icnt(j)), isoc(j)%magnitude(icnt(j)), isoc(j)%colour(icnt(j))
          end if
          old_teff=new_teff
          old_bol=new_bol
        end if
      end do each_line3

      deallocate(n_masses, log10_ages)
      close(2)

      !do j=1, 3
      !  write(30+j,*) isoc(j)%age
      !  write(30+j,*) trim(isofile)
      !  write(30+j,*) ' '
      !  do jage=1, isoc(j)%npts
      !    write(30+j,*) isoc(j)%colour(jage), isoc(j)%magnitude(jage)
      !  end do
      !end do
      !stop

      end subroutine read_padova_col

c ******************************************************************************
      subroutine read_newbaraffe_19(mage,lbol,teff,mass,nmod)
c ******************************************************************************

      implicit none
c import
      real mage          ! in Myr
c export
      integer nmod
      integer nmax1
      parameter (nmax1=38)
      real lbol(nmax1),teff(nmax1),mass(nmax1)   ! need to be in log units

c local
      integer nmax2,i,j,k,index
      parameter (nmax2=81)
      integer n(nmax1)         ! no. of points in each model
      
      character*32 junk

      real dlbol,dteff,age
      real logAge(nmax1,nmax2),logL(nmax1,nmax2),logT(nmax1,nmax2)
      real dmass(nmax1)
      real dummyage(nmax2),dummyl(nmax2),dummyt(nmax2)

2     format(1x,f5.3,1x,f8.5,1x,f6.1,8x,f5.2)

      data n/37,45,51,56,59,61,70,80,81,81,81,81,81,81,81,81,
     :       81,81,81,81,81,81,81,81,81,81,81,81,81,81,81,79,
     :       78,76,75,74,71,69/

      age=mage
      nmod=nmax1
      age = log10(age)+6.0

c      print*,'age=',age

c tracks are all in one big file
c idea is to read in a file for each mass and interpolate to get
c log L and log T at the appropriate log(age)

      open (1,file=trim(data_dir())//'/baraffe98alpha1.9',
     &form='formatted',access='sequential',status='old')

      read(1,'(A)') junk   
      read(1,'(A)') junk         
      read(1,'(A)') junk


      do i=1,nmod
c         print*,'n(i)=',n(i)

         read(1,2)
     :      (dmass(i),logAge(i,j),logT(i,j),logL(i,j), j=1,n(i))

         do j=1,n(i)

c convert to logs

            logAge(i,j)=log10(logAge(i,j))+9.0
            logT(i,j)=log10(logT(i,j))
         
         enddo

c          print*,dmass(i),logAge(i,1),logL(i,1),logT(i,1)
c          print*,dmass(i),logAge(i,40),logL(i,40),logT(i,40)
      enddo

c now have the models in the 2d grids
c for each row in the grid, need to find the index closest to the
c required age

      index=0

      do i=1,nmod

         do j=1,n(i)
            dummyage(j)=logAge(i,j)
c         if (i.eq.13) then
c            print*,j,dummyage(j)
c         endif
         enddo




c use numerical recipes subroutine - returns j, where age is between
c logAge(i,j) and logAge(i,j+1)

         call locate(dummyage,n(i),age,j)

c         if (j.gt.1.and.j.lt.n(i)-1) then
          if (j.gt.0.and.j.lt.n(i)) then
c piece of black magic from numerical recipes p.113

            index=index+1

            k = min(max(j-1,1),n(i)-2)

         
c now snip out the bits of the arrays we want into dummy arrays

            do j=1,3
               dummyage(j)=logAge(i,k+j-1)
               dummyl(j)=logL(i,k+j-1)
               dummyt(j)=logT(i,k+j-1)
            enddo

c now interpolate (parabolic so n=3), returns temp and dtemp (error)

            call polint(dummyage,dummyl,3,age,lbol(index),dlbol)
            call polint(dummyage,dummyt,3,age,teff(index),dteff)

            mass(index)=dmass(i)
c            print*,index,mass(index),lbol(index),teff(index)
   
         endif

      enddo

      nmod=index

      close(1)

      end subroutine read_newbaraffe_19


c ******************************
      real function bcbv(bv, flag)
c ******************************

c Original version by Rob to use Flower(1996)

      implicit none

      real, intent(in) :: bv
      integer, intent(out) :: flag
      
      integer j, k
      real dbcbv

      integer, parameter :: nmax_flower=44
      real, dimension(nmax_flower) :: flower_bc, flower_x

      integer, parameter :: nmax_bessell=75
      real, dimension(nmax_bessell) :: bessell_bcor, bessell_x

      integer :: nmax
      real, dimension(nmax_bessell) :: bc, x

      real, dimension(3) :: dumx, dumbc

      data flower_bc/
     :        -4.72, -3.234,-2.177,-1.525,-1.013,-0.614,-0.334,-0.155,
     :        -0.050, 0.004, 0.024, 0.032, 0.035,0.034,0.03,0.020,0.006,
     :        -0.012,-0.035,-0.061,-0.091,-0.124,-0.161,-0.200,-0.242,
     :        -0.285,-0.332,-0.382,-0.437,-0.494,-0.552,-0.614,-0.679,
     :        -0.752,-0.831,-0.920,-1.024,-1.150,-1.312,-1.539,-1.885,
     :        -2.460,-3.502,-5.535/
      data flower_x/
     :      -0.35,-0.3,-0.25,-0.2,-0.15,-0.1,-0.05,0.0,0.05,0.1,0.15,
     :       0.20,0.25,0.30,0.35,0.40,0.45,0.5,0.55,0.6,0.65,0.7,0.75,
     :       0.8,0.85,0.9,0.95,1.0,1.05,1.1,1.15,1.2,1.25,1.3,1.35,
     :       1.4,1.45,1.5,1.55,1.6,1.65,1.7,1.75,1.80/

      data bessell_bcor/
     :-4.27, -4.21, -4.14, -4.07, -4.00, -3.93, -3.86, -3.79, -3.71, 
     :-3.63, -3.56, -3.48, -3.40, -3.33, -3.25, -3.18, -3.12, -3.05,
     :-2.99, -2.92, -2.85, -2.77, -2.68, -2.60, -2.50, -2.41, -2.30,
     :-2.19, -2.08, -1.95, -1.82, -1.68, -1.53, -1.37, -1.20, -1.03, 
     :-0.85, -0.80, -0.75, -0.70, -0.66, -0.61, -0.56, -0.51, -0.46,
     :-0.41, -0.36, -0.31, -0.25, -0.21, -0.16, -0.11, -0.07, -0.03,
     : 0.00,  0.00,  0.01,  0.02,  0.02,  0.01,  0.01,  0.00,  0.00,
     :-0.02, -0.04, -0.07, -0.11, -0.18, -0.27, -0.40, -0.56, -0.77,
     :-1.00, -1.25, -1.62/
      data bessell_x/
     :-0.319, -0.318, -0.317, -0.316, -0.314, -0.313, -0.311, -0.310,
     :-0.308, -0.305, -0.303, -0.301, -0.299, -0.296, -0.294, -0.291,
     :-0.288, -0.283, -0.278, -0.272, -0.266, -0.259, -0.252, -0.246,
     :-0.240, -0.233, -0.226, -0.218, -0.209, -0.200, -0.189, -0.178,
     :-0.166, -0.153, -0.139, -0.124, -0.106, -0.101, -0.096, -0.090,
     :-0.084, -0.078, -0.071, -0.063, -0.053, -0.042, -0.030, -0.016,
     : 0.000,  0.017,  0.037,  0.058,  0.083,  0.102,  0.138,  0.176,
     : 0.215,  0.252,  0.290,  0.329,  0.372,  0.421,  0.477,  0.540,
     : 0.609,  0.684,  0.763,  0.844,  0.925,  1.012,  1.118,  1.240,
     : 1.348,  1.415,  1.447/  


      flag=0

      nmax=nmax_flower
      bc(1:nmax)=flower_bc(1:nmax)
      x(1:nmax)=flower_x(1:nmax)
      
      if (bv<x(1) .or. bv>x(nmax)) flag=1

      

c use numerical recipes subroutine - returns j, where b-v is between
c x(j) and x(j+1)

      call locate(x,nmax,bv,j)

c piece of black magic from numerical recipes p.113

      k = min(max(j-1,1),nmax-2)

         
c now snip out the bits of the arrays we want into dummy arrays

      do j=1,3
        dumx(j)=x(k+j-1)
        dumbc(j)=bc(k+j-1)
      enddo

c now interpolate (parabolic so n=3), returns bcbv and dbcbv (error)

      call polint(dumx,dumbc,3,bv,bcbv,dbcbv)

      end function bcbv


c ****************************
      real function bcvi(vi, logg, teff, flag)
c ****************************

c     Uses a polynomial fit to data (described in the Bols directory) for
c     the red stars, Bessell for the blue stars, and whichever is the
c     larger in the crossover region.
      
      real, intent(in) :: vi, logg, teff
      integer, intent(out) :: flag

      flag=0
      if (vi > 4.2) then
        flag=1
      else if (vi > 1.0) then  
        bcvi = -0.1216 + 2.0591*vi - 2.0137*vi**2.0
     :      +0.923752*vi**3.0 -0.205533*vi**4.0 +1.697281e-2*vi**5.0
     :      - vi
        if (vi < 1.5) bcvi=max(bcvi,bessell_bc_v(logg, teff, flag))
      else
        bcvi=bessell_bc_v(logg, teff, flag)
      end if

      end function bcvi


c ************************
      real function bcri(ri, flag)
c ************************

c     New function added by timn (see the README in his bols
c     directory April 2006.      
      
      real ri
      integer :: flag

      flag=0
      if (ri<-0.3 .or. ri>2.8) flag=1
                                                               
c     bcri=0.1738985+0.8821769*ri-0.4753564*ri**2
      
      bcri = -0.189001963  +5.58346367*ri -12.7131376*(ri**2.0)  
     :      +12.4155569*(ri**3.0) -5.44285679*(ri**4.0) 
     :      +0.862898767*(ri**5.0)
      
      
      end function bcri



c ************************
      real function bcik(ik, flag)
c ************************
      
c valid for 1.7 < I-K < 5.2     
      
      real ik
      integer :: flag

      flag=0
      if (ik<1.7 .or. ik>5.2) flag=1                                                        

      bcik=1.4316+0.6933*ik-6.189e-02*(ik**2.0)

      end function bcik


      end module colteff_subs
