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
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! CALCULATE DENSITY PROFILE MARCH 2013 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine calculate_density_profile(params,density_profile,osolve,mat, &
initialize,istr,jstr)
use definitions
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine calculates the density profile for a particular location on the |
! base of the model, using elemental material number to identify cells which |
! are/are not compactible. For compactible cells, an Athy curve is used to |
! determine density profile with depth.
!
! A temporary array of size matching the input density_str instance is used to |
! store the new density profile, at the same vertical devisions as the existing|
! profile. In cases where the new density is greater than |
! density_profile(i)%density (the previous stored current density) or where the|
! new density is 0 (when erosion has occurred), density_profile(i)%density is |
! replaced. |
!
! This routine works from the top of the string downwards, storing total |
! weight above a particular vertical location in the string, to inform Athy |
! calculation.
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
!include 'mpif.h'
type (parameters) params
type (string) density_profile
! How to specify a single instance of density_str, not the full structure?
! Note that we really only need density here, not densityp, so we could
! just past that, if easier?
! We also need material number for each element, and compaction parameters
! for compactible material types.
! We also need LSFs, to check if we're above/below the free surface.
! Suggest for now we consider only the free surface. This could be
! adapted to look for other LSFs, to better define transitions between
! material types.
type(octreesolve) osolve
type (material) mat(0:params%nmat)
integer :: initialize
double precision :: istr,jstr
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
!integer :: i,mat_maxi!,iproc,nproc,ierr
!logical :: mat_max_tie
type (string) density_temp
! How to specify that this is allocatable??
! We will need an allocatable array, density_temp, after the style of
! density_profile (an instance of density_str)
double precision :: cur_lsf,x,y,z,z_free_surf,
double precision :: dgrain,dfluid,spor,ccoef
integer :: matnum,i,j,k, free_surf_found
integer,dimension(:) :: size_str
double precision :: weight
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! For the given input instance of density_str (locally stored as density_profile):
! 1. Allocate density_temp with dimensions after density_profile.
! 2. Starting with the top (end?) of density_temp, move down the string, at each
! location:
! a. Is this location above the free surface? Is so density_temp%density=0.
! If not, continue.
! b. Is the element cut?
! - if so, use LSFs to determine what material number should be used
! for the present location on the string.
! - if not:
! i. Query which element we're in.
! ii. Determine the element's material number.
! c. Access compaction properties for this material number.
! d. Determine density at this location.
! e. Calculate weight contribution from this division of the string
! (density*string division length*surface area for one string).
! - add weight contribution to weight tally for this string.
! f. Compare denity_profile%density(i) with density_temp%density(i),
! overwrite density_profile%density(i) if:
! - density_temp%density(i)>density_profile%density(i) OR
! - density_temp%density(i)=0
! 3. Deallocate density_temp
!( Returns adjusted density_profile to compaction)
! Note: if using this routine to construct the original density profiles,
! could revise step f to always overprint density_profile, for all locations
! on the string.
eps = 1.d-10
size_str=size(density_profile%density)
x = 1/(2.d0**params%levelmax_oct)*(0.5+istr-1)
y = 1/(2.d0**params%levelmax_oct)*(0.5+jstr-1)
allocate (density_temp%density(size_str(1)),stat=err)
density_temp=0
!weight=0
! Loop over all locations in density_profile%density, starting at the
! uppermost (top of model) location.
! Will the array locations in density_profile%density start at 0 or 1?
! The following assumes 1; adjust if needed
free_surf_found=0
matnum=0
do i=0,size_str(1)-1
j=size_str(1)-i
! Calculate current vertical position
z=(1/(2.d0**params%levelmax_oct))/params%dstrnz*(0.5+j-1)
if (initialize.gt.0) then
! Use LSFs to determine material type
do k=1,osolve%nlsf
call octree_interpolate (osolve%octree,osolve%noctree,osolve%icon, &
osolve%nleaves,osolve%lsf(:,i),osolve%nnode, &
x,y,z,cur_lsf)
if (cur_lsf < eps) then
matnum = k
if (k=1.and.free_surf_found=0) then
z_free_surf=(1/(2.d0**params%levelmax_oct))/params%dstrnz*(0.5+j-1)
free_surf_found=1.
endif
endif
enddo
! Update material numbers for material changes
call mat_trans (params,mat,matnum,x,y,z)
endif
! Assign density for this vertical location
if (mat(matnum)%compactible) then
! Calculate density according to Athy curve.
! Note: for now this assumes that all compactible material is in a
! continuous stretch from the top of the string, and only one
! compactible material will be encountered.
dgrain=mat(matnum)%grain_density
dfliud=mat(matnum)%fluid%density
spor=mat(matnum)%surface%porosity
ccoef=mat(matnum)%compaction_coef
density_temp(j)=dgrain-(dgrain-dfluid)*spor*exp(-ccoef*(z_free_surf-z))
else
density_temp(j)=mat(matnum)%density
endif
! Increment weight above current position
!weight=weight+density_temp(j)*(1/(2.d0**params%levelmax_oct))/params%dstrnz*(1/(2.d0**params%levelmax_oct))**2
!Identify what element we're in
!if density_temp(j) is above the free surface then
!density_temp(j)=0
!else
!if the element we're in is cut then
! If we know what element we're in, can we check for cut element
! like done in make_cut?
!do i=1,nlsf
! prod=lsf(1,i)
! do k=2,params%mpe
! if (prod*lsf(k,i).le.0.d0) then
!if we're above the relevant LSF
!matnum=materialn---> for above this LSF
!else
!matnum=materialn---> for below this LSF
!endif
! endif
! enddo
!endo
!else
!matnum=material number for the element we're in
!endif
!if (compactible for this material) then
!density_temp(j)---> get from compaction properties for this material,
! weight of the column above us
!else
!density_temp(j)=density for this material
!endif
!endif
!weight=weight+density_temp(j)*STRING INCREMENT LENGTH*STRING AREA
enddo
if(initalize.gt.0)
do i=1,size_str(1)
density_profile(i)=density_temp(i)
enddo
!else
! do i=1,size_str(1)
! if (density_temp(i).gt.density_profile(i) then
! density_profile(i)=density_temp(i)
! elseif (density_temp(i)==0)
! density_profile(i)=density_temp(i)
! endif
! enddo
endif
! Need to first deallocate density_temp%density, then density_temp?
deallocate (density_temp)
end subroutine calculate_density_profile
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------