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

Updated define_bc_segmented_s_line_parabola to work properly in all three... 

Updated define_bc_segmented_s_line_parabola to work properly in all three convergence segments and to use velocity gradient smoothing for the BCs
parent a9c0ed33
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/define_bc_segmented_s_line_parabola.f90
+ 77
− 66
View file @ 6e7cd838
......@@ -83,18 +83,7 @@ dxy=1.d0/2**(params%levelmax_oct+1.d0)
alpha=alpha*pi/180
! allocate(x0(osolve%nnode))
! do i=1,osolve%nnode
! if (osolve%y(i).le.ystart) then
! x0(i)=xstart
! else if (osolve%y(i).le.yend) then
! x0(i)=(osolve%y(i)-ystart)*tan(theta)+xstart
! else
! x0(i)=xend
! endif
! enddo
flytt=(kink2-kink1)* tan(alpha)
flytt=(kink2-kink1) * tan(alpha)
a = (base-l) /((endp-startp)*(endp-startp))
d = startp
......@@ -129,66 +118,88 @@ do i=1,osolve%nnode
osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
osolve%kfixt(i)=1 ; osolve%temp(i)=1.d0
! Velocity in first margin-normal convergence segment
if (osolve%y(i) .lt. kink1) then
if (osolve%x(i) .le. startp) then
osolve%u(i)=1.d0
elseif (osolve%x(i) .gt. startp .AND. osolve%x(i) .lt. endp) then
dipangle=atan(2*a*(osolve%x(i)-d))
! Slab dip angle
dipangle = atan(2*a*(osolve%x(i)-startp))
! Velocity before transition to subduction
if (osolve%x(i) .le. startp-real(nelemz)*dxy) then
osolve%u(i) = vin
! Velocity in transition to subduction
elseif (osolve%x(i) .le. startp+real(nelemz)*dxy) then
! Velocity contribution from incoming horizontal velocity
! (diminishes to zero across transition)
osolve%u(i) = vin * (1.d0 - (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
! Velocity contribution from subduction region
! (increases from zero to full magnitude across transition)
osolve%u(i) = osolve%u(i) + vmag*cos(dipangle) * (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
osolve%w(i) = vmag*sin(dipangle) * (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
! Velocity in subduction region
elseif (osolve%x(i) .le. endp-real(nelemz)*dxy) then
osolve%u(i) = vmag*cos(dipangle)
osolve%w(i) = vmag*sin(dipangle)
else
osolve%u(i)=0.d0
osolve%v(i)=0.d0
osolve%w(i)=0.d0
! Velocity in transition from subduction region
elseif (osolve%x(i) .le. endp+real(nelemz)*dxy) then
! Velocity contribution from subduction region
! (decreases to zero across transition)
osolve%u(i) = vmag*cos(dipangle) * (1.d0 - (osolve%x(i)-(endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
osolve%w(i) = vmag*sin(dipangle) * (1.d0 - (osolve%x(i)-(endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
endif
endif
if ( osolve%y(i) .ge. kink1 .AND. osolve%y(i) .lt. kink2 ) then
if (osolve%x(i).lt. (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1)) + startp)) then
osolve%u(i)=1.d0
elseif (osolve%x(i) .gt. (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1)) + (startp)) .and. osolve%x(i) .lt. (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1))) + endp) then
osolve%u(i) = -(a*sqrt(8*(osolve%x(i)- (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1))) )*d)+0.2)*vin/1
osolve%w(i) = 10*2*a*sqrt( (osolve%x(i)- (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1)))) * (osolve%x(i)- (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1)))) - d*d )*vin/30
else
osolve%u(i)=0.d0
! Velocity in oblique convergence segment
elseif (osolve%y(i) .lt. kink2) then
! Slab dip angle
dipangle = atan(2*a*(osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp)))
! Velocity before transition to subduction
if (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy) then
osolve%u(i) = vin
! Velocity in transition to subduction
elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+startp+real(nelemz)*dxy) then
! Velocity contribution from incoming horizontal velocity
! (diminishes to zero across transition)
osolve%u(i)= vin * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
! Velocity contribution from subduction region
! (increases from zero to full magnitude across transition)
osolve%u(i)= osolve%u(i) + vmag*cos(dipangle) * (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
osolve%w(i) = vmag*sin(dipangle) * (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
! Velocity in subduction region
elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy) then
osolve%u(i) = vmag*cos(dipangle)
osolve%w(i) = vmag*sin(dipangle)
! Velocity in transition from subduction region
elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+endp+real(nelemz)*dxy) then
! Velocity contribution from subduction region
! (decreases to zero across transition)
osolve%u(i) = vmag*cos(dipangle) * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
osolve%w(i) = vmag*sin(dipangle) * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
endif
! Velocity in second margin-normal convergence segment
else
! Slab dip angle
dipangle = atan(2*a*(osolve%x(i)-(flytt+startp)))
! Velocity before transition to subduction
if (osolve%x(i) .le. flytt+startp-real(nelemz)*dxy) then
osolve%u(i) = vin
! Velocity in transition to subduction
elseif (osolve%x(i) .le. flytt+startp+real(nelemz)*dxy) then
! Velocity contribution from incoming horizontal velocity
! (diminishes to zero across transition)
osolve%u(i) = vin * (1.d0 - (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
! Velocity contribution from subduction region
! (increases from zero to full magnitude across transition)
osolve%u(i) = osolve%u(i) + vmag * cos(dipangle) * (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
osolve%w(i) = vmag * sin(dipangle) * (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
! Velocity in subduction region
elseif (osolve%x(i) .le. flytt+endp-real(nelemz)*dxy) then
osolve%u(i) = vmag * cos(dipangle)
osolve%w(i) = vmag * sin(dipangle)
! Velocity in transition from subduction region
elseif (osolve%x(i) .le. flytt+endp+real(nelemz)*dxy) then
! Velocity contribution from subduction region
! (decreases to zero across transition)
osolve%u(i) = vmag * cos(dipangle) * (1.d0 - (osolve%x(i)-((flytt+endp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
osolve%w(i) = vmag * sin(dipangle) * (1.d0 - (osolve%x(i)-((flytt+endp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
endif
endif
if ( osolve%y(i) .ge. kink2 .AND. osolve%y(i) .lt. (1.0-eps) ) then
if (osolve%x(i) .le. (startp + flytt)) then
osolve%u(i)=1.d0
! elseif (osolve%x(i) .gt. (startp+flytt) .AND. osolve%x(i) .lt. (startp+flytt+0.01)) then
!! osolve%w(i)= (-2*a*(base-d)*vin+(osolve%x(i)-flytt))
! osolve%w(i)= (-2*a*(base-d)*vin)/(2.9d0)
! osolve%u(i)= (10*vin*(osolve%x(i)-flytt))/(2.9d0)
! elseif (osolve%x(i) .gt. (startp+flytt+0.01) .AND. osolve%x(i) .lt. (startp+flytt+0.04)) then
! osolve%w(i)= (-2*a*(base-d)*vin)/(2.9d0)
! osolve%u(i)= (6*vin*(osolve%x(i)-flytt))/(2.9d0)
! elseif (osolve%x(i) .gt. (startp+flytt+0.04) .AND. osolve%x(i) .lt. (startp+flytt+0.05)) then
! osolve%w(i)= (-2*a*(base-d)*vin)/(2.9d0)
! osolve%u(i)= (4*vin*(osolve%x(i)-flytt))/(2.9d0)
! elseif (osolve%x(i) .gt. (startp+flytt+0.05) .AND. osolve%x(i) .lt. (endp+flytt-0.05)) then
! osolve%w(i)= (-2*a*(base-d)*vin)/(2.9d0)
! osolve%u(i)= (2*vin*(osolve%x(i)-flytt))/(2.9d0)
elseif (osolve%x(i) .gt. (startp+flytt) .AND. osolve%x(i) .lt. (endp+flytt)) then
! osolve%w(i)= (2*a*(osolve%x(i)-(d+flytt)))*vin
! osolve%u(i)= (2*vin*(osolve%x(i)-flytt))
osolve%u(i) = -(a*sqrt(8*(osolve%x(i)-flytt)*d)+0.2)*vin/1
osolve%w(i) = 10*2*a*sqrt( (osolve%x(i)-flytt)*(osolve%x(i)-flytt) - d*d )*vin/30
else
osolve%u(i)=0.d0
osolve%v(i)=0.d0
osolve%w(i)=0.d0
endif
endif
endif
if (osolve%z(i).gt.1.d0-eps) then
......
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