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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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,weight,eps
double precision :: dgrain,dfluid,spor,ccoef
integer :: matnum,i,j,k,free_surf_found,err,size_str
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
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! 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)
Dave Whipp
committed
x = 1.d0/(2.d0**dble(params%levelmax_oct))*(0.5d0+dble(istr-1))
y = 1.d0/(2.d0**dble(params%levelmax_oct))*(0.5d0+dble(jstr-1))
! JANICE: Made some changes below. Check them over in case they differ from what
! you intended to do here.
allocate (density_temp%density(size_str),stat=err)
!density_temp=0
density_temp%density=0.d0
!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
! Calculate current vertical position
Dave Whipp
committed
! Note we're now accounting for squishing here.
z=(1.d0/(2.d0**dble(params%levelmax_oct))*params%vex)/dble(params%dstrnz)*(0.5d0+dble(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
Dave Whipp
committed
! unsquished z into octree_interpolate
call octree_interpolate (osolve%octree,osolve%noctree,osolve%icon, &
Dave Whipp
committed
x,y,z/params%vex,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
dfluid=mat(matnum)%pore_fluid_density
spor=mat(matnum)%surface_porosity
ccoef=mat(matnum)%compaction_coef
! JANICE: More changes below
!density_temp(i)=dgrain-(dgrain-dfluid)*spor*exp(-ccoef*(z_free_surf-z))
density_temp%density(i)=dgrain-(dgrain-dfluid)*spor*exp(-ccoef*(z_free_surf-z))
else
! JANICE: More changes below
!density_temp(i)=mat(matnum)%density
density_temp%density(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 (initialize.gt.0) then
do i=1,size_str
! JANICE: More changes below
!density_profile(i)=density_temp(i)
density_profile%density(i)=density_temp%density(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
! JANICE: More changes below
!if (density_temp(i).gt.density_profile(i)) then
if (density_temp%density(i).gt.density_profile%density(i)) then
!density_profile(i)=density_temp(i)
density_profile%density(i)=density_temp%density(i)
! JANICE: I changed the line below from the original version you coded
! Did you want to test the value was exactly equal to zero, because this
! should only be done as coded for integer values. I presume
!elseif (density_temp(i)==0) then
elseif (density_temp%density(i) < eps) then
!density_profile(i)=density_temp(i)
density_profile%density(i)=density_temp%density(i)
endif
enddo
endif
! Need to first deallocate density_temp%density, then density_temp?
! JANICE: More changes below
!deallocate (density_temp)
deallocate (density_temp%density)
end subroutine calculate_density_profile
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------