      program msearch
      USE DFLIB
      use dfport
      character*80 line,moddir
      character*6  ziffer(2)
      character argv*60
      dimension nsel(4)
      logical loout
      integer*2 na
      common /mparm/knum,num(1000),teff(1000),glog(1000)
c      open (unit=6,file='msrchout.out')
      open (unit=6,file='msrchout.out',disp='delete')
c     ziffer is the key of 6 characters: ex. ap00k2 or am10k2
      na=1
      call getarg(na,argv)
      read(argv,*) teffs
      na=2
      call getarg(na,argv)
      read(argv,*) glogs
      na=3
      call getarg(na,argv)
      read(argv,*) xslog
c      if(xslog.ge.-5.0.and.xslog.lt.-4.5) then
      if(xslog.gt.-5.5.and.xslog.lt.-4.5) then
      ziffer(1)='am50k2'
      ziffer(2)='am45k2'
      xt=(xslog-(-5.0))/((-4.5)-(-5.0))
      else if (xslog.ge.-4.5.and.xslog.lt.-4.0) then
      ziffer(1)='am45k2'
      ziffer(2)='am40k2'
      xt=(xslog-(-4.5))/((-4.0)-(-4.5))
      else if (xslog.ge.-4.0.and.xslog.lt.-3.5) then
      ziffer(1)='am40k2'
      ziffer(2)='am35k2'
      xt=(xslog-(-4.0))/((-3.5)-(-4.0))
      else if (xslog.ge.-3.5.and.xslog.lt.-3.0) then
      ziffer(1)='am35k2'
      ziffer(2)='am30k2'
      xt=(xslog-(-3.5))/((-3.0)-(-3.5))
      else if (xslog.ge.-3.0.and.xslog.lt.-2.5) then
      ziffer(1)='am30k2'
      ziffer(2)='am25k2'
      xt=(xslog-(-3.0))/((-2.5)-(-3.0))
      else if (xslog.ge.-2.5.and.xslog.lt.-2.0) then
      ziffer(1)='am25k2'
      ziffer(2)='am20k2'
      xt=(xslog-(-2.5))/((-2.0)-(-2.5))
      else if (xslog.ge.-2.0.and.xslog.lt.-1.5) then
      ziffer(1)='am20k2'
      ziffer(2)='am15k2'
      xt=(xslog-(-2.0))/((-1.5)-(-2.0))
      else if (xslog.ge.-1.5.and.xslog.lt.-1.0) then
      ziffer(1)='am15k2'
      ziffer(2)='am10k2'
      xt=(xslog-(-1.5))/((-1.0)-(-1.5))
      else if (xslog.ge.-1.0.and.xslog.lt.-0.5) then
      ziffer(1)='am10k2'
      ziffer(2)='am05k2'
      xt=(xslog-(-1.0))/((-0.5)-(-1.0))
      else if (xslog.ge.-0.5.and.xslog.lt. 0.0) then
      ziffer(1)='am05k2'
      ziffer(2)='ap00k2'
      xt=(xslog-(-0.5))/(( 0.0)-(-0.5))
      else if (xslog.ge. 0.0.and.xslog.lt. 0.5) then
      ziffer(1)='ap00k2'
      ziffer(2)='ap05k2'
      xt=(xslog-( 0.0))/(( 0.5)-( 0.0))
      else if (xslog.ge. 0.5.and.xslog.lt. 1.0) then
      ziffer(1)='ap05k2'
      ziffer(2)='ap10k2'
      xt=(xslog-( 0.5))/(( 1.0)-( 0.5))
      else
      write(6,'( ''xslog='',f9.5,'' is out of range!'')') xslog
      stop
      end if
      write(6,*)
c      call getenv('MOD_DIR',moddir)
      call retenv('MOD_DIR',moddir)
c      open (88,file='passdir.tmp')
c      read(88,'(a)') linedir
c      read(88,'(a)') moddir
c      close(88)
c      na=4
c      call getarg(na,argv)
c      read(argv,'(a)') moddir
      call lablen(moddir,80,msdir,medir)

      write(6,*) '============ Output from msearch =============='
      open (unit=34,file='models1.tmp')
      rewind 4
      do 7777 nn=1,2
      write(6,7) ziffer(nn)
    7 format(1h ,'Metallicity ziffer: ',a6)
      write(6,8) teffs,glogs
    8 format(1h ,'Input Teff:',f7.0,1x,'& log(g):',f6.3)
      open (unit=32,file=moddir(msdir:medir)//'\'//ziffer(nn)//'.lst')
      rewind 32
      kn=0
      do k=1,1000
      read(32,22,end=98) num(k),teff(k),glog(k)
      kn=kn+1
   22 format(3x,i4,20x,f7.0,9x,f8.5)
      end do
   98 knum=kn
      call gfind(teffs,glogs,nsel,irt)
      if(irt.ne.0) then
      write(6,*) 'Will skip this parameter set!' 
      stop
      else
      write(6,*) 'Four original Kurucz models used for interpolation:'
      end if
      open (unit=33,file=moddir(msdir:medir)//'\'//ziffer(nn)//'.txt')
      rewind 33
      nc=0
      model=0
      do 1000 n=1,1000000
      read(33,'(a80)',end=999) line
      if(line(1:4).eq.'TEFF') then
      if(nc.eq.4) go to 999
      model=model+1
      loout=.false.
      do j=1,4
        if(nsel(j).eq.model) then
        loout=.true.
        nc=nc+1
        write(6,'(''<<<<< model'',i2,''>>>>>'')') nc
        go to 77
        end if
      end do
   77 continue
      end if
      if(loout) write(34,'(a80)') line
c      if(loout) write(6,'(1x,a80)') line
 1000 continue
  999 continue
      close(33)
 7777 continue
      write(34,'(a3)') 'END'
      close(34)
      close(36)
      stop
      end
      subroutine gfind(teffs,glogs,nsel,irt)
      common /mparm/knum,num(1000),teff(1000),glog(1000)
      dimension nsel(4),ks1(20),ks2(20)
      do j=1,4
      nsel(j)=0
      end do
      teff1=0.
      m1=0
      do 10 k=1,knum
      if(teff(k).gt.teffs) go to 11
      if(teff(k).gt.teff1) then 
      m1=k
      teff1=teff(k)
      end if
   10 continue
      write(6,'(''Assigned Teff='',f7.0,'' too large!'')') teffs
      irt=-1
      return
   11 m2=k
      teff2=teff(k)
       if(m1.eq.0) then
       write(6,'(''Assigned Teff='',f7.0,'' too small!'')') teffs
       irt=-2
       return
       end if
      k1=0
      do 120 k=m1,m1+20
      if(teff(k).ne.teff1) go to 130
      k1=k1+1
      ks1(k1)=k
  120 continue
  130 k1sum=k1
      k2=0
      do 220 k=m2,m2+20
      if(teff(k).ne.teff2) go to 230
      k2=k2+1
      ks2(k2)=k
  220 continue
  230 k2sum=k2
      do 300 k1=2,k1sum
      glog1=glog(ks1(k1))
      do 300 k2=2,k2sum
      glog2=glog(ks2(k2))
      if(glog1.ne.glog2) go to 300
      if(glog1.lt.glogs) go to 300
      n1p=ks1(k1)
      n1=n1p-1
      n2p=ks2(k2)
      n2=n2p-1
      gloga1=glog(n1)
      gloga2=glog(n2)
      if(gloga1.ne.gloga2) then
      write(6,*) 'check gloga1 & gloga2 -->',gloga1,gloga2
      irt=-999
      return
      end if
      if(gloga1.gt.glogs) go to 300
      go to 399
  300 continue
      write(6,'(''Not interpolate: log(g)='',f6.3,'' invalid range'')')
     # glogs
      irt=1
      return
  399 continue
      nsel(1)=n1
      nsel(2)=n1p
      nsel(3)=n2
      nsel(4)=n2p
      irt=0
c
      write(6,*)
      write(6,*) 'Arranged grid of model parameters:'
      write(6,444) teffs,glogs
  444 format(1h ,'<','Teff:',f7.0,3x,'log(g):',f6.2,'>')
      do i=1,4
      write(6,100) teff(nsel(i)),glog(nsel(i))
  100 format(1h ,1x,'Teff=',f7.0,3x,'log(g)=',f6.2)
      end do
c
      return
      end

      subroutine lablen(label,nlen,ns,ne)
      character*(*) label
      do n=nlen,1,-1
      if((label(n:n).ne.' ').and.(label(n:n).ne.'')) exit
      end do
      ne=n
      do n=1,nlen
      if((label(n:n).ne.' ').and.(label(n:n).ne.'')) exit
      end do
      ns=n
      return
      end

      subroutine retenv(varname,string)
      use dfport
      character*7 varname
      character*80 sptool,string
      call getenv('SPTOOL',sptool)
      call lablen(sptool,80,is,ie)
      call lablen(varname,7,ivs,ive)
      if(varname(ivs:ive).eq.'SPTOOL') then
          string=sptool(is:ie)
        else if(varname(ivs:ive).eq.'SPT_DIR') then
          string=sptool(is:ie)//'\exefiles'
        else if(varname(ivs:ive).eq.'MOD_DIR') then
          string=sptool(is:ie)//'\k2models'
        else if(varname(ivs:ive).eq.'KGF_DIR') then
          string=sptool(is:ie)//'\GF10'
        else if(varname(ivs:ive).eq.'SFX_DIR') then
          string=sptool(is:ie)//'\sflxkpno'
        else if(varname(ivs:ive).eq.'DAT_DIR') then
          string=sptool(is:ie)//'\miscdata'
        else
          write(*,*) 'I do not understand '//varname(ivs:ive)
          stop
      end if
      return
      end
