      use dflib
      use dfport
      parameter (MF1=65+235+2,MF2=13+15,MLN=MF1+MF2)
      real*8 wl,wla(MLN)
      dimension num(MLN),nc(MLN),ewsv(10,MLN),ew(MLN)
      character*80 linesv(10,MLN)
      character*3 infosv(10,MLN)
      character reg*1,ord*2
      common/pass/n1,n2,chid1(MF1),ewd1(MF1),abld1(MF1),
     & chid2(MF2),ewd2(MF2),abld2(MF2),mode1(MF1),mode2(MF2),
     & numb1(MF1),numb2(MF2)

      common/refdat/nrf,numr(MLN),ewr(MLN),abr(MLN),w_r(9999),a_r(9999)    

      logical lorex,loexist,loinit
      character*80 line,buf,line1(100),line2(100)
      common/labpas/ns,ne,ks,ke,head,refs,ffnam
      character*40 head,refs,ffnam
      character*10 labewcrt
      character*1 sgn,key,ch
      integer*2 narg1,narg2
      character*80 moddir,fnparams

c
c     read num and wla from params.dat
c
      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,*) c1dum
      read(2,*) c2dum
      read(2,*) c3dum
      read(2,*) ftoldum
      read(2,*,err=399,end=399) mzx
      do i=1,mzx
      	read(2,*) num(i),code,wla(i)
      end do
      go to 398
  399 mzx=0
  398 continue
      if(mzx.ne.MLN) then
      	write(*,*) 'mzx,MLN ',mzx,MLN
      	stop
      end if
      close (2)
c
      narg1=1
      call getarg(narg1,buf)
      read(buf,'(a)') ffnam
      narg2=2
      call getarg(narg2,buf)
      read(buf,'(a)',err=14,end=14) sgn
      if(sgn.ne.'i'.and.sgn.ne.'c'.and.sgn.ne.'w') go to 14
      go to 15
   14 write(*,*) 'Needed 2nd argument should be any of (i, c, w) !'
      stop
   15 continue

      call lablen(ffnam,40,nffs,nffe)
      nbar=index(ffnam,'/')
      if(nbar.gt.0) then
      read(ffnam(1:nbar-1),'(a)') head
      call lablen(head,40,ns,ne)
      read(ffnam(nbar+1:40),'(a)') refs
      call lablen(refs,40,ks,ke)
      else
      read(ffnam,'(a)') head
      call lablen(head,40,ns,ne)
      ks=0
      ke=0
      end if

      if(nbar.gt.0) then
c     read (refs).abd file usedv for comparison
c     check if the file exists
      inquire(file=refs(ks:ke)//'.abd',exist=lorex)
      if(.not.lorex) then
        write(*,*) 'R: why not exist ? --> '//refs(ks:ke)//'.abd'
        stop
      end if
      do ir=1,9999
        w_r(ir)=-999.9
        a_r(ir)=-999.9
      end do
      open (unit=9,file=refs(ks:ke)//'.abd')      
      nr=0
      do 791 n=1,999
      read(9,'(a)',end=799) line
      if(line(1:1).eq.'#') go to 791
      read(line,*)
     & 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
      close (9)
      end if

      if(sgn.eq.'i') loinit=.true.
      if(sgn.eq.'c'.or.sgn.eq.'w') loinit=.false.
      if(loinit) then
        go to 1000
      else
        go to 2000
      end if
c     initialize mode
 1000 continue
      open (unit=9,file=head(ns:ne)//'.wdt')
      do k=1,MLN
      	nc(k)=0
      end do
      do 99 n=1,99999
      read(9,'(a)',end=999) line
      if(line(2:6).ne.'26.00'.and.line(2:6).ne.'26.01') go to 99
      read(line,'(6x,f9.3,30x,f7.1,7x,a1,9x,a2)') wl,ewtmp,reg,ord
      do 20 k=1,MLN
      if(abs(wl-wla(k)).lt.0.02) go to 29
   20 continue
      go to 99
   29 km=k
      nc(km)=nc(km)+1
      ewsv(nc(km),km)=ewtmp
      linesv(nc(km),km)=line
      infosv(nc(km),km)=reg//ord
   99 continue
  999 continue

c     calculate average EW for use
      do 130 k=1,MLN
      ew(k)=0.
      if(nc(k).le.0) go to 130
      sum=0.
      nn=0
      do j=1,nc(k)
      if(ewsv(j,k).le.0.) cycle
      nn=nn+1
      sum=sum+ewsv(j,k)
      end do
      ew(k)=sum/nn
  130 continue
  605 write(*,*) 
     & 'ewcrt? (in mA) (just [Enter] --> infinity, no screening)'
      read(*,'(a)') labewcrt
      read(labewcrt,*,err=609,end=609) ewcrt
      if(ewcrt.le.0) ewcrt=+1.e30
      go to 608
  609 ewcrt=1.e30
  608 continue
      if(ewcrt.gt.1.e29) then
      write(*,*) 'EW threshould is set to infinity.' 
      else
      write(*,139) ewcrt
  139 format('OK, lines stronger than ',f8.1,' mA are not used.')
      end if
      write(*,*) 'Following lines will be rejected...'
      ndl=0
      do 331 k=1,MLN
      if(nc(k).le.0) go to 331
      if(ew(k).ge.ewcrt) then
        ndl=ndl+1
        key = '>'
        tew=-abs(ew(k))
        write(*,350) 'exceed  ',num(k),tew,linesv(1,k)(1:21),
     &  linesv(1,k)(22:27),linesv(1,k)(53:59),key,nc(k),
     &  (ewsv(j,k),j=1,nc(k))
        go to 331
      end if
      if((lorex).and.(w_r(num(k)).lt.0.)) then
        ndl=ndl+1
        key = '#'
        tew=-abs(ew(k))
        write(*,350) 'no refs ',num(k),tew,linesv(1,k)(1:21),
     &  linesv(1,k)(22:27),linesv(1,k)(53:59),key,nc(k),
     &  (ewsv(j,k),j=1,nc(k))
        go to 331
      end if
  350 format(a8,i4,f7.1,1x,a21,1x,a6,1x,a7,2x,a1,i4,1x,5f6.1)
  331 continue
      if (ndl.eq.0) write(*,*) '(In this case, no lines deleted.)'
      if(nrf.gt.0) then
      write(*,*)'May I update '//head(ns:ne)//'~'//refs(ks:ke)//'.obs ?'
      else
      write(*,*)'May I update '//head(ns:ne)//'.obs ?'
      end if
  619 write(*,
     & '(''go->"g", retry other ewcrt ->"r", stop-> "x"'')' ) 
      ch=getcharqq()
      if(ch.ne.'g'.and.ch.ne.'r'.and.ch.ne.'x') go to 619
      if(ch.eq.'x') then
        write(*,*) 'Aborted and stopted!'
        stop
      end if
      if(ch.eq.'r') go to 605
      if(nrf.gt.0) then
        open (unit=10,file=head(ns:ne)//'~'//refs(ks:ke)//'.obs')
      else
        open (unit=10,file=head(ns:ne)//'.obs')
      end if
      open (unit=13,file=head(ns:ne)//'.del')
      do 30 k=1,MLN
      if(nc(k).le.0) go to 30
      if(ew(k).ge.ewcrt) then
        key = '>'
        tew=-abs(ew(k))
        write(13,50) num(k),tew,linesv(1,k)(1:21),
     &  linesv(1,k)(22:27),linesv(1,k)(53:59),key,nc(k),
     &  (ewsv(j,k),infosv(j,k),j=1,nc(k))
        go to 30
      end if
      if((lorex).and.(w_r(num(k)).lt.0.)) then
        key = '#'
        tew=-abs(ew(k))
        write(13,50) num(k),tew,linesv(1,k)(1:21),
     &  linesv(1,k)(22:27),linesv(1,k)(53:59),key,nc(k),
     &  (ewsv(j,k),infosv(j,k),j=1,nc(k))
        go to 30
      end if

c      write(*,*)
c      do j=1,nc(k)
c      write(*,50) num(k),ewsv(j,k),linesv(j,k)(1:21),
c     & linesv(j,k)(22:27),linesv(j,k)(53:59),key
c      end do
c      if(nc(k).gt.1) write(*,60) tew
c   60 format(3x,'<',f7.1,'>')
  439 continue
      key=' '
      write(10,50) num(k),ew(k),linesv(1,k)(1:21),
     & linesv(1,k)(22:27),linesv(1,k)(53:59),key,nc(k),
     & (ewsv(j,k),infosv(j,k),j=1,nc(k))
   50 format(i4,f7.1,1x,a21,1x,a6,1x,a7,2x,a1,i4,1x,6(f6.1,1x,a3))
   30 continue
      close (10)
      close (13)
      write(*,*) 'file update finished !' 
      go to 9999

c     read current *.abd file if exists
 2000 continue
      inquire(file=head(ns:ne)//'.abd',exist=loexist)
      if(.not.loexist) then
        write(*,*) 'H: why not exist ? --> '//head(ns:ne)//'.abd'
        stop
      end if
      open (unit=11,file=head(ns:ne)//'.abd')      
      do n=1,MF1
        mode1(n)=1
      end do
      do n=1,MF2
        mode2(n)=2
      end do
C
      n1=0
      n2=0
      chimin=+1.e30
      chimax=-1.e30
      ewmin=+1.e30
      ewmax=-1.e30
      abmin=+1.e30
      abmax=-1.e30 
      do 91 n=1,999
      read(11,'(a)',end=9199) line
      if(line(1:1).eq.'#') go to 91
      read(line,*) 
     & num_,code_,wla_,chi_,gflog_,ew_,abl1_,abl2_
      if(chi_.gt.chimax) chimax=chi_
      if(chi_.lt.chimin) chimin=chi_
      if(ew_.gt.ewmax) ewmax=ew_
      if(ew_.lt.ewmin) ewmin=ew_
      if(code_.eq.26.00) then
      n1=n1+1
      line1(n1)=line
      numb1(n1)=num_
      chid1(n1)=chi_
      ewd1(n1)=ew_
      abld1(n1)=abl1_
      if(abl1_.gt.abmax) abmax=abl1_
      if(abl1_.lt.abmin) abmin=abl1_
      end if
      if(code_.eq.26.01) then
      n2=n2+1
      line2(n2)=line
      numb2(n2)=num_
      chid2(n2)=chi_
      ewd2(n2)=ew_
      abld2(n2)=abl2_
      if(abl2_.gt.abmax) abmax=abl2_
      if(abl2_.lt.abmin) abmin=abl2_
      end if
   91 continue
 9199 continue
      close (11)
c      write(*,*) 'ew range (raw): ',ewmin,ewmax
c      write(*,*) 'ab range (raw): ',abmin,abmax
      chimin=0.
      chimax=5.      
      ewmin=int(ewmin/50.)*50
      ewmax=(int(ewmax/50.)+1)*50
      abmin=int(abmin/0.5)*0.5
      abmax=(int(abmax/0.5)+1)*0.5
c      write(*,*) 'ew range (ajt): ',ewmin,ewmax
c      write(*,*) 'ab range (ajt): ',abmin,abmax
  566 continue
      if(sgn.eq.'w') call dataplot_w(ewmin,ewmax,abmin,abmax)
      if(sgn.eq.'c') call dataplot_c(chimin,chimax,abmin,abmax)
      write(*,*) 'You have escaped from the graphics mode...'
      write(*,*) 'List of manually rejected lines are as follows:'
      ndel=0
      do n=1,n1
        tew=abs(ewd1(n))
        if(mode1(n).le.0) then
          ndel=ndel+1
          tew=-abs(ewd1(n))
          write(*,505) numb1(n),tew,line1(n)(5:33)
          cycle
        end if
      end do
      do n=1,n2
        tew=abs(ewd2(n))
        if(mode2(n).le.0) then
          ndel=ndel+1
          tew=-abs(ewd2(n))
          write(*,505) numb2(n),tew,line2(n)(5:33)
          cycle
        end if
      end do
      write(*,
     & '(''('',i2,'' lines will be rejected)'')') ndel
      write(*,*) 'May I update '//head(ns:ne)//'.obs ?'
  519 write(*,
     & '('' GO->"g", redo->"r",stop -> "x"'')' ) 
      ch=getcharqq()
      if(ch.ne.'g'.and.ch.ne.'r'.and.ch.ne.'x') go to 519
      if(ch.eq.'x') then
        write(*,*) 'Aborted and stopted!'
        stop
      end if
      if(ch.eq.'r') then
      write(*,*) 'OK, then try again in the graphics window...'
      go to 566
      end if

      open (unit=10,file=head(ns:ne)//'.obs')
      open (unit=13,file=head(ns:ne)//'.del',access='append')
      ndel=0
      do n=1,n1
        tew=abs(ewd1(n))
        if(mode1(n).le.0) then
          ndel=ndel+1
          tew=-abs(ewd1(n))
          write(13,505) numb1(n),tew,line1(n)(5:33)
          cycle
        end if
        write(10,505) numb1(n),tew,line1(n)(5:33)
  505   format(i4,f7.1,1x,a29)
      end do
      do n=1,n2
        tew=abs(ewd2(n))
        if(mode2(n).le.0) then
          ndel=ndel+1
          tew=-abs(ewd2(n))
          write(13,505) numb2(n),tew,line2(n)(5:33)
          cycle
        end if
        write(10,505) numb2(n),tew,line2(n)(5:33)
      end do
      close (10)
      close (13)
      write(*,*) 'file update finished !' 
 9999 continue
      stop
      end

      function ifmode(num)
      parameter (MF1=65+235+2,MF2=13+15,MLN=MF1+MF2)
      common/pass/n1,n2,chid1(MF1),ewd1(MF1),abld1(MF1),
     & chid2(MF2),ewd2(MF2),abld2(MF2),mode1(MF1),mode2(MF2),
     & numb1(MF1),numb2(MF2)

      if(num.gt.200) go to 200     
  100 continue
c     Fe1 lines
      do n=1,n1
      if(num.eq.numb1(n)) go to 110
      end do
      write(*,*) 'Fe1--> Not found! num=',num
      write(*,*) (numb1(n),n=1,n1)
      stop
  110 continue
      ifmode=mode1(n)
      go to 900
  200 continue
c     Fe2 lines
      do n=1,n2
      if(num.eq.numb2(n)) go to 210
      end do
      write(*,*) 'Fe2--> Not found! num=',num
      write(*,*) (numb2(n),n=1,n2)
      stop
  210 continue
      ifmode=mode2(n)
  900 return
      end

      subroutine dataplot_w(ewmin,ewmax,abmin,abmax)
      use dflib
      use dfport
      parameter (MF1=65+235+2,MF2=13+15,MLN=MF1+MF2)
      common/labpas/ns,ne,ks,ke,head,refs,ffnam
      character*40 head,refs,ffnam
      common/pass/n1,n2,chid1(MF1),ewd1(MF1),abld1(MF1),
     & chid2(MF2),ewd2(MF2),abld2(MF2),mode1(MF1),mode2(MF2),
     & numb1(MF1),numb2(MF2)

      dimension x(2),yav(2),ysp(2),ysm(2),icod1(100),icod2(100)
      data ion1,ioff1/4,5/
      data ion2,ioff2/6,2/
      integer*2 ir
      character*1 ch
      character*4 cdev
c      data cdev/'/CGW'/
      data cdev/'/GW '/
      data itemp/0/
      if(itemp.gt.0) go to 700
      itemp=itemp+1
      ier=pgbeg(0,cdev,1,1)
c      ier=pgbeg(0,'?',1,1)
      call pgenv(ewmin,ewmax,abmin,abmax,0,-2)
      if(cdev.eq.'/GW ') then
      call pgscr(0,1.,1.,1.)
      call pgscr(1,0.,0.,0.)
      end if
      if(cdev.eq.'/CGW') then
      call pgscr(1,1.,1.,1.)
      call pgscr(0,0.,0.,0.)
      end if
      call pgsci(1)
      call pgbox('ABCTSN',0., 0, 'ABCTSN',0.,0)
      call pglabel('\frEquivalent width (mA)','\frFe Anundance', 
     & '\frA vs. W relation: '//head(ns:ne)) 
  700 continue
      loout=.false.
      do while (.not.loout)
      call avsig(av,sig,t_crt)
      call pgsci(1)
c
      x(1)=ewmin
      x(2)=ewmax
      yav(1)=av
      yav(2)=av
      call pgsls(1)
      call pgline(2,x,yav)
c
      ysp(1)=av + 2.*sig
      ysp(2)=av + 2.*sig
      call pgsls(4)
      call pgline(2,x,ysp)
c
      ysm(1)=av - 2.*sig
      ysm(2)=av - 2.*sig
      call pgsls(4)
      call pgline(2,x,ysm)
c
      do 11 n=1,n1
      if(mode1(n).gt.0) then
      icod1(n)=ion1
      else
      icod1(n)=ioff1
      end if
      CALL PGSCI(4)
      call pgpt1(ewd1(n),abld1(n),icod1(n))
   11 continue
      do 12 n=1,n2
      if(mode2(n).gt.0) then
      icod2(n)=ion2
      else
      icod2(n)=ioff2
      end if
      CALL PGSCI(3)
      call pgpt1(ewd2(n),abld2(n),icod2(n))
   12 continue
      ir=pgcurs(xm,ym,ch)
c      write(*,*) 'ir,xm,ym:',ir,xm,ym
c      write(*,'(''   ch: '',(a1))') ch
      if((xm.lt.ewmin.or.xm.gt.ewmax).or.
     % (ym.lt.abmin.or.ym.gt.abmax)) then
      loout=.true.
      if(ch.eq.'D') stop
      end if
      if(ichar(ch).eq.0) loout=.true.
      do 13 n=1,n1
      if(abs(xm-ewd1(n)).lt.2..and.abs(ym-abld1(n)).lt.0.02) then
      call pgsci(0)
      call pgpt1(ewd1(n),abld1(n),icod1(n))
c      call pgsci(1)      
      CALL PGSCI(4)
      if(ch.eq.'A'.and.mode1(n).gt.0) mode1(n)=-1
      if(ch.eq.'D'.and.mode1(n).lt.0) mode1(n)=+1
      end if
   13 continue
      do 14 n=1,n2
      if(abs(xm-ewd2(n)).lt.2..and.abs(ym-abld2(n)).lt.0.02) then
      call pgsci(0)
      call pgpt1(ewd2(n),abld2(n),icod2(n))
c      call pgsci(1)      
      CALL PGSCI(3)
      if(ch.eq.'A'.and.mode2(n).gt.0) mode2(n)=-2
      if(ch.eq.'D'.and.mode2(n).lt.0) mode2(n)=+2
      end if
   14 continue
      if(.not.loout) then
      call pgsci(0)
      call pgsls(1)
      call pgline(2,x,yav)
      call pgsls(4)
      call pgline(2,x,ysp)
      call pgsls(4)
      call pgline(2,x,ysm)
      call pgsci(1)
      end if
      end do
c  950 call gwquit(ir)
      return
      end

      subroutine dataplot_c(chimin,chimax,abmin,abmax)
      use dflib
      use dfport
      parameter (MF1=65+235+2,MF2=13+15,MLN=MF1+MF2)
      common/labpas/ns,ne,ks,ke,head,refs,ffnam
      character*40 head,refs,ffnam
      common/pass/n1,n2,chid1(MF1),ewd1(MF1),abld1(MF1),
     & chid2(MF2),ewd2(MF2),abld2(MF2),mode1(MF1),mode2(MF2),
     & numb1(MF1),numb2(MF2)

      dimension x(2),yav(2),ysp(2),ysm(2),icod1(100),icod2(100)
      data ion1,ioff1/4,5/
      data ion2,ioff2/6,2/
      integer*2 ir
      character*1 ch
      character*4 cdev
c      data cdev/'/CGW'/
      data cdev/'/GW '/
      data itemp/0/
      if(itemp.gt.0) go to 700
      itemp=itemp+1
      ier=pgbeg(0,cdev,1,1)
c      ier=pgbeg(0,'?',1,1)
      call pgenv(chimin,chimax,abmin,abmax,0,-2)
      if(cdev.eq.'/GW ') then
      call pgscr(0,1.,1.,1.)
      call pgscr(1,0.,0.,0.)
      end if
      if(cdev.eq.'/CGW') then
      call pgscr(1,1.,1.,1.)
      call pgscr(0,0.,0.,0.)
      end if
      call pgsci(1)
      call pgbox('ABCTSN',0., 0, 'ABCTSN',0.,0)
      call pglabel('\frExcitation potential (eV)','\frFe Anundance',
     & '\frA vs. chi relation: '//head(ns:ne)) 
  700 continue
      loout=.false.
      do while (.not.loout)
      call avsig(av,sig,t_crt)
      call pgsci(1)
c
      x(1)=chimin
      x(2)=chimax
      yav(1)=av
      yav(2)=av
      call pgsls(1)
      call pgline(2,x,yav)
c
      ysp(1)=av + 2.*sig
      ysp(2)=av + 2.*sig
      call pgsls(4)
      call pgline(2,x,ysp)
c
      ysm(1)=av - 2.*sig
      ysm(2)=av - 2.*sig
      call pgsls(4)
      call pgline(2,x,ysm)
c
      do 11 n=1,n1
      if(mode1(n).gt.0) then
      icod1(n)=ion1
      else
      icod1(n)=ioff1
      end if
      CALL PGSCI(4)
      call pgpt1(chid1(n),abld1(n),icod1(n))
   11 continue
      do 12 n=1,n2
      if(mode2(n).gt.0) then
      icod2(n)=ion2
      else
      icod2(n)=ioff2
      end if
      CALL PGSCI(3)
      call pgpt1(chid2(n),abld2(n),icod2(n))
   12 continue
      ir=pgcurs(xm,ym,ch)
c      write(*,*) 'ir,xm,ym:',ir,xm,ym
c      write(*,'(''   ch: '',(a1))') ch
      if((xm.lt.chimin.or.xm.gt.chimax).or.
     % (ym.lt.abmin.or.ym.gt.abmax)) then 
      loout=.true.
      if(ch.eq.'D') stop
      end if
      if(ichar(ch).eq.0) loout=.true.
      do 13 n=1,n1
      if(abs(xm-chid1(n)).lt.2..and.abs(ym-abld1(n)).lt.0.02) then
      call pgsci(0)
      call pgpt1(chid1(n),abld1(n),icod1(n))
c      call pgsci(1)      
      CALL PGSCI(4)
      if(ch.eq.'A'.and.mode1(n).gt.0) mode1(n)=-1
      if(ch.eq.'D'.and.mode1(n).lt.0) mode1(n)=+1
      end if
   13 continue
      do 14 n=1,n2
      if(abs(xm-chid2(n)).lt.2..and.abs(ym-abld2(n)).lt.0.02) then
      call pgsci(0)
      call pgpt1(chid2(n),abld2(n),icod2(n))
c      call pgsci(1)      
      CALL PGSCI(3)
      if(ch.eq.'A'.and.mode2(n).gt.0) mode2(n)=-2
      if(ch.eq.'D'.and.mode2(n).lt.0) mode2(n)=+2
      end if
   14 continue
      if(.not.loout) then
      call pgsci(0)
      call pgsls(1)
      call pgline(2,x,yav)
      call pgsls(4)
      call pgline(2,x,ysp)
      call pgsls(4)
      call pgline(2,x,ysm)
      call pgsci(1)
      end if
      end do
c  950 call gwquit(ir)
      return
      end

      subroutine avsig(av,sig,t_crt)
      parameter (MF1=65+235+2,MF2=13+15,MLN=MF1+MF2)
      common/pass/n1,n2,chid1(MF1),ewd1(MF1),abld1(MF1),
     & chid2(MF2),ewd2(MF2),abld2(MF2),mode1(MF1),mode2(MF2),
     & numb1(MF1),numb2(MF2)
c     calculate average
      sum=0.
      ic=0
      do n=1,n1
      if(mode1(n).gt.0) then
      ic=ic+1
      sum=sum+abld1(n)
      end if
      end do
      do n=1,n2
      if(mode2(n).gt.0) then
      ic=ic+1
      sum=sum+abld2(n)
      end if
      end do
      nc=ic
      av=sum/float(nc)
c     calculate sigma
      sumsq=0.
      ic=0
      do n=1,n1
      if(mode1(n).gt.0) then
      ic=ic+1
      sumsq=sumsq+(abld1(n)-av)**2
      end if
      end do
      do n=1,n2
      if(mode2(n).gt.0) then
      ic=ic+1
      sumsq=sumsq+(abld2(n)-av)**2
      end if
      end do
      nc=ic
      sig=sqrt(sumsq/float(nc))
      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


      function pp(t)
c     corresponding to the P function in table A1 of Taylor book
      pp=erf(t/1.41421356)
      return
      end

      FUNCTION erf(x)
      REAL erf,x
CU    USES gammp
      REAL gammp
      if(x.lt.0.)then
        erf=-gammp(.5,x**2)
      else
        erf=gammp(.5,x**2)
      endif
      return
      END

      FUNCTION gammp(a,x)
      REAL a,gammp,x
CU    USES gcf,gser
      REAL gammcf,gamser,gln
      if(x.lt.0..or.a.le.0.)pause 'bad arguments in gammp'
      if(x.lt.a+1.)then
        call gser(gamser,a,x,gln)
        gammp=gamser
      else
        call gcf(gammcf,a,x,gln)
        gammp=1.-gammcf
      endif
      return
      END

      SUBROUTINE gcf(gammcf,a,x,gln)
      INTEGER ITMAX
      REAL a,gammcf,gln,x,EPS,FPMIN
      PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
CU    USES gammln
      INTEGER i
      REAL an,b,c,d,del,h,gammln
      gln=gammln(a)
      b=x+1.-a
      c=1./FPMIN
      d=1./b
      h=d
      do 11 i=1,ITMAX
        an=-i*(i-a)
        b=b+2.
        d=an*d+b
        if(abs(d).lt.FPMIN)d=FPMIN
        c=b+an/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        del=d*c
        h=h*del
        if(abs(del-1.).lt.EPS)goto 1
11    continue
      pause 'a too large, ITMAX too small in gcf'

1     gammcf=exp(-x+a*log(x)-gln)*h
      return
      END

      SUBROUTINE gser(gamser,a,x,gln)
      INTEGER ITMAX
      REAL a,gamser,gln,x,EPS
      PARAMETER (ITMAX=100,EPS=3.e-7)
CU    USES gammln
      INTEGER n
      REAL ap,del,sum,gammln
      gln=gammln(a)
      if(x.le.0.)then
        if(x.lt.0.)pause 'x < 0 in gser'
        gamser=0.
        return
      endif
      ap=a
      sum=1./a
      del=sum
      do 11 n=1,ITMAX
        ap=ap+1.
        del=del*x/ap
        sum=sum+del
        if(abs(del).lt.abs(sum)*EPS)goto 1
11    continue
      pause 'a too large, ITMAX too small in gser'
1     gamser=sum*exp(-x+a*log(x)-gln)

      return
      END

      FUNCTION gammln(xx)
      REAL gammln,xx
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     *-.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)

      return
      END

      FUNCTION MAP1(XOLD,FOLD,NOLD,XNEW,FNEW,NNEW)
      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
