      program SubductionGeometry

      implicit none

      integer level,nx,ny,i,j,nx1,nxs,ny1,ny2
      integer i1,i2,i3,ij,k
      double precision x1,x2,x3,y1,y2,y3,z1,z2,z3
      double precision xne,yne,zne,xyzn
      integer nnode,nnode0,nxt,nxb,nnodemax,nelem,nelemmax
      double precision xl1,yl1,yl2,xs,dip,length,thick,pi,zfree
      double precision x,y,z,xsp,dist,depth
      double precision, dimension(:),allocatable::xp,yp,zp,xn,yn,zn
      integer,dimension(:,:),allocatable::icon

      pi=atan(1.)*4.

! xl1 is where the plate starts in the x-direction;
! xs is where the subduction starts; dip is the dip of the subducting
! slab; length is it length (along the dipping plane); thick is the
! plate thickness; yl1 is where the plate starts in the y direction and
! yl2 is where it ends

      xl1=.4376d0
      xs=.5624d0
      dip=60.d0
      dip=dip*pi/180.d0
      thick=.01d0
      depth=.1d0
      length=(depth-thick-.0001)/sin(dip)
      yl1=.4376d0
      yl2=.5624d0
      xsp=xs-thick*tan(dip/2.d0)
      zfree=1.d0

      level=8

      nx=2**level+1
      ny=2**level+1
      nnodemax=nx*ny*2
      nelemmax=nnodemax*3

      allocate (xp(nnodemax),yp(nnodemax),zp(nnodemax))
      allocate (xn(nnodemax),yn(nnodemax),zn(nnodemax))
      allocate (icon(3,nelemmax))

      nx1=int(xl1*nx)
      ny1=int(yl1*ny)
      ny2=int(yl2*ny)
      nxs=int(nx*xs)
      nxb=int(length*nx)
      nxt=nxb+1

      nnode=0
      nelem=0

      do j=1,ny1
! no plate
      y=float(j-1)/float(ny1-1)*yl1
      do i=1,nx1
      x=float(i-1)/float(nx1-1)*xl1
      z=zfree
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (j.ne.ny1) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
      enddo
      do i=nx1+1,nxs
      x=xl1+(xsp-xl1)*float(i-nx1)/float(nxs-nx1)
      z=zfree
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (j.ne.ny1) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
      enddo
      do i=nxs,nx
      x=xs+(1.-xs)*float(i-nxs)/float(nx-nxs)
      z=zfree
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (i.ne.nx.and.j.ne.ny1) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
      enddo
      enddo

      do j=ny1,ny2
! no plate
      do i=1,nx1
      x=float(i-1)/float(nx1-1)*xl1
      y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
      z=zfree
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
      enddo
! horizontal plate
      do i=nx1,nxs+1
      x=xl1+(xsp-xl1)*float(i-nx1)/float(nxs+1-nx1)
      y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
      z=zfree-thick
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
      if (i.ne.nxs+1.and.j.eq.ny1) call add_element (nelem,icon,nnode-2-nx,nnode-1-nx,nnode+1,nnode,nelemmax)
      if (i.ne.nxs+1.and.j.eq.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+4+nx+nxb+nxt,nnode+3+nx+nxb+nxt,nelemmax)
      enddo
! bottom dipping plate
      do i=1,nxb+1
      x=xsp+length*cos(dip)*float(i)/float(nxb+1)
      y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
      z=zfree-thick-length*sin(dip)*float(i)/float(nxb+1)
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
      if (j.eq.ny1) call add_element (nelem,icon,nnode+nxb+nxt-2*i+3,nnode+nxb+nxt-2*i+2,nnode,nnode-1,nelemmax)
      if (j.eq.ny2) call add_element (nelem,icon,nnode-1,nnode,nnode+nxb+nxt-2*i+2,nnode+nxb+nxt-2*i+3,nelemmax)
      enddo
! top of dipping plate
      do i=1,nxt
      x=xs+length*cos(dip)-length*cos(dip)*float(i-1)/float(nxt)
      y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
      z=zfree-length*sin(dip)+length*sin(dip)*float(i-1)/float(nxt)
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
      enddo
! no plate
      do i=nxs,nx
      x=xs+(1.-xs)*float(i-nxs)/float(nx-nxs)
      y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
      z=zfree
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (i.ne.nx.and.j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
      enddo
      enddo

      do j=ny2,ny
! no plate
      y=yl2+float(j-ny2)/float(ny-ny2)*(1.-yl2)
      do i=1,nx1
      x=float(i-1)/float(nx1-1)*xl1
      z=zfree
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (j.ne.ny) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
      enddo
      do i=nx1+1,nxs
      x=xl1+(xsp-xl1)*float(i-nx1)/float(nxs-nx1)
      z=zfree
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (j.ne.ny) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
      enddo
      do i=nxs,nx
      x=xs+(1.-xs)*float(i-nxs)/float(nx-nxs)
      z=zfree
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
      if (i.ne.nx.and.j.ne.ny) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
      enddo
      enddo

      print*,nnode,nelem

      i1=1
  111 continue
      do i=i1,nnode
      do j=i+1,nnode
      dist=(xp(i)-xp(j))**2+(yp(i)-yp(j))**2+(zp(i)-zp(j))**2
      if (dist.lt.1.d-10) then
              print*,'removing ',j
        do ij=j+1,nnode
        xp(ij-1)=xp(ij)
        yp(ij-1)=yp(ij)
        zp(ij-1)=zp(ij)
        enddo
        do ij=1,nelem
        do k=1,3
        if (icon(k,ij).eq.j) icon(k,ij)=i
        enddo
        enddo
        do ij=1,nelem
        do k=1,3
        if (icon(k,ij).gt.j) icon(k,ij)=icon(k,ij)-1
        enddo
        enddo
      nnode=nnode-1
      i1=i
      goto 111
      endif
      enddo
      enddo

      print*,nnode,nelem
      print*,minval(icon(:,1:nelem)),maxval(icon(:,1:nelem))

!add normals

      do i=1,nelem
      i1=icon(1,i)
      i2=icon(2,i)
      i3=icon(3,i)
      if (i1.eq.i2 .and. i2.eq.i3) print*,'element ',i,' collapsed: ',icon(:,i)
      x1=xp(i1)
      x2=xp(i2)
      x3=xp(i3)
      y1=yp(i1)
      y2=yp(i2)
      y3=yp(i3)
      z1=zp(i1) 
      z2=zp(i2) 
      z3=zp(i3) 
      xne=(y2-y1)*(z3-z1)-(y3-y1)*(z2-z1)
      yne=(z2-z1)*(x3-x1)-(z3-z1)*(x2-x1)
      zne=(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)
      xyzn=sqrt(xne**2+yne**2+zne**2)
      xne=xne/xyzn
      yne=yne/xyzn
      zne=zne/xyzn
        do j=1,3
        ij=icon(j,i)
        xn(ij)=xn(ij)+xne
        yn(ij)=yn(ij)+yne
        zn(ij)=zn(ij)+zne
        enddo
      enddo

      do i=1,nnode
      xyzn=sqrt(xn(i)**2+yn(i)**2+zn(i)**2)
      xn(i)=xn(i)/xyzn
      yn(i)=yn(i)/xyzn
      zn(i)=zn(i)/xyzn
      enddo

      open (71,file='surface-1.txt',status='unknown')
      write (71,*) nnode,nelem
        do i=1,nnode
        write (71,*) xp(i),yp(i),zp(i),xn(i),yn(i),zn(i)
        enddo
        do i=1,nelem
        write (71,*) (icon(j,i),j=1,3)
        enddo
      close (71)

      open(unit=30,file='Mesh.vtk')
      write(30,'(a)')'# vtk DataFile Version 3.0'
      write(30,'(a)')'surface'
      write(30,'(a)')'ASCII'
      write(30,'(a)')'DATASET UNSTRUCTURED_GRID'
      write(30,'(a7,i10,a6)')'POINTS ',nnode,' float'
      do i=1,nnode
      write(30,'(3f16.11)') xp(i),yp(i),zp(i)
      enddo
      write(30,'(A6, 2I10)') 'CELLS ',nelem,4*nelem
      do i=1,nelem
      write(30,'(9I10)')3,icon(1:3,i)-1
      enddo
      write(30,'(A11, I10)') 'CELL_TYPES ',nelem
      do i=1,nelem
      write(30,'(I2)')5 ! octree  (8 nodes)
      enddo
      write(30,'(a11,i10)')'POINT_DATA ',nnode
      write(30,'(a)')'VECTORS norm float'
      do i=1,nnode
      write(30,'(3f18.13)')xn(i),yn(i),zn(i)
      enddo
      close(30)

      deallocate (xp,yp,zp,icon)

      end

      !--------

      subroutine add_point (nnode,x,y,z,xp,yp,zp,nnodemax)

      double precision xp(nnodemax),yp(nnodemax),zp(nnodemax)
      double precision x,y,z

      nnode=nnode+1
      if (nnode.gt.nnodemax) then
      print*,nnode,nnodemax
      stop 'needs to increase nnodemax'
      endif

      xp(nnode)=x
      yp(nnode)=y
      zp(nnode)=z

      return
      end

      !----------

      subroutine add_element (nelem,icon,i1,i2,i3,i4,nelemmax)

      integer nelem,nelemmax,icon(3,nelemmax),i1,i2,i3,i4

      nelem=nelem+1
      if (nelem.gt.nelemmax) then
      print*,nelem,nelemmax
      stop 'nelemmax needs to be increased'
      endif
      icon(1,nelem)=i1
      icon(2,nelem)=i2
      icon(3,nelem)=i3

      nelem=nelem+1
      if (nelem.gt.nelemmax) then
      print*,nelem,nelemmax
      stop 'nelemmax needs to be increased'
      endif
      icon(1,nelem)=i1
      icon(2,nelem)=i3
      icon(3,nelem)=i4

      return
      end
