Skip to content
Snippets Groups Projects
Commit 8cbb4d96 authored by Dave Whipp's avatar Dave Whipp
Browse files

Added a second kink option for surface geometry 13

parent 1db72c22
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 96 additions and 34 deletions
src/create_surfaces.f90
+ 96
− 34
View file @ 8cbb4d96
...@@ -54,6 +54,7 @@ integer,dimension(:), allocatable :: add_tlist ...@@ -54,6 +54,7 @@ integer,dimension(:), allocatable :: add_tlist
integer,dimension(:,:),allocatable :: vertices,neighbour integer,dimension(:,:),allocatable :: vertices,neighbour
logical*1,dimension(:),allocatable :: lt_work,ln_work logical*1,dimension(:),allocatable :: lt_work,ln_work
double precision :: sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12 double precision :: sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12
double precision :: sp13,sp14
double precision :: delta,eps,epsil,delta_angle,delta_radius,angle double precision :: delta,eps,epsil,delta_angle,delta_radius,angle
double precision :: xx,yy,xini,xend,yini,yend,dx,dy,dist double precision :: xx,yy,xini,xend,yini,yend,dx,dy,dist
double precision,dimension(:), allocatable :: field double precision,dimension(:), allocatable :: field
...@@ -134,6 +135,8 @@ else ...@@ -134,6 +135,8 @@ else
sp10=surface%sp10 sp10=surface%sp10
sp11=surface%sp11 sp11=surface%sp11
sp12=surface%sp12 sp12=surface%sp12
sp13=surface%sp13
sp14=surface%sp14
surface_type=surface%surface_type surface_type=surface%surface_type
levelt=surface%levelt levelt=surface%levelt
...@@ -431,7 +434,8 @@ else ...@@ -431,7 +434,8 @@ else
!----------------------------------------------| !----------------------------------------------|
call zpoints (surface%nsurface,surface%x,surface%y,surface%z,surface_type, & call zpoints (surface%nsurface,surface%x,surface%y,surface%z,surface_type, &
sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12) sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12, &
sp13,sp14)
!----------------------------------------------| !----------------------------------------------|
!----computing the normals---------------------| !----computing the normals---------------------|
...@@ -452,8 +456,8 @@ end subroutine create_surf ...@@ -452,8 +456,8 @@ end subroutine create_surf
!------------------------------------------------------------------------------| !------------------------------------------------------------------------------|
!------------------------------------------------------------------------------| !------------------------------------------------------------------------------|
subroutine zpoints (ns,x,y,z,surface_type, & subroutine zpoints (ns,x,y,z,surface_type,sp01,sp02,sp03,sp04,sp05,sp06,sp07, &
sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12) sp08,sp09,sp10,sp11,sp12,sp13,sp14)
!------------------------------------------------------------------------------| !------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))) !(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
...@@ -463,7 +467,7 @@ subroutine zpoints (ns,x,y,z,surface_type, & ...@@ -463,7 +467,7 @@ subroutine zpoints (ns,x,y,z,surface_type, &
! ns is the number of points ! ns is the number of points
! x,y,z are the coordinates of the points ! x,y,z are the coordinates of the points
! surface_type is the type of surface under consideration ! surface_type is the type of surface under consideration
! sp01..sp12 are the surface parameters (not all are necessarily used) ! sp01..sp14 are the surface parameters (not all are necessarily used)
!------------------------------------------------------------------------------| !------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments )))))))))))))))))))) !(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------| !------------------------------------------------------------------------------|
...@@ -476,6 +480,7 @@ integer ns ...@@ -476,6 +480,7 @@ integer ns
double precision x(ns),y(ns),z(ns) double precision x(ns),y(ns),z(ns)
integer surface_type integer surface_type
double precision sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12 double precision sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12
double precision sp13,sp14
!------------------------------------------------------------------------------| !------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables ))))))))))))) !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
...@@ -669,16 +674,24 @@ select case(surface_type) ...@@ -669,16 +674,24 @@ select case(surface_type)
! sp08 is the angle of the first kink (theta1) ! sp08 is the angle of the first kink (theta1)
! sp09 is the y position of the second kink (y2) ! sp09 is the y position of the second kink (y2)
! sp10 is the angle of the second kink (theta2) ! sp10 is the angle of the second kink (theta2)
! sp11 is the reference elevation (elevation of the outer flat) ! sp11 is the reference elevation (elevation of the lower outer dip change)
! sp12 is the dip angle of the outer dip panel ! sp12 is the dip angle of the lower outer dip panel
! sp13 is the reference elevation (elevation of the upper outer flat)
! sp14 is the dip angle of the upper outer dip panel
dz1=sp02-sp01 dz1=sp02-sp01
!dz2=sp01-sp11
dz2=sp01-sp11 dz2=sp01-sp11
dz3=sp01-sp13
dz4=dz3-dz2
psi=sp03*pi/180.d0 psi=sp03*pi/180.d0
!psi2=sp12*pi/180.d0
psi2=sp12*pi/180.d0 psi2=sp12*pi/180.d0
psi3=sp14*pi/180.d0
theta1=sp08*pi/180.d0 theta1=sp08*pi/180.d0
theta2=sp10*pi/180.d0 theta2=sp10*pi/180.d0
psis=sign(tan(psi),dz1) psis=sign(tan(psi),dz1)
psi2s=sign(tan(psi2),dz2) psi2s=sign(tan(psi2),dz2)
psi3s=sign(tan(psi3,dz3)
theta1s=tan(theta1) theta1s=tan(theta1)
theta2s=tan(theta2) theta2s=tan(theta2)
y0s=1.d0 y0s=1.d0
...@@ -725,13 +738,20 @@ select case(surface_type) ...@@ -725,13 +738,20 @@ select case(surface_type)
z(i)=sp01 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.-y0s*x(i)+(y0+(x0-y0w)) .and. y(i).ge.y0s*x(i)+(y0-(x0+y0w))) then
if (y(i).ge.y0a) then if (y(i).ge.y0a) then
if (abs(dz2).gt.eps) then if (y(i).ge.y0a+(dz2/psi2s)) then
if ((dz2).gt.eps) then if (dz4.gt.eps) then
z(i)=max(sp01-psi2s*(y(i)-y0a),sp11) z(i)=max(sp01-dz2-psi3s*(y(i)-(y0a+(dz2/psi2s))),sp13)
else else
z(i)=min(sp01-psi2s*(y(i)-y0a),sp11) z(i)=min(sp01-dz2-psi3s*(y(i)-(y0a+(dz2/psi2s))),sp13)
endif endif
else
z(i)=sp01-psi2s*(y(i)-y0a)
endif endif
! 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
elseif (y(i).ge.y0b) then elseif (y(i).ge.y0b) then
z(i)=sp01+psis*(y0a-y(i)) z(i)=sp01+psis*(y0a-y(i))
elseif (y(i).lt.y0b) then elseif (y(i).lt.y0b) then
...@@ -739,13 +759,20 @@ select case(surface_type) ...@@ -739,13 +759,20 @@ select case(surface_type)
endif endif
elseif (y(i).ge.-y1s*x(i)+(y1s*(x0-y0w)+y1) .and. y(i).ge.y1s*x(i)+(-y1s*(x0+y0w)+y1)) then 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 (x(i).lt.((x0-y0w)-(dz1/(2.d0*psis)))) then
if (abs(dz2).gt.eps) then if (x(i).lt.((x0-y0w)-(dz1/(2.d0*psis)))-(dz2/psi2s)) then
if ((dz2).gt.eps) then if (dz4.gt.eps) then
z(i)=max(sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i)),sp11) z(i)=max(sp01-dz2-psi3s*(((x0-y0w)-(dz1/(2.d0*psis))-(dz2/psi2s))-x(i)),sp13)
else else
z(i)=min(sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i)),sp11) z(i)=min(sp01-dz2-psi3s*(((x0-y0w)-(dz1/(2.d0*psis))-(dz2/psi2s))-x(i)),sp13)
endif endif
else
z(i)=sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i))
endif endif
! 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
elseif (x(i).ge.((x0-y0w)-(dz1/(2.d0*psis))) .and. x(i).lt.((x0-y0w)+(dz1/(2.d0*psis)))) then 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)))) 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 elseif (x(i).ge.((x0-y0w)+(dz1/(2.d0*psis))) .and. x(i).lt.((x0+y0w)-(dz1/(2.d0*psis)))) then
...@@ -753,24 +780,38 @@ select case(surface_type) ...@@ -753,24 +780,38 @@ select case(surface_type)
elseif (x(i).ge.((x0+y0w)-(dz1/(2.d0*psis))) .and. x(i).lt.((x0+y0w)+(dz1/(2.d0*psis)))) then 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)))) z(i)=sp01+psis*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))))
else else
if (abs(dz2).gt.eps) then if (x(i).lt.((x0+y0w)+(dz1/(2.d0*psis)))+(dz2/psi2s)) then
if ((dz2).gt.eps) then z(i)=sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))))
z(i)=max(sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis)))),sp11) else
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))+(dz2/psi2s))),sp13)
else else
z(i)=min(sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis)))),sp11) z(i)=min(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))+(dz2/psi2s))),sp13)
endif endif
endif endif
! 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 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 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) dist=abs(a1*x(i)+b1*y(i)+c1)/sqrt(a1**2.d0+b1**2.d0)
if (abs(dz2).gt.eps) then if (dist.lt.(dz2/psi2s)) then
if ((dz2).gt.eps) then z(i)=sp01-psi2s*dist
z(i)=max(sp01-psi2s*dist,sp11) else
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else else
z(i)=min(sp01-psi2s*dist,sp11) z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif endif
endif endif
! if ((dz2).gt.eps) then
! z(i)=max(sp01-psi2s*dist,sp11)
! else
! z(i)=min(sp01-psi2s*dist,sp11)
! endif
elseif (y(i).ge.m2*x(i)+bb2 .and. y(i).lt.m1*x(i)+bb1) then 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) dist=abs(a1*x(i)+b1*y(i)+c1)/sqrt(a1**2.d0+b1**2.d0)
z(i)=sp01+dist*psis z(i)=sp01+dist*psis
...@@ -781,22 +822,36 @@ select case(surface_type) ...@@ -781,22 +822,36 @@ select case(surface_type)
z(i)=sp01+dist*psis z(i)=sp01+dist*psis
else else
dist=abs(a4*x(i)+b4*y(i)+c4)/sqrt(a4**2.d0+b4**2.d0) dist=abs(a4*x(i)+b4*y(i)+c4)/sqrt(a4**2.d0+b4**2.d0)
if (abs(dz2).gt.eps) then if (dist.lt.(dz2/psi2s)) then
if ((dz2).gt.eps) then z(i)=sp01-psi2s*dist
z(i)=max(sp01-psi2s*dist,sp11) else
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else else
z(i)=min(sp01-psi2s*dist,sp11) z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif endif
endif endif
! 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 elseif (x(i).lt.x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)) then
if (abs(dz2).gt.eps) then if (x(i).lt.x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)-(dz2/psi2s)) then
if ((dz2).gt.eps) then if (dz4.gt.eps) then
z(i)=max(sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i)),sp11) z(i)=max(sp01-dz2-psi3s*(((x0+y0w)+(y1-y2)*theta1s+(dz1/(2.d0*psis))+(dz2/psi2s))-x(i)),sp13)
else else
z(i)=min(sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i)),sp11) z(i)=min(sp01-dz2-psi3s*(((x0+y0w)+(y1-y2)*theta1s+(dz1/(2.d0*psis))+(dz2/psi2s))-x(i)),sp13)
endif endif
else
z(i)=sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i))
endif endif
! 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
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 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)))) 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 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
...@@ -804,13 +859,20 @@ select case(surface_type) ...@@ -804,13 +859,20 @@ select case(surface_type)
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 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)))) z(i)=sp01+psis*(x(i)-((x0+y0w+(y1-y2)*theta1s)+(dz1/(2.d0*psis))))
else else
if (abs(dz2).gt.eps) then if (x(i).gt.x0+y0w+(y1-y2)*theta1s-dz1/(2.d0*psis)+(dz2/psi2s)) then
if ((dz2).gt.eps) then if (dz4.gt.eps) then
z(i)=max(sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))),sp11) z(i)=max(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(y1-y2)*theta1s+(dz1/(2.d0*psis))+(dz2/psi2s))),sp13)
else else
z(i)=min(sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))),sp11) z(i)=min(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(y1-y2)*theta1s+(dz1/(2.d0*psis))+(dz2/psi2s))),sp13)
endif endif
else
z(i)=sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)))
endif endif
! 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 enddo
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment