c
c
c***********************************************************************
c
      subroutine bsread (b,s,mr,nr,nc,blab,ibin,idiv,ngain,llo,lup,
c
     1                    vlo,vup,forte,iovl,work,nloc2)
c
c-----------------------------------------------------------------------
c
      dimension b(nr),s(nr),blab(nr),work(nloc2)
c
      dimension ident(20),llo(nr),lup(nr),vlo(nr),vup(nr),forte(nr)
c
c----------------------------------------------------------------------
c
      if(ibin.le.0) go to 100 
c
c*******  read and bin the data
c
      nloc = nloc2 / 2
c
c------------------------------------------------------------------
c
      call towcbn ( b,s,nr,ident,idiv,ngain,llo,lup,vlo,vup,forte,
c
     1              work(1),work(nloc+1),nloc )
c
c-------------------------------------------------------------------
c
      scoff = 0.001 
      do 10 i=1,nr
      blab(i) = vlo(i)
      s(i) = amax1(s(i),scoff*b(i))
   10 continue
      go to 101
c
c********   read already binned data
c
  100 continue
       if (ibin.eq.(-1)) call matin (blab,mr,nr,1)
       if (ibin.eq.(-1)) call matin (b,mr,nr,2)
       if (ibin.eq.0) call matin (blab,mr,nr,3)
      write(6,601)
      do 11 i=1,nr
      if(b(i).le.0.0)  perr = 0.0
      if(b(i).gt.0.0)  perr = 100. * s(i)/b(i)
   11 write(6,602) i,blab(i),b(i),s(i),perr
c
  101 if(ngain.lt.2) go to 900
c
c*******  adjust overlapping multiple gain channels
c
      fudge=1.0
      i1 = iovl + 1 
      i2 = iovl + 2 
      i3 = iovl + 3 
      i4 = iovl + 4 
      i5 = iovl + 5 
      i6 = iovl + 6 
      i7 = iovl + 7 
      covh=b(iovl)+b(i2)+b(i4)+b(i6)+1.e-20
      sovh=sqrt(s(iovl)**2+s(i2)**2+s(i4)**2+s(i6)**2)
      rerrh=sovh/covh
      covl=b(i1)+b(i3)+b(i5)+b(i7)+1.e-20
      sovl=sqrt(s(i1)**2+s(i3)**2+s(i5)**2+s(i7)**2)
      rerrl=sovl/covl
      fudge=covl/covh
      relfu=abs(1.-fudge)
c     decide on strategy
   12 if(rerrh-relfu)13,13,14 
   13 if(rerrl-relfu)16,16,14 
   14 fudge=1.0
      go to 900
   16 fudgeh=1.+(covl-covh)/covh*sovl/(sovh+sovl) 
   17 fudgel=1.+(covh-covl)/covl*sovh/(sovh+sovl) 
      do 20 i=1,iovl
      b(i)=b(i)*fudgeh
   20 s(i)=s(i)*fudgeh
      do 30 i=i2,i6,2
      b(i)=b(i)*fudgeh
   30 s(i)=s(i)*fudgeh
      do 32 i=i7,nr 
      b(i)=b(i)*fudgel
   32 s(i)=s(i)*fudgel
      do 33 i=i1,i5,2
      b(i)=b(i)*fudgel
   33 s(i)=s(i)*fudgel
c
  900 continue
c
   34 write (6,90)scoff,fudge,relfu
   90 format(1h0, 45hcounts and standard deviation in grouped bins
     4          ,5x,7hscoff =,f8.4,5x,7hfudge =,f8.4,5x,7hrelfu =,f8.4)
  601 format(1h0,' bin',7x,'lo p.h.',10x,'counts',10x,'error',7x,
     1 '(-error'/)
  602 format(i4,1p3e16.4,0pf12.1)
      return
      end 
c
c
c***********************************************************************
c
      subroutine towcbn (gc,s,nr,ident,idiv,ngain,llo,lup,vlo,vup,forte,
c
     1                   c,cre,nloc)
c
c-----------------------------------------------------------------------
c
c  towcbn reads binning table, normalization = calibration data, and
c         channel data for 1 or 2 gain runs, fore- and back-ground.
c         returns gc(grouped counts) and s(errors) with pulse-height
c         boundaries vlo = vup.
c  modified by d.t.ingersoll 7/80 to allow a nonzero channel offset
c         and reduced storage requirements.
c
c---------------------------------------------------------------------
c
      character*4 forb(2)
c
      dimension c(nloc),cre(nloc)
c
      dimension  gc(nr), s(nr), llo(nr), lup(nr), vlo(nr), vup(nr),
c
     1           forte(nr)
c
      dimension uc(2,2),zi(2,2),gn(2,2),gnom(2)
c
      dimension  faf(19), ident(20), irel(3)
c
      double precision   title1(4)
c
      data forb/ 'fore', 'back'/
c
c------------------------------------------------------------------------
c
      read (5,909) faf
c
c  read binning table
c
      read (5,900) npoint,ngrps
      irel(2) = npoint
c
c----------------------------------------------------------------------
c
      call binrd ( ngrps, llo, lup, nr, nbins )
c
c----------------------------------------------------------------------
c
      write (6,901) 
c  read normalization and calibration data
      do 15 k=1,ngain
      read(5,902) uc(k,1),uc(k,2),zi(k,1),zi(k,2),gn(k,1),gn(k,2),
     1            gnom(k)
      write(6,903) k,uc(k,1),uc(k,2),zi(k,1),zi(k,2),gn(k,1),gn(k,2), 
     1             gnom(k)
   15 continue
c
c  read foreground then background data 
c
      do 20 i=1,nbins
      gc(i) = 0.0
   20  s(i) = 0.0
c
      do 2000 nfb=1,2
c
      write(6,911) forb(nfb)
c
c---------------------------------------------------------
c
      call cread ( c,nloc,title1,irel,nerr,npoint,ident )
c
c---------------------------------------------------------
c
c  shift to zero channel offset and change actual gains to gnom
c
      do 320 k=1,ngain
         kk = 1 + irel(k) - irel(1)
c
c----------------------------------------------------------------
c
         call zshift ( c(kk), cre(kk), nloc, npoint, zi(k,nfb) )
c
c----------------------------------------------------------------
c
         do 321 i=1,npoint
            vvlo=float(i-1)*gnom(k)
            vvup=float(i)  *gnom(k)
            j = kk+i-1
c
c-------------------------------------------------------------------
c
      cre(j) = rebin ( c(kk), nloc,  1, npoint, 0.0, gn(k,nfb), 1.0,
c
     1                 vvlo, vvup )
c
c-------------------------------------------------------------------
c
  321    continue
  320 continue
c
      do 323 k=1,ngain
         kk = 1 + irel(k) - irel(1)
         if (title1(k)) 325,326,325
  326    call equiv(title1(k),6hno bkg) 
  325    write (6,910) title1(k),forb(nfb),k
c
c----------------------------------------------
c
  323    call cprint ( cre(kk), nloc, npoint )
c
c----------------------------------------------
c
c  calculate p.h. boundaries and grouped counts + errors
      sgn = (-1.0)**(nfb+1)
      do 1999 i=1,nbins
         cnts = 0.0 
         var = 0.0
         l1=llo(i)
         l2=lup(i)
         if (ngain.eq.1)  ll = 1
         if (ngain.gt.1)  ll = 1 + l1/irel(2)
         lllo = l1 - irel(ll) 
         llup = l2 - irel(ll) 
         vlo(i) = float(lllo-1)*gnom(ll)
         vup(i) = float(llup)  *gnom(ll)
         if(lllo.gt.npoint) go to 1999
         if(llup.gt.npoint)  l2=npoint + irel(ll) 
      do 1998 l=l1,l2
         lk = l - irel(1)
         cnts = cnts + cre(lk)/uc(ll,nfb)
 1998    var = var + (cre(lk)+0.1)/uc(ll,nfb)**2
         gc(i) = gc(i) + sgn*cnts
         s(i) = s(i) + var
 1999    continue
c
 2000 continue
c
      do 30 i=1,nbins
      binw = 1.0
      if (idiv.eq.1)  binw = vup(i) - vlo(i)
      gc(i) = gc(i) / binw
   30 s(i) = sqrt(s(i)) / binw
c
c  adjust grouped counts and errors according to forte factors.  this 
c    allows for decreased efficiency of mca at low channels.
c
      do 2 i=1,nbins
        v=0.5*(vlo(i)+vup(i)) 
        if(v-0.10) 4,8,8
    4   if(v-0.01) 5,5,6
    5   forte(i) = faf(1)
        gc(i) =0.0
        s(i) = 0.0
        go to 2
    8   forte(i) = 1.0
        go to 7
    6   z = 200.0*(v-0.01)
        j=z
        z = z-float(j)
        forte(i)=(1.0-z)*faf(j+1)+z*faf(j+2)
    7   continue
        gc(i) = gc(i)/forte(i)
        s(i) = s(i)/forte(i) + 0.2*(1.0-forte(i))*gc(i)
    2 continue
      write (6,906) title1(1) 
      write (6,401) 
      do 50 i=1,nbins
   50 write (6,934) i,gc(i),s(i),llo(i),lup(i),vlo(i),vup(i),forte(i) 
  999 return
  401 format (1h ,6x,4hbins,6x,6hcounts,11x,5herror,6x,6hlochan,2x,
     16hupchan,6x,7hlo p.h.,9x,7hup p.h.,5x,6hf.a.f./)
  900 format(2i10)
  901 format(91h0  gain  fore norm   back norm   fore zero   back zero
     1 fore gain   back gain   noml gain /)
  902 format(7e10.4)
  903 format(1h ,i6,1p7e12.4) 
  906 format (1h1,10x,a6,11x,19hgrouped counts/unit,/,1h )
  909 format(19f4.3)
  910 format(1h1,a6,11x,13hrebin counts ,a4,13hground gain =,i2,/,1h )
  911 format(1h1,5x,11hraw counts ,a4,6hground)
  934 format (1h ,i8,1p2e16.4,2i8,2e16.4,0pf10.5) 
      end 
c
c
c******************************************************************** 
c
      subroutine binrd ( ngrps, llo, lup, nr, nbins )
c
c-------------------------------------------------------------------- 
c
c     reads bin card(s) with (klo(i),kup(i),kdel(i)),i=1,ngrps
c     and converts to expanded bin list (llo(j),lup(j)),j=1,nbins
c
c     needs ngrps from calling program
c     returns llo,lup,nbins
c
c-------------------------------------------------------------------- 
c
      dimension llo(nr),lup(nr),kup(40),klo(40),kdel(40)
c
c
c  read binning card(s)  10 max
c
      read (5,901) ( klo(i),kup(i),kdel(i), i=1,ngrps )
c
c  write out binning table on standard output tape
c
      write (6,902) 
      do 3 m=1,ngrps
      aup = kup(m)
      alo = klo(m)
      adel = kdel(m)
      bin = (aup-alo+1.0)/adel
      write (6,903) m,klo(m),kup(m),kdel(m),bin
    3 continue
c     convert table 
      nbins=0
      do 20 i=1,ngrps
           k1=klo(i)
           k2=kup(i)
           k3=kdel(i)
      k3 = max0(k3,1)
           do 10 k=k1,k2,k3
               nbins=nbins+1
               llo(nbins)=k
               lup(nbins)=k+k3-1
   10      continue 
   20 continue
c     write out computed no bins
      write (6,904) nbins
      return
  901 format(12i6)
  902 format(1h0,13hbinning table/6h0group,4x,6hloedge,4x,6hupedge,4x,8h
     1chan/bin,5x,4hbins)
  903 format (1h ,i3,6x,i5,5x,i5,8x,i3,7x,f4.1)
  904 format(43x,4h----/25x,14htotal no. bins,4x,i3,/////)
      end 
c
c
c**********************************************************************
c
      subroutine cread (c,n,title1,irel,nerr,npoint,ident)
c
c----------------------------------------------------------------------
c
c     to read standard bsf cards with no relocation, terminate with
c     blank card and preceed with  .......................
c
c     col 1-6      .....col 14-31....
c     dynmat       10 (1a6,i4,10f7.0)
c
c  d.t.i. -  irel now functions as intended to point to first channel 
c            of fore or back data.  irel(1) is internally moved to 0. 
c
c----------------------------------------------------------------------
c
      dimension c(n),card(20),fmt(20),ident(20)
c
      dimension  irel(3)
c
      double precision  title1(4)
c
c----------------------------------------------------------------------
c
      nerr = 0
      kk = 0
      do   30  k = 1,n
   30 c(k) = 0.0
    1 kk = kk + 1
      read (5,901) card(1),irel(kk),nword,(fmt(j),j=1,10)
      nchan = 0
      call equiv (x,4hdynm)
      if (card(1)/2.0-x/2.0) 3,2,3
    3 call equiv(x,4h    )
      if (card(1)/2.0-x/2.0) 2,88,2
    2 write (6,903) card(1),irel(kk),nword,(fmt(j),j=1,10)
      read (5,900) ident
      write (6,904) ident
   10 read(5,fmt) title1(kk),io,(card(j),j=1,nword)
      write(6,905) title1(kk), io, ( card(j), j=1,nword )
      io = io+irel(kk)-irel(1)
      do  17  i=1,nword
         nchan = nchan + 1
c                      *
c              .....index ok.....
c              no             yes
c       ....word zero....      *
c       no            yes    store
c        *             *       *
c      alarm           *       *
c        *             *       *
c       ............continue......
c                      *
           ii = io + i - 1
   12      if   (ii)             14,14,13
   13      if   (ii-n)           16,16,14
   14      if   (card(i))        99,17,99
   16      c(ii) = card(i)
   17      continue 
      if (nchan-npoint) 10,1,10
   88 return
   99 nerr = nerr+1 
      write (6,902) n,card
      go to 10
  900 format (20a4) 
  901 format (1a4,2x,i6,i3,1x,10a4)
  902 format ('0*** cread error, n =',i6,' card is **',20a4,2h**)
  903 format (//1h0,1a4,2x,i6,i3,1x,10a4)
  904 format (1h ,20a4)
  905 format (1h ,1a6,1i6,10f7.0,1a6)
      end 
c
c
c*************************************************************
c
      subroutine cprint ( c, nloc, np ) 
c
c-------------------------------------------------------------
c
      dimension  c(nloc), xbue(8), ibue(8)
c
c-------------------------------------------------------------
c
      tcount = 0.0
      npages = 1 + (np-1)/400 
      do   10   k=1,npages
      if   (k-1)   5,5,4
    4 write (6,904) 
    5 npoint = min0(400,np-(k-1)*400)
      nrows = min0(npoint,50) 
      do 7 i=1,nrows
           ncols=(npoint/50)+min0(1,(npoint-50*(npoint/50))/i)
      ncols = max0(1,ncols)
           do 8 j=1,ncols
               nc=i+(j-1)*50+(k-1)*400
               xbue(j)=c(nc)
      tcount = tcount+xbue(j) 
               ibue(j)=nc
    8          continue
      write (6,905) ( ibue(l),xbue(l), l=1,ncols )
    7      continue 
   10 continue
      write (6,906) tcount
      return
  904 format(1h1)
  905      format(1h ,8(i5,f9.1,2x))
  906 format (1h0/1h0,24htotal number of counts =,f12.1)
      end 
c
c
c*************************************************************************
c
      function rebin ( a, nloc, mr, n, rholo, rhodel, c, vlo, vup )
c
c------------------------------------------------------------------------
c
      dimension a(nloc)
c
c     the vector a gives the components of a histogram
c     the bottom edge of the lower bar is at rholo, and the width
c     of the bars is rhodel  this routine adds all the bins between
c     vlo and vup and splits any fractional bins linearly.
c     for greater flexibility, v = c*rho, where c is an input parameter.
c
c------------------------------------------------------------------------
c
    1 rlo=amax1(vlo/c,rholo)
    2 rup=amax1(vup/c,rholo)
c
   10 x1 = 2.0 + (rlo-rholo)/rhodel
      i1 = x1
      i1 = min0(i1,n+1)
      x1 = x1-float(i1)
      x1 = amin1(x1,1.0)
c
   15 x2 = 1.0 + (rup-rholo)/rhodel
      i2 = x2
      i2 = min0(i2,n)
      x2 = x2-float(i2)
      x2 = amin1(x2,1.0)
c
   20 l1 = ((i1-1) - 1)*mr + 1
      l2 = (i2-1)*mr + 1
      rebin = (1.0-x1)*a(l1) - (1.0-x2)*a(l2)
c
      if   (i2-i1+1)   40,40,30
   30 do   37   i = i1,i2
           l3 = (i-1)*mr + 1
           rebin = rebin + a(l3)
   37      continue 
c
   40 return
      end 
c
c
c*******************************************************************
c
      subroutine zshift ( a, b, nloc, n, zint )
c
c-------------------------------------------------------------------
c
c     purpose:  to correct a spectrum for zero intercept
c               ph = m*a  replaces  ph = m*(a-zint)
c               b is needed for scratch 
c
      dimension  a(nloc), b(nloc)
c
c------------------------------------------------------------------
c
      do 10 i=1,n
   10 b(i)=0.0
      fpart=zint - aint(zint) 
      if(zint) 18,24,11
c zint \ 0.0
   11 do 12 i=1,n
      k=i-1
      if(float(i-1).gt.zint)  go to 13
   12 continue
   13 do 14 i=1,n
      in = i
      if(k+i.gt.n)  go to 15
   14 b(i)=(1.0-fpart)*a(k+i-1) + fpart*a(k+i)
   15 b(in)=(1.0-fpart)*a(n)
      go to 22
c zint ) 0.0
   18 do 19 i=1,n
      k=i-1
      if(-float(i-1).lt.zint)  go to 20 
   19 continue
   20 b(k)=(1.0+fpart)*a(1)
      kk=k+1
      do 21 i=kk,n
   21 b(i)=(1.0+fpart)*a(i-k+1)-fpart*a(i-k)
c
   22 do 23 i=1,n
   23 a(i) = b(i)
   24 return
      end 
c
c
c***********************************************************
c
      subroutine equiv(x,y)
c
c---------------------------------------------------------- 
c
      x = y
c
      return
c
c
      end 
c
c
c******************************************************************** 
c
      subroutine wread (w, nc, nw, ntry, elab, plab, pctwid )
c
c-------------------------------------------------------------------- 
c
      dimension  w(nc), elab(nc), plab(nw), pctwid(nw)
c
c-------------------------------------------------------------------- 
c
c     the following is the ordinary definition of the gauss function
c
      gauss(x) = 0.3989422804014 * exp(-0.5*x**2) 
c
      sigma = (.01 * pctwid(ntry) * plab(ntry))/(2.0 * 1.17741002251) 
      do   10   i = 1,nc
           difos = abs(elab(i) - plab(ntry))/sigma
           if   (difos - 5.5)   3,2,2
    2           w(i) = 0.0
                     go to 10 
c*****    3           w(i) = 1.0e+6 * gauss(difos)/sigma
    3           w(i) = gauss(difos)/sigma
   10 continue
      return
      end 
c
c
c*************************************************************************
c
      subroutine matin ( a, mr, nr, nc )
c
c-------------------------------------------------------------------------
c
c     this program reads data into an array.
c     a format card must precede the data,
c       i.e.,format i  nn (...)
c            with  'format'  in col. 1-6,
c                  i = 'row' for row format, i = 'col' for column format.
c                  nn = no. data words per card,
c               (...) = format specification.
c     each data card must have two index numbers
c       indicating the matrix position of the first data word.
c     the reading of cards is completed by one blank card.
c     this routine has a provision for changing formats in midstream. 
c     a card with both indices negative implies that a format card follows.
c
c--------------------------------------------------------------------------
c
      character*3 irow,iroc
c
      character*4 icheck,iword, fmt
c
      dimension a(mr,nc),fmt(17),xlist(78)
c
      data icheck /'FORM'/
c
      data irow /'ROW'/
c
c-------------------------------------------------------------------------
c
      nerror = 0
      do  4  i = 1,mr
           do  4  j = 1,nc
    4           a(i,j) = 0.0
c  ***read format card
    5 read (5,900) iword,iroc,nn,fmt
  900 format(1a4,3x,1a3,1x,i2,1x,16a4,1a2)
      if (iword.eq.icheck) go to 10
    8 write (6,904) 
      go to 888
c  ***read data card
   10 read (5,fmt) i,j,(xlist(k),k=1,nn)
c      if(nn.eq.2)print 1000,i,j,(xlist(k),k=1,nn)
c 1000 format(' i=', i3,' j= ',i3,' data values ',2(e10.4,2x))
c  ***check for new format card
   32 if(i)50,11,11 
   50 if(j)5,721,701
c  ***check for blank card to end list
   11 if (i +j)  12,888,12
   12 i = max0(i,1) 
      j = max0(j,1) 
      if (iroc.ne.irow) go to 213
c *** row format ***
   13     if (i-nr)  14,14,701
  701 if(nerror-50) 702,702,10
  702 write (6,901) i,j,mr,nr,nc,(xlist(k1),k1=1,4)
      nerror=nerror+1
          go to 10
   14 do  20  jj = 1,nn
           jjj = j + jj - 1
           if (xlist(jj))16,20,16
   16         if (jjj - nc)  17,17,703
  703 if(nerror-50) 704,704,10
  704 write (6,902) i,jjj,mr,nr,nc,(xlist(k1),k1=1,4)
      nerror=nerror+1
                 go to 10
   17         if (a(i,jjj)) 705,18,705
  705 if(nerror-50) 706,706,20
  706 write (6,903) i,jjj,mr,nr,nc,(xlist(k1),k1=1,4)
      nerror=nerror+1
                 go to 20
   18         a(i,jjj) = xlist(jj)
c      if(nn.eq.2)print 1001,i,jjj,jj,a(i,jjj)
c 1001 format(' i,jjj,jj ',3(i3,2x),' a(i,jjj) ',e10.4)
   20 continue
   21 go to 10
c *** column format ***
  213     if (j-nc)  214,214,721
  721 if(nerror-50) 722,722,10
  722 write (6,902) i,j,mr,nr,nc,(xlist(k1),k1=1,4)
      nerror=nerror+1
          go to 10
  214 do  220  ii = 1,nn
           iii = i + ii - 1
           if (xlist(ii))  216,220,216
  216          if (iii-nr)  217,217,723 
  723 if(nerror-50) 724,724,10
  724 write (6,901) iii,j,mr,nr,nc,(xlist(k1),k1=1,4)
      nerror=nerror+1
                 go to 10
  217          if (a(iii,j))  725,218,725
  725 if(nerror-50) 726,726,220
  726 write (6,903) iii,j,mr,nr,nc,(xlist(k1),k1=1,4)
      nerror=nerror+1
               go to 220
  218      a(iii,j) = xlist(ii)
  220 continue
      go to 10
  888 return
  901 format (22h0illegal 1st index, i=,i3,2x,2hj=,i3,2x,3hmr=,i3,2x, 
     -        3hnr=,i3,2x,3hnc=,i3,5h   **,4e12.3,3h...)
  902 format (22h0illegal 2nd index, i=,i3,2x,2hj=,i3,2x,3hmr=,i3,2x, 
     -        3hnr=,i3,2x,3hnc=,i3,5h   **,4e12.3,3h...)
  903 format (22h0duplication error, i=,i3,2x,2hj=,i3,2x,3hmr=,i3,2x, 
     -        3hnr=,i3,2x,3hnc=,i3,5h   **,4e12.3,3h...)
  904 format(1h0,51h***  program expects format card at this point  ***)
      end 
c
c
c***********************************************************************
c
      subroutine wotac(x,ltbl,ncol,lg,top1,top2,top3,n6)
c
c----------------------------------------------------------------------
c
c   *** output writes 1,2, or 3-d arrays
c
c
      character*4 top1,top2,top3
c
      dimension  x(ltbl,ncol,lg)
c
c----------------------------------------------------------------------
c
      do 1 l=1,lg
      i02=0
      i03=(ncol+9)/10
      if(lg.gt.1)write(n6,10) top3,l
   10 format(/1h0,2x,a4,i5)
      do 2 i=1,i03
      i01=i02+1
      i02=min0(i01+9,ncol)
      write (n6,20) top1,(top2,j,j=i01,i02)
      do 3 k=1,ltbl 
      if(k.eq.1)go to 4
      do 6 j=i01,i02
    6 if(x(k,j,l).ne.x(k-1,j,l))go to 7 
      if(ke.eq.ks)go to 8
      ke=k
      if(k.eq.ltbl)go to 7
      go to 3
    8 if(k.eq.ltbl)go to 4
      do 9 j=i01,i02
    9 if(x(k,j,l).ne.x(k+1,j,l))go to 4 
      ke=ke+1
      go to 3
    7 if(ke.eq.ks)go to 4
                  write(n6,40) top1,ks,top1,ke
      if(k.eq.ltbl.and.ke.eq.k)  go to 3
   40 format(8x,a4,i5,6h thru ,a4,i5,14h same as above)
    4 write(n6,30) k,(x(k,j,l),j=i01,i02)
      ks=k+1
      ke=ks
    3 continue
    2 continue
    1 continue
   20 format('0 ',a4,3x,a4,i5,9(3x,a4,i5))
   30 format(i6,1p10e12.4)
      return
      end 
