-
Notifications
You must be signed in to change notification settings - Fork 0
/
fl_move.f90
299 lines (247 loc) · 8.72 KB
/
fl_move.f90
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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
! Move grid and adjust stresses due to rotation
subroutine fl_move
use arrays
include 'precision.inc'
include 'params.inc'
include 'arrays.inc'
! Move Grid
if (movegrid .eq. 0) return
! UPDATING COORDINATES
!$OMP parallel private(i)
!$OMP do
do i = 1,nx
! write(*,*) cord(j,i,1),cord(j,i,2),vel(j,i,1),vel(j,i,2),dt
cord(:,i,1) = cord(:,i,1) + vel(:,i,1)*dt
cord(:,i,2) = cord(:,i,2) + vel(:,i,2)*dt
! write(*,*) cord(j,i,1),cord(j,i,2)
enddo
!$OMP end do
!$OMP end parallel
! Diffuse topography
if( topo_kappa.gt.0.) call diff_topo
!$OMP parallel private(i,j,x1,y1,x2,y2,x3,y3,x4,y4, &
!$OMP vx1,vy1,vx2,vy2,vx3,vy3,vx4,vy4, &
!$OMP det,dw12,s11,s22,s12)
!$OMP do
!--- Adjusting Stresses And Updating Areas Of Elements
do i = 1,nx-1
do j = 1,nz-1
! Coordinates
x1 = cord (j ,i ,1)
y1 = cord (j ,i ,2)
x2 = cord (j+1,i ,1)
y2 = cord (j+1,i ,2)
x3 = cord (j ,i+1,1)
y3 = cord (j ,i+1,2)
x4 = cord (j+1,i+1,1)
y4 = cord (j+1,i+1,2)
! Velocities
vx1 = vel (j ,i ,1)
vy1 = vel (j ,i ,2)
vx2 = vel (j+1,i ,1)
vy2 = vel (j+1,i ,2)
vx3 = vel (j ,i+1,1)
vy3 = vel (j ,i+1,2)
vx4 = vel (j+1,i+1,1)
vy4 = vel (j+1,i+1,2)
! (1) Element A:
det=((x2*y3-y2*x3)-(x1*y3-y1*x3)+(x1*y2-y1*x2))
dvol(j,i,1) = det*area(j,i,1) - 1
area(j,i,1) = 1./det
! Adjusting stresses due to rotation
dw12 = 0.5*(vx1*(x3-x2)+vx2*(x1-x3)+vx3*(x2-x1) - &
vy1*(y2-y3)-vy2*(y3-y1)-vy3*(y1-y2))/det*dt
s11 = stress0(j,i,1,1)
s22 = stress0(j,i,2,1)
s12 = stress0(j,i,3,1)
stress0(j,i,1,1) = s11 + s12*2.*dw12
stress0(j,i,2,1) = s22 - s12*2.*dw12
stress0(j,i,3,1) = s12 + dw12*(s22-s11)
! rotate strains
s11 = strain(j,i,1)
s22 = strain(j,i,2)
s12 = strain(j,i,3)
strain(j,i,1) = s11 + s12*2.*dw12
strain(j,i,2) = s22 - s12*2.*dw12
strain(j,i,3) = s12 + dw12*(s22-s11)
! (2) Element B:
det=((x2*y4-y2*x4)-(x3*y4-y3*x4)+(x3*y2-y3*x2))
dvol(j,i,2) = det*area(j,i,2) - 1
area(j,i,2) = 1./det
! Adjusting stresses due to rotation
dw12 = 0.5*(vx3*(x4-x2)+vx2*(x3-x4)+vx4*(x2-x3) - &
vy3*(y2-y4)-vy2*(y4-y3)-vy4*(y3-y2))/det*dt
s11 = stress0(j,i,1,2)
s22 = stress0(j,i,2,2)
s12 = stress0(j,i,3,2)
stress0(j,i,1,2) = s11 + s12*2.*dw12
stress0(j,i,2,2) = s22 - s12*2.*dw12
stress0(j,i,3,2) = s12 + dw12*(s22-s11)
! (3) Element C:
det=((x2*y4-y2*x4)-(x1*y4-y1*x4)+(x1*y2-y1*x2))
dvol(j,i,3) = det*area(j,i,3) - 1
area(j,i,3) = 1./det
! Adjusting stresses due to rotation
dw12 = 0.5*(vx1*(x4-x2)+vx2*(x1-x4)+vx4*(x2-x1) - &
vy1*(y2-y4)-vy2*(y4-y1)-vy4*(y1-y2))/det*dt
s11 = stress0(j,i,1,3)
s22 = stress0(j,i,2,3)
s12 = stress0(j,i,3,3)
stress0(j,i,1,3) = s11 + s12*2.*dw12
stress0(j,i,2,3) = s22 - s12*2.*dw12
stress0(j,i,3,3) = s12 + dw12*(s22-s11)
! (4) Element D:
det=((x4*y3-y4*x3)-(x1*y3-y1*x3)+(x1*y4-y1*x4))
dvol(j,i,4) = det*area(j,i,4) - 1
area(j,i,4) = 1./det
! Adjusting stresses due to rotation
dw12 = 0.5*(vx1*(x3-x4)+vx4*(x1-x3)+vx3*(x4-x1) - &
vy1*(y4-y3)-vy4*(y3-y1)-vy3*(y1-y4))/det*dt
s11 = stress0(j,i,1,4)
s22 = stress0(j,i,2,4)
s12 = stress0(j,i,3,4)
stress0(j,i,1,4) = s11 + s12*2.*dw12
stress0(j,i,2,4) = s22 - s12*2.*dw12
stress0(j,i,3,4) = s12 + dw12*(s22-s11)
enddo
enddo
!$OMP end do
!$OMP end parallel
return
end subroutine fl_move
!============================================================
! Diffuse topography
!============================================================
subroutine diff_topo
use arrays
include 'precision.inc'
include 'params.inc'
include 'arrays.inc'
dimension tkappa(mnx+1)
!EROSION PROCESSES
if( topo_kappa .gt. 0. ) then
topomean = sum(cord(1,:,2)) / nx
tkappa = topo_kappa
! higher elevation has higher erosion rate
do i = 1, nx
if (cord(1,i,1) > topomean) tkappa(i) = topo_kappa * (1 + (cord(1,i,1) - topomean) * fac_kappa)
enddo
do i = 2, nx-1
snder = ( tkappa(i+1)*(cord(1,i+1,2)-cord(1,i ,2))/(cord(1,i+1,1)-cord(1,i ,1)) - &
tkappa(i-1)*(cord(1,i ,2)-cord(1,i-1,2))/(cord(1,i ,1)-cord(1,i-1,1)) ) / &
(cord(1,i+1,1)-cord(1,i-1,1))
dtopo(i) = dt * snder
end do
dtopo(1) = dtopo(2)
dtopo(nx) = dtopo(nx-1)
cord(1,1:nx,2) = cord(1,1:nx,2) + dtopo(1:nx)
! accumulated topo change since last resurface
dhacc(1:nx-1) = dhacc(1:nx-1) + 0.5 * (dtopo(1:nx-1) + dtopo(2:nx))
! adjust markers
if(mod(nloop, 100) .eq. 0) then
!!$ print *, 'max sed/erosion rate (m/yr):' &
!!$ , maxval(dtopo(1:nx)) * 3.16e7 / dt &
!!$ , minval(dtopo(1:nx)) * 3.16e7 / dt
call resurface
end if
endif
! magma extrusion
if ( .true. ) then
! grid spacing in x
extrusion(1:nx-1) = cord(1,2:nx,1) - cord(1,1:nx-1,1)
! height of extrusion = volume of extrusion / grid spacing in x
extrusion(1:nx-1) = andesitic_melt_vol(1:nx-1) / extrusion(1:nx-1)
extr_acc(1:nx-1) = extr_acc(1:nx-1) + extrusion(1:nx-1)
cord(1,1:nx-1,2) = cord(1,1:nx-1,2) + extrusion(1:nx-1)
cord(1,2:nx ,2) = cord(1,2:nx ,2) + extrusion(1:nx-1)
endif
return
end subroutine diff_topo
subroutine resurface
use marker_data
use arrays
include 'precision.inc'
include 'params.inc'
include 'arrays.inc'
include 'phases.inc'
dimension shp2(2,3,2)
do i = 1, nx-1
! averge thickness of this element
elz = 0.5 * (cord(1,i,2) - cord(2,i,2) + cord(1,i+1,2) - cord(2,i+1,2))
! change in topo
chgtopo = dhacc(i)
! # of markers in this element
kinc = sum(nphase_counter(:,1,i))
if (abs(chgtopo*kinc) >= elz) then
! add/remove markers if topo changed too much
if (chgtopo > 0.) then
! sedimentation, add a sediment marker
!print *, 'add sediment', i, chgtopo, elz
!call add_marker_at_top(i, 0.05d0, elz, ksed2, nmarkers)
else
! erosion, remove the top marker
!print *, 'erosion', i, chgtopo, elz
ymax = -1e30
nmax = 0
kmax = 0
call shape_functions(1,i,shp2)
do k = 1, ntopmarker(i)
n = itopmarker(k, i)
ntriag = mark(n)%ntriag
m = mod(ntriag,2) + 1
call bar2xy(mark(n)%a1, mark(n)%a2, shp2(:,:,m), x, y)
if(ymax < y) then
ymax = y
nmax = n
kmax = k
endif
end do
if (nmax .ne. 0) then
mark(nmax)%dead = 0
! replace topmarker k with last topmarker
itopmarker(kmax,i) = itopmarker(ntopmarker(i),i)
ntopmarker(i) = ntopmarker(i) - 1
endif
endif
dhacc(i) = 0
! recalculate phase ratio
kinc = sum(nphase_counter(:,1,i))
phase_ratio(1:nphase,1,i) = nphase_counter(1:nphase,1,i) / float(kinc)
else
! nothing to do
end if
! change in topo
chgtopo2 = extr_acc(i)
if (chgtopo2*kinc >= elz) then
! add/remove markers if topo changed too much
! extrusion, add an arc marker
n_to_add = ceiling(chgtopo2 / elz * kinc)
dz = chgtopo2 / elz / (n_to_add+1)
!print *, 'add arc', i, chgtopo2, elz, n_to_add, dz
do ii = 1, n_to_add
call add_marker_at_top(i, dz*ii, elz, karc1, nmarkers)
enddo
extr_acc(i) = 0
! recalculate phase ratio
kinc = sum(nphase_counter(:,1,i))
phase_ratio(1:nphase,1,i) = nphase_counter(1:nphase,1,i) / float(kinc)
end if
end do
end subroutine resurface
subroutine add_marker_at_top(i, dz_ratio, elz, kph, nmarkers)
use arrays
include 'precision.inc'
do while(.true.)
call random_number(rx)
xx = cord(1,i,1) + rx * (cord(1,i+1,1) - cord(1,i,1))
yy = cord(1,i,2) + rx * (cord(1,i+1,2) - cord(1,i,2)) - dz_ratio*elz
call add_marker(xx, yy, kph, time, nmarkers, 1, i, inc)
if(inc==1) exit
!write(333,*) 'add_marker_at_top failed: ', xx, yy, rx, elz, kph
!write(333,*) ' ', cord(1,i,:)
!write(333,*) ' ', cord(1,i+1,:)
!write(333,*) ' ', cord(2,i,:)
!write(333,*) ' ', cord(2,i+1,:)
!call SysMsg('Cannot add marker.')
end do
end subroutine add_marker_at_top