        character ifile*80,pre*76,id*56,pgm*8
        dimension x(2000),y(2000),z(2000),dip(2000),trend(2000)
        dimension out(2000)

        data ddval/1.e30/
        data dval/1.e38/
        d2r=3.1415927/180.
c
        print*,'enter filename prefix for derivatives:'
        read(5,100) pre
        len=len_trim(pre)
  100   format(a)
   10   ifile=pre(1:len)//'.x'
        open(11,file=ifile,form='unformatted',status='old')
        read(11) id,pgm,nc,nr,nz,xo,dx,yo,dy
        ifile=pre(1:len)//'.y'
        open(12,file=ifile,form='unformatted',status='old')
        read(12) id,pgm,nc1,nr1,nz,xo,dx,yo,dy
        if(nc1.ne.nc.or.nr1.ne.nr) stop 'grids are different sizes'
        ifile=pre(1:len)//'.z'
        open(13,file=ifile,form='unformatted',status='old')
        read(13) id,pgm,nc1,nr1,nz,xo,dx,yo,dy
        if(nc1.ne.nc.or.nr1.ne.nr) stop 'grids are different sizes'
        print*,'enter dip angle grid:'
        read(5,100) ifile
        open(15,file=ifile,form='unformatted',status='old')
        read(15) id,pgm,nc1,nr1,nz1,xo1,dx1,yo1,dy1
        if(nc1.ne.nc.or.nr1.ne.nr) stop 'grids are different sizes'
        print*,'enter dip trend grid:'
        read(5,100) ifile
        open(16,file=ifile,form='unformatted',status='old')
        read(16) id,pgm,nc1,nr1,nz1,xo1,dx1,yo1,dy1
        if(nc1.ne.nc.or.nr1.ne.nr) stop 'grids are different sizes'
        print*,'enter output grid:'
        read(5,100) ifile
        open(17,file=ifile,form='unformatted',status='unknown')
        write(17) id,pgm,nc,nr,nz,xo,dx,yo,dy
        print*,'enter inc,dec:'
        read(5,*) ainc,adec
        ainc=ainc*d2r
        adec=adec*d2r
c
        do 50 j=1,nr
        read(11) dlt,(x(i),i=1,nc)
        read(12) dlt,(y(i),i=1,nc)
        read(13) dlt,(z(i),i=1,nc)
        read(15) dlt,(dip(i),i=1,nc)
        read(16) dlt,(trend(i),i=1,nc)
        do 40 i=1,nc
        out(i)=dval
        if(trend(i).gt.ddval) go to 40
        if(dip(i).gt.ddval) go to 40
        if(x(i).gt.ddval) go to 40
        if(y(i).gt.ddval) go to 40
        if(z(i).gt.ddval) go to 40
        theta = trend(i)*d2r
        capa = theta - adec
        term = 2.*atan2(tan(ainc),cos(capa))
        adip = dip(i)*d2r
        phi = term - adip - 90.*d2r
        dh = x(i)*sin(theta) + y(i)*cos(theta)
        out(i) = cos(phi)*z(i) + sin(phi)*dh
   40   continue
        write(17) dlt,(out(i),i=1,nc)
   50   continue
        stop
        end

