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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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
! NOTE: as currently coded, this subroutine uses LSFs to determine material type.
! In the future, this should be done only when initializing the density profile.
! In subsequent calls, this subroutine should check the query the elemental
! material number.
! NOTE: Once the option to not use LSFs to query material number has been implemented,
! this routine should also check to see if the element corresponding to the current
! location on the string is cut. If so, use the LSF to determine material number,
! if not, query the elemental material number.
free_surf_found=0
matnum=0
do i=size_str(1),1,-1
! Calculate current vertical position
! Note we're not accounting for squishing here...
z=(1/(2.d0**params%levelmax_oct))/params%dstrnz*(0.5+i-1)
! If this call is initializing the density profiles, use LSFs to determine
! material number
! NOTE: for now this will calculate the temporary density profile based on the LSFs,
! whether this is an initialization or not.
! 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=z
free_surf_found=1
endif
endif
enddo
! Update material numbers for material changes
call mat_trans (params,mat,matnum,x,y,z)
! 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(i)=dgrain-(dgrain-dfluid)*spor*exp(-ccoef*(z_free_surf-z))
else
density_temp(i)=mat(matnum)%density
endif
! Increment weight above current position
!weight=weight+density_temp(i)*(1/(2.d0**params%levelmax_oct))/params%dstrnz*(1/(2.d0**params%levelmax_oct))**2
enddo
! If initializing the density strings, copy all locations from the temporary array to
! to the current density profile.
if(initalize.gt.0)
do i=1,size_str(1)
density_profile(i)=density_temp(i)
enddo
! If not initializing the density strings, check to see whether the density string should
! be updated at each location. Update only in cases where material has been eroded
! (density goes to 0), or where density at a given position on the string has increased.
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
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------