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