IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	dimension fe(10,100),dzxe(10,100),tinv(100)
	dimension c(10),af(100),zita(100)
	dimension t(100),f(100),zx(100),dtim(100),e(10),d0(10)
	character*9 tab1
	tab1=char(09)
	open(unit=7,file='temstep.in',status='old')
	open(unit=15,file='arr-me.in',status='old')
	open(unit=16,file='arr-me.out',status='unknown')
	open(unit=8,file='each-me.gr',status='unknown')
	open(unit=10,file='arr-me.dat',status='unknown')
	open(unit=1,file='each-me.dat',status='unknown')
	open(unit=11,file='logr-me.dat',status='unknown')
	pi=3.14159265
	ro=1.987E-3
	print *, 'spheres = 1, slabs = 2'
	read *, imp
	print *, 'imp=', imp
	b=6.
	acut = 0.85

	read(15,120)nsamp
	do 11 j=1, nsamp
	zita(j)=0.
   	read(15,140)e(j),d0(j),c(j)
	tc = tc + c(j)
	d0(j) = 10.**d0(j)
	if (imp.eq.1) goto 11
	b=8.
	d0(j)=d0(j)/4.
	print *, d0(j)
	acut=0.50
   11	continue
	if(abs(tc-1.).gt.1e-8) then
	print *,'Error  -  your concentration c(j) sum ',tc
	stop
	endif
	read(15,*)slop,ord
  120	format(I5)


	read(7,120)nstep
	do 10 nt=1,nstep
	f(nt)=0.
	read(7,*)t(nt),dtim(nt)
	dtim(nt)=60.*dtim(nt)
	do 20 j=1,nsamp
	fe(j,nt)=1.
	zita(j)=d0(j)*dtim(nt)*dexp(-e(j)/ro/t(nt))+zita(j)
	do 30 m=1,75000,imp
	xm=m
	xu=(pi*xm)**2*zita(j)
	if (xu.gt.99.)goto 20
	fe(j,nt)=fe(j,nt)-b*dexp(-xu)/(pi*xm)**2
   30	f(nt)=f(nt)-b*c(j)*dexp(-xu)/(pi*xm)**2
   20	continue
	nnt=nt
	f(nt)=1.+ f(nt)
	if (f(nt).gt..9999)then 
	nnt=nt-1
	goto 40
	endif
   10	continue
   40	continue

c	INVERSION OF THE DATA

 	do 50 k=1,nnt
        if (f(k).gt.acut) goto 60
	if (imp.eq.1) goto 22
	zx(k+1)=pi*(f(k)/4.)**2
	goto 50
   22   zx(k+1)=(2.-pi/3.*f(k)-2.*dsqrt(1.-pi/3.*f(k)))/pi
        goto 50
   60   zx(k+1)=-dlog(pi**2/b*(1.-f(k)))/pi**2
   50   continue
        zx(1)=0.
        do 80 k=1,nnt
        af(k)=(f(k)+f(k-1))/2.*100.
        dzx = dlog10((zx(k+1)-zx(k))/dtim(k)*imp**2)
        tinv(k)=1./t(k)*10000.
        xlogr=(ord-slop*tinv(k)-dzx)/2.
	write(10,110)tinv(k),tab1,dzx
	write(11,110)100.*f(k-1),tab1,xlogr
	write(11,110)100.*(f(k)-0.000001),tab1,xlogr 
  80   write(16,110)tinv(k),tab1,dzx,tab1,f(k)*100.,tab1,xlogr,tab1
     $  ,af(k)
 
	do 200 j=1,nsamp
	do 150 k=1,nnt
	if (fe(j,k).le.acut) then

	if (imp.eq.2) then
	zx(k+1)=pi*(fe(j,k)/4.)**2
	else
  	zx(k+1)=(2.-pi/3.*fe(j,k)-2.*dsqrt(1.-pi/3.*fe(j,k)))/pi 
        endif

	else

        if (fe(j,k).lt.1.)then
        zx(k+1)=-dlog(pi**2/b*(1.-fe(j,k)))/pi**2
	else
	zx(k+1) = 0.
	endif

	endif

  150   continue

        zx(1)=0.
        do 180 k=1,nnt
	if (zx(k+1).gt.0.) then
        dzxe(j,k) = dlog10((zx(k+1)-zx(k))*imp**2/dtim(k))
	endif
  180	continue
  200	continue
	do 195 j=1,nsamp
	do 190 k=1,nnt
	write(8,110)tinv(k),tab1,fe(j,k)
	if (dzxe(j,k).eq.0.)goto 190
        write(1,110)tinv(k),tab1,dzxe(j,k)
  190	continue
	write(8,130)
	write(1,130)
  195	continue

  110	format(1x,5(f12.8,a1))
  130   format(1x,'$')
  140	format(g20.8)
 	end