parameter (npzs=2,npz=2,npr=5,nply=2)
        complex rsp(5,npzs,npz,npr),rspcheck(5,npr)
        real gf(5)
        real r(npr),z(npz),zs(npzs)
	real beta(nply),alpha(nply),rho(nply),qb(nply),qa(nply),th(nply)
        open(8,file='wrt8.dat',status='old',form='unformatted')
        length=8*(5*npr)
        open(10,file='gf.dat',form='binary',status='unknown',
     #  recl=length,access='direct')
	read (8) NLY,NR,NZ,NZS,NFREQ,FI,DF,RINGLD,REFLAY,
     *   (R(I),I=1,NR),(Z(I),I=1,NZ),(ZS(I),I=1,NZS),
     *  (BETA(I),ALPHA(I),RHO(I),QB(I),QA(I),TH(I),I=1,NLY)     
        write (*,*) NLY,NR,NZ,NZS,NFREQ,FI,DF,RINGLD,REFLAY,
     *   (R(I),I=1,NR),(Z(I),I=1,NZ),(ZS(I),I=1,NZS),
     *  (BETA(I),ALPHA(I),RHO(I),QB(I),QA(I),TH(I),I=1,NLY)     

        read (8) rsp
        do i=1,nr
        do j=1,nz
        do k=1,nzs
        write (*,'(3i5)') i,j,k
        call writ(rsp(1,k,j,i))
        enddo
        enddo
        enddo

        do izs=1,2
        do iz=1,2
        irec=(izs-1)*nz+(iz)
        write (10,rec=irec) ((rsp(m,izs,iz,i),m=1,5),i=1,5)
        write (*,*) 'irec=',irec
        enddo
        enddo

        write (*,*) 'Enter ir,iz,izs to check:'
200     read (*,*) ir,iz,izs
        if(ir.eq.0)go to 300
        irec=(izs-1)*nz+(iz)
        write (*,*) 'Record Number:',irec
        read (10,rec=irec) ((rspcheck(m,i),m=1,5),i=1,5)
c----check static solution for accuracy----
        call mindlin(r(ir),z(iz),zs(izs),gf)

        do m=1,5
        write (*,201) rsp(m,izs,iz,ir),rspcheck(m,ir),gf(m)
201     format(2f11.4,2x,2f11.4,2x,f11.4)
        enddo
        go to 200
300     stop
        end
        subroutine writ(a)
        complex a(5)
        write (*,'(2f10.5)') a
        return
        end
c----
        subroutine mindlin(r,z,c,gf)
        real gf(5)
        real mu
        pi=3.1415926535
        theta=0.
c----mu is the Poisson's Ratio----
c----gf is to be multiplied by 1/(G*r)----
        mu=1./3.
        c1=(1.-mu)
        c2=(1.-2.*mu)
        c3=(3.-4.*mu)
        r1=sqrt(r**2+(z-c)**2)
        r2=sqrt(r**2+(z+c)**2)
        t1=(z-c)/r1**3
        t2=c3*(z-c)/r2**3
        t3=-4.*c1*c2/(r2*(r2+z+c))
        t4=6.*c*z*(z+c)/r2**5
        u=r*(t1+t2+t3+t4)/(16.*pi*c1)
        t1=c3/r1
        t2=(8.*c1*c1-c3)/r2
        t3=(z-c)**2/r1**3
        t4=(c3*(z+c)**2-2.*c*z)/r2**3
        t5=(6.*c*z*(z+c)**2)/r2**5
        w=(t1+t2+t3+t4+t5)/(16.*pi*c1)
        x=r*cos(theta)
        t1=c3/r1+1./r2+x**2/r1**3+c3*x**2/r2**3
        t2=(2.*c*z/r2**3)*(1.-3.*x**2/r2**2)
        t3=(4.*c1*c2/(r2+z+c))*(1.-x**2/(r2*(r2+z+c)))
        uh=(t1+t2+t3)/(16.*pi*c1)
        t1=(z-c)/r1**3+c3*(z-c)/r2**3
        t2=-6.*c*z*(z+c)/r2**5
        t3=4.*c1*c2/(r2*(r2+z+c))
        wh=(t1+t2+t3)*x/(16.*pi*c1)
        t1=1./r1**3 + c3/r2**3 - 6.*c*z/r2**5
        t2=-4.*c1*c2/(r2*(r2+z+c)**2)
        vh=(t1+t2)*r**2/(16.*pi*c1)
        gf(1)=u*r
        gf(2)=w*r
        gf(3)=uh*r
        gf(4)=(vh-uh)*r
        gf(5)=wh*r
        return
        end