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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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 (osolve%x(i).le.x0) then
! if (spu.le.1.d0-eps) then ! No change in B/Cs
if (osolve%u(i).gt.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
osolve%u(i)=sqrt(vinit**2-osolve%w(i)**2.d0) ! New x-velocity
endif
! endif
! else
! if (spu.ge.eps) then
! zdisp=osolve%u(i)*sin(atan(zsl)) ! Negative for negative slope
! osolve%w(i)=zdisp-osolve%w(i) ! New z-velocity
! osolve%u(i)=sqrt(vin**2-osolve%w(i)**2.d0) ! New x-velocity
! endif
! endif
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|