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

Merge up to r866 in svn: with new surface types and sedimentation options, for... 

Merge up to r866 in svn: with new surface types and sedimentation options, for minibasins and combined minibasins with progradation; clean up formatting
parent 5799fb2d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/create_surfaces.f90
+ 798
− 748
View file @ cd2fa96e
......@@ -168,7 +168,7 @@ else
select case(surface_type)
case (1,2,5,6,7,8,9,10,11,12,13,14,15,16,17)
case (1,2,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)
! if (surface%rand) then
indix=1
! convex hull
......@@ -694,576 +694,576 @@ select case(surface_type)
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 lower outer dip change)
! 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
!dz2=sp01-sp11
dz2=sp01-sp11
dz3=sp01-sp13
dz4=dz3-dz2
psi=sp03*pi/180.d0
!psi2=sp12*pi/180.d0
psi2=sp12*pi/180.d0
psi3=sp14*pi/180.d0
theta1=sp08*pi/180.d0
theta2=sp10*pi/180.d0
psis=sign(tan(psi),dz1)
psi2s=sign(tan(psi2),dz2)
psi3s=sign(tan(psi3),dz3)
theta1s=tan(theta1)
theta2s=tan(theta2)
y0s=1.d0
y1s=tan(theta1/2.d0)
y2s=tan((theta1+theta2)/2.d0)
x0=sp04
y0w=sp05+dz1/(2.d0*psis)
y0=sp06+dz1/(2.d0*psis)
y0a=y0+dz1/(2.d0*psis)
y0b=y0-dz1/(2.d0*psis)
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 (y(i).ge.y0a+(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(y(i)-(y0a+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(y(i)-(y0a+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(y(i)-y0a)
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
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 (x(i).lt.((x0-y0w)-(dz1/(2.d0*psis)))-(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(((x0-y0w)-(dz2/psi2s))-x(i)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(((x0-y0w)-(dz2/psi2s))-x(i)),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i))
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
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 (x(i).lt.((x0+y0w)+(dz1/(2.d0*psis)))+(dz2/psi2s)) then
z(i)=sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))))
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
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
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 (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
else
z(i)=sp13
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
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 (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
else
z(i)=sp13
endif
endif
! if ((dz2).gt.eps) then
! z(i)=max(sp01-psi2s*dist,sp11)
! else
! z(i)=min(sp01-psi2s*dist,sp11)
! endif
endif
elseif (x(i).lt.x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)) then
if (x(i).lt.x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)-(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(((x0-y0w)-(y1-y2)*theta1s-(dz2/psi2s))-x(i)),sp13)
!z(i)=sp13
else
z(i)=min(sp01-dz2-psi3s*(((x0-y0w)-(y1-y2)*theta1s-(dz2/psi2s))-x(i)),sp13)
!z(i)=sp13
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i))
!z(i)=sp13
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
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 (x(i).gt.x0+y0w+(y1-y2)*theta1s-dz1/(2.d0*psis)+(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-(x0+y0w+(y1-y2)*theta1s+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-(x0+y0w+(y1-y2)*theta1s+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
endif
else
!z(i)=sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)))
z(i)=sp01-psi2s*(x(i)-(x0+y0w+(y1-y2)*theta1s+dz1/(2.d0*psis)))
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
enddo
! 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 lower outer dip change)
! 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
!dz2=sp01-sp11
dz2=sp01-sp11
dz3=sp01-sp13
dz4=dz3-dz2
psi=sp03*pi/180.d0
!psi2=sp12*pi/180.d0
psi2=sp12*pi/180.d0
psi3=sp14*pi/180.d0
theta1=sp08*pi/180.d0
theta2=sp10*pi/180.d0
psis=sign(tan(psi),dz1)
psi2s=sign(tan(psi2),dz2)
psi3s=sign(tan(psi3),dz3)
theta1s=tan(theta1)
theta2s=tan(theta2)
y0s=1.d0
y1s=tan(theta1/2.d0)
y2s=tan((theta1+theta2)/2.d0)
x0=sp04
y0w=sp05+dz1/(2.d0*psis)
y0=sp06+dz1/(2.d0*psis)
y0a=y0+dz1/(2.d0*psis)
y0b=y0-dz1/(2.d0*psis)
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 (y(i).ge.y0a+(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(y(i)-(y0a+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(y(i)-(y0a+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(y(i)-y0a)
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
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 (x(i).lt.((x0-y0w)-(dz1/(2.d0*psis)))-(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(((x0-y0w)-(dz2/psi2s))-x(i)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(((x0-y0w)-(dz2/psi2s))-x(i)),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i))
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
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 (x(i).lt.((x0+y0w)+(dz1/(2.d0*psis)))+(dz2/psi2s)) then
z(i)=sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))))
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
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
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 (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
else
z(i)=sp13
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
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 (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
else
z(i)=sp13
endif
endif
! if ((dz2).gt.eps) then
! z(i)=max(sp01-psi2s*dist,sp11)
! else
! z(i)=min(sp01-psi2s*dist,sp11)
! endif
endif
elseif (x(i).lt.x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)) then
if (x(i).lt.x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)-(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(((x0-y0w)-(y1-y2)*theta1s-(dz2/psi2s))-x(i)),sp13)
!z(i)=sp13
else
z(i)=min(sp01-dz2-psi3s*(((x0-y0w)-(y1-y2)*theta1s-(dz2/psi2s))-x(i)),sp13)
!z(i)=sp13
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i))
!z(i)=sp13
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
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 (x(i).gt.x0+y0w+(y1-y2)*theta1s-dz1/(2.d0*psis)+(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-(x0+y0w+(y1-y2)*theta1s+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-(x0+y0w+(y1-y2)*theta1s+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
endif
else
!z(i)=sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)))
z(i)=sp01-psi2s*(x(i)-(x0+y0w+(y1-y2)*theta1s+dz1/(2.d0*psis)))
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
enddo
case (14)
! a 3D embankment of finite length and with rounded 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 (xsym)
! 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 lower outer dip change)
! 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
! sp15 is the radius used for rounding the second kink
dz1=sp02-sp01
dz2=sp01-sp11
dz3=sp01-sp13
dz4=dz3-dz2
psi1=sp03*pi/180.d0
psi2=sp12*pi/180.d0
psi3=sp14*pi/180.d0
theta1=sp08*pi/180.d0
theta2=sp10*pi/180.d0
y2rad=sp15
psi1s=sign(tan(psi1),dz1)
psi2s=sign(tan(psi2),dz2)
psi3s=sign(tan(psi3),dz3)
theta1s=tan(theta1)
theta2s=tan(theta2)
y0s=1.d0
y1s=tan(theta1/2.d0)
y2s=tan((theta1+theta2)/2.d0)
xsym=sp04
y0w=sp05
y0=sp06
y1=sp07
y2=sp09
y0u=y0
y0l=y0+dz1/psi1s
y1u=y1
y1l=y1+(dz1/psi1s)*y1s
y2u=y2
y2l=y2+(dz1/psi1s)*y2s
! if (iproc == 0) write (*,*) 'y0u: ',y0u
! if (iproc == 0) write (*,*) 'y0l: ',y0l
! if (iproc == 0) write (*,*) 'y1u: ',y1u
! if (iproc == 0) write (*,*) 'y1l: ',y1l
! if (iproc == 0) write (*,*) 'y2u: ',y2u
! if (iproc == 0) write (*,*) 'y2l: ',y2l
! if (iproc == 0) write (*,*) 'y2rad: ',y2rad
x0ut=xsym+y0w
x0lt=xsym+y0w+dz1/psi1s
x0ub=xsym-y0w
x0lb=xsym-y0w-dz1/psi1s
x1ut=x0ut
x1lt=x0lt
x1ub=x0ub
x1lb=x0lb
x2ut=x1ut+(y1-y2)*theta1s
x2lt=x1lt+(y1-y2)*theta1s
x2ub=x1ub-(y1-y2)*theta1s
x2lb=x1lb-(y1-y2)*theta1s
! if (iproc == 0) write (*,*) 'x0ut: ',x0ut
! if (iproc == 0) write (*,*) 'x0lt: ',x0lt
! if (iproc == 0) write (*,*) 'x0ub: ',x0ub
! if (iproc == 0) write (*,*) 'x0lb: ',x0lb
! if (iproc == 0) write (*,*) 'x1ut: ',x1ut
! if (iproc == 0) write (*,*) 'x1lt: ',x1lt
! if (iproc == 0) write (*,*) 'x1ub: ',x1ub
! if (iproc == 0) write (*,*) 'x1lb: ',x1lb
! if (iproc == 0) write (*,*) 'x2ut: ',x2ut
! if (iproc == 0) write (*,*) 'x2lt: ',x2lt
! if (iproc == 0) write (*,*) 'x2ub: ',x2ub
! if (iproc == 0) write (*,*) 'x2lb: ',x2lb
y12mt=(y2u-y1u)/(x2ut-x1ut)
y12mb=(y2u-y1u)/(x2ub-x1ub)
! if (iproc == 0) write (*,*) 'y12mt: ',y12mt
! if (iproc == 0) write (*,*) 'y12mb: ',y12mb
if (abs(theta1).gt.eps) then
y2as=tan(pi/2.d0-theta1)
else
y2as=1.d0
endif
if (abs(theta2).gt.eps) then
y2bs=tan(pi/2.d0-theta2)
else
y2bs=1.d0
endif
! if (iproc == 0) write (*,*) 'y2as: ',y2as
! if (iproc == 0) write (*,*) 'y2bs: ',y2bs
y2au=y2u+(y2rad*tan((theta1-theta2)/2.d0))*cos(theta1)
if (abs(theta1).gt.eps) then
y2al=y2u+(dz1/psi1s*y2as)+(y2rad*tan((theta1-theta2)/2.d0))*cos(theta1)
else
y2al=y2au
endif
y2bu=y2u-y2rad*tan((theta1-theta2)/2.d0)
! if (abs(theta2).gt.eps) then
! y2bl=y2u-(dz1/psi1s*y2as)-(y2rad*tan((theta1-theta2)/2.d0))
! else
y2bl=y2bu
! endif
! a 3D embankment of finite length and with rounded 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 (xsym)
! 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 lower outer dip change)
! 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
! sp15 is the radius used for rounding the second kink
dz1=sp02-sp01
dz2=sp01-sp11
dz3=sp01-sp13
dz4=dz3-dz2
psi1=sp03*pi/180.d0
psi2=sp12*pi/180.d0
psi3=sp14*pi/180.d0
theta1=sp08*pi/180.d0
theta2=sp10*pi/180.d0
y2rad=sp15
psi1s=sign(tan(psi1),dz1)
psi2s=sign(tan(psi2),dz2)
psi3s=sign(tan(psi3),dz3)
theta1s=tan(theta1)
theta2s=tan(theta2)
y0s=1.d0
y1s=tan(theta1/2.d0)
y2s=tan((theta1+theta2)/2.d0)
xsym=sp04
y0w=sp05
y0=sp06
y1=sp07
y2=sp09
y0u=y0
y0l=y0+dz1/psi1s
y1u=y1
y1l=y1+(dz1/psi1s)*y1s
y2u=y2
y2l=y2+(dz1/psi1s)*y2s
! if (iproc == 0) write (*,*) 'y0u: ',y0u
! if (iproc == 0) write (*,*) 'y0l: ',y0l
! if (iproc == 0) write (*,*) 'y1u: ',y1u
! if (iproc == 0) write (*,*) 'y1l: ',y1l
! if (iproc == 0) write (*,*) 'y2u: ',y2u
! if (iproc == 0) write (*,*) 'y2l: ',y2l
! if (iproc == 0) write (*,*) 'y2rad: ',y2rad
x0ut=xsym+y0w
x0lt=xsym+y0w+dz1/psi1s
x0ub=xsym-y0w
x0lb=xsym-y0w-dz1/psi1s
x1ut=x0ut
x1lt=x0lt
x1ub=x0ub
x1lb=x0lb
x2ut=x1ut+(y1-y2)*theta1s
x2lt=x1lt+(y1-y2)*theta1s
x2ub=x1ub-(y1-y2)*theta1s
x2lb=x1lb-(y1-y2)*theta1s
! if (iproc == 0) write (*,*) 'x0ut: ',x0ut
! if (iproc == 0) write (*,*) 'x0lt: ',x0lt
! if (iproc == 0) write (*,*) 'x0ub: ',x0ub
! if (iproc == 0) write (*,*) 'x0lb: ',x0lb
! if (iproc == 0) write (*,*) 'x1ut: ',x1ut
! if (iproc == 0) write (*,*) 'x1lt: ',x1lt
! if (iproc == 0) write (*,*) 'x1ub: ',x1ub
! if (iproc == 0) write (*,*) 'x1lb: ',x1lb
! if (iproc == 0) write (*,*) 'x2ut: ',x2ut
! if (iproc == 0) write (*,*) 'x2lt: ',x2lt
! if (iproc == 0) write (*,*) 'x2ub: ',x2ub
! if (iproc == 0) write (*,*) 'x2lb: ',x2lb
y12mt=(y2u-y1u)/(x2ut-x1ut)
y12mb=(y2u-y1u)/(x2ub-x1ub)
! if (iproc == 0) write (*,*) 'y12mt: ',y12mt
! if (iproc == 0) write (*,*) 'y12mb: ',y12mb
if (abs(theta1).gt.eps) then
y2as=tan(pi/2.d0-theta1)
else
y2as=1.d0
endif
if (abs(theta2).gt.eps) then
y2bs=tan(pi/2.d0-theta2)
else
y2bs=1.d0
endif
! if (iproc == 0) write (*,*) 'y2as: ',y2as
! if (iproc == 0) write (*,*) 'y2bs: ',y2bs
y2au=y2u+(y2rad*tan((theta1-theta2)/2.d0))*cos(theta1)
if (abs(theta1).gt.eps) then
y2al=y2u+(dz1/psi1s*y2as)+(y2rad*tan((theta1-theta2)/2.d0))*cos(theta1)
else
y2al=y2au
endif
y2bu=y2u-y2rad*tan((theta1-theta2)/2.d0)
! if (abs(theta2).gt.eps) then
! y2bl=y2u-(dz1/psi1s*y2as)-(y2rad*tan((theta1-theta2)/2.d0))
! else
y2bl=y2bu
! endif
! SHOULD THE y12mb's BELOW BE sin(theta1)'s INSTEAD?????????
x2aut=x2ut-(y2rad*tan((theta1-theta2)/2.d0))*y12mb
if (abs(theta1).gt.eps) then
x2alt=x2ut+(dz1/psi1s*(1.d0/y2as))-(y2rad*tan((theta1-theta2)/2.d0))*y12mb
else
x2alt=x2ut+dz1/psi1s-(y2rad*tan((theta1-theta2)/2.d0))*y12mb
endif
x2but=x2ut
x2blt=x2lt
x2aub=x2ub+(y2rad*tan((theta1-theta2)/2.d0))*y12mb
if (abs(theta2).gt.eps) then
x2alb=x2ub-(dz1/psi1s*(1.d0/y2as))+(y2rad*tan((theta1-theta2)/2.d0))*y12mb
else
x2alb=x2ub-dz1/psi1s+(y2rad*tan((theta1-theta2)/2.d0))*y12mb
endif
x2bub=x2ub
x2blb=x2lb
! if (iproc == 0) write (*,*) 'y2au: ',y2au
! if (iproc == 0) write (*,*) 'y2al: ',y2al
! if (iproc == 0) write (*,*) 'y2bu: ',y2bu
! if (iproc == 0) write (*,*) 'y2bl: ',y2bl
! if (iproc == 0) write (*,*) 'x2aut: ',x2aut
! if (iproc == 0) write (*,*) 'x2alt: ',x2alt
! if (iproc == 0) write (*,*) 'x2but: ',x2but
! if (iproc == 0) write (*,*) 'x2blt: ',x2blt
! if (iproc == 0) write (*,*) 'x2aub: ',x2aub
! if (iproc == 0) write (*,*) 'x2alb: ',x2alb
! if (iproc == 0) write (*,*) 'x2bub: ',x2bub
! if (iproc == 0) write (*,*) 'x2blb: ',x2blb
m1=(y2l-y1l)/(x2lb-x1lb)
!bb1=(y1l+(dz1/psi1s)*y1s)-m1*(x1lb-(dz1/psi1s))
bb1=y1l-m1*x1lb
a1=-m1
b1=1
c1=-bb1
m2=(y2u-y1u)/(x2ub-x1ub)
bb2=y1u-m2*x1ub
a2=-m2
b2=1
c2=-bb2
m3=(y2u-y1u)/(x2ut-x1ut)
bb3=y1u-m3*x1ut
a3=-m3
b3=1
c3=-bb3
m4=(y2l-y1l)/(x2lt-x1lt)
bb4=(y1l+(dz1/psi1s)*y1s)-m4*(x1lt+dz1/psi1s)
a4=-m4
b4=1
c4=-bb4
do i=1,ns
z(i)=sp01
if (y(i).ge.-y0s*x(i)+(y0u+x0ub) .and. y(i).ge.y0s*x(i)+(y0u-x0ut)) then
if (y(i).ge.y0l) then
if (y(i).ge.y0l+(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(y(i)-(y0l+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(y(i)-(y0l+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(y(i)-y0l)
endif
elseif (y(i).ge.y0u) then
z(i)=sp01+psi1s*(y0l-y(i))
elseif (y(i).lt.y0u) then
z(i)=sp02
endif
elseif (y(i).ge.-y1s*x(i)+(y1s*x1ub+y1u) .and. y(i).ge.y1s*x(i)+(-y1s*x1ut+y1u)) then
if (x(i).lt.x1lb) then
if (x(i).lt.x1lb-dz2/psi2s) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x1ub-dz2/psi2s)-x(i),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x1ub-dz2/psi2s)-x(i),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(x1lb-x(i))
endif
elseif (x(i).ge.x1lb .and. x(i).lt.x1ub) then
z(i)=sp01+psi1s*(x(i)-x1lb)
elseif (x(i).ge.x1ub .and. x(i).lt.x1ut) then
z(i)=sp02
elseif (x(i).ge.x1ut .and. x(i).lt.x1lt) then
z(i)=sp01-psi1s*(x(i)-x1lt)
else
if (x(i).lt.x1lt+(dz2/psi2s)) then
z(i)=sp01-psi2s*(x(i)-x1lt)
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-x1lt),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-x1lt),sp13)
endif
else
z(i)=sp13
endif
endif
endif
elseif (y(i).ge.-y2as*x(i)+(y2as*x2aub+y2au) .and. y(i).ge.y2as*x(i)+(-y2as*x2aut+y2au)) then
! elseif (y(i).ge.-y2s*x(i)+(y2s*x2ub+y2u) .and. y(i).ge.y2s*x(i)+(-y2s*x2ut+y2u)) then
! z(i)=0.13d0
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 (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
else
z(i)=sp13
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*psi1s
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*psi1s
else
dist=abs(a4*x(i)+b4*y(i)+c4)/sqrt(a4**2.d0+b4**2.d0)
if (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
else
z(i)=sp13
endif
endif
endif
!elseif (y(i).ge.-y2as*x(i)+(y2as*x2aub+y2au) .and. y(i).ge.y2as*x(i)+(-y2as*x2aut+y2au)) then
elseif (y(i).ge.y2bu) then
dist=sqrt((x(i)-(x2bub+y2rad))**2.d0+(y(i)-y2bu)**2.d0)
if (dist <= y2rad) then
z(i)=sp02
elseif (dist <= y2rad+dz1/psi1s) then
z(i)=sp01+((y2rad+dz1/psi1s)-dist)*psi1s
elseif (dist <= y2rad+dz1/psi1s+dz2/psi2s) then
z(i)=sp01-psi2s*(dist-(y2rad+dz1/psi1s))
else
if (abs(dz4) > eps) then
if (dz4 > eps) then
z(i)=max(sp01-dz2-psi3s*((y2rad+dz1/psi1s+dz2/psi2s+dz3/psi3s)-dist),sp13)
else
z(i)=min(sp01-dz2-psi3s*((y2rad+dz1/psi1s+dz2/psi2s+dz3/psi3s)-dist),sp13)
endif
else
z(i)=sp13
endif
endif
elseif (x(i) < x2lb) then
if (x(i) < x2lb-(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x2lb-(dz2/psi2s)-x(i)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x2lb-(dz2/psi2s)-x(i)),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(x2lb-x(i))
endif
elseif (x(i).ge.x2lb .and. x(i).lt.x2ub) then
z(i)=sp01+psi1s*(x(i)-x2lb)
elseif (x(i).ge.x2ub .and. x(i).lt.x2ut) then
z(i)=sp02
elseif (x(i).ge.x2ut .and. x(i).lt.x2lt) then
z(i)=sp01+psi1s*(x2lb-x(i))
else
if (x(i) > x2ut+(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-(x2ut+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-(x2ut+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(x(i)-x2lt)
endif
endif
enddo
x2aut=x2ut-(y2rad*tan((theta1-theta2)/2.d0))*y12mb
if (abs(theta1).gt.eps) then
x2alt=x2ut+(dz1/psi1s*(1.d0/y2as))-(y2rad*tan((theta1-theta2)/2.d0))*y12mb
else
x2alt=x2ut+dz1/psi1s-(y2rad*tan((theta1-theta2)/2.d0))*y12mb
endif
x2but=x2ut
x2blt=x2lt
x2aub=x2ub+(y2rad*tan((theta1-theta2)/2.d0))*y12mb
if (abs(theta2).gt.eps) then
x2alb=x2ub-(dz1/psi1s*(1.d0/y2as))+(y2rad*tan((theta1-theta2)/2.d0))*y12mb
else
x2alb=x2ub-dz1/psi1s+(y2rad*tan((theta1-theta2)/2.d0))*y12mb
endif
x2bub=x2ub
x2blb=x2lb
! if (iproc == 0) write (*,*) 'y2au: ',y2au
! if (iproc == 0) write (*,*) 'y2al: ',y2al
! if (iproc == 0) write (*,*) 'y2bu: ',y2bu
! if (iproc == 0) write (*,*) 'y2bl: ',y2bl
! if (iproc == 0) write (*,*) 'x2aut: ',x2aut
! if (iproc == 0) write (*,*) 'x2alt: ',x2alt
! if (iproc == 0) write (*,*) 'x2but: ',x2but
! if (iproc == 0) write (*,*) 'x2blt: ',x2blt
! if (iproc == 0) write (*,*) 'x2aub: ',x2aub
! if (iproc == 0) write (*,*) 'x2alb: ',x2alb
! if (iproc == 0) write (*,*) 'x2bub: ',x2bub
! if (iproc == 0) write (*,*) 'x2blb: ',x2blb
m1=(y2l-y1l)/(x2lb-x1lb)
!bb1=(y1l+(dz1/psi1s)*y1s)-m1*(x1lb-(dz1/psi1s))
bb1=y1l-m1*x1lb
a1=-m1
b1=1
c1=-bb1
m2=(y2u-y1u)/(x2ub-x1ub)
bb2=y1u-m2*x1ub
a2=-m2
b2=1
c2=-bb2
m3=(y2u-y1u)/(x2ut-x1ut)
bb3=y1u-m3*x1ut
a3=-m3
b3=1
c3=-bb3
m4=(y2l-y1l)/(x2lt-x1lt)
bb4=(y1l+(dz1/psi1s)*y1s)-m4*(x1lt+dz1/psi1s)
a4=-m4
b4=1
c4=-bb4
do i=1,ns
z(i)=sp01
if (y(i).ge.-y0s*x(i)+(y0u+x0ub) .and. y(i).ge.y0s*x(i)+(y0u-x0ut)) then
if (y(i).ge.y0l) then
if (y(i).ge.y0l+(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(y(i)-(y0l+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(y(i)-(y0l+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(y(i)-y0l)
endif
elseif (y(i).ge.y0u) then
z(i)=sp01+psi1s*(y0l-y(i))
elseif (y(i).lt.y0u) then
z(i)=sp02
endif
elseif (y(i).ge.-y1s*x(i)+(y1s*x1ub+y1u) .and. y(i).ge.y1s*x(i)+(-y1s*x1ut+y1u)) then
if (x(i).lt.x1lb) then
if (x(i).lt.x1lb-dz2/psi2s) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x1ub-dz2/psi2s)-x(i),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x1ub-dz2/psi2s)-x(i),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(x1lb-x(i))
endif
elseif (x(i).ge.x1lb .and. x(i).lt.x1ub) then
z(i)=sp01+psi1s*(x(i)-x1lb)
elseif (x(i).ge.x1ub .and. x(i).lt.x1ut) then
z(i)=sp02
elseif (x(i).ge.x1ut .and. x(i).lt.x1lt) then
z(i)=sp01-psi1s*(x(i)-x1lt)
else
if (x(i).lt.x1lt+(dz2/psi2s)) then
z(i)=sp01-psi2s*(x(i)-x1lt)
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-x1lt),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-x1lt),sp13)
endif
else
z(i)=sp13
endif
endif
endif
elseif (y(i).ge.-y2as*x(i)+(y2as*x2aub+y2au) .and. y(i).ge.y2as*x(i)+(-y2as*x2aut+y2au)) then
! elseif (y(i).ge.-y2s*x(i)+(y2s*x2ub+y2u) .and. y(i).ge.y2s*x(i)+(-y2s*x2ut+y2u)) then
! z(i)=0.13d0
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 (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
else
z(i)=sp13
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*psi1s
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*psi1s
else
dist=abs(a4*x(i)+b4*y(i)+c4)/sqrt(a4**2.d0+b4**2.d0)
if (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
else
z(i)=sp13
endif
endif
endif
!elseif (y(i).ge.-y2as*x(i)+(y2as*x2aub+y2au) .and. y(i).ge.y2as*x(i)+(-y2as*x2aut+y2au)) then
elseif (y(i).ge.y2bu) then
dist=sqrt((x(i)-(x2bub+y2rad))**2.d0+(y(i)-y2bu)**2.d0)
if (dist <= y2rad) then
z(i)=sp02
elseif (dist <= y2rad+dz1/psi1s) then
z(i)=sp01+((y2rad+dz1/psi1s)-dist)*psi1s
elseif (dist <= y2rad+dz1/psi1s+dz2/psi2s) then
z(i)=sp01-psi2s*(dist-(y2rad+dz1/psi1s))
else
if (abs(dz4) > eps) then
if (dz4 > eps) then
z(i)=max(sp01-dz2-psi3s*((y2rad+dz1/psi1s+dz2/psi2s+dz3/psi3s)-dist),sp13)
else
z(i)=min(sp01-dz2-psi3s*((y2rad+dz1/psi1s+dz2/psi2s+dz3/psi3s)-dist),sp13)
endif
else
z(i)=sp13
endif
endif
elseif (x(i) < x2lb) then
if (x(i) < x2lb-(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x2lb-(dz2/psi2s)-x(i)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x2lb-(dz2/psi2s)-x(i)),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(x2lb-x(i))
endif
elseif (x(i).ge.x2lb .and. x(i).lt.x2ub) then
z(i)=sp01+psi1s*(x(i)-x2lb)
elseif (x(i).ge.x2ub .and. x(i).lt.x2ut) then
z(i)=sp02
elseif (x(i).ge.x2ut .and. x(i).lt.x2lt) then
z(i)=sp01+psi1s*(x2lb-x(i))
else
if (x(i) > x2ut+(dz2/psi2s)) then
if (abs(dz4).gt.eps) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-(x2ut+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-(x2ut+(dz2/psi2s))),sp13)
endif
else
z(i)=sp13
endif
else
z(i)=sp01-psi2s*(x(i)-x2lt)
endif
endif
enddo
case (15)
! sp01 is the z base level (z0)
......@@ -1342,191 +1342,241 @@ select case(surface_type)
end if
enddo
case (16)
! sp01 = z1, left hand side z level
! sp02 = z2, right hand side z level
! sp03 = x0, x coord of corner point
! sp04 = y0, y coord of corner point
! sp05 = alpha1, 0 < alpha1 < 90deg
! sp06 = alpha2
! sp07 = theta, 0 <= theta < 90deg
! sp08 = theta_switch, <0 = theta is perpendicular to strike
! >0 = theta is perpendicular to y axis
sp05r = sp05*pi/180.0
sp06r = sp06*pi/180.0
sp07r = sp07*pi/180.0
if (sp02.lt.sp01) then
dzsign = -1.0
! NB! Currently does not work
! TODO: Fix to allow z_left > z_right
case (16)
! sp01 = z1, left hand side z level
! sp02 = z2, right hand side z level
! sp03 = x0, x coord of corner point
! sp04 = y0, y coord of corner point
! sp05 = alpha1, 0 < alpha1 < 90deg
! sp06 = alpha2
! sp07 = theta, 0 <= theta < 90deg
! sp08 = theta_switch, <0 = theta is perpendicular to strike
! >0 = theta is perpendicular to y axis
sp05r = sp05*pi/180.0
sp06r = sp06*pi/180.0
sp07r = sp07*pi/180.0
if (sp02.lt.sp01) then
dzsign = -1.0
! NB! Currently does not work
! TODO: Fix to allow z_left > z_right
else
dzsign = 1.0
endif
do i=1,ns
if (y(i).le.sp04) then
! lower y region
dysign = 1.0
alpha = sp06
alphar = sp06r
alphar_con = sp05r
else
! y(i) > sp04=y0,
! upper y region
alpha = sp05
alphar = sp05r
dysign = -1.0
alphar_con = sp06r
endif
dy = dysign * (sp04 - y(i))
x_upper_line = sp03-dy*tan(alphar) ! 'upper line' === right hand line
x_upper_line_con = sp03-dysign*dy*tan(alphar_con)
if (sp08.lt.0.0) then
! theta is given as perpendicular to strike
! shortest horizontal distance of the upper line from the lower line
dist_upper_lower = dzsign*(sp02-sp01)/tan(sp07r)
! location x of the lower line at y=y(i):
! first, distance of the lower line from the upper line, along x axis
dx_upper_lower = dist_upper_lower / cos(alphar)
dx_upper_lower_con = dist_upper_lower / cos(alphar_con)
! then the absolute location
x_lower_line = x_upper_line - dx_upper_lower
x_lower_line_con = x_upper_line_con - dx_upper_lower_con
if (x(i).le.x_upper_line) then
! left of upper line
!!if (x(i).le.x_lower_line) then
!! ! also left of lower line
!! z(i) = sp01
!!else
! left of upper line, right of lower line, i.e. on the slope:
! calc effective angle at the slope,
! the angle thetaeff, in a XZ-section.
! (we only need tangent of that)
thetar_eff_tan = (sp02-sp01)/(dx_upper_lower*dzsign)
thetar_eff_tan_con = (sp02-sp01)/(dx_upper_lower_con*dzsign)
! calc z
dx = x(i) - x_lower_line
dx_con = x(i) - x_lower_line_con
zi = sp01 + dx * thetar_eff_tan
zi_con = sp01 + dx_con * thetar_eff_tan_con
if (dzsign.lt.0.0) then
z(i) = min(sp01,zi_con,zi) ! near the groove by the corner
! the other side (cf dysign) might overlap
else
z(i) = max(sp01,zi_con,zi) ! near the groove by the corner
! the other side (cf dysign) might overlap
endif
!!endif
else
! right of upper line
z(i) = sp02
endif
else
! theta is given assuming strike is perpendicular to x axis
! shortest horizontal distance of the upper line from the lower line
dist_upper_lower = dzsign*(sp02-sp01)/tan(sp07r)
! ... which, by assumption, is also the distance in x direction
dx_upper_lower = dist_upper_lower
x_lower_line = x_upper_line - dx_upper_lower
if (x(i).le.x_upper_line) then
! left of upper line
if (x(i).le.x_lower_line) then
! also left of lower line
z(i) = sp01
else
! left of upper line, right of lower line, i.e. on the slope:
dx = x(i) - x_lower_line
z(i) = sp01 + dx * tan(sp07r)
endif
else
! right of upper line
z(i) = sp02
endif
endif
! take care of the corner ...
if (sp08.lt.0.0) then
endif
if (z(i).lt.0.0) then
call stop_run ('z less than zero$')
elseif (z(i).gt.1.0) then
call stop_run ('z larger than one$')
endif
enddo
case (17)
! basin with differing height basement steps, geometry constant in y
! sp01 is the z level at the bottom of the basin
! sp02 and 03 are x1,x2
! sp04 is not used (this could specify a 2nd slope later)
! sp05 is the distal basement step thickness
! sp06 is the proximal basement step thickness
! sp07 is the slope
m=tan(sp07*pi/180.d0)
m=sign(m,-sp06)
x1a=sp02+(sp06/(2.d0*m))
x1b=sp02-(sp06/(2.d0*m))
x2a=sp03+(sp05/(2.d0*m))
x2b=sp03-(sp05/(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
z(i)=sp01+sp06
elseif (x(i).le.x1b) then
z(i)=(sp01+sp06)-(((x(i)-x1a)/(-sp06/m))*sp06)
elseif (x(i).le.x2a) then
z(i)=sp01
elseif (x(i).le.x2b) then
z(i)=(sp01+sp05)-((1.d0-((x(i)-x2a)/(-sp05/m)))*sp05)
!zx=sp01-((1.d0-((x(i)-x2a)/(sp06/m)))*sp06)
else
dzsign = 1.0
z(i)=sp01+sp05
endif
do i=1,ns
if (y(i).le.sp04) then
! lower y region
dysign = 1.0
alpha = sp06
alphar = sp06r
alphar_con = sp05r
else
! y(i) > sp04=y0,
! upper y region
alpha = sp05
alphar = sp05r
dysign = -1.0
alphar_con = sp06r
endif
dy = dysign * (sp04 - y(i))
x_upper_line = sp03-dy*tan(alphar) ! 'upper line' === right hand line
x_upper_line_con = sp03-dysign*dy*tan(alphar_con)
if (sp08.lt.0.0) then
! theta is given as perpendicular to strike
! shortest horizontal distance of the upper line from the lower line
dist_upper_lower = dzsign*(sp02-sp01)/tan(sp07r)
! location x of the lower line at y=y(i):
! first, distance of the lower line from the upper line, along x axis
dx_upper_lower = dist_upper_lower / cos(alphar)
dx_upper_lower_con = dist_upper_lower / cos(alphar_con)
! then the absolute location
x_lower_line = x_upper_line - dx_upper_lower
x_lower_line_con = x_upper_line_con - dx_upper_lower_con
if (x(i).le.x_upper_line) then
! left of upper line
!!if (x(i).le.x_lower_line) then
!! ! also left of lower line
!! z(i) = sp01
!!else
! left of upper line, right of lower line, i.e. on the slope:
! calc effective angle at the slope,
! the angle thetaeff, in a XZ-section.
! (we only need tangent of that)
thetar_eff_tan = (sp02-sp01)/(dx_upper_lower*dzsign)
thetar_eff_tan_con = (sp02-sp01)/(dx_upper_lower_con*dzsign)
! calc z
dx = x(i) - x_lower_line
dx_con = x(i) - x_lower_line_con
zi = sp01 + dx * thetar_eff_tan
zi_con = sp01 + dx_con * thetar_eff_tan_con
if (dzsign.lt.0.0) then
z(i) = min(sp01,zi_con,zi) ! near the groove by the corner
! the other side (cf dysign) might overlap
else
z(i) = max(sp01,zi_con,zi) ! near the groove by the corner
! the other side (cf dysign) might overlap
endif
!!endif
else
! right of upper line
z(i) = sp02
endif
else
! theta is given assuming strike is perpendicular to x axis
!! There is a smarter way to do this, but...
!if (sp06.gt.eps) then
! if (y(i).le.y1a) then
! z(i)=sp01
! elseif (y(i).le.y1b) then
! z(i)=max(sp01-(((y(i)-y1a)/(-sp06/m))*sp06),zx)
! elseif (y(i).le.y2a) then
! z(i)=max(sp01-sp06,zx)
! elseif (y(i).le.y2b) then
! z(i)=max(sp01-((1.d0-((y(i)-y2a)/(-sp06/m)))*sp06),zx)
! else
! z(i)=sp01
! endif
!else
! 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
!endif
end do
! shortest horizontal distance of the upper line from the lower line
dist_upper_lower = dzsign*(sp02-sp01)/tan(sp07r)
! ... which, by assumption, is also the distance in x direction
dx_upper_lower = dist_upper_lower
x_lower_line = x_upper_line - dx_upper_lower
if (x(i).le.x_upper_line) then
! left of upper line
if (x(i).le.x_lower_line) then
! also left of lower line
z(i) = sp01
else
! left of upper line, right of lower line, i.e. on the slope:
dx = x(i) - x_lower_line
z(i) = sp01 + dx * tan(sp07r)
endif
else
! right of upper line
z(i) = sp02
endif
endif
! take care of the corner ...
if (sp08.lt.0.0) then
endif
if (z(i).lt.0.0) then
call stop_run ('z less than zero$')
elseif (z(i).gt.1.0) then
call stop_run ('z larger than one$')
endif
enddo
case (17)
! basin with differing height basement steps, geometry constant in y
! sp01 is the z level at the bottom of the basin
! sp02 and 03 are x1,x2
! sp04 is not used (this could specify a 2nd slope later)
! sp05 is the distal basement step thickness
! sp06 is the proximal basement step thickness
! sp07 is the slope
m=tan(sp07*pi/180.d0)
m=sign(m,-sp06)
x1a=sp02+(sp06/(2.d0*m))
x1b=sp02-(sp06/(2.d0*m))
x2a=sp03+(sp05/(2.d0*m))
x2b=sp03-(sp05/(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
z(i)=sp01+sp06
elseif (x(i).le.x1b) then
z(i)=(sp01+sp06)-(((x(i)-x1a)/(-sp06/m))*sp06)
elseif (x(i).le.x2a) then
z(i)=sp01
elseif (x(i).le.x2b) then
z(i)=(sp01+sp05)-((1.d0-((x(i)-x2a)/(-sp05/m)))*sp05)
!zx=sp01-((1.d0-((x(i)-x2a)/(sp06/m)))*sp06)
else
z(i)=sp01+sp05
endif
!! There is a smarter way to do this, but...
!if (sp06.gt.eps) then
! if (y(i).le.y1a) then
! z(i)=sp01
! elseif (y(i).le.y1b) then
! z(i)=max(sp01-(((y(i)-y1a)/(-sp06/m))*sp06),zx)
! elseif (y(i).le.y2a) then
! z(i)=max(sp01-sp06,zx)
! elseif (y(i).le.y2b) then
! z(i)=max(sp01-((1.d0-((y(i)-y2a)/(-sp06/m)))*sp06),zx)
! else
! z(i)=sp01
! endif
!else
! 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
!endif
end do
case (18)
! JA 050415. Sinusoidally-perturbed surface, to match aggradation options
! sp01 is the z level at the base of the sinusoid
! sp02 and 03 are x1,x2
! sp04 05 are y1, y2
! sp06 is the wavelength in x
! sp07 is the amplitude x
! sp08 is the wavelength in y
! sp09 is the amplitude in y
do i=1,ns
if (x(i).ge.sp02.and.x(i).le.sp03.and.y(i).ge.sp04.and.y(i).le.sp05) then
if (sp09.lt.eps) then
z(i)=sp01+sp07*(1+sin((x(i)-sp02)*3.14159/(sp06/2)-3.14159/2))
else
z(i)=sp01+sp07*(1+sin((x(i)-sp02)*3.14159/(sp06/2)-3.14159/2))*sp09*(1+sin((y(i)-sp04)*3.14159/(sp08/2)-3.14159/2))
endif
endif
end do
case (19)
! JA 070415. half-Gaussian surface, with the option of a Gaussian deformation
! in y. Matches progradation options
! sp01 and 02 are the z levels at the landward and seaward edges of the progradation profile
! sp03 and 04 are x1,x2
! sp05 is the half width of the Gaussian in x
! sp06 sets whether there will be a Gaussian in y (if 1)
! sp07 is the centre position of the Gaussian in y
! sp08 is the half width of the Gaussian in y
do i=1,ns
if (x(i).le.sp03) then
if(sp06.gt.0) then
z(i)=sp02+(sp01-sp02)*exp(-((y(i)-sp07)**2)/(2*sp08**2))
else
z(i)=sp01
endif
elseif(x(i).ge.sp03.and.x(i).le.sp04) then
if(sp06.gt.0) then
z(i)=sp02+(sp01-sp02)*exp(-((x(i)-sp03)**2)/(2*sp05**2))*exp(-((y(i)-sp07)**2)/(2*sp08**2))
else
z(i)=sp02+(sp01-sp02)*exp(-((x(i)-sp03)**2)/(2*sp05**2))
endif
elseif (x(i).ge.sp04) then
z(i)=sp02
endif
end do
case default
call stop_run ('surface type not defined$')
......
src/sedimentation.f90
+ 7
− 4
View file @ cd2fa96e
......@@ -126,17 +126,20 @@ baselevelx1=params%baselevelx1
baselevely0=params%baselevely0
baselevely1=params%baselevely1
! JA 070415. Amplitudes now account for squishing factor, which will match the initial
! surface definition. Aggradation rate also includes squishing factor
sedimentation_type=params%sedimentation_type
zaggrade_init=params%zaggrade_init
x_agg_start=params%x_agg_start
x_agg_end=params%x_agg_end
y_agg_start=params%y_agg_start
y_agg_end=params%y_agg_end
x_agg_sinus_amp=params%x_agg_sinus_amp
x_agg_sinus_amp=params%x_agg_sinus_amp*params%vex
x_agg_sinus_wavelth=params%x_agg_sinus_wavelth
y_agg_sinus_amp=params%y_agg_sinus_amp
y_agg_sinus_amp=params%y_agg_sinus_amp*params%vex
y_agg_sinus_wavelth=params%y_agg_sinus_wavelth
aggrade_rate=params%aggrade_rate
aggrade_rate=params%aggrade_rate*params%vex
z_prog_init=params%z_prog_init
z_prog_fin=params%z_prog_fin
x_prog_start=params%x_prog_start
......@@ -293,7 +296,7 @@ select case(sedimentation_type)
case (2)
i=1
do i=1,surface%nsurface
if (surface%x(i).gt.x_agg_start.and.surface%x(i).lt.x_agg_end.and. surface%y(i).gt.y_agg_start.and.surface%y(i).lt.y_agg_end) then
if (surface%x(i).ge.x_agg_start.and.surface%x(i).le.x_agg_end.and. surface%y(i).ge.y_agg_start.and.surface%y(i).le.y_agg_end) then
if (y_agg_sinus_amp.lt.eps) then
z_compare=zaggrade_current+x_agg_sinus_amp*(1+sin((surface%x(i)-x_agg_start)*3.14159/(x_agg_sinus_wavelth/2)-3.14159/2))
else
......
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