crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

added projection on normal scores, minor fixes, code clean up

Browse source
This commit is contained in:
Guillaume 2017-06-16 23:03:27 -04:00
parent df12cb6c1d
commit 21c3d6083c
9 changed files with 120 additions and 252 deletions
Download patch file Download diff file Expand all files Collapse all files

5
src/constants.F90
View file

@ -305,7 +305,7 @@ module constants
EVENT_ABSORB = 2
! Tally score type
integer, parameter :: N_SCORE_TYPES = 25
integer, parameter :: N_SCORE_TYPES = 26
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
@ -331,7 +331,8 @@ module constants
SCORE_FISS_Q_PROMPT = -22, & ! prompt fission Q-value
SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value
SCORE_DECAY_RATE = -24, & ! delayed neutron precursor decay rate
SCORE_CELL_TO_CELL_TYPE = -25 ! cell to cell partial current
SCORE_CELL_TO_CELL = -25, & ! cell to cell partial current
SCORE_CELL_TO_CELL_NORMAL_PROJECTION = -26 ! cell to cell partial current projected on surface normal
! Maximum scattering order supported
integer, parameter :: MAX_ANG_ORDER = 10

4
src/endf.F90
View file

@ -66,6 +66,10 @@ contains
string = "fission-q-prompt"
case (SCORE_FISS_Q_RECOV)
string = "fission-q-recoverable"
case (SCORE_CELL_TO_CELL)
string = "cell-to-cell-partial-current"
case (SCORE_CELL_TO_CELL_NORMAL_PROJECTION)
string = "cell-to-cell-partial-current-projection-on-normal"
! Normal ENDF-based reactions
case (TOTAL_XS)

41
src/geometry.F90
View file

@ -12,7 +12,7 @@ module geometry
use stl_vector, only: VectorInt
use string, only: to_str
use tally, only: score_surface_current, &
&score_cell_to_cell
&score_partial_current
implicit none
@ -399,8 +399,6 @@ contains
logical :: rotational ! if rotational periodic BC applied
logical :: found ! particle found in universe?
class(Surface), pointer :: surf
!print *, "cross surface called"
i_surface = abs(p % surface)
surf => surfaces(i_surface)%obj
@ -613,23 +611,33 @@ contains
! cells on the positive side
call find_cell(p, found, surf%neighbor_pos)
! /CHANGE/ Score cell_to_cell
! print *, "Positive side, calling score_cell_to_cell"
call score_cell_to_cell(p, last_cell)
! Save cosine of angle between surface normal and particle direction
p % normal_proj = dot_product(p % coord(1) % uvw, &
& surf % normal(p % coord(1) % xyz)) / &
& sqrt(dot_product(surf % normal(p % coord(1) % xyz)&
&, surf % normal(p % coord(1) % xyz)))
! Score cell to cell partial currents
call score_partial_current(p)
if (found) return
elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then
! If coming from positive side of surface, search all the neighboring
! cells on the negative side
! Save cosine of angle between surface normal and particle direction
p % normal_proj = dot_product(p % coord(1) % uvw, &
& surf % normal(p % coord(1) % xyz)) / &
& sqrt(dot_product(surf % normal(p % coord(1) % xyz)&
&, surf % normal(p % coord(1) % xyz)))
call find_cell(p, found, surf%neighbor_neg)
! /CHANGE/ Score cell_to_cell
! print *, "Negative side, calling score_cell_to_cell"
call score_cell_to_cell(p, last_cell)
! Score cell to cell partial currents
call score_partial_current(p)
if (found) return
end if
@ -662,10 +670,9 @@ contains
return
end if
end if
! /CHANGE/ Score cell_to_cell
! print *, "Searched all cells then calling score_cell_to_cell"
call score_cell_to_cell(p, last_cell)
! Score cell to cell partial currents
call score_partial_current(p)
end subroutine cross_surface

26
src/input_xml.F90
View file

@ -3004,7 +3004,7 @@ contains
end if
n_words = node_word_count(node_filt, "bins")
case ("mesh", "universe", "material", "cell", "distribcell", &
"cellborn", "surface", "delayedgroup")
"cellborn", "cellto", "cellfrom", "surface", "delayedgroup")
if (.not. check_for_node(node_filt, "bins")) then
call fatal_error("Bins not set in filter " // trim(to_str(filter_id)))
end if
@ -3035,7 +3035,7 @@ contains
allocate(filt % cells(n_words))
call get_node_array(node_filt, "bins", filt % cells)
end select
case ('cellfrom')
! Allocate and declare the filter type
allocate(CellFromFilter :: f % obj)
@ -3046,7 +3046,7 @@ contains
allocate(filt % cells(n_words))
call get_node_array(node_filt, "bins", filt % cells)
end select
case ('cellto')
! Allocate and declare the filter type
allocate(CellToFilter :: f % obj)
@ -3057,7 +3057,7 @@ contains
allocate(filt % cells(n_words))
call get_node_array(node_filt, "bins", filt % cells)
end select
case ('cellborn')
! Allocate and declare the filter type
allocate(CellbornFilter :: f % obj)
@ -3414,17 +3414,10 @@ contains
t % find_filter(FILTER_DISTRIBCELL) = j
type is (CellFilter)
t % find_filter(FILTER_CELL) = j
type is (CellFromFilter)
t % find_filter(FILTER_CELL_TO_CELL) = j
t % estimator = ESTIMATOR_ANALOG
type is (CellToFilter)
t % find_filter(FILTER_CELL_TO_CELL) = j
t % estimator = ESTIMATOR_ANALOG
type is (CellbornFilter)
t % find_filter(FILTER_CELLBORN) = j
type is (MaterialFilter)
@ -3927,13 +3920,14 @@ contains
t % find_filter(FILTER_SURFACE) = n_filter
t % filter(n_filter) = i_filt
case ('events')
t % score_bins(j) = SCORE_EVENTS
case ('partial_current')
t % type = TALLY_CELL_TO_CELL
t % score_bins(j) = SCORE_CELL_TO_CELL_TYPE
t % score_bins(j) = SCORE_CELL_TO_CELL
case ('projected_partial_current')
t % type = TALLY_CELL_TO_CELL
t % score_bins(j) = SCORE_CELL_TO_CELL_NORMAL_PROJECTION
case ('events')
t % score_bins(j) = SCORE_EVENTS
case ('elastic', '(n,elastic)')
t % score_bins(j) = ELASTIC
case ('(n,2nd)')

3
src/output.F90
View file

@ -733,6 +733,9 @@ contains
score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity"
score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power"
score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power"
score_names(abs(SCORE_CELL_TO_CELL)) = "Partial current"
score_names(abs(SCORE_CELL_TO_CELL_NORMAL_PROJECTION)) = "Partial &
&current projected on normal of surface"
! Create filename for tally output
filename = trim(path_output) // "tallies.out"

20
src/particle_header.F90
View file

@ -15,9 +15,7 @@ module particle_header
!===============================================================================
type LocalCoord
integer :: last_cell = NONE
! Indices in various arrays for this level
integer :: cell = NONE
integer :: universe = NONE
@ -60,13 +58,11 @@ module particle_header
real(8) :: wgt ! particle weight
real(8) :: mu ! angle of scatter
logical :: alive ! is particle alive?
! Previous cell before crossing surface
integer :: last_cell
! Particle coordinates
integer :: last_n_coord ! number of current coordinates
! type(LocalCoord) :: last_coord(MAX_COORD) ! coordinates for all levels
! Crossing surface data
integer :: last_cell ! last cell the particle was in before crossing
real(8) :: normal_proj ! cos of angle to surface normal when crossing
! Pre-collision physical data
real(8) :: last_xyz_current(3) ! coordinates of the last collision or
! reflective/periodic surface crossing
@ -136,6 +132,7 @@ contains
! clear attributes
this % surface = NONE
this % last_cell = NONE
this % cell_born = NONE
this % material = NONE
this % last_material = NONE
@ -143,6 +140,7 @@ contains
this % wgt = ONE
this % last_wgt = ONE
this % absorb_wgt = ZERO
this % normal_proj = 1
this % n_bank = 0
this % wgt_bank = ZERO
this % sqrtkT = ERROR_REAL
@ -181,8 +179,6 @@ contains
elemental subroutine reset_coord(this)
class(LocalCoord), intent(inout) :: this
this % last_cell = NONE
this % cell = NONE
this % universe = NONE
this % lattice = NONE

79
src/tally.F90
View file

@ -116,12 +116,8 @@ contains
!#########################################################################
! Determine appropirate scoring value.
! print *, "Score bin", score_bin
select case(score_bin)
! use default case for cell_to_from tallies, with t % estimator == ESTIMATOR_ANALOG
case (SCORE_FLUX, SCORE_FLUX_YN)
if (t % estimator == ESTIMATOR_ANALOG) then
@ -1180,16 +1176,13 @@ contains
end if
end if
case(SCORE_CELL_TO_CELL_TYPE)
! Not using default because score_bin and event_MT thing
! score_bin is -25, and event_MT seems to match the cell id
! print*, p % event_MT, score_bin
! print *, "Scoring", p % last_wgt, flux
case(SCORE_CELL_TO_CELL)
score = p % last_wgt * flux
case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION)
score = p % last_wgt * flux * p % normal_proj
case default
!print *, "Found right score type"
!print * , t % estimator
if (t % estimator == ESTIMATOR_ANALOG) then
! Any other score is assumed to be a MT number. Thus, we just need
! to check if it matches the MT number of the event
@ -2091,14 +2084,13 @@ contains
case (SCORE_EVENTS)
! Simply count number of scoring events
score = ONE
case(SCORE_CELL_TO_CELL_TYPE)
! Not using default because score_bin and event_MT thing
! score_bin is -25, and event_MT seems to match the cell id
! print*, p % event_MT, score_bin
! print *, "Scoring", p % last_wgt, flux
case(SCORE_CELL_TO_CELL)
score = p % last_wgt * flux
case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION)
score = p % last_wgt * flux * p % normal_proj
end select
!#########################################################################
@ -2130,8 +2122,6 @@ contains
integer :: num_nm ! Number of N,M orders in harmonic
integer :: n ! Moment loop index
real(8) :: uvw(3)
! print *, "Expand and score, score bin", score_bin, "score", score, "filter_index", filter_index
select case(score_bin)
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
@ -2215,7 +2205,6 @@ contains
case default
! print *, "Adding to total in expand-and-score", score
!$omp atomic
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score
@ -3160,18 +3149,17 @@ contains
end subroutine score_collision_tally
!===============================================================================
! SCORE_CELL_TO_CELL tallies partial currents from cell_to to cell_from
! score_partial_current tallies partial currents from cell_to to cell_from
!===============================================================================
subroutine score_cell_to_cell(p, last_cell)
subroutine score_partial_current(p)
type(Particle), intent(in) :: p
integer, intent(in) :: last_cell
! if particle has leaked
! DF could preserve leakage aiight
!if p % alive = .false. then
! DF could preserve leakage
! if p % alive = .false. then
!end if
integer :: i
integer :: i_tally
@ -3183,9 +3171,8 @@ contains
real(8) :: flux ! collision estimate of flux
real(8) :: filter_weight ! combined weight of all filters
logical :: finished ! found all valid bin combinations
real(8) :: dummy_real8
type(TallyObject), pointer :: t
! Determine collision estimate of flux
if (survival_biasing) then
! We need to account for the fact that some weight was already absorbed
@ -3193,13 +3180,15 @@ contains
else
flux = p % last_wgt
end if
! A loop over all tallies is necessary because we need to simultaneously
! determine different filter bins for the same tally in order to score to it
! However, we can restrict the list to cell_to_cell tallies since we know we
! are crossing a surface
TALLY_LOOP: do i = 1, active_collision_tallies % size()
TALLY_LOOP: do i = 1, active_cell_to_cell_tallies % size()
! Get index of tally and pointer to tally
i_tally = active_collision_tallies % data(i)
i_tally = active_cell_to_cell_tallies % data(i)
t => tallies(i_tally)
! Find all valid bins in each filter if they have not already been found
@ -3245,11 +3234,10 @@ contains
filter_weight = filter_weight * filter_matches(i_filt) % weights % &
data(i_bin)
end do
call score_general(p, t, 0, filter_index, &
0, dummy_real8, filter_weight)
0, 0d0, flux * filter_weight)
! ======================================================================
! Filter logic
@ -3282,15 +3270,15 @@ contains
if (assume_separate) exit TALLY_LOOP
end do TALLY_LOOP
! Reset filter matches flag
filter_matches(:) % bins_present = .false.
! Reset tally map positioning
position = 0
end subroutine score_cell_to_cell
end subroutine score_partial_current
!===============================================================================
! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually
! determining which mesh surfaces were crossed
@ -4584,9 +4572,8 @@ contains
call active_current_tallies % push_back(i_user_tallies + i)
elseif (user_tallies(i) % type == TALLY_CELL_TO_CELL) then
call active_cell_to_cell_tallies % push_back(i_user_tallies + i)
call active_current_tallies % push_back(i_user_tallies + i)
end if
end do
call active_tallies % shrink_to_fit()
@ -4594,6 +4581,7 @@ contains
call active_tracklength_tallies % shrink_to_fit()
call active_collision_tallies % shrink_to_fit()
call active_current_tallies % shrink_to_fit()
call active_cell_to_cell_tallies % shrink_to_fit()
end subroutine setup_active_usertallies
@ -4617,6 +4605,9 @@ contains
else if (active_current_tallies % size() > 0) then
call fatal_error("Active current tallies should not exist before CMFD &
&tallies!")
else if (active_cell_to_cell_tallies % size() > 0) then
call fatal_error("Active cell to cell tallies should not exist before &
&CMFD tallies!")
end if
do i = 1, n_cmfd_tallies
@ -4627,6 +4618,7 @@ contains
if (cmfd_tallies(i) % type == TALLY_VOLUME) then
if (cmfd_tallies(i) % estimator == ESTIMATOR_ANALOG) then
call active_analog_tallies % push_back(i_cmfd_tallies + i)
call active_cell_to_cell_tallies % push_back(i_cmfd_tallies + 1)
elseif (cmfd_tallies(i) % estimator == ESTIMATOR_TRACKLENGTH) then
call active_tracklength_tallies % push_back(i_cmfd_tallies + i)
end if
@ -4640,6 +4632,7 @@ contains
call active_tracklength_tallies % shrink_to_fit()
call active_collision_tallies % shrink_to_fit()
call active_current_tallies % shrink_to_fit()
call active_cell_to_cell_tallies % shrink_to_fit()
end subroutine setup_active_cmfdtallies

183
src/tally_filter.F90
View file

@ -70,20 +70,7 @@ module tally_filter
procedure :: text_label => text_label_cell
procedure :: initialize => initialize_cell
end type CellFilter
!===============================================================================
! CELLFROMTOFILTER specifies which geometric cells particles exit.
!===============================================================================
type, extends(TallyFilter) :: CellToCellFilter
integer, allocatable :: cells(:)
type(DictIntInt) :: map
contains
procedure :: get_next_bin => get_next_bin_cell_to_cell
procedure :: to_statepoint => to_statepoint_cell_to_cell
procedure :: text_label => text_label_cell_to_cell
procedure :: initialize => initialize_cell_to_cell
end type CellToCellFilter
!===============================================================================
! CELLFROMFILTER specifies which geometric cells particles exit.
!===============================================================================
@ -96,7 +83,7 @@ module tally_filter
procedure :: text_label => text_label_cell_from
procedure :: initialize => initialize_cell_from
end type CellFromFilter
!===============================================================================
! CELLTOFILTER specifies which geometric cells particles enter.
!===============================================================================
@ -757,99 +744,10 @@ contains
end function text_label_cell
!===============================================================================
! Cell_to_cell_Filter methods
! CellFromFilter methods
!===============================================================================
subroutine get_next_bin_cell_to_cell(this, p, estimator, current_bin, &
subroutine get_next_bin_cell_from(this, p, estimator, current_bin, &
next_bin, weight)
class(CellToCellFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
integer, value, intent(in) :: current_bin
integer, intent(out) :: next_bin
real(8), intent(out) :: weight
integer :: i, start
! Find the coordinate level of the last bin we found.
if (current_bin == NO_BIN_FOUND) then
start = 1
else
do i = 1, p % n_coord
if (p % coord(i) % cell == this % cells(current_bin)) then
start = i + 1
exit
end if
end do
end if
! Starting one coordinate level deeper, find the next bin. !!!!!!!!!!!!!!!!!! OLD !!!!!!!!!!!!!!!!!!!!!!!!!!
next_bin = NO_BIN_FOUND
weight = ERROR_REAL
do i = start, p % n_coord
if (this % map % has_key(p % coord(i) % cell) .and. &
this % map % has_key(p % last_cell)) then
next_bin = this % map % get_key(p % coord(i) % cell)
weight = ONE
exit
end if
end do
end subroutine get_next_bin_cell_to_cell
subroutine to_statepoint_cell_to_cell(this, filter_group)
class(CellToCellFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
integer :: i
integer, allocatable :: cell_ids(:)
call write_dataset(filter_group, "type", "cell")
call write_dataset(filter_group, "n_bins", this % n_bins)
allocate(cell_ids(size(this % cells)))
do i = 1, size(this % cells)
cell_ids(i) = cells(this % cells(i)) % id
end do
call write_dataset(filter_group, "bins", cell_ids)
end subroutine to_statepoint_cell_to_cell
subroutine initialize_cell_to_cell(this)
class(CellToCellFilter), intent(inout) :: this
integer :: i, id
! Convert ids to indices.
do i = 1, this % n_bins
id = this % cells(i)
! print *,this % n_bins, this % cells(i), cell_dict % get_key(id)
if (cell_dict % has_key(id)) then
this % cells(i) = cell_dict % get_key(id)
else
call fatal_error("Could not find cell " // trim(to_str(id)) &
&// " specified on tally filter.")
end if
end do
! Generate mapping from cell indices to filter bins.
do i = 1, this % n_bins
call this % map % add_key(this % cells(i), i)
end do
end subroutine initialize_cell_to_cell
function text_label_cell_to_cell(this, bin) result(label)
class(CellToCellFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
label = "Cell to from" // to_str(cells(this % cells(bin)) % id)
end function text_label_cell_to_cell
!===============================================================================
! Cell_From_Filter methods
!===============================================================================
subroutine get_next_bin_cell_from(this, p, estimator, current_bin, &
next_bin, weight)
class(CellFromFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
@ -859,45 +757,24 @@ contains
integer :: i, start
! Find the coordinate level of the last bin we found.
! if (current_bin == NO_BIN_FOUND) then
! start = 1
! else
! do i = 1, p % last_n_coord !shouldnt be done this way, modify later!!
! ! print *, p % last_cell, this % cells(current_bin)
! if (p % last_cell == this % cells(current_bin)) then
! start = i + 1
! exit
! end if
! end do
! end if
! Starting one coordinate level deeper, find the next bin.
! Particle can only have one cell_from
! If current_bin is already a bin, then no need to look for another bin
next_bin = NO_BIN_FOUND
weight = ERROR_REAL
weight = ONE
if (current_bin == NO_BIN_FOUND) then
! print *, "Looking in map for cell at", p % last_cell," :",this % map % get_key(p % last_cell)
if (this % map % has_key(p % last_cell)) then
! print *, "Found in map for cell at", p % last_cell," :",this % map % get_key(p % last_cell)
next_bin = this % map % get_key(p % last_cell)
weight = ONE
end if
else
next_bin = NO_BIN_FOUND
end if
! do i = start, p % last_n_coord
! if (this % map % has_key(p % last_cell)) then
! next_bin = this % map % get_key(p % coord(i) % last_cell) !need to use coord(i) to dig one level deeper
! weight = ONE
! exit
! end if
! end do
end subroutine get_next_bin_cell_from
subroutine to_statepoint_cell_from(this, filter_group)
class(CellFromFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
integer(HID_T), intent(in) :: filter_group
integer :: i
integer, allocatable :: cell_ids(:)
@ -920,7 +797,6 @@ contains
! Convert ids to indices.
do i = 1, this % n_bins
id = this % cells(i)
! print *,this % n_bins, this % cells(i), cell_dict % get_key(id)
if (cell_dict % has_key(id)) then
this % cells(i) = cell_dict % get_key(id)
else
@ -931,31 +807,29 @@ contains
! Generate mapping from cell indices to filter bins.
do i = 1, this % n_bins
! print *, "Filling bin/cell map : cell",this % cells(i), "bin", i
call this % map % add_key(this % cells(i), i)
end do
end subroutine initialize_cell_from
function text_label_cell_from(this, bin) result(label)
class(CellFromFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
label = "Cell to from" // to_str(cells(this % cells(bin)) % id)
label = "Cell from " // to_str(cells(this % cells(bin)) % id)
end function text_label_cell_from
!===============================================================================
! Cell_To_Filter methods
! CellToFilter methods
!===============================================================================
subroutine get_next_bin_cell_to(this, p, estimator, current_bin, &
subroutine get_next_bin_cell_to(this, p, estimator, current_bin, &
next_bin, weight)
class(CellToFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
integer, value, intent(in) :: current_bin
integer, intent(out) :: next_bin
real(8), intent(out) :: weight
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
integer, value, intent(in) :: current_bin
integer, intent(out) :: next_bin
real(8), intent(out) :: weight
integer :: i, start
@ -985,7 +859,7 @@ contains
subroutine to_statepoint_cell_to(this, filter_group)
class(CellToFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
integer(HID_T), intent(in) :: filter_group
integer :: i
integer, allocatable :: cell_ids(:)
@ -1008,7 +882,7 @@ contains
! Convert ids to indices.
do i = 1, this % n_bins
id = this % cells(i)
! print *,this % n_bins, this % cells(i), cell_dict % get_key(id)
if (cell_dict % has_key(id)) then
this % cells(i) = cell_dict % get_key(id)
else
@ -1022,16 +896,15 @@ contains
call this % map % add_key(this % cells(i), i)
end do
end subroutine initialize_cell_to
function text_label_cell_to(this, bin) result(label)
class(CellToFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
label = "Cell to to" // to_str(cells(this % cells(bin)) % id)
label = "Cell to " // to_str(cells(this % cells(bin)) % id)
end function text_label_cell_to
!===============================================================================
! DistribcellFilter methods
!===============================================================================

11
src/tracking.F90
View file

@ -149,11 +149,10 @@ contains
! PARTICLE CROSSES SURFACE
if (next_level > 0) p % n_coord = next_level
! Saving previous cell data
last_cell = p % coord(p % n_coord) % cell
last_n_coord = p % n_coord
p % coord(p % n_coord) % cell = NONE
if (any(lattice_translation /= 0)) then
! Particle crosses lattice boundary
@ -163,12 +162,10 @@ contains
else
! Particle crosses surface
p % surface = surface_crossed
! /CHANGE/ Saving stuff on particle
p % last_cell = last_cell
p % last_n_coord = last_n_coord
p % coord(next_level) % last_cell = last_cell
call cross_surface(p, last_cell)
p % event = EVENT_SURFACE
end if