Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DEFINE_ISOSTASY_BC OCT. 2009 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine define_isostasy_bc(params,osolve,vo,zisodisp,firstcall,l,x0,spu,spv,&
ldisp)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! New subroutine to apply the modified isostasy velocity boundary conditions
! - This routine uses the basal displacement array defined in isostasy.f90 to
! modify the basal velocity boundary conditions to have the incoming flux
! velocity tangent to the moho. This might not be exactly correct, but works
! for now.
!
! dwhipp (11/09)
!
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
implicit none
type (parameters) params
type (octreesolve) osolve
type (void) vo
double precision zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1)
double precision x0(osolve%nnode),ldisp(osolve%nnode)
double precision l,spu,spv
logical firstcall
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i,j,iproc,nproc,ierr,nb,xdisp,ydisp,xnow,ynow,ie,jp
double precision eps,pi,zsl,dxy,zdisp,vinit
double precision,dimension(:,:),allocatable::zisoslx,zisosly
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
! General variable definition
eps=1.d-10 ! Tiny number
pi=atan(1.d0)*4.d0 ! Pi
nb=2**params%levelmax_oct ! Number of nodes along one side of model
dxy=1./real(nb) ! Non-dimensional node spacing
if (firstcall) then
do i=1,osolve%nnode
! Shift s-line depth
xdisp=nint(x0(i)*nb)+1
ydisp=nint(osolve%y(i)*nb)+1
ldisp(i)=l+zisodisp(xdisp,ydisp)
enddo
else
allocate(zisoslx(nb+1,nb+1))!,zisosly(nb+1,nb+1))
do j=1,nb+1
do i=1,nb+1
if (i==1) then
zisoslx(i,j)=(zisodisp(i+1,j)-zisodisp(i,j))/dxy ! Displacement between first and second nodes along x-axis
elseif (i==nb+1) then
zisoslx(i,j)=(zisodisp(i,j)-zisodisp(i-1,j))/dxy ! Displacement between second-to-last and last nodes along x-axis
else
zisoslx(i,j)=((zisodisp(i,j)-zisodisp(i-1,j))/dxy+& ! Average displacement of nodes on either side of current node, along x-axis
(zisodisp(i+1,j)-zisodisp(i,j))/dxy)/2.d0
endif
!if (j==1) then
! zisosly(i,j)=(zisodisp(i,j+1)-zisodisp(i,j))/dxy ! Displacement between first and second nodes along y-axis
!elseif (j==nb+1) then
! zisosly(i,j)=(zisodisp(i,j)-zisodisp(i,j-1))/dxy ! Displacement between second-to-last and last nodes along y-axis
!else
! zisosly(i,j)=((zisodisp(i,j)-zisodisp(i,j-1))/dxy+& ! Average displacement of nodes on either side of current node, along y-axis
! (zisodisp(i,j+1)-zisodisp(i,j))/dxy)/2.d0
!endif
enddo
enddo
do i=1,osolve%nnode
xnow=nint(osolve%x(i)*nb)+1
ynow=nint(osolve%y(i)*nb)+1
zsl=zisoslx(xnow,ynow)
! Modify velocities
if (osolve%z(i).lt.eps) then
if (spu.le.1.d0-eps) then ! No change in B/Cs
zdisp=osolve%u(i)*sin(atan(zsl)) ! Negative for negative slope
vinit=(osolve%u(i)**2.d0+osolve%w(i)**2.d0)**0.5d0 ! Incoming velocity vector magnitude
osolve%w(i)=osolve%w(i)+zdisp ! New z-velocity
if (osolve%u(i).ge.0.d0) then
osolve%u(i)=sqrt(vinit**2-osolve%w(i)**2.d0) ! New x-velocity
else
osolve%u(i)=-sqrt(vinit**2-osolve%w(i)**2.d0) ! New x-velocity
endif
else
if (spu.ge.eps) then
zdisp=osolve%u(i)*sin(atan(zsl)) ! Negative for negative slope
vinit=(osolve%u(i)**2.d0+osolve%w(i)**2.d0)**0.5d0 ! Incoming velocity vector magnitude
osolve%w(i)=osolve%w(i)+zdisp ! New z-velocity
if (osolve%u(i).ge.0.d0) then
osolve%u(i)=sqrt(vinit**2-osolve%w(i)**2.d0) ! New x-velocity
else
osolve%u(i)=-sqrt(vinit**2-osolve%w(i)**2.d0) ! New x-velocity
endif
endif
endif
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
endif
enddo
deallocate(zisoslx)
endif
! here we need to impose that the bc satisfy the bad faces constrains
!do ie=1,osolve%nface
! do j=1,4
! jp=1+mod(j,4)
! i=osolve%iface(j+4,ie)
! if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=(osolve%u(osolve%iface(j,ie))+&
! osolve%u(osolve%iface(jp,ie)))/2.d0
! if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=(osolve%v(osolve%iface(j,ie))+&
! osolve%v(osolve%iface(jp,ie)))/2.d0
! if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=(osolve%w(osolve%iface(j,ie))+&
! osolve%w(osolve%iface(jp,ie)))/2.d0
! if (osolve%kfixt(i).eq.1) osolve%temp(i)=(osolve%temp(osolve%iface(j,ie))+&
! osolve%temp(osolve%iface(jp,ie)))/2.d0
! enddo
! i=osolve%iface(9,ie)
! if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=&
! sum(osolve%u(osolve%iface(1:4,ie)))/4.d0
! if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=&
! sum(osolve%v(osolve%iface(1:4,ie)))/4.d0
! if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=&
! sum(osolve%w(osolve%iface(1:4,ie)))/4.d0
! if (osolve%kfixt(i).eq.1) osolve%temp(i)=&
! sum(osolve%temp(osolve%iface(1:4,ie)))/4.d0
!enddo
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|