Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! COMPACTION MARCH 2013 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine compaction(params,density_str,osolve,ov,mat,dtcur)
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
use definitions
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine compares density profiles, for each possible element centre |
! location at the base of the model, between the start of the time step (past |
! density) and the current mechanical solution (current density). Compaction |
! is called after erosion and sedimentation have been applied. This subroutine|
! determines the velocity field necessary to account for increased density in |
! compactible materials, compared to the start of time step density |
! configuration. Velocities are stores on ____, to be added to the nodal |
! velocities used to move surfaces and advect the cloud. Current density is |
! updated with each call to compaction.
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
!include 'mpif.h'
type (parameters) params
type (string) :: density_str(2**params%levelmax_oc,2**params%levelmax_oct)
! In how much detail do we need to specify density_str? We will not be
! allocating/deallocating and of the density_str entires in compaction or
! calculate_density_profile.
! Note we are passing in the full structure here.
! 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(octreev) ov
type (material) mat(0:params%nmat)
double precision :: dtcur
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
!integer :: i,mat_maxi!,iproc,nproc,ierr
!logical :: mat_max_tie
double precision,dimension(:),allocatable :: str_velo,countnode
double precision :: vtot,xstr,ystr,zstr,x0,y0,z0,dxyz,zmax_elem,str_velo
double precision :: str_length_inc
integer :: i,j,k,size_str(2),leaf,level,loc
logical :: next_elem
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! Loop over all instances of density_str (i,j 1:2**params%levelmax_oct), for each
! string:
! 1. Call calculate_density_profile(density_str(i,j), ??)
! - this will also need material number for each element, and compaction
! parameters for all materials (some of which may not be compactible)
! 2. Compare density_str(i,j)%density and %densityp
! - look at each point on density_str(i,j), moving from base of model
! to intersection with free surface (found when density==0+/-smallnumber)
! 3. If density>densityp, calculate velocity required to account for compaction
! 4. Interpolate velocities from density_str locations onto element nodes.
! - easier to first interpolate onto element centers, looking over all
! elements?
! - will need to account for the number of velocity contributions to each
! node
! Loop through all instances of density_str (all possible element center
! locations at the base of the model
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
! Zero out existing osolve%compaction_density
osolve%compaction_density=0.d0
! Could calculate this before calc_density_profile and pass in
size_str=size(density_str(1,1)%density)
! Allocate str_velo
allocate(str_velo(size_str(1)))
! Allocate countnode
allocate (countnode(osolve%nnode))
countnode=0
do i=1,osolve%nleaves
do j=5,8
k = osolve%icon(j,i)
countnode(k) = countnode(k) + 1
enddo
enddo
! String length increment
str_length_inc = 1/(2.d0**dble(params%levelmax_oct))/dble(params%dstrnz)
do j=1,2**params%levelmax_oct
do i=1,2**params%levelmax_oct
!Update current density, density_str(i,j)%density
call calculate_density_profile(params,density_str(i,j),osolve,mat,0,i,j)
! Could calculate this before calc_density_profile and pass in
xstr = 1/(2.d0**params%levelmax_oct)*(0.5+i-1)
ystr = 1/(2.d0**params%levelmax_oct)*(0.5+j-1)
next_elem = .true.
str_velo = 0.d0
do k=1,size_str(1)
! Note we're not accounting for squishing here...
zstr = 1/(2.d0**params%levelmax_oct)*(0.5+k-1)/params%dstrnz
if (next_elem) then
call octree_find_leaf (osolve%octree,osolve%noctree,xstr,ystr,zstr, &
leaf,level,loc,x0,y0,z0,dxyz)
nstr_per_elem = int(2.d0**(2.d0*dble(params%levelmax_oct - level)))
nstr_pts_per_elem = int(2.d0**(dble(params%levelmax_oct - level)) * &
dble(nstr_per_elem) * dble(params%dstrnz))
zmax_elem = z0 + dxyz / 2.d0
next_elem = .false.
endif
osolve%compaction_density(leaf) = osolve%compaction_density(leaf) + density_str(i,j)%density(k)/dble(nstr_pts_per_elem)
str_velo(k) = str_length_inc * ((density_str(i,j)%densityp(k)/density_str(i,j)%density(k))-1.d0)/dtcur
if (k > 1) str_velo(k) = str_velo(k) + str_velo(k-1)
if (zmax_elem - zstr .le. str_length_inc) then
do inode=5,8
! This is OK for a uniform octree at octreelvl_max, but the velocites
! should be scaled by distance to the nodes for elements with multiple
! strings!
ov%wnodecompact(ov%icon(inode,leaf)) = ov%wnodecompact(ov%icon(inode,leaf)) + str_velo(k) / dble(nstr_per_elem)
enddo
next_elem = .true.
endif
enddo
enddo
enddo
do i = 1,osolve%nnode
ov%wnodecompact(i) = ov%wnodecompact(i) / dble(countnode(i))
ov%wnode(i)=ov%wnode(i)+ov%wnodecompact(i)
enddo
!Note: if storing velocity at a differently sized array, could compare
! density, densityp outside the above loop over all instances of density_str
end subroutine compaction
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------