!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              CREATE_SURF    May. 2008                                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine create_surf(surface,is,debug)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! for a given surfaces passed as argument, this routine computes
! the 8 fields x,y,z,xn,yn,zn,r,s for each point of the surface
! according to its type and levelt, and the triangulation

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use constants
use definitions
!use mpi
use randomodule

implicit none

include 'mpif.h'

type (sheet) :: surface
integer is
integer debug

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer :: levelt,surface_type,nproc,iproc,ierr,i,j,seed,indix
integer :: err,ntmax,nhmax,npmax,nmax,nmode,nvmax,nnpnmax,nh,nohalt_hull,loc
integer :: dmode,nt,nhull,icircles,irect
integer,dimension(:),  allocatable :: hulltriangles,vis_tlist,vis_elist
integer,dimension(:),  allocatable :: add_tlist
integer,dimension(:,:),allocatable :: vertices,neighbour
logical*1,dimension(:),allocatable :: lt_work,ln_work
double precision :: sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12
double precision :: delta,eps,epsil,delta_angle,delta_radius,angle
double precision :: xx,yy,xini,xend,yini,yend,dx,dy,dist
double precision,dimension(:),  allocatable :: field
double precision,dimension(:,:),allocatable :: points,centres
logical clockwise
character ch
character*72  :: shift
integer nnn(1),nnlist(1),ntrilist(1)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

shift=' '

if (surface%surface_type.lt.0) then
   !---------------------------------
   ! reading surface from ascii file
   !---------------------------------
   if (iproc.eq.0) then
      write (ch,'(i1)') -surface%surface_type
      write(*,*) shift//'surface-'//ch//'.txt'
      open(71,file='surface-'//ch//'.txt',status='old')
      read(71,*) surface%nsurface,surface%nt
      write(*,*) shift//'-> nb of points   :',surface%nsurface
      write(*,*) shift//'-> nb of triangles:',surface%nt
   end if
   call mpi_bcast (surface%nsurface,1,mpi_integer,0,mpi_comm_world,ierr)
   call mpi_bcast (surface%nt,1,mpi_integer,0,mpi_comm_world,ierr)
   deallocate (surface%x,surface%y,surface%z)
   deallocate (surface%xn,surface%yn,surface%zn)
   deallocate (surface%r,surface%s)
   deallocate (surface%icon)
   allocate (surface%x(surface%nsurface),surface%y(surface%nsurface),surface%z(surface%nsurface))
   allocate (surface%xn(surface%nsurface),surface%yn(surface%nsurface),surface%zn(surface%nsurface))
   allocate (surface%r(surface%nsurface),surface%s(surface%nsurface))
   allocate (surface%icon(3,surface%nt))
   if (iproc.eq.0) then
      write(*,*) surface%nsurface,surface%nt
      do i=1,surface%nsurface
         read (71,*) surface%x(i),surface%y(i),surface%z(i),surface%xn(i),surface%yn(i),surface%zn(i)
         surface%r(i)=surface%x(i)
         surface%s(i)=surface%y(i)
         if (surface%x(i)<0 .or. surface%x(i)>1 .or. surface%y(i)<0 .or. surface%y(i)>1 .or. surface%z(i)<0 .or. surface%z(i)>1 &
              .and. iproc.eq.0)  &
            print *,surface%x(i),surface%y(i),surface%z(i) 
      enddo
      do i=1,surface%nt
         read (71,*) (surface%icon(j,i),j=1,3)
      enddo
      close (71)
   endif
   call mpi_bcast(surface%x,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%y,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%z,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%xn,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%yn,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%zn,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%r,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%s,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%icon,surface%nt*3,mpi_integer,0,mpi_comm_world,ierr)
else
   !----------------
   ! create surface
   !----------------

   seed=1234567
   sp01=surface%sp01
   sp02=surface%sp02
   sp03=surface%sp03
   sp04=surface%sp04
   sp05=surface%sp05
   sp06=surface%sp06
   sp07=surface%sp07
   sp08=surface%sp08
   sp09=surface%sp09
   sp10=surface%sp10
   sp11=surface%sp11
   sp12=surface%sp12
   surface_type=surface%surface_type
   levelt=surface%levelt

!   delta=2.d0**(-levelt) ! spacing between points on hull

   delta=1.d0 / (2.d0**levelt +1) ! opla

   epsil=(surface%stretch-1)/2.d0

   if (iproc.eq.0) write(8,*) 'create_surfaces: delta=',delta
   if (iproc.eq.0) write(8,*) 'create_surfaces: epsilon=',epsil

   !-----building the surface-------------
   ! type 1  : flat surface
   ! type 2  : rectangular emboss
   ! type 3  : convex spherical emboss 
   ! type 4  : concave spherical emboss
   ! type 5  : double rectangular emboss
   ! type 6  : sinus 
   ! type 7  : noisy surface 
   ! type 8  : double sinus 
   ! type 9  : cosinus sinus
   ! type 10 : slope
   ! type 11 : 2D hill
   ! type 12 : rectangular emboss with specified slope
   ! type 13 : slope with 2 kinks

   select case(surface_type)

      case (1,2,5,6,7,8,9,10,11,12,13) 
!           if (surface%rand) then
              indix=1
              ! convex hull
              do i=1,2**levelt   +1
                 surface%x(indix)=dble(i-1)*delta
                 surface%y(indix)=0.d0
                 indix=indix+1
              end do
              do i=1,2**levelt   +1
                 surface%x(indix)=1.d0
                 surface%y(indix)=dble(i-1)*delta
                 indix=indix+1
              end do
              do i=1,2**levelt    +1
                 surface%x(indix)=1.d0-dble(i-1)*delta
                 surface%y(indix)=1.d0
                 indix=indix+1
              end do
              do i=1,2**levelt   +1
                 surface%x(indix)=0.d0
                 surface%y(indix)=1.d0-dble(i-1)*delta
                 indix=indix+1
              end do
              !interior nodes
              if(surface%rand) then
                 call random_number (surface%x(indix:surface%nsurface))
                 call random_number (surface%y(indix:surface%nsurface))
              else
                 surface%x(indix:surface%nsurface)=0.5d0
                 surface%y(indix:surface%nsurface)=0.5d0
              end if
              do j=2,2**levelt +1
                 do i=2,2**levelt  +1
                    surface%x(indix)=(((surface%x(indix)-0.5d0)*2.d0)*epsil+dble(i-1))*delta
                    surface%y(indix)=(((surface%y(indix)-0.5d0)*2.d0)*epsil+dble(j-1))*delta
                    indix=indix+1
                 enddo
              enddo
!           else
!              ! regular grid of points
!              indix=1
!              do j=1,2**levelt+1
!                 do i=1,2**levelt+1
!                    surface%x(indix)=float(i-1)*delta
!                    surface%y(indix)=float(j-1)*delta
!                    indix=indix+1 
!                 end do
!              end do
!           end if

      case (22)
           if (surface%rand) then
              indix=1
              ! convex hull
              do i=1,2**levelt
                 surface%x(indix)=(i-1)*delta
                 surface%y(indix)=0.d0
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=1.d0
                 surface%y(indix)=(i-1)*delta
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=1.d0-(i-1)*delta
                 surface%y(indix)=1.d0
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=0.d0
                 surface%y(indix)=1.d0-(i-1)*delta
                 indix=indix+1
              end do
              if (iproc.eq.0) write(8,*) 'nhull=',indix-1,'for surface ',is
              ! placing rows of points outside, on and outside the rectangle
              eps=0.001
              do irect=-5,3,2
                 xini=sp02+irect*eps
                 xend=sp03-irect*eps
                 yini=sp04+irect*eps
                 yend=sp05-irect*eps
                 dx=(xend-xini)/2**levelt
                 dy=(yend-yini)/2**levelt
                 do i=1,2**levelt
                    xx=xini+(i-1)*dx
                    yy=yini
                    if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
                       surface%x(indix)=xx
                       surface%y(indix)=yy
                       indix=indix+1
                    end if
                 end do
                 do i=1,2**levelt
                    xx=xend
                    yy=yini+(i-1)*dy
                    if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
                       surface%x(indix)=xx
                       surface%y(indix)=yy
                       indix=indix+1
                    end if
                 end do
                 do i=1,2**levelt
                    xx=xend-(i-1)*dx
                    yy=yend
                    if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
                       surface%x(indix)=xx
                       surface%y(indix)=yy
                       indix=indix+1
                    end if
                 end do
                 do i=1,2**levelt
                    xx=xini
                    yy=yend-(i-1)*dy
                    if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
                       surface%x(indix)=xx
                       surface%y(indix)=yy
                       indix=indix+1
                    end if 
                 end do
              end do
              ! throwing remaining points at random
              do while (indix.le.surface%nsurface)
                 xx=ran(seed)
                 yy=ran(seed)
                 if ((xx<sp02-10*eps .or. xx>sp03+10*eps .or. yy<sp04-10*eps .or. yy>sp05+10*eps) .or. &
                     (xx>sp02+10*eps .and. xx<sp03-10*eps .and. yy>sp04+10*eps .and. yy<sp05-10*eps)) then 
                    surface%x(indix)=xx
                    surface%y(indix)=yy
                    indix=indix+1
                 end if
              end do
           else
              ! regular grid of points
              indix=1
              do j=1,2**levelt+1
                 do i=1,2**levelt+1
                 surface%x(indix)=float(i-1)*delta
                 surface%y(indix)=float(j-1)*delta
                 indix=indix+1 
                 end do
              end do
           end if   

      case (3,4) 
           if (surface%rand) then
              ! convex hull
              indix=1
              do i=1,2**levelt
                 surface%x(indix)=(i-1)*delta
                 surface%y(indix)=0.d0
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=1.d0
                 surface%y(indix)=(i-1)*delta
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=1.d0-(i-1)*delta
                 surface%y(indix)=1.d0
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=0.d0
                 surface%y(indix)=1.d0-(i-1)*delta
                 indix=indix+1
              end do
              nhull=indix-1 
              if (iproc.eq.0) write(8,*) 'nhull=',nhull,'on surface ',is
              ! placing rows of nhull points outside, in, and inside the circle
              delta_angle=2.d0*pi/nhull
              delta_radius=0.0025d0*sp04
              do icircles=-4,2
                 do i=0,nhull-1
                    angle=i*delta_angle
                    surface%x(indix)=sp02+(sp04+icircles*delta_radius)*cos(angle)
                    surface%y(indix)=sp03+(sp04+icircles*delta_radius)*sin(angle)
                    indix=indix+1
                 end do
              end do
              ! throwing the rest of the points at random provided they are not
              ! too close to the concentric circles so that it does not disturb
              ! the regular triangulation
              do while (indix.le.surface%nsurface)
                 xx=ran(seed)
                 yy=ran(seed)
                 dist=sqrt((xx-sp02)**2+(yy-sp03)**2)
                 if (dist>(sp04+10*delta_radius) .or. dist<(sp04-10*delta_radius)) then
                    surface%x(indix)=xx
                    surface%y(indix)=yy
                    indix=indix+1
                 end if
              end do
           else
              ! regular grid of points
              indix=1
              do j=1,2**levelt+1
                 do i=1,2**levelt+1
                    surface%x(indix)=float(i-1)*delta
                    surface%y(indix)=float(j-1)*delta
                    indix=indix+1 
                 end do
              end do
           end if

      case default
           call stop_run ('surface type not defined$')
   end select

   !-------------------------------------------|
   !-----building icon array-------------------|
   !-------------------------------------------|
   ! we use the NN library to compute the triangulation 
   ! of points for the surface
   ntmax=surface%nsurface*3
   nhmax=surface%nsurface
   npmax=surface%nsurface
   nnpnmax=100
   nmax=3*ntmax+npmax
   nvmax=ntmax
   allocate (points(2,surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc points in create_surface$')
   allocate (field(surface%nsurface),stat=err)    ; if (err.ne.0) call stop_run ('Error alloc field in create_surface$')
   allocate (vertices(3,ntmax),stat=err)          ; if (err.ne.0) call stop_run ('Error alloc vertices in create_surface$')
   allocate (centres(3,ntmax),stat=err)           ; if (err.ne.0) call stop_run ('Error alloc centress in create_surface$')
   allocate (neighbour(3,ntmax),stat=err)         ; if (err.ne.0) call stop_run ('Error alloc neighbour in create_surface$')
   allocate (hulltriangles(nhmax),stat=err)       ; if (err.ne.0) call stop_run ('Error alloc hulltriangles in create_surface$')
   allocate (vis_tlist(nvmax),stat=err)           ; if (err.ne.0) call stop_run ('Error alloc vis_tlist in create_surface$')
   allocate (vis_elist(nvmax),stat=err)           ; if (err.ne.0) call stop_run ('Error alloc vis_elist in create_surface$')
   allocate (add_tlist(nvmax),stat=err)           ; if (err.ne.0) call stop_run ('Error alloc add_tlist in icreate_surface$')
   allocate (lt_work(ntmax),stat=err)             ; if (err.ne.0) call stop_run ('Error alloc lt_work in create_surface$')
   allocate (ln_work(surface%nsurface),stat=err)  ; if (err.ne.0) call stop_run ('Error alloc ln_work in create_surface$')
   points(1,:)=surface%x
   points(2,:)=surface%y
   dmode=-2
   nmode=0
   loc=1
   nohalt_hull=0
   clockwise=.true.
   eps=1.d-8
   call nn2d_setup (surface%nsurface,ntmax,nhmax,npmax,nnpnmax,nmax, &
                    points,dmode,nmode,clockwise,field,nt,vertices, &
                    centres,neighbour,nh,hulltriangles,nohalt_hull, &
                    loc,nnn,nnlist,ntrilist, &
                    eps,nvmax,vis_tlist,vis_elist,add_tlist, &
                    lt_work,ln_work)
   if (iproc.eq.0) write (8,*) 'surface ',is,' has ',nt,' triangles'

   surface%icon(1:3,1:surface%nt)=vertices(1:3,1:nt)

   deallocate (points)
   deallocate (field)
   deallocate (vertices)
   deallocate (centres)
   deallocate (neighbour)
   deallocate (hulltriangles)
   deallocate (vis_tlist)
   deallocate (vis_elist)
   deallocate (add_tlist)
   deallocate (lt_work)
   deallocate (ln_work)

   !----------------------------------------------|
   !----computing z coordinate--------------------|
   !----------------------------------------------|

   call zpoints (surface%nsurface,surface%x,surface%y,surface%z,surface_type, &
                 sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12)

   !----------------------------------------------|
   !----computing the normals---------------------|  
   !----------------------------------------------|

   call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
                         surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)

   if (debug>=2) call output_surf (surface,is,'init',-1,-1) 

end if

end subroutine create_surf




!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine zpoints (ns,x,y,z,surface_type, &
                    sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! knowing the two-dimensional spatial distribution of points on the plane (x,y)
! this routine computes the z-coordinate accordingly to the surface type.
! ns is the number of points
! x,y,z are the coordinates of the points
! surface_type is the type of surface under consideration
! sp01..sp12 are the surface parameters (not all are necessarily used)
!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use constants

implicit none

integer ns
double precision x(ns),y(ns),z(ns)
integer surface_type
double precision sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer i
double precision dist,ran,m,x1a,x1b,x2a,x2b,y1a,y1b,y2a,y2b,zx
double precision psi,theta1,theta2,theta1s,theta2s,psis,y0s,y1s,y2s,x0,y0w
double precision y0a,y0b,y0,y1,y2,a1,b1,bb1,c1,m1,a2,b2,bb2,c2,m2,a3,b3,bb3,c3
double precision m3,a4,b4,bb4,c4,m4,eps,psi2,psi2s,dz1,dz2

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

eps=1.d-10

select case(surface_type)
   case (1)
        ! a flat surface, 
        ! sp01 is the z level
        z=sp01
   case (2)
        ! rectangular emboss, 
        ! sp01 is the z level
        ! sp02 and 03 are x1,x2 
        ! sp04 and 05 are y1,y2 
        ! sp06 is the thickness 
        do i=1,ns
           if ((x(i)-sp02)*(x(i)-sp03).le.0.d0 .and. (y(i)-sp04)*(y(i)-sp05).le.0.d0) then
              z(i)=sp01-sp06
           else
              z(i)=sp01
           end if
        end do
   case (3)
        ! convex spherical emboss, 
        ! sp01 is the z level
        ! sp02 and 03 are x0,y0 
        ! sp04 is the radius 
        do i=1,ns
           dist=sqrt((x(i)-sp02)**2+(y(i)-sp03)**2)
           if (dist.le.sp04) then
              z(i)=sp01+sp04*sin(acos(dist/sp04))
           else
              z(i)=sp01
           endif
       end do
   case (4)
        ! convex spherical emboss, 
        ! sp01 is the z level
        ! sp02 and 03 are x0,y0 
        ! sp04 is the radius 
        do i=1,ns
           dist=sqrt((x(i)-sp02)**2+(y(i)-sp03)**2)
           if (dist.le.sp04) then
              z(i)=sp01-sp04*sin(acos(dist/sp04))
           else
              z(i)=sp01
           endif
       end do
   case (5)
        ! double rectangular emboss, 
        ! sp01 is the z level
        ! sp02 and 03 are x1,x2 
        ! sp04 and 05 are x3,x4 
        ! sp06 and 07 are y1,y2 
        ! sp08 and 09 are y3,y4 
        ! sp10 is the thickness 
         z=sp01
         do i=1,ns
            if ((x(i)-sp02)*(x(i)-sp03).le.0.d0 .and. (y(i)-sp06)*(y(i)-sp07).le.0.d0) then
               z(i)=z(i)-sp10
            endif
            if ((x(i)-sp04)*(x(i)-sp05).le.0.d0 .and. (y(i)-sp08)*(y(i)-sp09).le.0.d0) then
              z(i)=z(i)-sp10
           endif
        end do
   case (6)
        ! a sinus, 
        ! sp01 is the z level
        ! sp02 is the wavelength 
        ! sp03 is the amplitude
        ! sp04 is the phase
        do i=1,ns
           z(i)=sp01+sp03*sin(2.d0*pi*(x(i)/sp02+sp04))
        end do
   case (7)
        ! a noisy surface, 
        ! sp01 is the z level
        ! sp02 is the noise amplitude
        do i=1,ns
           call random_number (ran)
           z(i)=sp01+ran*sp02
        end do
   case (8)
        ! a double sinus.
        ! sp01 is the z level
        ! sp02 is the x-wavelength 
        ! sp03 is the x-amplitude 
        ! sp04 is the y-wavelength 
        ! sp05 is the y-amplitude
        do i=1,ns 
           z(i)=sp01+sp03*sin(x(i)*2.*pi/sp02)+sp05*sin(y(i)*2.*pi/sp04)
        end do
   case (9)
        ! a cosinus, 
        ! sp01 is the z level
        ! sp02 is the wavelength 
        ! sp03 is the amplitude
        ! sp04 is the phase
        do i=1,ns
           z(i)=sp01+sp03*cos(2.d0*pi*(x(i)/sp02+sp04))
        end do
   case (10)
        ! a 2D embankment
        ! sp01 is z0
        ! sp02 is y0
        ! sp03 is psi
        ! sp04 is the thickness
        do i=1,ns 
           if (y(i)<sp02) then
              z(i)=sp01
           else
              z(i)=sp01+min((y(i)-sp02)*tan(sp03*pi/180.d0),sp04)
           end if
        end do
   case (11)
       ! a 2D hill
       ! sp01 is z level
       ! sp02 height
       ! sp03 is y0
       ! sp04 is width
        do i=1,ns 
           if (y(i)>sp03-sp04/2.d0 .and. y(i)<sp03+sp04/2.d0) then
              z(i)=sp01+sp02*(sin(pi*((y(i)-sp03)/sp04+0.5d0)))**2
           else
              z(i)=sp01
           end if
        end do
   case (12)
        ! rectangular emboss with specified slope
        ! sp01 is the z level
        ! sp02 and 03 are x1,x2 
        ! sp04 and 05 are y1,y2 
        ! sp06 is the thickness
        ! sp07 is the slope
        m=tan(sp07*pi/180.d0)
        x1a=sp02+(sp06/(2.d0*m))
        x1b=sp02-(sp06/(2.d0*m))
        x2a=sp03+(sp06/(2.d0*m))
        x2b=sp03-(sp06/(2.d0*m))
        y1a=sp04+(sp06/(2.d0*m))
        y1b=sp04-(sp06/(2.d0*m))
        y2a=sp05+(sp06/(2.d0*m))
        y2b=sp05-(sp06/(2.d0*m))
        do i=1,ns
           if (x(i).le.x1a) then
             zx=sp01
           elseif (x(i).le.x1b) then
             zx=sp01-(((x(i)-x1a)/(-sp06/m))*sp06)
           elseif (x(i).le.x2a) then
             zx=sp01-sp06
           elseif (x(i).le.x2b) then
             zx=sp01-((1.d0-((x(i)-x2a)/(-sp06/m)))*sp06)
             !zx=sp01-((1.d0-((x(i)-x2a)/(sp06/m)))*sp06)
           else
             zx=sp01
           endif
           if (y(i).le.y1a) then
             z(i)=sp01
           elseif (y(i).le.y1b) then
             z(i)=min(sp01-(((y(i)-y1a)/(-sp06/m))*sp06),zx)
           elseif (y(i).le.y2a) then
             z(i)=min(sp01-sp06,zx)
           elseif (y(i).le.y2b) then
             z(i)=min(sp01-((1.d0-((y(i)-y2a)/(-sp06/m)))*sp06),zx)
           else
             z(i)=sp01
           endif           
        end do

   case (13)
     ! a 3D embankment of finite length and with 2 kinks
     ! sp01 is the z level
     ! sp02 is the z level 2
     ! sp03 is slope angle (psi)
     ! sp04 is x value for the center of the slope (x0)
     ! sp05 is the x width of the narrow end of the plateau (y0w)
     ! sp06 is the y position for the start of the slope (y0)
     ! sp07 is the y position of the first kink (y1)
     ! sp08 is the angle of the first kink (theta1)
     ! sp09 is the y position of the second kink (y2)
     ! sp10 is the angle of the second kink (theta2)
     ! sp11 is the reference elevation (elevation of the outer flat)
     ! sp12 is the dip angle of the outer dip panel
     dz1=sp02-sp01
     dz2=sp01-sp11
     psi=sp03*pi/180.d0
     psi2=sp12*pi/180.d0
     theta1=sp08*pi/180.d0
     theta2=sp10*pi/180.d0
     psis=sign(tan(psi),dz1)
     psi2s=sign(tan(psi2),dz2)
     theta1s=tan(theta1)
     theta2s=tan(theta2)
     y0s=1.d0
     y1s=tan(theta1/2.d0)
     y2s=tan((theta1+theta2)/2.d0)
     x0=sp04

     !y0w=sp05
     y0w=sp05+dz1/(2.d0*psis)
     y0=sp06
     !y0a=y0+sp02/(2.d0*psis)
     !y0b=y0-sp02/(2.d0*psis)
     y0a=y0+dz1/psis
     y0b=y0
     y1=sp07+(dz1/(2.d0*psis))*y1s
     y2=sp09+(dz1/(2.d0*psis))*y2s

     m1=((y2+(dz1/(2.d0*psis))*y1s)-(y1+(dz1/(2.d0*psis))*y1s))/&
        ((x0-y0w-(y1-y2)*theta1s-(dz1/(2.d0*psis)))-(x0-y0w-(dz1/(2.d0*psis))))
     bb1=(y1+(dz1/(2.d0*psis))*y1s)-m1*(x0-y0w-(dz1/(2.d0*psis)))
     a1=-m1
     b1=1
     c1=-bb1

     m2=((y2-(dz1/(2.d0*psis))*y1s)-(y1-(dz1/(2.d0*psis))*y1s))/&
        ((x0-y0w-(y1-y2)*theta1s+(dz1/(2.d0*psis)))-(x0-y0w+(dz1/(2.d0*psis))))
     bb2=(y1-(dz1/(2.d0*psis))*y1s)-m1*(x0-y0w+(dz1/(2.d0*psis)))
     a2=-m2
     b2=1
     c2=-bb2

     m3=((y2-(dz1/(2.d0*psis))*y1s)-(y1-(dz1/(2.d0*psis))*y1s))/&
        ((x0+y0w+(y1-y2)*theta1s-(dz1/(2.d0*psis)))-(x0+y0w-(dz1/(2.d0*psis))))
     bb3=(y1-(dz1/(2.d0*psis))*y1s)-m3*(x0+y0w-(dz1/(2.d0*psis)))
     a3=-m3
     b3=1
     c3=-bb3

     m4=((y2+(dz1/(2.d0*psis))*y1s)-(y1+(dz1/(2.d0*psis))*y1s))/&
        ((x0+y0w+(y1-y2)*theta1s+(dz1/(2.d0*psis)))-(x0+y0w+(dz1/(2.d0*psis))))
     bb4=(y1+(dz1/(2.d0*psis))*y1s)-m3*(x0+y0w+(dz1/(2.d0*psis)))
     a4=-m4
     b4=1
     c4=-bb4

     do i=1,ns
       z(i)=sp01
       if (y(i).ge.-y0s*x(i)+(y0+(x0-y0w)) .and. y(i).ge.y0s*x(i)+(y0-(x0+y0w))) then
         if (y(i).ge.y0a) then
           if (abs(dz2).gt.eps) then
             if ((dz2).gt.eps) then
               z(i)=max(sp01-psi2s*(y(i)-y0a),sp11)
             else
               z(i)=min(sp01-psi2s*(y(i)-y0a),sp11)
             endif
           endif
         elseif (y(i).ge.y0b) then
           z(i)=sp01+psis*(y0a-y(i))
         elseif (y(i).lt.y0b) then
           z(i)=sp02
         endif
       elseif (y(i).ge.-y1s*x(i)+(y1s*(x0-y0w)+y1) .and. y(i).ge.y1s*x(i)+(-y1s*(x0+y0w)+y1)) then
         if (x(i).lt.((x0-y0w)-(dz1/(2.d0*psis)))) then
           if (abs(dz2).gt.eps) then
             if ((dz2).gt.eps) then
               z(i)=max(sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i)),sp11)
             else
               z(i)=min(sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i)),sp11)
             endif
           endif
         elseif (x(i).ge.((x0-y0w)-(dz1/(2.d0*psis))) .and. x(i).lt.((x0-y0w)+(dz1/(2.d0*psis)))) then
           z(i)=sp01+psis*(x(i)-((x0-y0w)-(dz1/(2.d0*psis))))
         elseif (x(i).ge.((x0-y0w)+(dz1/(2.d0*psis))) .and. x(i).lt.((x0+y0w)-(dz1/(2.d0*psis)))) then
           z(i)=sp02
         elseif (x(i).ge.((x0+y0w)-(dz1/(2.d0*psis))) .and. x(i).lt.((x0+y0w)+(dz1/(2.d0*psis)))) then
           z(i)=sp01+psis*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))))
         else
           if (abs(dz2).gt.eps) then
             if ((dz2).gt.eps) then
               z(i)=max(sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis)))),sp11)
             else
               z(i)=min(sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis)))),sp11)
             endif
           endif
         endif
       elseif (y(i).ge.-y2s*x(i)+(y2s*(x0-y0w-(y1-y2)*theta1s)+y2) .and. y(i).ge.y2s*x(i)+(-y2s*(x0+y0w+(y1-y2)*theta1s)+y2)) then
         if (y(i).ge.m1*x(i)+bb1) then
           dist=abs(a1*x(i)+b1*y(i)+c1)/sqrt(a1**2.d0+b1**2.d0)
           if (abs(dz2).gt.eps) then
             if ((dz2).gt.eps) then
               z(i)=max(sp01-psi2s*dist,sp11)
             else
               z(i)=min(sp01-psi2s*dist,sp11)
             endif
           endif
         elseif (y(i).ge.m2*x(i)+bb2 .and. y(i).lt.m1*x(i)+bb1) then
           dist=abs(a1*x(i)+b1*y(i)+c1)/sqrt(a1**2.d0+b1**2.d0)
           z(i)=sp01+dist*psis
         elseif (y(i).lt.m2*x(i)+bb2 .and. y(i).lt.m3*x(i)+bb3) then
           z(i)=sp02
         elseif (y(i).ge.m3*x(i)+bb3 .and. y(i).lt.m4*x(i)+bb4) then
           dist=abs(a4*x(i)+b4*y(i)+c4)/sqrt(a4**2.d0+b4**2.d0)
           z(i)=sp01+dist*psis
         else
           dist=abs(a4*x(i)+b4*y(i)+c4)/sqrt(a4**2.d0+b4**2.d0)
           if (abs(dz2).gt.eps) then
             if ((dz2).gt.eps) then
               z(i)=max(sp01-psi2s*dist,sp11)
             else
               z(i)=min(sp01-psi2s*dist,sp11)
             endif
           endif
         endif
       elseif (x(i).lt.x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)) then
         if (abs(dz2).gt.eps) then
           if ((dz2).gt.eps) then
             z(i)=max(sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i)),sp11)
           else
             z(i)=min(sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i)),sp11)
           endif
         endif
       elseif (x(i).ge.((x0-y0w-(y1-y2)*theta1s)-(dz1/(2.d0*psis))) .and. x(i).lt.((x0-y0w-(y1-y2)*theta1s)+(dz1/(2.d0*psis)))) then
         z(i)=sp01+psis*(x(i)-((x0-y0w-(y1-y2)*theta1s)-(dz1/(2.d0*psis))))
       elseif (x(i).ge.((x0-y0w-(y1-y2)*theta1s)+(abs(dz1)/(2.d0*psis))) .and. x(i).lt.((x0+y0w+(y1-y2)*theta1s)-(abs(dz1)/(2.d0*psis)))) then
         z(i)=sp02
       elseif (x(i).ge.((x0+y0w+(y1-y2)*theta1s)-(abs(dz1)/(2.d0*psis))) .and. x(i).lt.((x0+y0w+(y1-y2)*theta1s)+(abs(dz1)/(2.d0*psis)))) then
         z(i)=sp01+psis*(x(i)-((x0+y0w+(y1-y2)*theta1s)+(dz1/(2.d0*psis))))
       else
         if (abs(dz2).gt.eps) then
           if ((dz2).gt.eps) then
             z(i)=max(sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))),sp11)
           else
             z(i)=min(sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))),sp11)
           endif
         endif
       endif
     enddo

   case default
         call stop_run ('surface type not defined$')
end select

return
end
