      program tgv
      use dflib
      use dfport
      implicit real*4 (a-h,o-z)
      parameter (MLN=330,MIT=41)
      parameter(mp=4,np=3,msp=50)
      dimension p(mp,np),x(np),y(mp)
      external fcn1
      external fcn2
      external fcn3
      common/linedata/num(MLN),code(MLN),wla(MLN),chi(MLN),gflog(MLN)
      common/numdata/mf1,mf2,mz,mt,mg,mv,ma,ml
      real*8 wla,wlax
      character*80 dfname,title,line
      character*80 fnobs
      character*3 cvf(0:1)
      common/abdata/ablcal(MLN)
      common/ewdata/ifuse(MLN),ewobslg(MLN),ewobs(MLN)
      common/tgvpass/t0,g0,v0,ift,ifg,ifv,ndim,nmod
      common/initry/t1,t2,g1,g2,v1,v2
c
c     dt,dg,dv are the interval for evaluating the initial vertices values
c     from the input t0,g0,v0
      data dt/300./,dg/0.5/,dv/0.5/
c
c .......  deltat=t_gsp/msp, deltag=g_gsp/msp, deltav=v_gsp/msp .......
      data t_gsp/250./, g_gsp/0.5/, v_gsp/0.5/
c
      dimension ap1(-msp:msp,3),sp1(-msp:msp,3)
      dimension ap2(-msp:msp,3),sp2(-msp:msp,3)
      dimension ac(-msp:msp,3),bc(-msp:msp,3)
      dimension aw(-msp:msp,3),bw(-msp:msp,3)
      dimension ifoverc(-msp:msp,3),ifoverw(-msp:msp,3)
      dimension ifovera(-msp:msp,3)
      dimension dtdim(-msp:msp),dgdim(-msp:msp),dvdim(-msp:msp)
      dimension del(3)

      common/lc2num/lctonum(9999)
      common/avsig/a1,s1,a2,s2
      equivalence (a1,ave1), (s1,sig1), (a2,ave2), (s2,sig2)
      common/metallicity/xcur,met
      character*3 met
      character*3 metdum
      character*5 metal
      character*3 mlab(6)
      data mlab/'m15','m10','m05','p00','p05','xxx'/
      common/params/c1,c2,c3,ftol
      common/nrfpas/nrf
      common/refdat/numr(MLN),ewr(MLN),abr(MLN),w_r(9999),a_r(9999)    
      common/refpar/afer,xxt,xxg,xxv,xxf,xxa1,xxs1,xxa2,xxs2,nxx1,nxx2

      character*80 labrsv(4),fnparams,moddir

      integer*2 narg1,narg2,narg3
      character*80 buf,option,head,str,optdim(MIT)
      logical loexist,lorex
      common/mzxpass/mzx
      common/ldatax/numx(MLN),codex(MLN),wlax(MLN),chix(MLN),
     & gflogx(MLN),depsl(MLN)

      mzx=0
      do n=1,MLN
      	depsl(n)=0.
      end do

      inquire(file='.\params.dat',exist=loexist)
      if(.not.loexist) then
      call getenv('MOD_DIR',moddir)
      call lablen(moddir,80,msdir,medir)
      fnparams=moddir(msdir:medir)//'\params.dat'
      else
      fnparams='.\params.dat'
      end if
      open (unit=2,file=fnparams)
      read(2,*) c1
      read(2,*) c2
      read(2,*) c3
      read(2,*) ftol
      read(2,*,err=399,end=399) mzx
      do i=1,mzx
      read(2,*) numx(i),codex(i),wlax(i),chix(i),gflogx(i)
      end do
      go to 398
  399 mzx=0
  398 continue
      close (2)

      narg1=1
      call getarg(narg1,buf)
c      read(buf,'(a)') option
      read(buf,*) ncyc

      if(ncyc*2+1.gt.MIT) then
      write(*,*) 'ncyc*2+1, MIT = ',ncyc*2+1,MIT
      stop
      end if
      do j=1,ncyc*2
      if(mod(j,2).eq.1) then
      optdim(j)='-tg'
      if(j.eq.1) optdim(j)='+tg'
      else
      optdim(j)='-v'
      end if
      end do
      optdim(ncyc*2+1)='-'
      

      narg2=2
      call getarg(narg2,buf)
      read(buf,'(a)') dfname
c      npex=index(dfname,'.obs')
c      if(npex.le.0) then
c      write(*,*) 'Data file name should end with *.obs !'
c      stop
c      end if
c      call lablen(dfname(1:npex-1),npex-1,ns,ne)
c      nsobs=ns
c      neobs=ne
      call lablen(dfname,80,ns,ne)
      nsla=index(dfname,'/')
      if(nsla.gt.0) then
      call lablen(dfname(nsla+1:80),80-nsla,nsmm,nemm)
      ns=nsmm+nsla
      ne=nemm+nsla
      end if

      nbar=index(dfname,'~')
      if(nbar.gt.0) then
        nmax=80
        if(nsla.gt.0) nmax=nsla-1
        call lablen(dfname(1:nbar-1),nbar-1,ns1,ne1)
        call lablen(dfname(nbar+1:nmax),nmax-nbar,nspp,nepp)
        ns2=nspp+nbar
        ne2=nepp+nbar
      else
        nmax=80
        if(nsla.gt.0) nmax=nsla-1
        call lablen(dfname(1:nmax),nmax,ns1,ne1)
      end if
      nsobs=ns1
      neobs=ne1
      nsref=ns2
      neref=ne2
c
c     set nrf=0 and afer=0. for non-reference mode
      nrf=0
      afer=0.
      if(nbar.gt.0) then
c     read (refs).abd file used for comparison
c     check if the file exists
      inquire(file=dfname(nsref:neref)//'.abd',exist=lorex)
      if(.not.lorex) then
        write(*,*) 
     &  'R: why not exist ? --> '//dfname(nsref:neref)//'.abd'
        stop
      end if
      do ir=1,9999
        w_r(ir)=-999.9
        a_r(ir)=-999.9
      end do
      open (unit=9,file=dfname(nsref:neref)//'.abd')      
      do n=1,4
        read(9,'(a)') labrsv(n)
      end do
      read(labrsv(3)(2:80),*) 
     & xxt,xxg,xxv,xxf,xxa1,xxs1,xxa2,xxs2,nxx1,nxx2
      afer=0.5*(xxa1+xxa2)
      nr=0
      do 791 n=1,999
      read(9,*,end=799)
     & num__,code__,wla__,chi__,gflog__,ew__,abl1__,abl2__
      nr=nr+1
      numr(nr)=num__
      ewr(nr)=ew__
      abr(nr)=abs(abl1__)
      w_r(num__)= ewr(nr)
      a_r(num__)= abr(nr)
  791 continue
  799 nrf=nr
      if(nrf.ne.nxx1+nxx2) then
        write(*,*) 'nrf=',nrf,'  nxx1=',nxx1,'  nxx2=',nxx2
        stop
      end if
      close (9)
      end if
c

cc      narg3=3
cc      call getarg(narg3,buf)
cc      read(buf,'(a)',err=904,end=904) met
cc      call lablen(met,3,mms,mme)
cc      if(mms.gt.mme) go to 904
ccc      write(*,'(''met has been set to be '',a3)') met
cc      go to 905
cc  904 continue
      met='xxx'
c      write(*,'(''The default '',a3,'' is used for met.'')')met
cc  905 continue    

      do im=1,6
      if(met.eq.mlab(im)) go to 3
      end do
      write (*,*) 
     & 'Invalid metallicity label -->: ',met
      stop
    3 continue

c     output file string coordination
      if(met.eq.'xxx') then
        ml=1
        metal='.'
      else
        ml=5
        metal='_'//met//'.'
      end if
cc*iter
      do 8900 itnum=1,ncyc*2+1
      option=optdim(itnum)
      write(*,*) 'itnum=',itnum
      write(*,'(a)') option(1:10)
cc*iter


c     find which parameter is to be varied
      ift=0
      ifg=0
      ifv=0
      if1=0
      if(index(option,'t').gt.0) ift=1
      if(index(option,'g').gt.0) ifg=1
      if(index(option,'v').gt.0) ifv=1
      if(index(option,'1').gt.0) if1=1

c     find whether initialing or continuation      
      if(index(option,'+').gt.0) then
      ifinit=1
      else
      ifinit=0
      end if
      if(index(option,'-').gt.0) then
      ifcont=1
      else
      ifcont=0
      end if
      if(ifinit+ifcont.ne.1) then
      write(*,'(''I do not understand this option -->'',a)') option
      stop
      end if

      if(ifinit.eq.1) then
      	cvf(0)='fix'
      	cvf(1)='var'
c      	write(*,'(a)') 'Your specified mode is: Teff --> '//cvf(ift)//
c     &	',  glog --> '//cvf(ifg)//', vt --> '//cvf(ifv) 
      	write(*,*) 
     &	'Please input starting Teff, log g, and vt:'
c     &	'Please input starting (or fixed)  Teff, log g, and vt:'
      	write(*,*) 
     &	'[If OK with default values (5800,4.0,1.0), press enter.]'
      	read(*,'(a)',err=549,end=549) str
        call lablen(str,80,nstr_s,nstr_e)
        if(nstr_s.ge.nstr_e) go to 549
        read(str,*) t0,g0,v0
        go to 548
  549   continue
        t0=5800.
        g0=4.0
        v0=1.
  548   continue
c
c       initialize xcur
        if(met.eq.'m15') then 
          xcur = -1.5
        else if (met.eq.'m10') then 
          xcur = -1.0 
        else if (met.eq.'m05') then 
          xcur = -0.5 
        else if (met.eq.'p00') then 
          xcur =  0.0 
        else if (met.eq.'p05') then 
          xcur = +0.5 
        else if (met.eq.'xxx') then 
          xcur =  0.0
        else 
          write(*,*) met,'<-- ???'
          stop
        end if
c
c
cc     	open (unit=9,file=dfname(ns:ne)//metal(1:ml)//'tmp')
cc     	write(9,'(f6.3,1x,a3,1x,6e15.7)') 
cc     &  xcur,met,t0,g0,v0,-1.e30,-1.e30,-1.e30
cc     	close (9)
cc      else
cc     	inquire(file=dfname(ns:ne)//metal(1:ml)//'tmp',exist=loexist)
cc     	if(.not.loexist) then
cc     	write(*,*) 'I cannot find the file --> ',
cc     &   dfname(ns:ne)//metal(1:ml)//'tmp'
cc      	stop
cc      end if
      end if
c     reads log(W) data for the first time
      ilzero=0
      tdum=0.
      gdum=0.
      vdum=0.
      edum=0.
      call abint(ilzero,tdum,gdum,vdum,edum,adum)
c
cc     now reading the current solution from file 9
cc      open (unit=9,file=dfname(ns:ne)//metal(1:ml)//'tmp')
cc      read(9,'(f6.3,1x,a3,1x,6e15.7)') 
cc     & xcur,metdum,t0,g0,v0,dumdt,dumdg,dumdv
cc      close (9)

c     numbering setup
      do j=1,9999
      lctonum(j)=-999
      end do
      do i=1,mz
      lctonum(num(i))=i
      end do
c
c     open file 3 containing the EW data
      if(nbar.gt.0) then
        fnobs=dfname(nsobs:neobs)//'~'//dfname(nsref:neref)//'.obs'
      else
        fnobs=dfname(nsobs:neobs)//'.obs'
      end if
      inquire(file=fnobs,exist=loexist)
      if(.not.loexist) then
      write(*,*)'file does not exist! -->'//fnobs
      stop
      end if
      
      open (unit=3,file=fnobs)

      ndim=ift+ifg+ifv
      nmod=100*ift+10*ifg+ifv

      t1=max(5000.,t0-dt)
      t2=min(6500.,t0+dt)
      g1=max(1.5,g0-dg)
      g2=min(4.5,g0+dg)
      v1=max(0.5,v0-dv)
      v2=min(5.0,v0+dv)

c
c     initializing
      do 13 i=1,MLN
   13 ifuse(i)=0
      n1=0
      n2=0
c     now begins to read data
      do 88 k=1,99999
      read(3,'(a)',end=99) line
      do j=1,80
      if(line(j:j).ne.' ') go to 98
      end do
      go to 88
   98 js=j
      if(ifnum(line(js:js)).eq.1) then
      	read(line,*) lcd,ewma
       	  il=lctonum(lcd)
      	if(ewma.le.0.) go to 88
          if(if1.eq.1.and.il.gt.mf1) go to 88
           ewobs(il)=ewma
      	  ewobslg(il)=log10(ewma)
      	  ifuse(il)=1
          if(il.le.mf1) then
          	n1=n1+1
           else
          	n2=n2+1
          end if
      end if
   88 continue
   99 continue
      close (3)
c
c
      if(ndim.eq.0.and.ifinit.eq.1) then
c     Case of non-iteration
      write(*,*) 
     &'Non-iteration mode: abundance calculation for given (t,g,v)'
      tf=t0
      gf=g0
      vf=v0
      call avsigcal(t0,g0,v0,ave1,sig1,ave2,sig2)
      y(1)=0.
      go to 999
      end if
c

c
c
      if(ndim.eq.0.and.ifcont.eq.1) then
c     Perturbation mode (non-iteration)
      write(*,*) 
     &'<<Perturbation mode around the solution for error estimation>>'
      deltat=t_gsp/msp
      deltag=g_gsp/msp
      deltav=v_gsp/msp
      do 310 ip=-msp,+msp
      dtdim(ip)=deltat*ip
      dgdim(ip)=deltag*ip
      dvdim(ip)=deltav*ip
      call avsigcal
     & (t0+dtdim(ip),g0,v0,ap1(ip,1),sp1(ip,1),ap2(ip,1),sp2(ip,1))
      call coeff(ac(ip,1),bc(ip,1),aw(ip,1),bw(ip,1))
      call avsigcal
     & (t0,g0+dgdim(ip),v0,ap1(ip,2),sp1(ip,2),ap2(ip,2),sp2(ip,2))
      call coeff(ac(ip,2),bc(ip,2),aw(ip,2),bw(ip,2))
      call avsigcal
     & (t0,g0,v0+dvdim(ip),ap1(ip,3),sp1(ip,3),ap2(ip,3),sp2(ip,3))
      call coeff(ac(ip,3),bc(ip,3),aw(ip,3),bw(ip,3))
  310 continue
      do 311 ih=1,3
      do 311 ip=-msp,msp
      ap1(ip,ih)=max(0.,ap1(ip,ih))
      sp1(ip,ih)=max(0.,sp1(ip,ih))
      ap2(ip,ih)=max(0.,ap2(ip,ih))
      sp2(ip,ih)=max(0.,sp2(ip,ih))
  311 continue
      cmin=+1.e30
      cmax=-1.e30
      wmin=+1.e30
      wmax=-1.e30
      do 312 i=1,mz
      if(ifuse(i).eq.0) go to 312
      if(chi(i).gt.cmax) cmax=chi(i)
      if(chi(i).lt.cmin) cmin=chi(i)
      if(ewobs(i).gt.wmax) wmax=ewobs(i)
      if(ewobs(i).lt.wmin) wmin=ewobs(i)
  312 continue

c     Now going to estimate the uncertainties of the solutions
      chispan=cmax-cmin
      ewspan=wmax-wmin
c     evaluate the Teff range within the critical limit
      fdampc=1.
      fdampw=1.
      fdampa=1.
      sigcrt_c=sp1(0,1)*fdampc
      sigcrt_w=sp1(0,1)*fdampw
      if(ap2(0,1).lt.0.) then
      sigcrt_a=sp1(0,1)*fdampa
      else
      sigcrt_a=
     & (sp1(0,1)/sqrt(float(n1))+sp2(0,1)/sqrt(float(n2)))*fdampa
      end if
      do 370 ih=1,3
      do 380 ip=-msp,msp
      ifoverc(ip,ih)=0
      ifoverw(ip,ih)=0
  380 ifovera(ip,ih)=0
      do 381 ip=-msp,msp
      if(abs(bc(ip,ih))*chispan.gt.sigcrt_c) ifoverc(ip,ih)=1 
      if(abs(bw(ip,ih))*ewspan.gt.sigcrt_w)  ifoverw(ip,ih)=1
      if(abs(ap1(ip,ih)-ap2(ip,ih)).gt.sigcrt_a) ifovera(ip,ih)=1
  381 continue
  370 continue
      dtu=0.
      dgu=0.
      dvu=0.
      do 390 ip=msp,1,-1
      itsum=ifoverc(ip,1)+ifoverw(ip,1)+ifovera(ip,1)
      igsum=ifoverc(ip,2)+ifoverw(ip,2)+ifovera(ip,2)
      ivsum=ifoverc(ip,3)+ifoverw(ip,3)+ifovera(ip,3)
      if(itsum.eq.0.and.dtu.eq.0.) dtu=dtdim(ip)
      if(igsum.eq.0.and.dgu.eq.0.) dgu=dgdim(ip)
      if(ivsum.eq.0.and.dvu.eq.0.) dvu=dvdim(ip)
  390 continue
      idum=map1(dtdim(-msp),ap1(-msp,1),2*msp+1,dtu,at1u,1)
      idum=map1(dgdim(-msp),ap1(-msp,2),2*msp+1,dgu,ag1u,1)
      idum=map1(dvdim(-msp),ap1(-msp,3),2*msp+1,dvu,av1u,1)
      if(ap2(0,1).lt.0.) then
      at2u=0.
      ag2u=0.
      av2u=0.
      else
      idum=map1(dtdim(-msp),ap2(-msp,1),2*msp+1,dtu,at2u,1)
      idum=map1(dgdim(-msp),ap2(-msp,2),2*msp+1,dgu,ag2u,1)
      idum=map1(dvdim(-msp),ap2(-msp,3),2*msp+1,dvu,av2u,1)
      end if

      dtl=0.
      dgl=0.
      dvl=0.
      do 391 ip=-msp,-1,1
      itsum=ifoverc(ip,1)+ifoverw(ip,1)+ifovera(ip,1)
      igsum=ifoverc(ip,2)+ifoverw(ip,2)+ifovera(ip,2)
      ivsum=ifoverc(ip,3)+ifoverw(ip,3)+ifovera(ip,3)
      if(itsum.eq.0.and.dtl.eq.0.) dtl=dtdim(ip)
      if(igsum.eq.0.and.dgl.eq.0.) dgl=dgdim(ip)
      if(ivsum.eq.0.and.dvl.eq.0.) dvl=dvdim(ip)
  391 continue

      idum=map1(dtdim(-msp),ap1(-msp,1),2*msp+1,dtl,at1l,1)
      idum=map1(dgdim(-msp),ap1(-msp,2),2*msp+1,dgl,ag1l,1)
      idum=map1(dvdim(-msp),ap1(-msp,3),2*msp+1,dvl,av1l,1)
      if(ap2(0,1).lt.0.) then
      at2l=0.
      ag2l=0.
      av2l=0.
      else
      idum=map1(dtdim(-msp),ap2(-msp,1),2*msp+1,dtl,at2l,1)
      idum=map1(dgdim(-msp),ap2(-msp,2),2*msp+1,dgl,ag2l,1)
      idum=map1(dvdim(-msp),ap2(-msp,3),2*msp+1,dvl,av2l,1)
      end if

      errt=(dtu-dtl)*0.5
      errg=(dgu-dgl)*0.5
      errv=(dvu-dvl)*0.5

      ea1t=abs(at1l-at1u)*0.5
      ea1g=abs(ag1l-ag1u)*0.5
      ea1v=abs(av1l-av1u)*0.5
      erra1=sqrt(ea1t**2+ea1g**2+ea1v**2)

      ea2t=abs(at2l-at2u)*0.5
      ea2g=abs(ag2l-ag2u)*0.5
      ea2v=abs(av2l-av2u)*0.5
      erra2=sqrt(ea2t**2+ea2g**2+ea2v**2)


      open (unit=29,file=dfname(ns:ne)//metal(1:ml)//'ptb')
      write(29,354) xcur,met,c1,c2,c3,ftol
  354 format('# xcur:',f6.3,' met: ',a3,
     & '  c1:',0pf9.5,'  c2:',0pf9.5,'  c3:',0pf9.5,'  ftol=',1pe12.5)
      write(29,355) 
  355 format('# ','   Teff(f) ','  logg(f) ','   vt(f) ',11x,
     & ' A(Fe1) ','  sig1  ',' A(Fe2) ',' sig2 ',2x,' cmin ',' cmax ',
     & ' ewmin',' ewmax','  n1','  n2')
      write(29,360) t0,g0,v0,ap1(0,1),sp1(0,1),ap2(0,1),sp2(0,1),
     & cmin,cmax,wmin,wmax,n1,n2
  360 format('# ',f10.1,f10.3,f9.2,10x,4f8.3,2x,2f6.3,2f6.1,2i4)
      write(29,359) errt,errg,errv,erra1,erra2
  359 format('# ',' +/- ',f5.1,' +/- ',f5.3,' +/- ',f4.2,
     & 9x,'+/- ',f5.3,7x,'+/- ',f5.3)
      if(n2.eq.0) then
      write(29,357)
  357 format('#',' CAUTION ! ---> ^^^^^ <--- This error in log g is NOT'
     & ' reliable because Fe II lines are lacking !!!')
      end if 
      write(29,361) deltat,deltag,deltav
  361 format('# perturbation with steps of ',f4.1,'K, ',f5.3,' and ',
     & f5.3,'km/s')
      write(29,362)
  362 format('#',3x,
     & 24x,'<<<<tvar>>>>',24x,
     & 24x,'<<<<gvar>>>>',24x,
     & 24x,'<<<<vvar>>>>')
      write(29,365)
  365 format('# ip ',3(2x,'  del ','  A1  ','  A2  ','  sig1 ',
     & '  sig2 ','   ac   ','   bc   ','   aw   ','bw*1000 '))
      do 348 ip=-msp,msp
      del(1)=dtdim(ip)
      del(2)=dgdim(ip)
      del(3)=dvdim(ip)
      write(29,369) ip, (del(ih),ap1(ip,ih),ap2(ip,ih),sp1(ip,ih),
     & sp2(ip,ih),ac(ip,ih),bc(ip,ih),aw(ip,ih),bw(ip,ih)*1000.,ih=1,3)
  369 format(i4,2x,f6.1,2f6.3,2f7.4,4f8.4,2x,f6.3,2f6.3,0p2f7.4,4f8.4,
     & 2x,f6.3,2f6.3,2f7.4,4f8.4)
  348 continue
      write(29,393)
  393 format('#  ip',2x,
     & 2x,'  dt  ',' ic iw ia',
     & 2x,'  dg  ',' ic iw ia',
     & 2x,'  dv  ',' ic iw ia')
      do 397 ip=-msp,msp
      write(29,394) ip,
     & dtdim(ip),ifoverc(ip,1),ifoverw(ip,1),ifovera(ip,1),
     & dgdim(ip),ifoverc(ip,2),ifoverw(ip,2),ifovera(ip,2),
     & dvdim(ip),ifoverc(ip,3),ifoverw(ip,3),ifovera(ip,3)
  394  format('#',1x,i3,2x, 2x,f6.1,3i3, 2x,f6.3,3i3, 2x,f6.3,3i3)
  397 continue
      if(nrf.gt.0) then
        do i=1,4
          write(29,'(a)') labrsv(i)
        end do
      end if
      close (29)
      stop
      end if
c

      go to (1001,1002,1003),ndim
 1003 continue
c     3D case (ndim=3)
      p(1,1)=t1
      p(1,2)=g1
      p(1,3)=v1
      p(2,1)=t2
      p(2,2)=g2
      p(2,3)=v2
      p(3,1)=t1
      p(3,2)=g2
      p(3,3)=v1
      p(4,1)=t2
      p(4,2)=g1
      p(4,3)=v2
      do 123 i=1,4
      do 113 j=1,3
  113 x(j)=p(i,j)
      y(i)=fcn3(x)
  123 continue   
      call amoeba(p,y,mp,np,ndim,ftol,fcn3,iter)
      write(*,'(1x,a,i3)') 'Number of iterations: ',iter
      write(*,'(1x,a)') 'Vertices of final 3-D simplex and'
      write(*,'(1x,a)') 'function values at the vertices:'
      write(*,'(3x,a,t11,a,t23,a,t35,a,t45,a)') 'I',
     *     'X(I)','Y(I)','Z(I)','FUNCTION'
      do i=1,ndim+1
        write(*,'(1x,i3,0p3f12.5,1pe12.5)') i,(p(i,j),j=1,ndim),y(i)
      end do
      tf=p(ndim+1,1)
      gf=p(ndim+1,2)
      vf=p(ndim+1,3)
      go to 999
 1002 continue
c     2D case (ndim=2)
      if(nmod.eq.101) then
      	xx1=t1
      	xx2=t2
      	yy1=v1
      	yy2=v2
      end if
      if(nmod.eq.011) then
      	xx1=g1
      	xx2=g2
      	yy1=v1
      	yy2=v2
      end if
      if(nmod.eq.110) then
      	xx1=t1
      	xx2=t2
      	yy1=g1
      	yy2=g2
      end if
      p(1,1)=xx1
      p(1,2)=yy1
      p(2,1)=xx2
      p(2,2)=yy1
      p(3,1)=xx1
      p(3,2)=yy2
      do 122 i=1,3
      do 112 j=1,2
  112 x(j)=p(i,j)
      y(i)=fcn2(x)
  122 continue   
      call amoeba(p,y,mp,np,ndim,ftol,fcn2,iter)
      write(*,'(1x,a,i3)') 'Number of iterations: ',iter
      write(*,'(1x,a)') 'Vertices of final 2-D simplex and'
      write(*,'(1x,a)') 'function values at the vertices:'
      write(*,'(3x,a,t11,a,t23,a,t35,a)') 'I',
     *     'X(I)','Y(I)','FUNCTION'
      do i=1,ndim+1
        write(*,'(1x,i3,0p2f12.5,1pe12.5)') i,(p(i,j),j=1,ndim),y(i)
      end do
      if(nmod.eq.101) then
      	tf=p(ndim+1,1)
      	gf=g0
      	vf=p(ndim+1,2)
      end if
      if(nmod.eq.011) then
      	tf=t0
      	gf=p(ndim+1,1)
      	vf=p(ndim+1,2)
      end if
      if(nmod.eq.110) then
      	tf=p(ndim+1,1)
      	gf=p(ndim+1,2)
      	vf=v0
      end if
      go to 999
 1001 continue
c     1D case (ndim=1)
      if(nmod.eq.100) then
      	xx1=t1
      	xx2=t2
      end if
      if(nmod.eq.010) then
      	xx1=g1
      	xx2=g2
      end if
      if(nmod.eq.001) then
      	xx1=v1
      	xx2=v2
      end if
      p(1,1)=xx1
      p(2,1)=xx2
      do 121 i=1,2
      do 111 j=1,1
  111 x(j)=p(i,j)
      y(i)=fcn1(x)
  121 continue   
      call amoeba(p,y,mp,np,ndim,ftol,fcn1,iter)
      write(*,'(1x,a,i3)') 'Number of iterations: ',iter
      write(*,'(1x,a)') 'Vertices of final 1-D simplex and'
      write(*,'(1x,a)') 'function values at the vertices:'
      write(*,'(3x,a,t11,a,t23,a)') 'I',
     *     'X(I)','FUNCTION'
      do i=1,ndim+1
        write(*,'(1x,i3,0pf12.5,1pe12.5)') i,(p(i,j),j=1,ndim),y(i)
      end do
      if(nmod.eq.100) then
      	tf=p(ndim+1,1)
      	gf=g0
      	vf=v0
      end if
      if(nmod.eq.010) then
      	tf=t0
      	gf=p(ndim+1,1)
      	vf=v0
      end if
      if(nmod.eq.001) then
      	tf=t0
      	gf=g0
      	vf=p(ndim+1,1)
      end if
  999 continue
        delt=tf-t0
        delg=gf-g0
        delv=vf-v0
cc      	open (unit=9,file=dfname(ns:ne)//metal(1:ml)//'tmp')
cc        xcurnew = 0.5*(a1+a2)-7.5
cc      	write(9,'(f6.3,1x,a3,1x,6e15.7)') 
cc     &  xcurnew,met,tf,gf,vf,delt,delg,delv
cc        close (9)
cc
cc      revise xcur
        xcur=0.5*(a1+a2)-7.5
cc      replace (t0,g0,v0) by (tf,gf,vf)
        t0=tf
        g0=gf
        v0=vf
cc
        if(ifinit.eq.1) then
      	open (unit=19,file=dfname(ns:ne)//metal(1:ml)//'log')
        rewind 19
        else
      	open (unit=19,file=dfname(ns:ne)//metal(1:ml)//'log',
     &   access='append')
        end if
      	write(19,'(f6.3,1x,a3,1x,f8.1,f7.3,f6.2,1pe12.5,
     &   0p4f6.3,3i3,f8.1,f7.3,f6.2)') 
     %   xcur,met,tf,gf,vf,y(ndim+1),a1,s1,a2,s2,ift,ifg,ifv,
     %   delt,delg,delv
        close (19)
      open (unit=21,file=dfname(ns:ne)//metal(1:ml)//'abd')
      write(21,254) xcur,met,c1,c2,c3,ftol
  254 format('# xcur:',f6.3,' met: ',a3,
     & '  c1:',0pf9.5,'  c2:',0pf9.5,'  c3:',0pf9.5,'  ftol=',1pe12.5)
      write(21,255) 
  255 format('# ',' Teff(f)',' glog(f)',' vt(f)','    func    ',
     & 'A(Fe1)',' sig1 ','A(Fe2)',' sig2 ','  n1','  n2')
      write(21,260) tf,gf,vf,y(ndim+1),a1,s1,a2,s2,n1,n2
  260 format('# ',0pf8.1,0pf8.3,0pf6.2,1pe12.5,0p4f6.3,2i4)
      write(21,265)
  265 format('#No.','  code ','   wla   ','  chi  ',' gflog ',
     & '  ew  ',' epslog ',' epslog ')
      do 301 i=1,mf1
      if(ifuse(i).eq.0) go to 301
      write(21,270) num(i),code(i),wla(i),chi(i),gflog(i)-depsl(i),
     & ewobs(i),ablcal(i),-ablcal(i)
  270 format(i4,f6.2,f9.3,f7.3,f7.3,f6.1,f8.3,f8.3)
  301 continue
      do 302 i=mf1+1,mz
      if(ifuse(i).eq.0) go to 302
      write(21,270) num(i),code(i),wla(i),chi(i),gflog(i)-depsl(i),
     & ewobs(i),-ablcal(i),ablcal(i)
  302 continue
      if(nrf.gt.0) then
        do i=1,4
          write(21,'(a)') labrsv(i)
        end do
      end if
      close (21)

cc*iter
 8900 continue
cc*iter

      stop
      end 

      function fcn3(x)
      implicit real*4 (a-h,o-z)
      dimension x(3)
      common/tgvpass/t0,g0,v0,ift,ifg,ifv,ndim,nmod
      common/avsig/a1,s1,a2,s2
      equivalence (a1,ave1), (s1,sig1), (a2,ave2), (s2,sig2)
      common/params/c1,c2,c3,ftol
c     x(1)--> teff, x(2)--> glog, x(3)--> vt
      call avsigcal(x(1),x(2),x(3),ave1,sig1,ave2,sig2)
      if(ave2.lt.0.) then
      fcn3=sig1**2
      else
      fcn3=(sig1**2+c1*sig2**2)+c2*(ave1-ave2+c3)**2
      end if
      return
      end

      function fcn2(x)
      implicit real*4 (a-h,o-z)
      dimension x(2)
      common/tgvpass/t0,g0,v0,ift,ifg,ifv,ndim,nmod
      common/avsig/a1,s1,a2,s2
      equivalence (a1,ave1), (s1,sig1), (a2,ave2), (s2,sig2)
      common/params/c1,c2,c3,ftol
c     if nmod=110, then  x(1)--> teff and x(2)--> glog
c     if nmod=101, then  x(1)--> teff and x(2)--> vt
c     if nmod=011, then  x(1)--> glog and x(2)--> vt
      if (nmod.eq.110) call avsigcal(x(1),x(2),v0,ave1,sig1,ave2,sig2)
      if (nmod.eq.101) call avsigcal(x(1),g0,x(2),ave1,sig1,ave2,sig2)
      if (nmod.eq. 11) call avsigcal(t0,x(1),x(2),ave1,sig1,ave2,sig2)
      if(ave2.lt.0.) then
      fcn2=sig1**2
      else
      fcn2=(sig1**2+c1*sig2**2)+c2*(ave1-ave2+c3)**2
      end if
      return
      end

      function fcn1(x)
      implicit real*4 (a-h,o-z)
      dimension x(1)
      common/tgvpass/t0,g0,v0,ift,ifg,ifv,ndim,nmod
      common/avsig/a1,s1,a2,s2
      equivalence (a1,ave1), (s1,sig1), (a2,ave2), (s2,sig2)
      common/params/c1,c2,c3,ftol
c     if nmod=100, then  x(1)--> teff
c     if nmod= 10, then  x(1)--> glog
c     if nmod=  1, then  x(1)--> vt
      if (nmod.eq.100) call avsigcal(x(1),g0,v0,ave1,sig1,ave2,sig2)
      if (nmod.eq. 10) call avsigcal(t0,x(1),v0,ave1,sig1,ave2,sig2)
      if (nmod.eq.  1) call avsigcal(t0,g0,x(1),ave1,sig1,ave2,sig2)
      if(ave2.lt.0.) then
      fcn1=sig1**2
      else
      fcn1=(sig1**2+c1*sig2**2)+c2*(ave1-ave2+c3)**2
      end if
      return
      end
      
      subroutine coeff(ac,bc,aw,bw)
      implicit real*4 (a-h,o-z)
      parameter (MLN=330)
      common/ewdata/ifuse(MLN),ewobslg(MLN),ewobs(MLN)
      common/abdata/ablcal(MLN)
      common/linedata/num(MLN),code(MLN),wla(MLN),chi(MLN),gflog(MLN)
      common/numdata/mf1,mf2,mz,mt,mg,mv,ma
      real*8 wla,wlax
      dimension ifsel(MLN),x(MLN),y(MLN)
      do 10 i=1,mz
      ifsel(i)=ifuse(i)
      x(i)=chi(i)
      y(i)=ablcal(i)
   10 continue
      call linreg(x,y,ifsel,mz,ac,bc)
      do 20 i=1,mz
      ifsel(i)=ifuse(i)
      x(i)=ewobs(i)
      y(i)=ablcal(i)
   20 continue
      call linreg(x,y,ifsel,mz,aw,bw)
      return
      end
      
      subroutine linreg (x,y,ifsel,n,a,b)
      implicit real*4 (a-h,o-z)
      dimension x(n),y(n),ifsel(n)
      sx=0.
      sy=0.
      sxx=0.
      sxy=0.
      nc=0
      do j=1,n
        if(ifsel(j).eq.0) cycle
        nc=nc+1
        sx=sx+x(j)
        sy=sy+y(j)
        sxx=sxx+x(j)*x(j)
        sxy=sxy+x(j)*y(j)
      end do
      bumbo=sx**2-nc*sxx
      bunshia=sx*sxy-sxx*sy
      bunshib=sx*sy-nc*sxy
      a=bunshia/bumbo
      b=bunshib/bumbo
      return
      end

      subroutine avsigcal(teff,glog,vt,ave1,sig1,ave2,sig2)
      implicit real*4 (a-h,o-z)
      parameter (MLN=330)
      common/ewdata/ifuse(MLN),ewobslg(MLN),ewobs(MLN)
      common/abdata/ablcal(MLN)
      common/numdata/mf1,mf2,mz,mt,mg,mv,ma
      dimension nc(2),aveab(2),sigab(2),suma(2),sumd(2)
      dimension ks(2),ke(2)
      ks(1)=1
      ke(1)=mf1
      ks(2)=mf1+1
      ke(2)=mz
      do 5 i=1,mz
    5 ablcal(i)=0.
      do 10 m=1,2 
      nc(m)=0
      suma(m)=0.
      sumd(m)=0.
      do 20 i=ks(m),ke(m)
      if(ifuse(i).eq.0) go to 20
      nc(m)=nc(m)+1
      call abint(i,teff,glog,vt,ewobslg(i),ablcal(i))
      suma(m)=suma(m)+ablcal(i)
   20 continue
      if(nc(m).le.0) then
      aveab(m)=-9.99
      else
      aveab(m)=suma(m)/nc(m)
      end if
      do 30 i=ks(m),ke(m)
      if(ifuse(i).eq.0) go to 30
      sumd(m)=sumd(m)+(ablcal(i)-aveab(m))**2
   30 continue
      if(nc(m).le.0) then
      sigab(m)=-9.99
      else
      sigab(m)=sqrt(sumd(m)/nc(m))
      end if
   10 continue
      ave1=aveab(1)
      ave2=aveab(2)
      sig1=sigab(1)
      sig2=sigab(2)
      return
      end

      subroutine abint(iline,teff,glog,vt,ewlobs,abl)
      USE DFLIB
      use dfport
      implicit real*4 (a-h,o-z)
      parameter (MLN=330)
      parameter (MMT=9,MMG=9,MMV=12,MMA=18)
      common/linedata/num(MLN),code(MLN),wla(MLN),chi(MLN),gflog(MLN)
      common/numdata/mf1,mf2,mz,mt,mg,mv,ma
      real*8 wla,wlax
      dimension ewlog(MMT,MMG,MMV,MMA,MLN)
      dimension ewlog1(MMT,MMG,MMV,MMA,MLN),ewlog2(MMT,MMG,MMV,MMA,MLN)
      dimension gablog(MMA),gvt(MMV),gteff(MMT),gglog(MMG)
      dimension ewltemp(MMV),gewl(MMA)
      common/metallicity/xcur,met
      character*3 met
      character*80 label,moddir
      character*3 tem(2)
      dimension wt(2),ewdum(MMA)
      common/nrfpas/nrf
      common/refdat/numr(MLN),ewr(MLN),abr(MLN),w_r(9999),a_r(9999)    
      common/refpar/afer,xxt,xxg,xxv,xxf,xxa1,xxs1,xxa2,xxs2,nxx1,nxx2
      common/lc2num/lctonum(9999)
      common/mzxpass/mzx
      common/ldatax/numx(MLN),codex(MLN),wlax(MLN),chix(MLN),
     & gflogx(MLN),depsl(MLN)
      data nnxsv/999/

      if(iline.gt.0) go to 555

c-----beginning of initialize part

c     initialize section (defined met case)
c
      call getenv("MOD_DIR",moddir)
      call lablen(moddir,80,msdir,medir)
c
      if(met.eq.'xxx') go to 444

      inquire(file=moddir(msdir:medir)//'\ewmes'//met//'.asc',
     & exist=loexist)
      if(.not.loexist) then
      write(*,*) 'file does not exist! -->'//moddir(msdir:medir)//
     & '\ewdat'//met//'.asc'
      stop
      end if
      open (unit=88,
     & file=moddir(msdir:medir)//'\ewmes'//met//'.asc')
      read(88,*)
      read(88,*) mz,mf1,mf2
      if(mz.ne.MLN) then
      	write(*,*) 'mz,MLN= ',mz,MLN
      	if(mz.gt.MLN) stop
      end if
      if(mzx.gt.0) call ncheck(mz,mzx,'mz|mzx',0)
      do 100 i=1,mz
      read(88,*) 
     & num(i),code(i),wla(i),chi(i),gflog(i)
      if(mzx.gt.0) then
      call ncheck(num(i),numx(i),'num|numx',i)
      call check(code(i),codex(i),'code|codex',i)
      call check(wla(i),wlax(i),'wla|wlax',i)
      call check(chi(i),chix(i),'chi|chix',i)
      end if
  100 continue
      read(88,*) mt,mg,mv,ma
      read(88,*) (gteff(it),it=1,mt)
      read(88,*) (gglog(ig),ig=1,mg)
      read(88,*) (gvt(iv),iv=1,mv)
      read(88,*) (gablog(ia),ia=1,ma)
      do 200 il=1,mz
      do 200 it=1,mt
      do 200 ig=1,mg
      do 200 iv=1,mv
      read(88,*) (ewlog(it,ig,iv,ia,il),ia=1,ma)
  200 continue
      close (88)
      if(mzx.gt.0) then
      	do i=1,mz
      	  depsl(i)=gflog(i)-gflogx(i)
      	end do
      end if
      return

  444 continue

c     initialize section (undefined met case)
c
c     check if necessary
c
      if(xcur.lt.-1.5) then
      write(*,*) 'xcur=',xcur, ' is out of range!'
      stop
      else if (xcur.ge.-1.5.and.xcur.lt.-1.0) then
      nnx=-3
      tem(1)='m15'
      tem(2)='m10'
      wt(1)=(-1.0-xcur)/0.5
      wt(2)=1.-wt(1)
      else if (xcur.ge.-1.0.and.xcur.lt.-0.5) then
      nnx=-2
      tem(1)='m10'
      tem(2)='m05'
      wt(1)=(-0.5-xcur)/0.5
      wt(2)=1.-wt(1)
      else if (xcur.ge.-0.5.and.xcur.lt. 0.0) then
      nnx=-1
      tem(1)='m05'
      tem(2)='p00'
      wt(1)=( 0.0-xcur)/0.5
      wt(2)=1.-wt(1)
      else if (xcur.ge. 0.0.and.xcur.lt.+0.5) then
      nnx=0
      tem(1)='p00'
      tem(2)='p05'
      wt(1)=(+0.5-xcur)/0.5
      wt(2)=1.-wt(1)
      else if (xcur.ge.+0.5) then
      write(*,*) 'xcur=',xcur, ' is out of range!'
      stop
      end if
c
c     if unnecessary, skip to read "ewmes*.asc" files
      if(nnx.eq.nnxsv) go to 520
c
      do 499 k=1,2
c      if(k.eq.1.and.nnx.eq.nnxsv+1) then
c        do 401 il=1,mz
c        do 401 ia=1,ma
c        do 401 iv=1,mv
c        do 401 ig=1,mg
c        do 401 it=1,mt
c        ewlog1(it,ig,iv,ia,il)=ewlog2(it,ig,iv,ia,il)
c  401   continue
c        go to 499
c      end if
c      if(k.eq.2.and.nnx.eq.nnxsv-1) then
c        do 402 il=1,mz
c        do 402 ia=1,ma
c        do 402 iv=1,mv
c        do 402 ig=1,mg
c        do 402 it=1,mt
c        ewlog2(it,ig,iv,ia,il)=ewlog1(it,ig,iv,ia,il)
c  402   continue
c        go to 499
c      end if
c
c     read the ewmes*.asc file
c      
      inquire(file=moddir(msdir:medir)//'\ewmes'//tem(k)//'.asc',
     & exist=loexist)
      if(.not.loexist) then
      write(*,*) 'file does not exist! -->'//moddir(msdir:medir)//
     & '\ewdat'//tem(k)//'.asc'
      stop
      end if
      open (unit=88,
     & file=moddir(msdir:medir)//'\ewmes'//tem(k)//'.asc')
      read(88,*)
      read(88,*) mz,mf1,mf2
      if(mz.ne.MLN) then
      	write(*,*) 'mz,MLN= ',mz,MLN
      	if(mz.gt.MLN) stop
      end if
      if(mzx.gt.0) call ncheck(mz,mzx,'mz|mzx',0)
      do 104 i=1,mz
      read(88,*) 
     & num(i),code(i),wla(i),chi(i),gflog(i)
      if(mzx.gt.0) then
      call ncheck(num(i),numx(i),'num|numx',i)
      call check(code(i),codex(i),'code|codex',i)
      call check(wla(i),wlax(i),'wla|wlax',i)
      call check(chi(i),chix(i),'chi|chix',i)
      end if
  104 continue
      read(88,*) mt,mg,mv,ma
      read(88,*) (gteff(it),it=1,mt)
      read(88,*) (gglog(ig),ig=1,mg)
      read(88,*) (gvt(iv),iv=1,mv)
      read(88,*) (gablog(ia),ia=1,ma)
      go to (1201,1202),k
 1201 continue
      do 204 il=1,mz
      do 204 it=1,mt
      do 204 ig=1,mg
      do 204 iv=1,mv
      read(88,*) (ewlog1(it,ig,iv,ia,il),ia=1,ma)
  204 continue
      go to 1209
 1202 continue
      do 205 il=1,mz
      do 205 it=1,mt
      do 205 ig=1,mg
      do 205 iv=1,mv
      read(88,*) (ewlog2(it,ig,iv,ia,il),ia=1,ma)
  205 continue
 1209 continue
      close (88)

      if(mzx.gt.0) then
      	do i=1,mz
      	  depsl(i)=gflog(i)-gflogx(i)
      	end do
      end if

  499 continue

  520 continue

c     setting-up ewlog by interpolation for the current xcur  
      do 9204 il=1,mz
      do 9204 ia=1,ma
      do 9204 iv=1,mv
      do 9204 ig=1,mg
      do 9204 it=1,mt
      ewlog(it,ig,iv,ia,il)=wt(1)*ewlog1(it,ig,iv,ia,il)+
     %                      wt(2)*ewlog2(it,ig,iv,ia,il)
 9204 continue
c
c     reset nnxsv parameter
      nnxsv=nnx
c
      return
c----- end of initialize section

c
c----- Now computing abl for given ewlobs
  555 continue
      if(iline.lt.1.or.iline.gt.mz) then
      write(*,*) 'Invalid iline=',iline
      stop
      end if
      do 120 ia=1,ma
      do 123 iv=1,mv
      call zint(gteff,gglog,mt,mg,teff,glog,ewlog(1,1,iv,ia,iline),
     & ewltemp(iv))      
  123 continue
      idum=map1(gvt,ewltemp,mv,vt,gewl(ia),1)
  120 continue
      idum=map1(gewl,gablog,ma,ewlobs,abl,1)
c
c     abundance correction according to the gf difference
c     (for the case of mzx > 0)
c
      if(mzx.gt.0) then
        abl=abl+depsl(iline)
      end if
c
c     calculate differential abundance in case of ref mode
      if(nrf.gt.0) then
        if(a_r(num(iline)).lt.0.) then
          write(*,*) 'ref set incompatible! num(iline)=',num(iline)
          stop
        end if
c
        abl=(abl-a_r(num(iline)))+afer
c        abl=abl-a_r(num(iline))
c
      end if
      return
      end

      subroutine zint(gx,gy,nx,ny,x,y,gz,z)
      implicit real*4 (a-h,o-z)
      dimension gx(nx),gy(ny),gz(nx,ny)
c     x ordinate
      if(x.lt.gx(1)) then
      i1=1
      go to 10
      end if
      if(x.ge.gx(nx-1)) then
      i1=nx-1
      go to 10
      end if
      do i=1,nx-1
      if(x.ge.gx(i).and.x.lt.gx(i+1)) then
      i1=i
      go to 10
      end if
      end do
   10 i2=i1+1
c     y ordinate
      if(y.lt.gy(1)) then
      j1=1
      go to 20
      end if
      if(y.ge.gy(ny-1)) then
      j1=ny-1
      go to 20
      end if
      do j=1,ny-1
      if(y.ge.gy(j).and.y.lt.gy(j+1)) then
      j1=j
      go to 20
      end if
      end do
   20 j2=j1+1
      dx=gx(i2)-gx(i1)
      dy=gy(j2)-gy(j1)
      p=(x-gx(i1))/dx
      q=(y-gy(j1))/dy
      z0=+gz(i1,j1)
      z1=-gz(i1,j1)+gz(i2,j1)
      z2=-gz(i1,j1)+gz(i1,j2)
      z3=+gz(i1,j1)-gz(i1,j2)-gz(i2,j1)+gz(i2,j2)
      z=z0+z1*p+z2*q+z3*p*q
      return
      end

      FUNCTION MAP1(XOLD,FOLD,NOLD,XNEW,FNEW,NNEW)
      implicit real*4 (a-h,o-z)
      DIMENSION XOLD(1),FOLD(1),XNEW(1),FNEW(1)
      L=2
      LL=0
      DO 50 K=1,NNEW
   10 IF(XNEW(K).LT.XOLD(L))GO TO 20
      L=L+1
      IF(L.GT.NOLD)GO TO 30
      GO TO 10
   20 IF(L.EQ.LL)GO TO 50
      IF(L.EQ.2)GO TO 30
      IF(L.EQ.3)GO TO 30
      L1=L-1
      IF(L.GT.LL+1.OR.L.EQ.3)GO TO 21
      IF(L.GT.LL+1.OR.L.EQ.4)GO TO 21
      CBAC=CFOR
      BBAC=BFOR
      ABAC=AFOR
      IF(L.EQ.NOLD)GO TO 22
      GO TO 25
   21 L2=L-2
      D=(FOLD(L1)-FOLD(L2))/(XOLD(L1)-XOLD(L2))
      CBAC=FOLD(L)/((XOLD(L)-XOLD(L1))*(XOLD(L)-XOLD(L2)))+
     1(FOLD(L2)/(XOLD(L)-XOLD(L2))-FOLD(L1)/(XOLD(L)-XOLD(L1)))/
     2(XOLD(L1)-XOLD(L2))
      BBAC=D-(XOLD(L1)+XOLD(L2))*CBAC
      ABAC=FOLD(L2)-XOLD(L2)*D+XOLD(L1)*XOLD(L2)*CBAC
      IF(L.LT.NOLD)GO TO 25
   22 C=CBAC
      B=BBAC
      A=ABAC
      LL=L
      GO TO 50
   25 D=(FOLD(L)-FOLD(L1))/(XOLD(L)-XOLD(L1))
      CFOR=FOLD(L+1)/((XOLD(L+1)-XOLD(L))*(XOLD(L+1)-XOLD(L1)))+
     1(FOLD(L1)/(XOLD(L+1)-XOLD(L1))-FOLD(L)/(XOLD(L+1)-XOLD(L)))/
     2(XOLD(L)-XOLD(L1))
      BFOR=D-(XOLD(L)+XOLD(L1))*CFOR
      AFOR=FOLD(L1)-XOLD(L1)*D+XOLD(L)*XOLD(L1)*CFOR
      WT=0.
      IF(ABS(CFOR).NE.0.)WT=ABS(CFOR)/(ABS(CFOR)+ABS(CBAC))
      A=AFOR+WT*(ABAC-AFOR)
      B=BFOR+WT*(BBAC-BFOR)
      C=CFOR+WT*(CBAC-CFOR)
      LL=L
      GO TO 50
   30 IF(L.EQ.LL)GO TO 50
      L=MIN0(NOLD,L)
      C=0.
      B=(FOLD(L)-FOLD(L-1))/(XOLD(L)-XOLD(L-1))
      A=FOLD(L)-XOLD(L)*B
      LL=L
   50 FNEW(K)=A+(B+C*XNEW(K))*XNEW(K)
      MAP1=LL-1
      RETURN
      END
      
      SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter)
      implicit real*4 (a-h,o-z)
      INTEGER iter,mp,ndim,np,NMAX,ITMAX
      REAL ftol,p(mp,np),y(mp),funk
      PARAMETER (NMAX=20,ITMAX=5000)
      EXTERNAL funk
CU    USES amotry,funk
      INTEGER i,ihi,ilo,inhi,j,m,n
      REAL rtol,sum,swap,ysave,ytry,psum(NMAX),amotry
      iter=0
1     do 12 n=1,ndim
        sum=0.
        do 11 m=1,ndim+1
          sum=sum+p(m,n)
11      continue
        psum(n)=sum
12    continue
2     ilo=1
      if (y(1).gt.y(2)) then
        ihi=1
        inhi=2
      else
        ihi=2

        inhi=1
      endif
      do 13 i=1,ndim+1
        if(y(i).le.y(ilo)) ilo=i
        if(y(i).gt.y(ihi)) then
          inhi=ihi
          ihi=i
        else if(y(i).gt.y(inhi)) then
          if(i.ne.ihi) inhi=i
        endif
13    continue
      rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
      if (rtol.lt.ftol) then
        swap=y(1)
        y(1)=y(ilo)
        y(ilo)=swap
        do 14 n=1,ndim
          swap=p(1,n)
          p(1,n)=p(ilo,n)
          p(ilo,n)=swap
14      continue

        return
      endif
      if (iter.ge.ITMAX) pause 'ITMAX exceeded in amoeba'
      iter=iter+2
      ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.0)
      if (ytry.le.y(ilo)) then
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.0)
      else if (ytry.ge.y(inhi)) then
        ysave=y(ihi)
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.5)
        if (ytry.ge.ysave) then
          do 16 i=1,ndim+1
            if(i.ne.ilo)then
              do 15 j=1,ndim
                psum(j)=0.5*(p(i,j)+p(ilo,j))
                p(i,j)=psum(j)
15            continue
              y(i)=funk(psum)
            endif
16        continue
          iter=iter+ndim
          goto 1

        endif
      else
        iter=iter-1
      endif
      goto 2
      END

      FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac)
      implicit real*4 (a-h,o-z)
      INTEGER ihi,mp,ndim,np,NMAX
      REAL amotry,fac,p(mp,np),psum(np),y(mp),funk
      PARAMETER (NMAX=20)
      EXTERNAL funk
CU    USES funk
      INTEGER j
      REAL fac1,fac2,ytry,ptry(NMAX)
      fac1=(1.-fac)/ndim
      fac2=fac1-fac
      do 11 j=1,ndim
        ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
11    continue
      ytry=funk(ptry)
      if (ytry.lt.y(ihi)) then
        y(ihi)=ytry
        do 12 j=1,ndim
          psum(j)=psum(j)-p(ihi,j)+ptry(j)
          p(ihi,j)=ptry(j)

12      continue
      endif
      amotry=ytry
      return
      END

      subroutine lablen(label,nlen,ns,ne)
      implicit real*4 (a-h,o-z)
      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

      function ifnum(ch)
      implicit real*4 (a-h,o-z)
      character*1 ch,num(10)
      data num/'1','2','3','4','5','6','7','8','9','0'/
      do 100 i=1,10
      if(ch.eq.num(i)) go to 99
  100 continue
      ifnum=0
      return
   99 ifnum=1
      return
      end

      subroutine ncheck(n1,n2,str,ndum)
      implicit real*4 (a-h,o-z)
      character*8 str
      if(n1.ne.n2) then
      write(*,'(a16,2i6,a8,i6)') 'ncheck disorder!',n1,n2,str,ndum
      end if
      return
      end

      subroutine check(x1,x2,str,xdum)
      implicit real*4 (a-h,o-z)
      character*8 str
      if(abs(x1-x2).gt.1.e-6) then
      write(*,'(a16,2e12.5,a8,e12.5)') ' check disorder!',x1,x2,str,xdum
      end if
      return
      end

