calculate_density_profile.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!                     CALCULATE DENSITY PROFILE    MARCH 2013                  |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine calculate_density_profile(params,density_profile,osolve,mat,        &
                                     initialize,istr,jstr,calc_comp_velo,      &
                                     write_compp)

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, write_compp
logical :: calc_comp_velo(size(density_profile%density))


!------------------------------------------------------------------------------|
!(((((((((((((((( 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,eps
double precision :: dgrain,dfluid,spor,ccoef
integer :: matnum,i,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



 ! 100414 JCA 
    ! If this position of the string is current in a compactible material, and also
    ! this position was compactible at the last time step, set calc_comp_velo=T
      
    if (mat(matnum)%compactible.and.density_profile%compactiblep(i)) then
      !write (*,*) 'size(calc_comp_velo): ',size(calc_comp_velo)
      !write (*,*) 'i: ',i
      !write (*,*) 'size_str: ',size_str
      calc_comp_velo(i)=.true.
    endif
         
      ! 100414 JCA
      ! If this calculation of compaction is at the end of the timestep, overwrite
      ! compactiblep for this location with the current status (this will now store
      ! compactibility, for use in the next time step's calculations
    
    if (write_compp.gt.eps.or.initialize.gt.eps) then
      density_profile%compactiblep(i)=mat(matnum)%compactible
    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

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------