      module colteff_subs

      use bessell
      use ken_hart
      use phoenix

      implicit none

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

      ! These two structures are slowly evolving to be the same.
      type isochron
        real :: age
        integer :: npts
        real, dimension(mmod) :: col, mag, mass, teff
        real, dimension(mmod) :: lbol, logg
      end type isochron

      type chrone 
        ! Now the specification of the iscohrone.  Its age, colours and model.
        real :: age
        integer :: colflag, modflag, npts
        real, dimension(:), allocatable :: col, mag, mass, teff
        character(len=50), dimension(:), allocatable :: 
     &  colflg, magflg, massflg, teffflg
      end type chrone

      !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_z019_new.dat'
      !character(len=50), parameter :: isofile='/Padova/isoc_z030.dat
      !character(len=50),parameter:: isofile='/Padova_HST/isoc_z008.dat'

      ! 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=70
      character(len=12), dimension(n_names), parameter, public :: 
     &col_name=(/
     &'V-I        ','R-I        ','I-K        ','B-V        ',
     &'B-V        ','           ','           ','           ',
     &'           ','           ',

     &'V-I        ','           ','           ','B-V        ',
     &'B-V        ','           ','           ','           ',
     &'           ','           ',

     &'V-I        ','R-I        ','I-K        ','B-V        ',
     &'B-V        ','F555W-F814W','           ','           ',
     &'           ','           ',

     &'           ','           ','           ','B-V        ',
     &'           ','           ','           ','           ',
     &'           ','           ',

     &'V-I        ','           ','           ','           ',
     &'           ','           ','J-H        ','           ',
     &'           ','           ',

     &'g-i        ','i-z        ','           ','           ',
     &'           ','           ','           ','           ',
     &'           ','           ',

     &'V-I        ','           ','           ','B-V        ',
     &'           ','V-I        ','           ','           ',
     &'           ','           '/)

      character(len=12), dimension(n_names), public, parameter :: 
     &mag_name=(/
     &'V    ','I    ','I    ','V    ','U-B  ',
     &'     ','     ','     ','     ','     ',

     &'V    ','     ','     ','V    ','U-B  ',
     &'     ','     ','     ','     ','     ',

     &'V    ','I    ','I    ','V    ','U-B  ',
     &'F555W','     ','     ','     ','     ',

     &'     ','     ','     ','V    ','     ',
     &'     ','     ','     ','     ','     ',

     &'V    ','     ','     ','     ','     ',
     &'     ','J    ','     ','     ','     ',

     &'g    ','i    ','     ','     ','     ',
     &'     ','     ','     ','     ','     ',

     &'V    ','     ','     ','V    ','     ',
     &'I    ','     ','     ','     ','     '/)

      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)               ', 
     &'Bessell (solar abundance)               ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
c 
     &'As supplied with isochrone.             ', 
     &'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).               ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'Kenyon & Hartmann (1995).               ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ',
c 
     &'Phoenix ugriz                           ', 
     &'Phoenix ugriz                           ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ',
c
     &'Phoenix UBVRI                           ', 
     &'                                        ', 
     &'                                        ', 
     &'Phoenix UBVRI                           ', 
     &'                                        ', 
     &'Phoenix UBVRI                           ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        ', 
     &'                                        '/)
      ! 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        ', '             ',
     &'             ', 'siessz02     ', 'newbaraffe   ', 
     &'newbaraffe_19', 'siessz01     ', 'padova       ',
     &'geneva       ', 'user         ', 'palla_stahler', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             ', '             ', '             ', 
     &'             '/)
      
      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, 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_pms(trim(data_dir())//'/DAM97/dam97.iso', 
     &  age, lbol, teff, mass, logg, nmod)
      else if (modflag.eq.4) then
        call read_pms(trim(data_dir())//'/Siess/siessz02.iso', 
     &  age, lbol, teff, mass, logg, nmod)
      else if (modflag.eq.5) then
        call read_pms(trim(data_dir())//'/baraffe98alpha1.0', 
     &  age, lbol, teff, mass, logg, nmod)
      else if (modflag.eq.6) then
        call read_pms(trim(data_dir())//'/baraffe98alpha1.9', 
     &  age, lbol, teff, mass, logg, nmod)
      else if (modflag.eq.7) then
        call read_pms(trim(data_dir())//'/Siess/siessz01.iso', 
     &  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
        character(len=30) :: wrkflg

        type(isochron), dimension(1) :: isoc
      

                                ! Read the appropriate models
        
        if (modflag.eq.1) then
          call read_pms(trim(data_dir())//'/DAM97/dam97.iso', 
     &    cage, lbol, teff, mass, logg, nmod)
        else if (modflag.eq.4) then
          call read_pms(trim(data_dir())//'/Siess/siessz02.iso', 
     &    cage, lbol, teff, mass, logg, nmod)
        else if (modflag.eq.5) then
          call read_pms(trim(data_dir())//'/baraffe98alpha1.0', 
     &    cage, lbol, teff, mass, logg, nmod)
        else if (modflag.eq.6) then
          call read_pms(trim(data_dir())//'/baraffe98alpha1.9', 
     &    cage, lbol, teff, mass, logg, nmod)
        else if (modflag.eq.7) then
          call read_pms(trim(data_dir())//'/Siess/siessz01.iso', 
     &    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_00.000.dat'
          write(ifname(6:10), '(f5.2)') cage
          do i=6, 10
            if (ifname(i:i) == ' ') ifname(i:i)='0'
          end do
          call read_padova_col(ifname, cage, colflag, isoc)
          nmod=isoc(1)%npts
          lbol=isoc(1)%lbol
          teff=isoc(1)%teff
          mass=isoc(1)%mass
          logg=isoc(1)%logg
        else if (modflag.eq.11) then
          call read_pms(trim(data_dir())//'/palla_stahler.dat', 
     &    cage, lbol, teff, mass, logg, nmod)
        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), 
     &      colflg(imod))
          else if (colflag==14 .or. colflag==15) then
            isocol(imod)=bessell_col_bv(logg(imod), teff(imod), 
     &      colflg(imod))
          else if (colflag == 34) then
            isocol(imod)=bessell_col_bv(logg(imod), teff(imod), wrkflg)
            call tycho_col(isocol(imod), colflg(imod))
            if (trim(wrkflg) /= 'OK') colflg(imod) = wrkflg
         else if (colflag == 41) then
            isocol(imod)=ken_hart_col_vi(teff(imod), colflg(imod))
         else if (colflag == 47) then
            isocol(imod)=ken_hart_col_jh(teff(imod), colflg(imod))
          else if (colflag == 51) then
            isocol(imod)=phoenix_sloan_gi(logg(imod), teff(imod), 
     &      colflg(imod))
          else if (colflag == 52) then
            isocol(imod)=phoenix_sloan_iz(logg(imod), teff(imod), 
     &      colflg(imod))
          else if (colflag == 61 .or. colflag == 66) then
            isocol(imod)=phoenix_bess_vi(logg(imod), teff(imod),
     &      colflg(imod))
          else if (colflag == 64) then
            isocol(imod)=phoenix_bess_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), 
     &      magflg(imod))
          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), magflg(imod))
            else if (colflag==34) then
              bc=bessell_bc_v(logg(imod), teff(imod), magflg(imod))
              bc=bc+tycho_bc(isocol(imod), magflg(imod))
              if (trim(wrkflg) /= 'OK') magflg(imod) = wrkflg
            else if (colflag==41) then
              bc=ken_hart_bc_v(teff(imod), bcflag)
            else if (colflag==47) then
              bc=ken_hart_bc_j(teff(imod), bcflag)
            else if (colflag==51) then
              bc=phoenix_sloan_g(logg(imod), teff(imod), magflg(imod))
            else if (colflag==52) then
              bc=phoenix_sloan_i(logg(imod), teff(imod), magflg(imod))
            else if (colflag == 61 .or. colflag == 64) then
              bc=phoenix_bess_v(logg(imod), teff(imod), magflg(imod))
            else if (colflag == 66) then
              bc=phoenix_bess_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<11) 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

!          print*, imod, isomag(imod), trim(colflg(imod)), 
!     &    isocol(imod), trim(magflg(imod)), 10.0**teff(imod),
!     &    lbol(imod), mass(imod), bc
          
        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 == 47) then
           ! Vector from Nick Wright, who got it from Schlegel.
	   ! A(J)=2.76*E(J-H)
	   mag = mag + 2.76*red
        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  


      subroutine read_pms(ifname, age, lbol, teff, mass, logg, nmod)

      ! Modelled on Rob's Baraffe routine, this is a more generalisable
      ! routine for reading and interpolating PMS isochrones.

      ! At this stage not very general, designed to read the
      ! Palla and Stahler files Cameron created, but should work for
      ! Baraffe.

      character(len=*) :: ifname

      ! Required age in Myr.
      real, intent(in) :: age

      ! Need to be in log units
      real, dimension(:), intent(out) :: lbol, teff, mass, logg
      integer, intent(out) :: nmod

      ! Locals
      integer :: n_masses, i_masses, max_n_ages, i_ages
      integer, allocatable, dimension(:) :: n_ages
      
      real, allocatable, dimension(:,:) :: logAge, logL, logT, log_G
      real, allocatable, dimension(:) :: dmass

      real :: old_mass, new_mass
      integer :: iostat
      character(len=30) :: flag


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

      open (1,file=ifname,
     &form='formatted',access='sequential',status='old')

      
      ! Start by counting the number of masses.
      read(1,'(/,/)') 
      n_masses=0
      old_mass=-1
      do
        read(1,*,iostat=iostat) new_mass
        if (iostat < 0) exit
        if (new_mass > nearest(old_mass,1.0)) then
          ! If the new mass is bigger than the old one.
          n_masses=n_masses+1
          old_mass=new_mass
        end if
      end do

      ! Now count the number of ages at each mass.
      allocate(n_ages(n_masses))
      rewind(1)
      read(1,'(/,/)') 
      old_mass=-1.0
      i_masses=0
      each_mass: do
        read(1,*,iostat=iostat) new_mass
        if (iostat < 0) exit each_mass
        if (new_mass > nearest(old_mass,1.0)) then
          ! If the new mass is bigger than the old one.
          old_mass=new_mass
          i_masses=i_masses+1
          n_ages(i_masses)=1
        else
          n_ages(i_masses)=n_ages(i_masses)+1
        end if
      end do each_mass
      max_n_ages=maxval(n_ages,1)

      ! Now read the data in.
      rewind(1)
      read(1,'(/,/)') 
      allocate(dmass(n_masses), logAge(n_masses, max_n_ages))
      allocate(logT(n_masses, max_n_ages), logL(n_masses, max_n_ages))
      allocate(log_G(n_masses, max_n_ages))
      do i_masses=1, n_masses

        do i_ages=1, n_ages(i_masses)
          read(1,*) dmass(i_masses),logAge(i_masses,i_ages),
     :    logT(i_masses,i_ages), log_G(i_masses,i_ages), 
     :    logL(i_masses,i_ages)
        end do

        ! Convert to logs
        do i_ages=1,n_ages(i_masses)
          logAge(i_masses,i_ages)=log10(logAge(i_masses,i_ages))+9.0
          logT(i_masses,i_ages)=log10(logT(i_masses,i_ages))
        end do

      end do

      close(1)

      ! Now have the models in the 2d grids
      ! for each row in the grid, need to find the index closest to the
      ! required age.

      nmod=0

      do i_masses=1, n_masses

        lbol(nmod+1)=quadintflg(logAge(i_masses,1:n_ages(i_masses)), 
     &  logL(i_masses,1:n_ages(i_masses)), age, flag)
        if (trim(flag) /= 'OK') cycle

        teff(nmod+1)=quadintflg(logAge(i_masses,1:n_ages(i_masses)), 
     &  logT(i_masses,1:n_ages(i_masses)), age, flag)
        if (trim(flag) /= 'OK') cycle
        
        logg(nmod+1)=quadintflg(logAge(i_masses,1:n_ages(i_masses)), 
     &  log_G(i_masses,1:n_ages(i_masses)), age, flag)
        if (trim(flag) /= 'OK') cycle
        
        nmod=nmod+1
        mass(nmod)=dmass(i_masses)

        !print*, nmod, mass(nmod), lbol(nmod), teff(nmod)

      end do
            
      if (nmod == 0) then
            
        print*, 'Could not find appropriate age.'
        print*, 'Required log10(age in years) is ', age
        stop    
      endif

      end subroutine read_pms


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=30) :: flag

      if (modflag == 8) then
        call read_padova_col(trim(data_dir())//isofile, 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)%teff(imass)
        lbol_in(2)=isoc(2)%lbol(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)=linint(isoc(iage)%mass(1:isoc(iage)%npts), 
     &    isoc(iage)%teff(1:isoc(iage)%npts), mass(nmod), flag)
          lbol_in(iage)=linint(isoc(iage)%mass(1:isoc(iage)%npts), 
     &    isoc(iage)%lbol(1:isoc(iage)%npts), mass(nmod), flag)
          logg_in(iage)=linint(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)=linint(age_in, teff_in, mage, flag)
        lbol(nmod)=linint(age_in, lbol_in, mage, flag)
        logg(nmod)=linint(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, returning the three isochrones
      ! closest to the requested age which can then be interpolated.

      ! 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)
      jage=minloc(abs(log10_ages-age),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) < age) jage=2
      if (jage == 1) then
        print*, 'Youngest isochrone available is ', log10_ages(1)
        print*, 'So cannot do an age of ', age
        stop
      else if (jage > nages) then
        print*, 'Oldest isochrone available is ', log10_ages(nages)
        print*, 'So cannot do an age of ', age
        stop
      end if
      
      isoc%age=log10_ages(jage-1:jage+1)

      do iage=1, 3
        call read_one_geneva(trim(data_dir())//'/'//trim(isodir)//'/'//
     &  trim(isofiles(jage-2+iage)), colflag, isoc(iage))

      end do

      end subroutine read_geneva_col

      subroutine read_one_geneva(filnam, colflag, isoc)

      ! Reads a single Geneva isochrone.

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

      type(isochron), intent(inout) :: isoc

      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, isoc%npts
      read(2,*)
      read(2,*)
      read(2,*)
      if (isoc%npts > mmod) then
         print*, 'mmod needs expanding to ', isoc%npts
         stop
      end if
      do i=1, isoc%npts
         read(2,*) junk, isoc%mass(i), junk, isoc%teff(i), 
     &        isoc%logg(i), isoc%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
            isoc%mag(i)=v
            isoc%col(i)=v_i
         else if (colflag == 24) then
            isoc%mag(i)=v
            isoc%col(i)=b_v
         else if (colflag == 25) then
            isoc%mag(i)=u_b
            isoc%col(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=isoc%npts
      do i=kage, 1, -1
         if (isoc%teff(i) < -98.0) isoc%npts=kage-1
         if (isoc%lbol(i) < -98.0) isoc%npts=kage-1
      end do
      close(2)

      end subroutine read_one_geneva 


      subroutine read_padova_col(ifname, age, colflag, isoc)

      ! A very general isochrone reader.  If isoc is of dimension one
      ! the reader expects to find an isochrone at exactly the right 
      ! age.  If it is if dimension three it will return the three
      ! nearest in age.

      character(len=*) :: ifname
      real, intent(in) :: age
      ! The colours required.  Set to zero for lbol and Teff.
      integer, intent(in) :: colflag
      type(isochron), dimension(:), intent(out) :: isoc

      character(len=150) :: test_str
      integer :: n_ages, iostat, i_ages, i, jage, j
      real :: old_mass, new_mass, old_age, new_age, junk
      integer, allocatable, dimension(:) :: n_masses
      real, allocatable, dimension(:) :: log10_ages
      integer, allocatable, dimension(:) :: icnt
      real :: teff, logg, mag_bol, lbol
      real :: mag_u, mag_b, mag_v, mag_r, mag_i, mag_j, mag_h, mag_k
      integer :: icol_plus, icol_minus, imag_plus, imag_minus
      integer :: imass, igrav, iteff, ilbol
      character(len=10) :: col_plus, col_minus, mag_plus, mag_minus
      character(len=20), dimension(20) :: col_head
      real, dimension(20) :: data
      integer :: ihed, nhed, ipos
      character(len=20) work_str

      logical :: first=.true.

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

      open (unit=2, file=ifname, status='old')
      ! First get a feel for the ages available.
      n_ages=0
      each_line1: do
        read(2,'(a50)',iostat=iostat) test_str
        test_str=adjustl(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
        test_str=adjustl(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, new_mass
        if (i_ages == 0) then
          i_ages=1
          old_age =new_age
          old_mass=new_mass
          n_masses(i_ages)=n_masses(i_ages)+1
          log10_ages(i_ages)=new_age
        else if (abs(new_age-old_age) > 2.0*tiny(old_age)) then
          i_ages=i_ages+1
          old_age =new_age
          old_mass=new_mass
          log10_ages(i_ages)=new_age
        else
          if (abs(new_mass-old_mass) > 2.0*tiny(old_mass)) then
            old_age =new_age
            old_mass=new_mass
            n_masses(i_ages)=n_masses(i_ages)+1
          end if
        end if
      end do each_line2

      ! Find the nearest age.
      jage=minloc(abs(log10_ages-age),1)
!      print*, 'Nearest age is ', log10_ages(jage), 
!     &'Myr number ', jage
      if (size(isoc) == 1) then
        ! Do we have an exact match?
        if (abs(log10_ages(jage)-age) < 0.0001) then
          isoc(1)%age=age
          if (size(isoc) > 1) isoc(2)%age=0.0
          if (size(isoc) > 2) isoc(3)%age=0.0
          isoc(1)%npts=n_masses(jage)
        else 
          print*, 'read_padova_col called for one isochrone when there '
          print*, 'is not an exact match in the file.'
          stop
        end if
      else
        ! 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) < age) jage=2
        if (jage  == 1) then
          print*, 'Youngest isochrone available is ', log10_ages(1)
          print*, 'So cannot do an age of ', age
          stop
        else if (jage > n_ages) then
          print*, 'Oldest isochrone available is ', log10_ages(n_ages)
          print*, 'So cannot do an age of ', age
          stop
        end if
        isoc%age=log10_ages(jage-1:jage+1)
        isoc%npts=n_masses(jage-1:jage+1)
      end if

      ! Now get the names of the column heads.
      rewind(2)
      col_head=' '
      each_line4: do 
        read(2,'(a150)',iostat=iostat) test_str
        test_str=adjustl(test_str)
        ! Remove the tabs.
        do ihed=1, 150
          if (test_str(ihed:ihed) == char(9)) test_str(ihed:ihed)=' '
        end do
        if (test_str(1:1) /= '#') exit each_line4
        work_str='#'
        each_header: do ihed=1, 20
          ipos=index(test_str, trim(work_str))+len_trim(work_str)
          test_str=test_str(ipos:len(test_str))
          test_str=adjustl(test_str)
          read(test_str(1:index(test_str, ' ')) ,'(20a)',
     &    iostat=iostat) col_head(ihed)
          if (iostat<0 .or. len_trim(col_head(ihed))<1) exit
          work_str=col_head(ihed)
          if (col_head(ihed) == 'Mu') col_head(ihed)='U'
          if (col_head(ihed) == 'Mb') col_head(ihed)='B'
          if (col_head(ihed) == 'Mv') col_head(ihed)='V'
          if (col_head(ihed) == 'Mr') col_head(ihed)='R'
          if (col_head(ihed) == 'Mi') col_head(ihed)='I'
          if (col_head(ihed) == 'Mj') col_head(ihed)='J'
          if (col_head(ihed) == 'Mh') col_head(ihed)='H'
          if (col_head(ihed) == 'Mk') col_head(ihed)='K'
          nhed=ihed
        end do each_header
      end do each_line4

      mag_plus='NOT SET'
      mag_minus='NOT SET'
      col_plus='NOT SET'
      col_minus='NOT SET'
      if (colflag == 0) then
        mag_plus='logL/Lo'
        col_plus='logTe'
        icol_minus=-1
        imag_minus=-1
      else
        ipos=index(mag_name(colflag), '-')
        if (ipos > 0) then
          mag_plus =mag_name(colflag)(1:ipos-1)
          mag_minus=mag_name(colflag)(ipos:len(mag_name(colflag)))
          imag_minus=0
        else
          mag_plus = mag_name(colflag)
          imag_minus=-1
        end if
        ipos=index(col_name(colflag), '-')
        if (ipos > 0) then
          col_plus =col_name(colflag)(1:ipos-1)
          col_minus=col_name(colflag)(ipos+1:len(col_name(colflag)))
          icol_minus=0
        else
          col_plus =col_name(colflag)
          icol_minus=-1
        end if
      end if

      imag_plus=-1
      icol_plus=-1
      do ihed=1, nhed
        if (trim(col_head(ihed)) == 'M_ini') imass=ihed
        if (trim(col_head(ihed)) == 'logG') igrav=ihed
        if (trim(col_head(ihed)) == 'logL/Lo') ilbol=ihed
        if (trim(col_head(ihed)) == 'logTe') iteff=ihed
        if (trim(col_head(ihed)) == mag_plus) imag_plus=ihed
        if (trim(col_head(ihed)) == col_plus) icol_plus=ihed
        if (trim(col_head(ihed)) == mag_minus .and. imag_minus>=0) 
     &  imag_minus=ihed
        if (trim(col_head(ihed)) == col_minus .and. icol_minus>=0) 
     &  icol_minus=ihed
      end do

      if (imag_plus == -1) then
        print*, 'Could not find the colour ', trim(mag_plus) 
        print*, 'amoungst ', col_head(1:nhed)
        stop
      end if
      if (icol_plus == -1) then
        print*, 'Could not find the colour ', trim(col_plus) 
        print*, 'amoungst ', col_head(1:nhed)
        stop
      end if

      ! Now get the three ages we want.
      rewind(2)
      allocate(icnt(size(isoc)))
      icnt=0
      old_mass=0.0
      each_line3: do
        read(2,'(a150)',iostat=iostat) test_str
        test_str=adjustl(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, new_mass
        !j=minval(minloc(abs(new_age-(log10(isoc%age)+6.0))))
        j=minloc(abs(new_age-isoc%age),1)
        if (abs(new_age-isoc(j)%age) < 2.0*tiny(new_age)) then
          if (abs(old_mass-new_mass) > 2.0*tiny(new_mass)) then 
            icnt(j)=icnt(j)+1
            read(test_str,*) data(1:nhed)
            isoc(j)%lbol(icnt(j))=data(ilbol)
            isoc(j)%teff(icnt(j))=data(iteff)
            isoc(j)%mass(icnt(j))=data(imass)
            isoc(j)%logg(icnt(j))=data(igrav)
            if (icol_minus > 0) then
              isoc(j)%col(icnt(j))=data(icol_plus)-data(icol_minus)
            else
              isoc(j)%col(icnt(j))=data(icol_plus)
            end if
            if (imag_minus > 0) then
              isoc(j)%mag(icnt(j))=data(imag_plus)-data(imag_minus)
            else
              isoc(j)%mag(icnt(j))=data(imag_plus)
            end if
!           print*, new_age, icnt(j), isoc(j)%mass(icnt(j)), 
!     &     isoc(j)%mag(icnt(j)), isoc(j)%col(icnt(j))
          end if
          old_mass=new_mass
        end if
      end do each_line3
      deallocate(icnt)

!      if (colflag == 0) then
!        print*, 'Searched for logL/Lo and logTe amoungst'
!      else
!        print*, 'Searched for ', mag_name(colflag), ' and ', 
!     &  col_name(colflag), ' amoungst '
!      end if
!      print*, col_head(1:nhed)
!      print*, 'Needed columns ', imag_plus, imag_minus, icol_plus, 
!     &icol_minus
!      print*, 'Corresponding to ', mag_plus, mag_minus, col_plus, 
!     &col_minus

      deallocate(n_masses, log10_ages)
      close(2)

      !do j=1, 3
      !  write(30+j,*) isoc(j)%age
      !  write(30+j,*) trim(ifname)
      !  write(30+j,*) ' '
      !  do jage=1, isoc(j)%npts
      !    write(30+j,*) isoc(j)%col(jage), isoc(j)%mag(jage)
      !  end do
      !end do
      !stop

      end subroutine read_padova_col



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

      character(len=30) :: cflag
      real :: work1, work2

      flag=0
      if (vi > 4.2) then
        flag=1
      else if (vi > 1.0) then  
        work1 = -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) then
          cflag='OK'
          work2=bessell_bc_v(logg, teff, cflag)
          if (trim(cflag) =='OK') work1=max(work1,work2)
        end if
      else
        cflag='OK'
        work1=bessell_bc_v(logg, teff, cflag)
        if (trim(cflag) /='OK') flag=1
      end if

      bcvi=work1

      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
