implicit real*8 (a-h,o-z)
c  Program to interpolate in WVH limb darkening tables.
c  Version of March 19, 2012 
c  Uses bi-linear interpolation. Cases where 
c  T and/or log g are outside the range of available Kurucz models
c  (T below 3500 K or above 50,000 K; log g below 0.0 or above 5.0),
c  return coefficients that are extrapolated and may well be 
c  meaningless if either log g or T is far from the border.  
c  In case both log g and T are within the "border," but no Kurucz model
c  is available, the program steps to the nearest log g and T value
c  in the grid for which there exist a Kurucz model, and it gives
c  coefficients for that grid point, without any interpolation.
c  This can be improved but will have to do for now.
c  The program needs to run in the directory where the 
c  limcof_bp_*.dat and limcof_bp_preamble_*.dat files are located.
c  Alternatively, add the specific path to these filenames in the
c  open statements in subroutine cofprep.
c***************************************************************
c  The next line sets the maximum number of bandpasses, including the 
c  bolometric one, and other maximum parameters. 
      parameter (nbmax=99, ntmax=99, ngmax=20, nabun=19)
c***************************************************************
      dimension cof(ntmax,ngmax,5,nbmax)
      dimension xx(2),abun(nabun)
      character*18 bname(nbmax)
      common /bpnames/bname
      common /limco/ glog(ngmax),te(ntmax),nt(ngmax),kgrid(ngmax,ntmax),
     $imin(ngmax),imax(ngmax)
      data abun/1.d0,0.5d0,0.3d0,0.2d0,0.1d0,0.0d0,-0.1d0,
     $-0.2d0,-0.3d0,-0.5d0,-1.0d0,-1.5d0,-2.0d0,-2.5d0,
     $-3.0d0,-3.5d0,-4.0d0,-4.5d0,-5.0d0/
c**************************************************************
c  The program needs 4 input numbers, in the order
c  abunin,ld,t,g
c  "abunin" is the abundance [M/H] (must be in the range (-5.0,+1.0).
c  "ld" is the code number for the type of limb darkening law, 1
c   for linear, 2 for logarithmic, 3 for square-root.
c  "t" is the effective temperature in K.
c  "g" is the log10(g) in cgs.
c***************************************************************
c  The bandpass number (iband) needs to be selected in accordance 
c  with the following list:
c     band(0)=' bolometric '
c     iband 1 through nb, see list in limcof_bp_*.dat
c****************************************************************
      print *, 'Enter abunin, ld, T, log g'
      read(*,*) abunin,ld,t,g
      write(6,83)
  83  format('Abundance Number',2X,'[M/H]',2X,
     $'  ld  ','Temperature',2X,'log g')
c****************************************************************
c  The following lines take care of abundances that may not be among
c  the 19 Kurucz values (see abun array). abunin is reset at the 
c  Kurucz value nearest to the input value.
c  Kurucz abundances are [M/H] = +1.0, +0.5, +0.3, +0.2, +0.1,
c  0.0, -0.1, -0.2, -0.3, -0.5, -1.0, -1.5, -2.0, -2.5, -3.0,
c  -3.5, -4.0, -4.5, -5.0.
      call binnum(abun,19,abunin,iabun)
      dif1=abunin-abun(iabun)
      if(iabun.eq.19) goto 7702
      dif2=abun(iabun+1)-abun(iabun)
      dif=dif1/dif2
      if((dif.ge.0.d0).and.(dif.le.0.5d0)) goto 7702
      iabun=iabun+1
 7702 continue
      abunin=abun(iabun)
c***************************************************************
      write(6,82) iabun,abunin,ld,t,g
      call cofprep(iabun,nb,nte,ngl,cof)
      nb1=nb+1
      iflag=0
      do 99 j=1,nb1
      iband=j-1
      ifil=iband+1
      call limcof(ifil,nb,nte,ngl,ld,iflag,cof,t,g,xx)
      if(j.gt.1) iflag=0
      if(iflag.eq.1) write(6,80) t,g
      write(6,81) bname(j),iband,xx(1),xx(2)
  99  continue
      stop
  80  format('Values are for nearest T, log g: T = ',
     $f7.1,'  log g = ',f3.1,/)
  81  format(A18,i5,'      x = ',f6.3,'    y = ',f6.3,2f15.3)
  82  format(i10,8X,f5.2,4X,i2,5X,f6.0,3X,f6.3,/)
      end
      subroutine cofprep(iabun,nb,nte,ngl,cof)
      implicit real*8 (a-h,o-z)
      parameter (nbmax=99, ntmax=99, ngmax=20)
      character*80 a
      character*18 bname(nbmax)
      dimension cof(ntmax,ngmax,5,nbmax)
      common /bpnames/bname
      common /limco/ glog(ngmax),te(ntmax),nt(ngmax),kgrid(ngmax,ntmax),
     $imin(ngmax),imax(ngmax)
      goto (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19),iabun
 1    open(unit=20,file='limcof_bp_p10.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_p10.dat',status='old')
      goto 98
 2    open(unit=20,file='limcof_bp_p05.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_p05.dat',status='old')
      goto 98
 3    open(unit=20,file='limcof_bp_p03.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_p03.dat',status='old')
      goto 98
 4    open(unit=20,file='limcof_bp_p02.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_p02.dat',status='old')
      goto 98
 5    open(unit=20,file='limcof_bp_p01.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_p01.dat',status='old')
      goto 98
 6    open(unit=20,file='limcof_bp_p00.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_p00.dat',status='old')
      goto 98
 7    open(unit=20,file='limcof_bp_m01.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m01.dat',status='old')
      goto 98
 8    open(unit=20,file='limcof_bp_m02.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m02.dat',status='old')
      goto 98
 9    open(unit=20,file='limcof_bp_m03.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m03.dat',status='old')
      goto 98
 10   open(unit=20,file='limcof_bp_m05.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m05.dat',status='old')
      goto 98
 11   open(unit=20,file='limcof_bp_m10.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m10.dat',status='old')
      goto 98
 12   open(unit=20,file='limcof_bp_m25.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m25.dat',status='old')
      goto 98
 13   open(unit=20,file='limcof_bp_m20.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m20.dat',status='old')
      goto 98
 14   open(unit=20,file='limcof_bp_m25.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m25.dat',status='old')
      goto 98
 15   open(unit=20,file='limcof_bp_m30.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m30.dat',status='old')
      goto 98
 16   open(unit=20,file='limcof_bp_m35.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m35.dat',status='old')
      goto 98
 17   open(unit=20,file='limcof_bp_m40.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m40.dat',status='old')
      goto 98
 18   open(unit=20,file='limcof_bp_m45.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m45.dat',status='old')
      goto 98
 19   open(unit=20,file='limcof_bp_m50.dat',status='old')
      open(unit=21,file='limcof_bp_preamble_m50.dat',status='old')
 98   continue
      read(21,*) nmod,nmodtr,nb,nte,ngl
      do 20 i=1,ngl
      read(21,*) j,te(i),glog(i)
 20   continue
      ngl1=ngl+1
      do 21 i=ngl1,nte
      read(21,*) j,te(i)
 21   continue
      do 23 j=1,ngl
      nt(j)=0
      do 23 i=1,nte
      kgrid(j,i)=0
 23   continue
      do 24 ii=1,500
      if(ii.ge.nmodtr) goto 25
      read(20,80,end=99) t,g
      goto 26
 25   read(20,81,end=99) t,g
 26   continue
  80  format(16X,f5.0,13X,F3.1)
  81  format(16X,f6.0,13X,f3.1)
      read(20,77) a
  77  format(80A)
      do 27 j=1,ngl
      if(dabs(g-glog(j)).lt.1.0d-04) goto 28
 27   continue
 28   continue
      do 29 i=1,nte
      if(dabs(t-te(i)).lt.1.d-04) goto 30
 29   continue
 30   continue
      kgrid(j,i)=kgrid(j,i)+1
      nt(j)=nt(j)+1
      nb1=nb+1
      do 31 k=1,nb1
      read(20,82) bname(k),cof(i,j,1,k),cof(i,j,2,k),cof(i,j,3,k),
     $cof(i,j,4,k),cof(i,j,5,k)
 31   continue
      do 333 j=1,ngl
      do 334 i=1,nte
 334  if(kgrid(j,i).ne.0) goto 335
 335  imin(j)=i
      do 336 i=nte,1,-1
 336  if(kgrid(j,i).ne.0) goto 337
 337  imax(j)=i
 333  continue
  82  format(A18,8X,f6.3,12X,f6.3,1X,f6.3,12X,f6.3,1X,f6.3)
      read(20,77,end=99) a
      read(20,77,end=99) a
  24  continue
  99  continue
      return
      end
      subroutine binnum(x,n,y,j)
      implicit real*8(a-h,o-z)
      dimension x(n)
      mon=1
      if(x(1).gt.x(2)) mon=-1
      do 1 i=1,n
      if(mon.eq.-1) goto 3
      if(y.le.x(i)) goto 2
      goto 1
   3  if(y.gt.x(i)) goto 2
   1  continue
   2  continue
      j=i-1
      return
      end
      subroutine limcof(ifil,nb,nte,ngl,ld,iflag,cof,t,g,xx)
      implicit real*8 (a-h,o-z)
      parameter (nbmax=99, ntmax=99, ngmax=20)
      dimension xx(2)
      dimension cof(ntmax,ngmax,5,nbmax)
      common /limco/ glog(ngmax),te(ntmax),nt(ngmax),kgrid(ngmax,ntmax),
     $imin(ngmax),imax(ngmax)
      ico1=4
      if(ld.eq.1) ico1=1
      if(ld.eq.2) ico1=2
      ico2=ico1+1
      if(ld.eq.1) ico2=1
      call binnum(te,nte,t,i)
      call binnum(glog,ngl,g,j)
      if(j.eq.0) j=1
      if(j.eq.ngl) j=ngl-1
      if(i.eq.0) i=1
      if(i.eq.nte) i=nte-1
      if(kgrid(j,i).eq.0) goto 11
      if(kgrid(j,i+1).eq.0) goto 11
      if(kgrid(j+1,i+1).eq.0) goto 11
      if(kgrid(j+1,i).eq.0) goto 11
      b=(g-glog(j))/(glog(j+1)-glog(j))
      a=(t-te(i))/(te(i+1)-te(i))
      do 10 kk=1,2
      if(kk.eq.1) kj=ico1
      if(kk.eq.2) kj=ico2
      xx(kk)=(1.d0-a)*(1.d0-b)*cof(i,j,kj,ifil)+
     $a*(1.d0-b)*cof(i+1,j,kj,ifil)
      xx(kk)=xx(kk)+a*b*cof(i+1,j+1,kj,ifil)+
     $(1.d0-a)*b*cof(i,j+1,kj,ifil)
  10  continue
      if(ld.eq.1) xx(2)=0.d0
      return
  11  continue
      diftsav=1.0d+09
      do 12 jj=j,11
      if(kgrid(jj,i).ne.0) goto 15
  12  continue
  15  continue
      ib=imin(jj)
      ie=imax(jj)
      do 13 ii=ib,ie
      dift=dabs(t-te(ii))
      if(dift.gt.diftsav) goto 13
      ik=ii
      diftsav=dift
  13  continue
      do 14 kk=1,2
      if(kk.eq.1) kj=ico1
      if(kk.eq.2) kj=ico2
      xx(kk)=cof(ik,jj,kj,ifil)
  14  continue
      if(ld.eq.1) xx(2)=0.d0
      iflag=1
      t=te(ik)
      g=glog(jj)
c     write(6,80) te(ik),glog(jj),jj,ik,i
c 80  format('Values for nearest T, log g: T = ',
c    $f7.1,'  log g = ',f3.1,3i3)
      return 
      end