Added sp11 and sp12 as input options for the defining the surface geometry

src/create_surfaces.f90

@@ -50,7 +50,7 @@ integer,dimension(:),  allocatable :: hulltriangles,vis_tlist,vis_elist
 integer,dimension(:),  allocatable :: add_tlist,nb
 integer,dimension(:,:),allocatable :: vertices,neighbour,nn
 logical*1,dimension(:),allocatable :: lt_work,ln_work
-double precision :: sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10
+double precision :: sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12
 double precision :: x1,x2,x3,y1,y2,y3,z1,z2,z3,xne,yne,zne,xyzn
 double precision :: delta,eps,epsil,xtemp,ytemp,delta_angle,delta_radius,angle
 double precision :: xx,yy,xini,xend,yini,yend,dx,dy,dist
@@ -133,6 +133,8 @@ else
    sp08=surface%sp08
    sp09=surface%sp09
    sp10=surface%sp10
+   sp11=surface%sp11
+   sp12=surface%sp12
    surface_type=surface%surface_type
    levelt=surface%levelt
 
@@ -430,7 +432,7 @@ else
    !----------------------------------------------|
 
    call zpoints (surface%nsurface,surface%x,surface%y,surface%z,surface_type, &
-                 sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10)
+                 sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12)
 
    !----------------------------------------------|
    !----computing the normals---------------------|  
@@ -452,7 +454,7 @@ end subroutine create_surf
 !------------------------------------------------------------------------------|
 
 subroutine zpoints (ns,x,y,z,surface_type, &
-                    sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10)
+                    sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12)
 
 !------------------------------------------------------------------------------|
 !(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
@@ -462,7 +464,7 @@ subroutine zpoints (ns,x,y,z,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..sp10 are the surface parameters (not all are necessarily used)
+! sp01..sp12 are the surface parameters (not all are necessarily used)
 !------------------------------------------------------------------------------|
 !((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
 !------------------------------------------------------------------------------|
@@ -474,7 +476,7 @@ 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
+double precision sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12
 
 !------------------------------------------------------------------------------|
 !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
@@ -665,6 +667,7 @@ select case(surface_type)
      ! 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 dip angle of the outer 'flat' region
      psi=sp03*pi/180.d0
      theta1=sp08*pi/180.d0
      theta2=sp10*pi/180.d0

src/module_definitions.f90

@@ -27,8 +27,8 @@ module definitions
       integer nsurface,nt
       double precision,dimension(:),pointer::r,s,x,y,z,xn,yn,zn,u,v,w
       integer levelt,itype,material,surface_type
-      double precision :: sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10
-      double precision activation_time
+      double precision :: sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11
+      double precision :: sp12,activation_time
       integer,dimension(:,:),pointer::icon
       logical rand
       integer leveloct
@@ -330,7 +330,7 @@ module definitions
 
       end interface
 
-!=====[ZISODISP================================================================
+!=====[ZISODISP]================================================================
 ! type ziso is used to store the isostatic deflection of the basal boundary
 ! due to isostasy
 

src/read_input_file.f90

@@ -422,6 +422,10 @@ do i=1,params%ns
    if (iproc==0) call scanfile (params%infile,'surface_param_09_'//cm(il:3),surface(i)%sp09,ires)
    surface(i)%sp10 = 0.d0
    if (iproc==0) call scanfile (params%infile,'surface_param_10_'//cm(il:3),surface(i)%sp10,ires)
+   surface(i)%sp11 = 0.d0
+   if (iproc==0) call scanfile (params%infile,'surface_param_11_'//cm(il:3),surface(i)%sp11,ires)
+   surface(i)%sp12 = 0.d0
+   if (iproc==0) call scanfile (params%infile,'surface_param_12_'//cm(il:3),surface(i)%sp12,ires)
    
    call mpi_bcast(surface(i)%sp01,1,mpi_double_precision,0,mpi_comm_world,ierr)
    call mpi_bcast(surface(i)%sp02,1,mpi_double_precision,0,mpi_comm_world,ierr)
@@ -433,6 +437,8 @@ do i=1,params%ns
    call mpi_bcast(surface(i)%sp08,1,mpi_double_precision,0,mpi_comm_world,ierr)
    call mpi_bcast(surface(i)%sp09,1,mpi_double_precision,0,mpi_comm_world,ierr)
    call mpi_bcast(surface(i)%sp10,1,mpi_double_precision,0,mpi_comm_world,ierr)
+   call mpi_bcast(surface(i)%sp11,1,mpi_double_precision,0,mpi_comm_world,ierr)
+   call mpi_bcast(surface(i)%sp12,1,mpi_double_precision,0,mpi_comm_world,ierr)
 
    write(threadinfo%Logunit,*) 'surface ',i,'sp01',surface(i)%sp01
    write(threadinfo%Logunit,*) 'surface ',i,'sp02',surface(i)%sp02
@@ -444,6 +450,8 @@ do i=1,params%ns
    write(threadinfo%Logunit,*) 'surface ',i,'sp08',surface(i)%sp08
    write(threadinfo%Logunit,*) 'surface ',i,'sp09',surface(i)%sp09
    write(threadinfo%Logunit,*) 'surface ',i,'sp10',surface(i)%sp10
+   write(threadinfo%Logunit,*) 'surface ',i,'sp11',surface(i)%sp11
+   write(threadinfo%Logunit,*) 'surface ',i,'sp12',surface(i)%sp12
 
    surface(i)%activation_time=-1.d0
    if (iproc==0) call scanfile (params%infile,'activation_time_'//cm(il:3),surface(i)%activation_time,ires)