cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Input a grid of z values (representing a surface) and a thickness -
c      produce xyz values for a surface offset by that constant thickness
c
c     saltus 16sep99 for aleut project

      character title*80,pgm*8,ingrid*80,outxyz*80
      real z1(10000)
      real z2(10000)
      real nx,ny,nz
c     parameter (dval='ffff7fff'x)
      parameter (dval=1.0e38)
      parameter (dvaltest=1.0e20)

      call pfinit('gdthicklayer')

      print *,
     & 'I calculate xyz positions of the other edge of a thick layer',
     & 'you provide a surface (grid) and constant thickness',
     & '(use a negative thickness if you are giving the top)'

      ival=iaquest('Input grid (topo in km)',ingrid,'(a80)',0)
      call gopen(21,ingrid,'old','read',ierr)
      ival=irquest('Thickness',thick,'(f15.5)',1)
      if (thick.lt.0.) then
        outxyz='bottom.xyz'
      else
        outxyz='top.xyz'
      end if
      ival=iaquest('Output xyz file',outxyz,'(a80)',1)
      open(11,file=outxyz,form='formatted',status='unknown')
c     call getfile(21,'*GIVE IN GRID:','in','unformatted')
c     call getfile(22,'*GIVE OUT GRID:','out','unformatted')

      call gh1f4('r',21,title,ncol,nrow,xo,dx,yo,dy,ierr)
c     call rdheader(21,id,pgm,ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl)
c     ival=iaquest('Title',title,'(a80)',1)
c     pgm='gdtmoho '
c     title='Crustal thickness from topography'
c     call wrheader(22,id,pgm,ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl)
c     call gh1f4('w',22,title,ncol,nrow,xo,dx,yo,dy,ierr)

      j=1
      call growf4('r',21,j,z1,ncol,ierr)
      do 50 j=2,nrow
      call growf4('r',21,j,z2,ncol,ierr)
       do 40 i=1,ncol-1
         if ((z1(i).lt.dvaltest).and.
     &       (z1(i+1).lt.dvaltest).and.
     &       (z2(i).lt.dvaltest)) then
         bx=(i+1)*dx
         by=(j-1)*dy
         bz=z1(i+1)
         cx=i*dx
         cy=j*dy
         cz=z2(i)
         ax=cx
         ay=by
         az=z1(i)
         nx=by*cz-bz*cy
         ny=bz*cx-bx*cz
         nz=bx*cy-by*cx
c        if (nz.lt.0.) then
c          nx=-nx
c          ny=-ny
c          nz=-nz
c        end if
         rlen=sqrt(nx**2+ny**2+nz**2)
         x=xo+ax+thick*nx/rlen
         y=yo+ay+thick*ny/rlen
         z=az+thick*nz/rlen
         write(11,*)x,y,z,xo+ax,yo+ay,az,z-az
         end if
         z1(i)=z2(i)
 40    continue
      z1(ncol)=z2(ncol)
c     call growf4('w',22,j,z,ncol,ierr)
c     call putrow(22,ncol,z)
 50   continue
 
      call gclose(21,'keep')
      close(11)
c     call gclose(22,'keep')

      print *,'Have a thick day...'

      end