C *** This program compares the interpolated values of the bolometric
C *** corrections with those given in the files created by executing
C *** gettestbcs.for (which reproduce the BC values from the synthetic
C *** magnitudes as derived directly from the MARCS model atmosphere grids.
C
C *** Unit 33 contains the tabulated BCs that are interpolated: the name of
C *** the file is "bc_std.data", "bc_p04.data", "bc_p00.data", or
C *** "bc_m04.data", depending on whether the value of ialf in selectbc.data
C *** is 1, 2, 3, or 4, respectively.  Output is sent to unit 6, which is
C *** assumed to be attached to the monitor.
C 
C *** Very small adjustments are made to the input gv and Teff values
C *** to ensure that they stay within the set boundaries of gravity and
C *** temperature within which interpolations are performed.  Otherwise
C *** round-off or integer-to-real conversions causes input entries to
C *** be outside of those boundaries sometimes.
C
C *** the total number of interpolations (n) should equal the total
C *** number of data points in the input file minus number of entries
C *** for 2500 K, for [Fe/H] = -5.0, and for Teff values > the maximum
C *** value considered at a given value of log g.
C
      character*13 fnme
      fnme(1:13)='selectbc.data'
      call iofile(9,'in',fnme,13)
      read(9,*) ialf
      read(9,*) nfil
      if(ialf.gt.1) go to 5
      call intrp_s
      go to 10
    5 call intrp_a(ialf)
   10 stop
      end
      subroutine intrp_s
      real*4 g(13),bcin(5),bc(5),dif(5),bcx(5),dfx(5)
      integer*4 itm(13),num(13)
      character*13 fnme,fzero
      character*6 filter(41),gravity(13),fgrav,fil(5),fff(5),ff0
      character*5 ggg(5),gg1,gg0,ftype
      data g/-0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5/,
     1  itm/4250,5000,5500,6000,6500,6750,7750,8000,8000,8000,8000,
     2  8000,5000/,num/57,161,237,278,318,332,400,411,408,418,412,
     3  415,169/,gravity/'sphm05','sph00 ','sph05 ','sph10 ','sph15 ',
     4  'sph20 ','sph25 ','sph30 ','sph35 ','ppl40 ','ppl45 ','ppl50 ',
     5  'ppl55 '/,ftype/'.data'/,fzero/'             '/,gg1/'delta'/,
     6  filter/'  J   ','  H   ','  K   ','acs435','acs475','acs555',
     7  'acs606','acs814',' f218w',' f225w',' f275w',' f336w','f350lp',
     8  ' f390m',' f390w',' f438w',' f475w',' f547m',' f555w',' f606w',
     9  ' f625w',' f775w',' f814w','f850lp',' ir98m','ir110w','ir125w',
     *  'ir140w','ir160w','sdss_u','sdss_g','sdss_r','sdss_i','sdss_z',
     *  '  U   ','  B   ','  V   ','  R   ','  I   ','  UX  ','  BX  '/,
     *  ff0/'      '/,gg0/'     '/
 1000 format(i4,2f6.2,5f8.4)
 1001 format(' bc(i),bcin(i):',10f8.4)
 1002 format(/,8x,'input/',/,' log g',1x,'checked',2x,'Teff  [Fe/H]',
     1  5(2x,a6,3x,a5))
 1003 format(9x,f5.2)
 1004 format(i3)
 1005 format(f5.1,2x,i3,'/',i3,2x,i4,f7.2,1x,10f8.4)
c *** read the input bolometric correction tables (= bc_std.data)
      fe=0.0
      call getbc_s(0,fe,gv,teff,ebv,bc,fil,nbc)
      do 5 i=1,nbc
      fff(i)=ff0
      ggg(i)=gg0
    5 continue
      do 8 i=1,nbc
      read(9,*) isys,ifil
      fff(i)=filter(ifil)
      ggg(i)=gg1
    8 continue
      nf1=9
      write(6,1002) (fff(i),ggg(i),i=1,nbc)
      do 45 j=1,13
      nf1=nf1+1
      fnme=fzero
      fgrav=gravity(j)
      fnme(1:6)=fgrav(1:6)
      nch=6
      if(j.eq.1) go to 10
      nch=5
   10 fnme(nch+1:nch+5)=ftype(1:5)
      nch=nch+5
      call iofile(nf1,'in',fnme,nch)
      gv=g(j)
      itx=itm(j)
      n=0
      dmax=0.0
      read(nf1,1003) redden
      read(nf1,1004) npt
      do 40 m=1,npt
      read(nf1,1000) ite,fe,afe,(bcin(k),k=1,nbc)
      temp=float(ite)
      teff=alog10(temp)
      if(fe.lt.-4.0) go to 40
      if(ite.lt.2600.or.ite.gt.itx) go to 40
      call getbc_s(-1,fe,gv,teff,ebv,bc,fil,nbc)
      call getbc_s(1,fe,gv,teff,ebv,bc,fil,nbc)
      difmax=0.0
      do 25 l=1,nbc
      dif(l)=abs(bc(l)-bcin(l))
      difmax=amax1(dif(l),difmax)
      if(difmax.lt.dmax) go to 20
      dmax=difmax
      itmx=ite
      femx=fe
      afex=afe
      do 15 k=1,nbc
      dfx(k)=dif(k)
      bcx(k)=bc(k)
   15 continue
   20 if(abs(bc(l)-bcin(l)).gt.0.001) go to 30
   25 continue
      go to 35
   30 write(6,1001) (bc(k),bcin(k),k=1,nbc)
      go to 40
   35 n=n+1
   40 continue
      write(6,1005) g(j),num(j),n,itmx,femx,(bcx(k),dfx(k),k=1,nbc)
   45 continue
      stop
      end
      subroutine intrp_a(ialf)
      real*4 g(12),bcin(5),bc(5),dif(5),bcx(5),dfx(5)
      integer*4 itm(12),num(12),nump04(12),nump00(12),numm04(12)
      character*13 fnme,fzero
      character*6 filter(41),fil(5),fff(5),ff0
      character*5 gravity(12),ggg(5),gg1,gg0,fgrav,ftype
      data g/0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5/,
     1  itm/5000,5250,5500,5750,6250,7750,8000,8000,8000,8000,8000,
     2  5000/,nump04/80,108,122,132,139,173,179,175,170,171,178,80/,
     3  nump00/87,136,146,151,158,192,210,207,211,210,212,88/,
     4  numm04/133,177,186,193,202,240,251,245,256,258,255,134/,
     5  gravity/'sph00','sph05','sph10','sph15','sph20','sph25',
     6  'sph30','sph35','ppl40','ppl45','ppl50','ppl55'/,gg1/'delta'/,
     7  ftype/'.data'/,fzero/'             '/,filter/'  J   ','  H   ',
     8  '  K   ','acs435','acs475','acs555','acs606','acs814',' f218w',
     9  ' f225w',' f275w',' f336w','f350lp',' f390m',' f390w',' f438w',
     *  ' f475w',' f547m',' f555w',' f606w',' f625w',' f775w',' f814w',
     *  'f850lp',' ir98m','ir110w','ir125w','ir140w','ir160w','sdss_u',
     *  'sdss_g','sdss_r','sdss_i','sdss_z','  U   ','  B   ','  V   ',
     *  '  R   ','  I   ','  UX  ','  BX  '/,ff0/'      '/,gg0/'     '/
 1000 format(i4,2f6.2,5f8.4)
 1001 format(' bc(i),bcin(i):',10f8.4)
 1002 format(/,8x,'input/',/,' log g',1x,'checked',2x,'Teff  [Fe/H]',
     1  5(2x,a6,3x,a5))
 1003 format(9x,f5.2)
 1004 format(i3)
 1005 format(f5.1,2x,i3,'/',i3,2x,i4,f7.2,1x,10f8.4)
      if(ialf.gt.2) go to 10
      do 5 i=1,12
      num(i)=nump04(i)
    5 continue
      go to 30
   10 if(ialf.gt.3) go to 20
      do 15 i=1,12
      num(i)=nump00(i)
   15 continue
      go to 30
   20 do 25 i=1,12
      num(i)=numm04(i)
   25 continue
c *** read the input bolometric correction tables (e.g., bc_p04.data)
   30 fe=0.0
      call getbc_a(0,ialf,fe,gv,teff,ebv,bc,fil,nbc)
      do 35 i=1,nbc
      fff(i)=ff0
      ggg(i)=gg0
   35 continue
      do 40 i=1,nbc
      read(9,*) isys,ifil
      fff(i)=filter(ifil)
      ggg(i)=gg1
   40 continue
      nf1=9
      write(6,1002) (fff(i),ggg(i),i=1,nbc)
      do 75 j=1,12
      nf1=nf1+1
      fnme=fzero
      fgrav=gravity(j)
      fnme(1:5)=fgrav(1:5)
      fnme(6:10)=ftype(1:5)
      nch=10
      call iofile(nf1,'in',fnme,nch)
      gv=g(j)
      itx=itm(j)
      n=0
      dmax=0.0
      read(nf1,1003) redden
      read(nf1,1004) npt
      do 70 m=1,npt
      read(nf1,1000) ite,fe,afe,(bcin(k),k=1,nbc)
      temp=float(ite)
      teff=alog10(temp)
      if(fe.lt.-2.5.or.fe.gt.0.5) go to 70
      if(ite.lt.2600.or.ite.gt.itx) go to 70
      call getbc_a(-1,ialf,fe,gv,teff,ebv,bc,fil,nbc)
      call getbc_a(1,ialf,fe,gv,teff,ebv,bc,fil,nbc)
      difmax=0.0
      do 55 l=1,nbc
      dif(l)=abs(bc(l)-bcin(l))
      difmax=amax1(dif(l),difmax)
      if(difmax.lt.dmax) go to 50
      dmax=difmax
      itmx=ite
      femx=fe
      afex=afe
      do 45 k=1,nbc
      dfx(k)=dif(k)
      bcx(k)=bc(k)
   45 continue
   50 if(abs(bc(l)-bcin(l)).gt.0.001) go to 60
   55 continue
      go to 65
   60 write(6,1001) (bc(k),bcin(k),k=1,nbc)
      go to 70
   65 n=n+1
   70 continue
      write(6,1005) g(j),num(j),n,itmx,femx,(bcx(k),dfx(k),k=1,nbc)
   75 continue
      stop
      end
      Subroutine getbc_s(mode,fe,gv,teff,ebv,bc,fil,nbc)
C ----------------------------------------------------------------------
C *** This routine interpolates in the tables of bolometric corrections
C     derived from MARCS model atmospheres for up to 3 filter bandpasses
C     of choice.  The standard variation of [alpha/Fe] with [Fe/H] is
C     assumed; i.e., [alpha/Fe] = +0.4 at [Fe/H] <= -1.0, 0.0 at [Fe/H]
C     >= 0.0, and [alpha/Fe] = -0.4[Fe/H] at intermediate iron abundances.
C     The input tables, in a file named `bc_std.data', are assigned to
C     unit 30.
C
C *** If MODE=0, the BC tables are read and interpolated to the desired
C     value of [Fe/H] (= fe in the argument list).  MODE must be set to
C     zero the first time that this subroutine is called.  A message is
C     sent to unit 6 (which is assumed to be attached to the monitor) to
C     specify the [Fe/H] value and the magnitudes (e.g., B, I, F606W) for
C     which bolometric corrections are provided in the input tables.
C
C *** If MODE=1, interpolations are carried out to obtain the bolometric
C     corrections at the input values of log g and log Te (gv and teff in
C     the argument list) and the specified [Fe/H] value.  It is only for
C     this value of MODE that interpolations are performed to obtain BC
C     information.
C
C *** If MODE=-1, the input tables (which were previously read by setting
C     MODE=0) are re-interpolated for a new value of [Fe/H].  If desired,
C     the subroutine may be modified (after the statement `45 continue')
C     so that a message is sent to unit 6 each time an [Fe/H] interpolation
C     is performed.
C
C *** The interpolation code will issue a STOP code if (i) the input [Fe/H]
C     value is outside the range of the input tables, (ii) the input log g
C     or log Teff values are outside the tabulated ranges, or (ii) the
C     input file (`bc_std.data') does not exist.
C -------------------------------------------------------------------7--
      real*4 ttt(8),bc(5)
      character*29 alfe(15)
      character*10 namfe
      character*6 fil(5)
      common /bcstd/ a(31,13,5,15),b(31,13,5),c(4,5),te(31),t(31),g(13),
     1  feh(15),ttm(8),r(4),af(4),ag(4),at(4),maxt(13),nt,ng,nfe,ndx
      data ttt/4250.,5000.,5500.,6000.,6500.,6750.,7750.,8000./
 1000 format(/,i3,14x,i2,14x,i2,14x,i2,21x,f6.3)
 1001 format(13f6.0)
 1002 format(20f4.1)
 1003 format(10f8.4)
 1004 format(' BCs for E(B-V) =',f6.3,' interpolated to [Fe/H] =',f6.2,
     1  /,' filters: ',5(a6,2x))
 1005 format(' log g =',f9.6,'  or Teff =',f9.3,' are outside the',
     1  1x,'range of the input tables')
 1006 format(a10,f5.2,a29,a6)
 1007 format(' input  fe =',f5.2,' is outside the range of table',i2)
 1008 format(20i4)
      if(mode.gt.0) go to 50
      if(mode.lt.0) go to 25
C *** The input BC tables are read only if MODE=0.
      open(unit=33,iostat=ierr,file='bc_std.data',status='old')
      if(ierr.ne.0) go to 150
      read(33,1000) nt,ng,nfe,ndx,ebv
      read(33,1001) (te(i),i=1,nt)
      read(33,1002) (g(i),i=1,ng)
      read(33,1008) (maxt(i),i=1,ng)
      nbc=ndx
      do 5 i=1,8
      ttm(i)=alog10(ttt(i))
    5 continue
      do 10 i=1,nt
      t(i)=alog10(te(i))
   10 continue
      do 20 k=1,ndx
      do 20 l=1,nfe
      read(33,1006) namfe,feh(l),alfe(l),fil(k)
      do 15 j=1,ng
      nend=maxt(j)
      read(33,1003) (a(i,j,k,l),i=1,nend)
   15 continue
   20 continue
      close(unit=33,status='keep')
C *** When MODE=0 OR MODE=-1, the input BC table is interpolated to
C *** create a table of BC values for the input value of [Fe/H) (=fe).
   25 if(fe.lt.feh(1).or.fe.gt.feh(nfe)) go to 35
      nfem1=nfe-1
      do 30 m=3,nfem1
      l=m-2
      if(fe.le.feh(m)) go to 40
   30 continue
      go to 40
   35 write(6,1007) fe
      stop ' input [Fe/H] value is outside the range of the tables'
   40 r(1)=feh(l)
      r(2)=feh(l+1)
      r(3)=feh(l+2)
      r(4)=feh(l+3)
      call lgr4(r,af,fe)
      do 45 k=1,ndx
      do 45 j=1,ng
      nend=maxt(j)
      do 45 i=1,nend
      b(i,j,k)=af(1)*a(i,j,k,l)+af(2)*a(i,j,k,l+1)+af(3)*a(i,j,k,l+2)+
     1  af(4)*a(i,j,k,l+3)
   45 continue
C *** If desired, the following statement may be commented out so that
C *** the write statement is executed after each [Fe/H] interpolation.
      if(mode.ne.0) go to 155
      write(6,1004) ebv,fe,(fil(k),k=1,ndx)
C *** When MODE=0 or MODE=-1, control returns to the calling program
C *** once the [Fe/H] interpolations have been completed; i.e., no
C *** interpolations are performed within the tables.
      go to 155
C *** MODE=1: interpolations are performed for the input values of
C *** log g (= gv) and log Teff (= teff).
   50 nbc=ndx
      temp=10.**teff
      if(abs(teff-t(nt)).ge.1.e-5) go to 55
      teff=t(nt)
   55 if(abs(teff-t(1)).ge.1.e-5) go to 60 
      teff=t(1)
   60 if(teff.le.t(nt).and.teff.ge.t(1)) go to 70
   65 write(6,1005) gv,temp
      stop ' range of the tables exceeded'
   70 if(gv.lt.-0.5.or.gv.gt.5.5) go to 65
C *** Note that ngmin=1 only if the input value of gv is < 0.0.
      ngmin=2
      if(gv.ge.0.0) go to 75
      ngmin=1
C *** This loop determines the index ngm where ng(ngm) is the lowest
C *** value of log g for which data are available for the input Teff.
   75 do 80 i=ngmin,8
      ngm=i
      if(teff.le.ttm(i)) go to 85
   80 continue
   85 ngl=ngm+2
C *** Note that ngmax=13 only if the input value of gv is > 5.0.
      ngmax=13
      if(gv.gt.5.0) go to 90
      ngmax=12
   90 ngm=ngmax-1
C *** This loop determines the index mg where g(mg) is the lowest 
C *** of the grid values of log g that are used in the interpolations.
      do 95 i=ngl,ngm
      mg=i-2
      if(gv.le.g(i)) go to 105
   95 continue
      mg=ngmax-3
  105 mmx=maxt(mg)
C *** This section increases the value of mg if the maximum value of
C *** Teff at the initial value of mg is less than the input Teff.
      if(teff.le.t(mmx)) go to 110
      mg=mg+1
      mmx=maxt(mg)
C *** This section tests whether the maximum value of Teff at the
C *** revised value of mg is >= the input Teff value.  If this condition
C *** is not satisfied, the execution terminates.
      if(teff.gt.t(mmx)) go to 65
  110 r(1)=g(mg)
      r(2)=g(mg+1)
      r(3)=g(mg+2)
      r(4)=g(mg+3)
C *** The Lagrangian coefficients for the log g interpolation are computed.
      call lgr4(r,ag,gv)
      ntm=nt-1
C *** This loop determines the value of mt where t(mt) is the lowest of
C *** the grid values of log Teff that are used in the interpolations.
      do 115 i=3,ntm
      mt=i-2
      if(teff.le.t(i)) go to 120
  115 continue
      mt=nt-3
C *** This section checks that the input T is <= the maximum value for
C *** which tabulated data are available at the lowest value of mg 
C *** in the interpolations.
  120 limt=maxt(mg)
      if(gv.le.5.0) go to 125
      limt=19
      if(teff.gt.t(limt)) go to 65
  125 if((mt+3).le.limt) go to 130
      mt=limt-3
  130 r(1)=t(mt)
      r(2)=t(mt+1)
      r(3)=t(mt+2)
      r(4)=t(mt+3)
C *** The Lagrangian coefficients for the log Teff interpolations are computed.
      call lgr4(r,at,teff)
      l=mg-1
C *** The desired bolometric corrections are calculated.
      do 140 i=1,4
      l=l+1
      do 135 k=1,ndx
      c(i,k)=at(1)*b(mt,l,k)+at(2)*b(mt+1,l,k)+at(3)*b(mt+2,l,k)+
     1   at(4)*b(mt+3,l,k)
  135 continue
  140 continue
      do 145 i=1,ndx
      bc(i)=ag(1)*c(1,i)+ag(2)*c(2,i)+ag(3)*c(3,i)+ag(4)*c(4,i)
  145 continue
      go to 155
  150 stop ' input file bc_std.data does not exist'
  155 return
      end
      subroutine iofile(nunit,fstat,fnme,nch)
      integer*4 nunit,nch
      character*13 fnme
      character*2 fstat
 1000 format(' unit',i3,' input file does not exist',/,'name = ',a52)
      if(fstat.eq.'in') go to 10
c *** the named file is apparently an output file
      open(unit=nunit,iostat=ierr,file=fnme(1:nch),status='new')
      if(ierr.eq.0) go to 15
c *** if the named file already exists it is deleted and opened as a
c *** new file; i.e., the existing file will be overwritten
      open(unit=nunit,file=fnme(1:nch),status='old')
      close(unit=nunit,status='delete')
      open(unit=nunit,file=fnme(1:nch),status='new')
      go to 15
c *** the named file is an input file 
   10 open(unit=nunit,iostat=ierr,file=fnme(1:nch),status='old')
      if(ierr.eq.0) go to 15
c *** the named input file apparently does not exist: execution terminated
      write(6,1000) nunit,fnme
      stop ' specified input file does not exist'
   15 return
      end
      Subroutine getbc_a(mode,ialf,fe,gv,teff,ebv,bc,fil,nbc)
C ----------------------------------------------------------------------
C *** This routine interpolates in the tables of bolometric corrections
C     derived from MARCS model atmospheres for up to 3 filter bandpasses
C     of choice.  The input tables, in a file named `bc_p04.data' which is
C     assigned to unit 33, assume that [alpha/Fe] = +0.4 for all [Fe/H]
C     values.
C
C *** If MODE=0, the BC tables are read and interpolated to the desired
C     value of [Fe/H] (= fe in the argument list).  MODE must be set to
C     zero the first time that this subroutine is called.  A message is
C     sent to unit 6 (which is assumed to be attached to the monitor) to
C     specify the [Fe/H] value and the magnitudes (e.g., B, I, F606W) for
C     which bolometric corrections are provided in the input tables.  
C
C *** If MODE=1, interpolations are carried out to obtain the bolometric
C     corrections at the input values of log g and log Te (gv and teff in
C     the argument list) and the specified [Fe/H] value.  It is only for
C     this value of MODE that interpolations are performed to obtain BC
C     information.
C
C *** If MODE=-1, the input tables (which were previously read by setting
C     MODE=0) are re-interpolated for a new value of [Fe/H].  If desired,
C     the subroutine may be modified (after the statement `45 continue')
C     so that a message is sent to unit 6 each time an [Fe/H] interpolation
C     is performed.
C
C *** The interpolation code will issue a STOP code if (i) the input [Fe/H]
C     value is outside the range of the input tables, (ii) the input log g
C     or log Teff values are outside the tabulated ranges, or (ii) the
C     input file (`bc_p04.data') does not exist.
C -------------------------------------------------------------------7--
      real*4 ttt(7),bc(5)
      character*29 alfe(15)
      character*11 filen(4),fnme
      character*10 namfe
      character*6 fil(5)
      common /bcp04/ a(31,12,5,15),b(31,12,5),c(4,5),te(31),t(31),g(12),
     1  feh(15),ttm(8),r(4),af(4),ag(4),at(4),maxt(12),nt,ng,nfe,ndx
      data ttt/5000.,5250.,5500.,5750.,6250.,7750.,8000./,
     1  filen/'bc_std.data','bc_p04.data','bc_p00.data','bc_m04.data'/
 1000 format(/,i3,14x,i2,14x,i2,14x,i2,21x,f6.3)
 1001 format(13f6.0)
 1002 format(20f4.1)
 1003 format(10f8.4)
 1004 format(' BCs for E(B-V) =',f6.3,' interpolated to [Fe/H] =',f6.2,
     1  /,' filters: ',5(a6,2x))
 1005 format(' log g =',f9.6,'  or Teff =',f9.3,' are outside the',
     1  1x,'range of the input tables')
 1006 format(a10,f5.2,a29,a6)
 1007 format(' input  fe =',f5.2,' is outside the range of table',i2)
 1008 format(20i4)
      if(mode.gt.0) go to 50
      if(mode.lt.0) go to 25
C *** The input BC tables are read only if MODE=0.
      fnme=filen(ialf)
      open(unit=33,iostat=ierr,file=fnme,status='old')
      if(ierr.ne.0) go to 150
      read(33,1000) nt,ng,nfe,ndx,ebv
      read(33,1001) (te(i),i=1,nt)
      read(33,1002) (g(i),i=1,ng)
      read(33,1008) (maxt(i),i=1,ng)
      nbc=ndx
      do 5 i=1,7
      ttm(i)=alog10(ttt(i))
    5 continue
      do 10 i=1,nt
      t(i)=alog10(te(i))
   10 continue
      do 20 k=1,ndx
      do 20 l=1,nfe
      read(33,1006) namfe,feh(l),alfe(l),fil(k)
      do 15 j=1,ng
      nend=maxt(j)
      read(33,1003) (a(i,j,k,l),i=1,nend)
   15 continue
   20 continue
      close(unit=33,status='keep')
C *** When MODE=0 OR MODE=-1, the input BC table is interpolated to
C *** create a table of BC values for the input value of [Fe/H) (=fe).
   25 if(fe.lt.feh(1).or.fe.gt.feh(nfe)) go to 35
      nfem1=nfe-1
      do 30 m=3,nfem1
      l=m-2
      if(fe.le.feh(m)) go to 40
   30 continue
      go to 40
   35 write(6,1007) fe
      stop ' input [Fe/H] value is outside the range of the tables'
   40 r(1)=feh(l)
      r(2)=feh(l+1)
      r(3)=feh(l+2)
      r(4)=feh(l+3)
      call lgr4(r,af,fe)
      do 45 k=1,ndx
      do 45 j=1,ng
      nend=maxt(j)
      do 45 i=1,nend
      b(i,j,k)=af(1)*a(i,j,k,l)+af(2)*a(i,j,k,l+1)+af(3)*a(i,j,k,l+2)+
     1  af(4)*a(i,j,k,l+3)
   45 continue
C *** If desired, the following statement may be commented out so that
C *** the write statement is executed after each [Fe/H] interpolation.
      if(mode.ne.0) go to 155
      write(6,1004) ebv,fe,(fil(k),k=1,ndx)
C *** When MODE=0 or MODE=-1, control returns to the calling program
C *** once the [Fe/H] interpolations have been completed; i.e., no
C *** interpolations are performed within the tables.
      go to 155
C *** MODE=1: interpolations are performed for the input values of
C *** log g (= gv) and log Teff (= teff).
   50 nbc=ndx
      temp=10.**teff
      if(abs(teff-t(nt)).ge.1.e-5) go to 55
      teff=t(nt)
   55 if(abs(teff-t(1)).ge.1.e-5) go to 60 
      teff=t(1)
   60 if(teff.le.t(nt).and.teff.ge.t(1)) go to 70
   65 write(6,1005) gv,temp
      stop ' range of the tables exceeded'
   70 if(gv.lt.0.0.or.gv.gt.5.5) go to 65
C *** Note that ngmin=1 only if the input value of gv is < 0.5.
      ngmin=2
      if(gv.ge.0.5) go to 75
      ngmin=1
C *** This loop determines the index ngm where ng(ngm) is the lowest
C *** value of log g for which data are available for the input Teff.
   75 do 80 i=ngmin,7
      ngm=i
      if(teff.le.ttm(i)) go to 85
   80 continue
   85 ngl=ngm+2
C *** Note that ngmax=12 only if the input value of gv is > 5.0.
      ngmax=12
      if(gv.gt.5.0) go to 90
      ngmax=11
   90 ngm=ngmax-1
C *** This loop determines the index mg where g(mg) is the lowest 
C *** of the grid values of log g that are used in the interpolations.
      do 95 i=ngl,ngm
      mg=i-2
      if(gv.le.g(i)) go to 105
   95 continue
      mg=ngmax-3
  105 mmx=maxt(mg)
C *** This section increases the value of mg if the maximum value of
C *** Teff at the initial value of mg is less than the input Teff.
      if(teff.le.t(mmx)) go to 110
      mg=mg+1
      mmx=maxt(mg)
C *** This section tests whether the maximum value of Teff at the
C *** revised value of mg is >= the input Teff value.  If this condition
C *** is not satisfied, the execution terminates.
      if(teff.gt.t(mmx)) go to 65
  110 r(1)=g(mg)
      r(2)=g(mg+1)
      r(3)=g(mg+2)
      r(4)=g(mg+3)
C *** The Lagrangian coefficients for the log g interpolation are computed.
      call lgr4(r,ag,gv)
      ntm=nt-1
C *** This loop determines the value of mt where t(mt) is the lowest of
C *** the grid values of log Teff that are used in the interpolations.
      do 115 i=3,ntm
      mt=i-2
      if(teff.le.t(i)) go to 120
  115 continue
      mt=nt-3
C *** This section checks that the input T is <= the maximum value for
C *** which tabulated data are available at the lowest value of mg 
C *** in the interpolations.
  120 limt=maxt(mg)
      if(gv.le.5.0) go to 125
      limt=19
      if(teff.gt.t(limt)) go to 65
  125 if((mt+3).le.limt) go to 130
      mt=limt-3
  130 r(1)=t(mt)
      r(2)=t(mt+1)
      r(3)=t(mt+2)
      r(4)=t(mt+3)
C *** The Lagrangian coefficients for the log Teff interpolations are computed.
      call lgr4(r,at,teff)
      l=mg-1
C *** The desired bolometric corrections are calculated.
      do 140 i=1,4
      l=l+1
      do 135 k=1,ndx
      c(i,k)=at(1)*b(mt,l,k)+at(2)*b(mt+1,l,k)+at(3)*b(mt+2,l,k)+
     1   at(4)*b(mt+3,l,k)
  135 continue
  140 continue
      do 145 i=1,ndx
      bc(i)=ag(1)*c(1,i)+ag(2)*c(2,i)+ag(3)*c(3,i)+ag(4)*c(4,i)
  145 continue
      go to 155
  150 stop ' input file bc_p04.data does not exist'
  155 return
      end
      subroutine lgr4(x,a,xx)
C *** a 4-point Lagrangian interpolation subroutine
      real*4 x(4),a(4)
      r1=(x(1)-x(2))*(x(1)-x(3))*(x(1)-x(4))
      r2=(x(2)-x(1))*(x(2)-x(3))*(x(2)-x(4))
      r3=(x(3)-x(1))*(x(3)-x(2))*(x(3)-x(4))
      r4=(x(4)-x(1))*(x(4)-x(2))*(x(4)-x(3))
      a(1)=((xx-x(2))*(xx-x(3))*(xx-x(4)))/r1
      a(2)=((xx-x(1))*(xx-x(3))*(xx-x(4)))/r2
      a(3)=((xx-x(1))*(xx-x(2))*(xx-x(4)))/r3
      a(4)=((xx-x(1))*(xx-x(2))*(xx-x(3)))/r4
      return
      end
