!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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,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
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! 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.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))
allocate (density_temp%density(size_str),stat=err)
density_temp%density=0.d0
! 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
! Calculate current vertical position
! 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
! JANICE2: Commented out debug output below
!if (jstr.eq.2**(params%levelmax_oct-1)) then
! if (istr.eq.3 .or. istr.eq.2**(params%levelmax_oct-1) .or. istr.eq.2**(params%levelmax_oct)-3) then
! !if (iproc.eq.0) then
! write (*,*) 'density string output (',istr,',',jstr,')'
! write (*,*) 'x, y, z :',x,', ',y,', ',z
! !endif
! endif
!endif
! Use LSFs to determine material type
cur_lsf=1.d0
do k=1,osolve%nlsf
! unsquished z into octree_interpolate
call octree_interpolate (osolve%octree,osolve%noctree,osolve%icon, &
osolve%nleaves,osolve%lsf(:,k),osolve%nnode, &
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
! JANICE2: Commented out debug output below
!if (jstr.eq.2**(params%levelmax_oct-1)) then
! if (istr.eq.3 .or. istr.eq.2**(params%levelmax_oct-1) .or. istr.eq.2**(params%levelmax_oct)-3) then
!
! !if (iproc.eq.0) then
! write (*,*) 'k = ', k, '; matnum = ', matnum
! !endif
!
! endif
!endif
enddo
! Update material numbers for material changes
call mat_trans (params,mat,matnum,x,y,z)
!endif
! JANICE2: Commented out debug output below
!if (jstr.eq.2**(params%levelmax_oct-1)) then
! if (istr.eq.3 .or. istr.eq.2**(params%levelmax_oct-1) .or. istr.eq.2**(params%levelmax_oct)-3) then
!
!
! !if (iproc.eq.0) then
! write (*,*) 'at location ',i,' z = ',z,'mat_num =',matnum
! !endif
!
! endif
!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
dfluid=mat(matnum)%pore_fluid_density
spor=mat(matnum)%surface_porosity
ccoef=mat(matnum)%compaction_coef
density_temp%density(i)=-(dgrain-(dgrain-dfluid)*spor*exp(-ccoef*(z_free_surf-z)))
! JANICE2: Commented out debug output below
! if (jstr.eq.2**(params%levelmax_oct-1)) then
! if (istr.eq.3 .or. istr.eq.2**(params%levelmax_oct-1) .or. istr.eq.2**(params%levelmax_oct)-3) then
!
!
! !if (iproc.eq.0) then
! write (*,*) 'Compacting density'
! write (*,*) 'density_temp%density(',i,'): ',density_temp%density(i)
! !endif
!
!
! endif
! endif
else
density_temp%density(i)=mat(matnum)%density
! JANICE2: Commented out debug output below
! if (jstr.eq.2**(params%levelmax_oct-1)) then
! if (istr.eq.3 .or. istr.eq.2**(params%levelmax_oct-1) .or. istr.eq.2**(params%levelmax_oct)-3) then
!
!
! !if (iproc.eq.0) then
! write (*,*) 'Non-compacting density'
! write (*,*) 'density_temp%density(',i,'): ',density_temp%density(i)
! !endif
!
!
! endif
! endif
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
! JANICE2: Commented out debug output below
! if (jstr.eq.2**(params%levelmax_oct-1)) then
! if (istr.eq.3 .or. istr.eq.2**(params%levelmax_oct-1) .or. istr.eq.2**(params%levelmax_oct)-3) then
!
!
! !if (iproc.eq.0) then
! write (*,*) 'Initializing density'
! !endif
!
!
! endif
! endif
density_profile%density=density_temp%density
! JANICE2: Commented out debug output below
! if (jstr.eq.2**(params%levelmax_oct-1)) then
! if (istr.eq.3 .or. istr.eq.2**(params%levelmax_oct-1) .or. istr.eq.2**(params%levelmax_oct)-3) then
!
!
! !if (iproc.eq.0) then
! write (*,*) 'density_profile%density: ',density_profile%density
! !endif
!
!
! endif
! endif
! 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
! JANICE2: Commented out debug output below
! if (jstr.eq.2**(params%levelmax_oct-1)) then
! if (istr.eq.3 .or. istr.eq.2**(params%levelmax_oct-1) .or. istr.eq.2**(params%levelmax_oct)-3) then
!
!
! !if (iproc.eq.0) then
! write (*,*) 'Updating density'
! !endif
!
!
! endif
! endif
do i=1,size_str
if (abs(density_temp%density(i)).gt.abs(density_profile%density(i))) then
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 (abs(density_temp%density(i)) < eps) then
density_profile%density(i)=density_temp%density(i)
endif
! JANICE2: Commented out debug output below
! if (jstr.eq.2**(params%levelmax_oct-1)) then
! if (istr.eq.3 .or. istr.eq.2**(params%levelmax_oct-1) .or. istr.eq.2**(params%levelmax_oct)-3) then
!
!
! !if (iproc.eq.0) then
! write (*,*) 'density_profile%density(',i,'): ',density_profile%density(i)
! !endif
!
!
! endif
! endif
enddo
endif
deallocate (density_temp%density)
end subroutine calculate_density_profile
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------