From e6af8d35442b5d690d0fd306bf57e6a17691e6bf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 14 Feb 2011 05:19:45 +0000 Subject: [PATCH] Sampling of reactions, several important bugs fixed. --- ChangeLog | 19 ++++++ src/Makefile | 4 +- src/ace.f90 | 58 +++++++++++++---- src/endf.f90 | 147 +++++++++++++++++++++++++++++++++++++++++ src/geometry.f90 | 74 +++++++++++---------- src/global.f90 | 16 +++-- src/input_sample | 13 ++-- src/main.f90 | 8 +-- src/physics.f90 | 165 ++++++++++++++++++++++++++++++++++++++--------- src/types.f90 | 21 +++--- 10 files changed, 422 insertions(+), 103 deletions(-) create mode 100644 src/endf.f90 diff --git a/ChangeLog b/ChangeLog index 5cbfb09c1..e6a35fe30 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,22 @@ +2011-02-13 Paul Romano + + * ace.f90: Fixed xs_continuous allocation again -- cross sections + were not being assigned correctly to materials. Added elastic + scattering in list of reactions. Changed name of rxn % + energy_index to IE. + * endf.f90: New module. Includes subroutine reaction_name which + returns the name of a reaction for a given ENDF MT value. + * geometry.f90: Fixed bug in dist_to_boundary (the 'w' in uvw was + actually being set to 'z', not 'w'). Changed 0's and 1's to ZEROs + and ONEs. + * global.f90: Added ZERO and ONE. + * physics.f90: Can now sample reactions once nuclide is + selected. Added subroutines elastic_scatter, level_inelastic, and + n_gamma (level_inelastic is just a stub), and significantly + changed elastic scattering (currently based on stationary target). + * types.f90: Added n_reaction attribute on AceContinuous, changed + AceReaction % energy_index to IE. + 2011-02-09 Paul Romano * ace.f90: Fixed bug in the way xs_continuous was allocated. If, diff --git a/src/Makefile b/src/Makefile index 341826cf1..ce1816838 100644 --- a/src/Makefile +++ b/src/Makefile @@ -5,6 +5,7 @@ test = unittest.f90 modules = ace.f90 \ cross_section.f90 \ data_structures.f90 \ + endf.f90 \ energy_grid.f90 \ fileio.f90 \ geometry.f90 \ @@ -50,6 +51,7 @@ unittest: $(test_objects) ace.o: global.o output.o string.o fileio.o string.o cross_section.o: global.o string.o data_structures.o output.o data_structures.o: global.o +endf.o: global.o energy_grid.o: global.o output.o data_structures.o fileio.o: types.o global.o string.o output.o data_structures.o geometry.o: types.o global.o output.o string.o data_structures.o @@ -59,7 +61,7 @@ main.o: global.o fileio.o output.o geometry.o mcnp_random.o \ ace.o energy_grid.o output.o: global.o physics.o: types.o global.o mcnp_random.o geometry.o output.o \ - search.o + search.o endf.o search.o: output.o source.o: global.o mcnp_random.o string.o: global.o output.o diff --git a/src/ace.f90 b/src/ace.f90 index 7a5bce9e3..911f495c5 100644 --- a/src/ace.f90 +++ b/src/ace.f90 @@ -6,7 +6,7 @@ module ace use fileio, only: read_line, read_data, skip_lines use string, only: split_string, str_to_real use data_structures, only: dict_create, dict_add_key, dict_has_key, & - & dict_get_key + & dict_get_key, dict_delete integer :: NXS(16) integer :: JXS(32) @@ -37,7 +37,7 @@ contains integer :: n integer :: index_continuous integer :: index_thermal - type(ListKeyValueCI), pointer :: elem + type(DictionaryCI), pointer :: temp_dict call dict_create(ace_dict) @@ -78,14 +78,31 @@ contains allocate(xs_thermal(n_thermal)) ! loop over all nuclides in xsdata + call dict_create(temp_dict) + index_continuous = 0 - do i = 1, size(xsdatas) - key = xsdatas(i)%alias - if (dict_has_key(ace_dict, key)) then - index_continuous = index_continuous + 1 - call read_ACE_continuous(index_continuous, i) - end if + index_thermal = 0 + do i = 1, n_materials + mat => materials(i) + do j = 1, mat%n_isotopes + index = mat%isotopes(j) + key = xsdatas(index)%id + n = len_trim(key) + call lower_case(key) + select case (key(n:n)) + case ('c') + if (.not. dict_has_key(temp_dict, key)) then + index_continuous = index_continuous + 1 + call read_ACE_continuous(index_continuous, index) + end if + case ('t') + n_thermal = n_thermal + 1 + end select + end do end do + + ! delete dictionary + call dict_delete(temp_dict) end subroutine read_xs @@ -222,7 +239,11 @@ contains allocate(table%sigma_el(NE)) allocate(table%heating(NE)) - ! read data from XSS + ! read data from XSS -- right now the total, absorption and + ! elastic scattering are read in to these special arrays, but in + ! reality, it should be necessary to only store elastic scattering + ! and possible total cross-section for total material xs + ! generation. XSS_index = 1 table%energy = get_real(NE) table%sigma_t = get_real(NE) @@ -259,11 +280,22 @@ contains JXS7 = JXS(7) NMT = NXS(4) - ! allocate array of reactions - allocate(table%reactions(NMT)) + ! allocate array of reactions. Add one since we need to include an + ! elastic scattering channel + table%n_reaction = NMT + 1 + allocate(table%reactions(NMT+1)) + + ! Store elastic scattering cross-section on reaction one + rxn => table%reactions(1) + rxn%MT = 2 + rxn%Q_value = 0.0_8 + rxn%TY = 1 + rxn%IE = 1 + allocate(rxn%sigma(table%n_grid)) + rxn%sigma = table%sigma_el do i = 1, NMT - rxn => table%reactions(i) + rxn => table%reactions(i+1) ! read MT number, Q-value, and neutrons produced rxn%MT = XSS(LMT+i-1) @@ -272,7 +304,7 @@ contains ! read cross section values LOCA = XSS(LXS+i-1) - rxn%energy_index = XSS(JXS7 + LOCA - 1) + rxn%IE = XSS(JXS7 + LOCA - 1) NE = XSS(JXS7 + LOCA) allocate(rxn%sigma(NE)) XSS_index = JXS7 + LOCA + 1 diff --git a/src/endf.f90 b/src/endf.f90 new file mode 100644 index 000000000..a76fea201 --- /dev/null +++ b/src/endf.f90 @@ -0,0 +1,147 @@ +module endf + + use global, only: int_to_str + +contains + +!===================================================================== +! REACTION_NAME gives the name of the reaction for a given MT value +!===================================================================== + + function reaction_name(MT) result(string) + + integer, intent(in) :: MT + character(20) :: string + + select case (MT) + case (1) + string = '(n,total)' + case (2) + string = '(elastic)' + case (11) + string = '(n,2nd)' + case (16) + string = '(n,2n)' + case (17) + string = '(n,3n)' + case (18) + string = '(n,fission)' + case (19) + string = '(n,f)' + case (20) + string = '(n,2nf)' + case (21) + string = '(n,2nf)' + case (22) + string = '(n,na)' + case (23) + string = '(n,n3a)' + case (24) + string = '(n,2na)' + case (25) + string = '(n,3na)' + case (28) + string = '(n,np)' + case (29) + string = '(n,n2a)' + case (30) + string = '(n,2n2a)' + case (32) + string = '(n,nd)' + case (33) + string = '(n,nt)' + case (34) + string = '(n,n He-3)' + case (35) + string = '(n,nd3a)' + case (36) + string = '(n,nt2a)' + case (37) + string = '(n,4n)' + case (38) + string = '(n,3nf)' + case (41) + string = '(n,2np)' + case (42) + string = '(n,3np)' + case (44) + string = '(n,2np)' + case (45) + string = '(n,npa)' + case (51 : 90) + string = '(n,n' // trim(int_to_str(MT-50)) // ')' + case (91) + string = '(n,nc)' + case (102) + string = '(n,gamma)' + case (103) + string = '(n,p)' + case (104) + string = '(n,d)' + case (105) + string = '(n,t)' + case (106) + string = '(n,3He)' + case (107) + string = '(n,a)' + case (108) + string = '(n,2a)' + case (109) + string = '(n,3a)' + case (111) + string = '(n,2p)' + case (112) + string = '(n,pa)' + case (113) + string = '(n,t2a)' + case (114) + string = '(n,d2a)' + case (115) + string = '(n,pd)' + case (116) + string = '(n,pt)' + case (117) + string = '(n,da)' + case (201) + string = '(n,Xn)' + case (202) + string = '(n,Xgamma)' + case (203) + string = '(n,Xp)' + case (204) + string = '(n,Xd)' + case (205) + string = '(n,Xt)' + case (206) + string = '(n,X3He)' + case (207) + string = '(n,Xa)' + case (444) + string = '(damage)' + case (600 : 648) + string = '(n,p' // trim(int_to_str(MT-600)) // ')' + case (649) + string = '(n,pc)' + case (650 : 698) + string = '(n,d' // trim(int_to_str(MT-650)) // ')' + case (699) + string = '(n,dc)' + case (700 : 748) + string = '(n,t' // trim(int_to_str(MT-600)) // ')' + case (749) + string = '(n,tc)' + case (750 : 798) + string = '(n,3He' // trim(int_to_str(MT-650)) // ')' + case (799) + string = '(n,3Hec)' + case (800 : 848) + string = '(n,a' // trim(int_to_str(MT-800)) // ')' + case (849) + string = '(n,tc)' + case default + string = 'MT=' // trim(int_to_str(MT)) + end select + + end function reaction_name + +end module endf diff --git a/src/geometry.f90 b/src/geometry.f90 index 4aea337be..7fcc49e44 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -87,7 +87,7 @@ contains subroutine find_cell(neut) - type(Neutron), pointer, intent(inout) :: neut + type(Neutron), pointer :: neut type(Cell), pointer :: this_cell logical :: found_cell @@ -123,7 +123,7 @@ contains subroutine cross_boundary(neut) - type(Neutron), pointer, intent(in) :: neut + type(Neutron), pointer :: neut type(Surface), pointer :: surf type(Cell), pointer :: c @@ -155,6 +155,8 @@ contains c => cells(index_cell) if (cell_contains(c, neut)) then neut%cell = index_cell + cCell => cells(index_cell) + cMaterial => materials(cCell%material) return end if end do @@ -166,6 +168,8 @@ contains c => cells(index_cell) if (cell_contains(c, neut)) then neut%cell = index_cell + cCell => cells(index_cell) + cMaterial => materials(cCell%material) return end if end do @@ -177,6 +181,8 @@ contains c => cells(i) if (cell_contains(c, neut)) then neut%cell = i + cCell => cells(i) + cMaterial => materials(cCell%material) return end if end do @@ -196,7 +202,7 @@ contains subroutine dist_to_boundary(neut, dist, surf, other_cell) - type(Neutron), intent(in) :: neut + type(Neutron), pointer :: neut real(8), intent(out) :: dist integer, intent(out) :: surf integer, optional, intent(in) :: other_cell @@ -224,18 +230,18 @@ contains current_surf = neut%surface - x = neut%xyz(1) - y = neut%xyz(2) - z = neut%xyz(3) - u = neut%uvw(1) - v = neut%uvw(2) - z = neut%uvw(3) - dist = INFINITY n_boundaries = size(cell_p%boundary_list) allocate(expression(n_boundaries)) expression = cell_p%boundary_list do i = 1, n_boundaries + x = neut%xyz(1) + y = neut%xyz(2) + z = neut%xyz(3) + u = neut%uvw(1) + v = neut%uvw(2) + w = neut%uvw(3) + ! check for coincident surfaec index_surf = expression(i) if (index_surf == current_surf) cycle @@ -247,30 +253,30 @@ contains surf_p => surfaces(index_surf) select case (surf_p%type) case (SURF_PX) - if (u == 0.0) then + if (u == ZERO) then d = INFINITY else x0 = surf_p%coeffs(1) d = (x0 - x)/u - if (d < 0) d = INFINITY + if (d < ZERO) d = INFINITY end if case (SURF_PY) - if (v == 0.0) then + if (v == ZERO) then d = INFINITY else y0 = surf_p%coeffs(1) d = (y0 - y)/v - if (d < 0) d = INFINITY + if (d < ZERO) d = INFINITY end if case (SURF_PZ) - if (w == 0.0) then + if (w == ZERO) then d = INFINITY else z0 = surf_p%coeffs(1) d = (z0 - z)/w - if (d < 0.0) d = INFINITY + if (d < ZERO) d = INFINITY end if case (SURF_PLANE) @@ -280,16 +286,16 @@ contains D = surf_p%coeffs(4) tmp = A*u + B*v + C*w - if (tmp == 0.0) then + if (tmp == ZERO) then d = INFINITY else d = -(A*x + B*y + C*w - D)/tmp - if (d < 0.0) d = INFINITY + if (d < ZERO) d = INFINITY end if case (SURF_CYL_X) - a = 1.0 - u**2 ! v^2 + w^2 - if (a == 0.0) then + a = ONE - u**2 ! v^2 + w^2 + if (a == ZERO) then d = INFINITY else y0 = surf_p%coeffs(1) @@ -302,7 +308,7 @@ contains c = y**2 + z**2 - r**2 quad = k**2 - a*c - if (c < 0) then + if (c < ZERO) then ! particle is inside the cylinder, thus one distance ! must be negative and one must be positive. The ! positive distance will be the one with negative sign @@ -317,14 +323,14 @@ contains ! positive sign on sqrt(quad) d = -(k + sqrt(quad))/a - if (d < 0) d = INFINITY + if (d < ZERO) d = INFINITY end if end if case (SURF_CYL_Y) - a = 1.0 - v**2 ! u^2 + w^2 - if (a == 0.0) then + a = ONE - v**2 ! u^2 + w^2 + if (a == ZERO) then d = INFINITY else x0 = surf_p%coeffs(1) @@ -337,7 +343,7 @@ contains c = x**2 + z**2 - r**2 quad = k**2 - a*c - if (c < 0) then + if (c < ZERO) then ! particle is inside the cylinder, thus one distance ! must be negative and one must be positive. The ! positive distance will be the one with negative sign @@ -352,14 +358,14 @@ contains ! positive sign on sqrt(quad) d = -(k + sqrt(quad))/a - if (d < 0) d = INFINITY + if (d < ZERO) d = INFINITY end if end if case (SURF_CYL_Z) - a = 1.0 - w**2 ! u^2 + v^2 - if (a == 0.0) then + a = ONE - w**2 ! u^2 + v^2 + if (a == ZERO) then d = INFINITY else x0 = surf_p%coeffs(1) @@ -372,12 +378,12 @@ contains c = x**2 + y**2 - r**2 quad = k**2 - a*c - if (quad < 0) then + if (quad < ZERO) then ! no intersection with cylinder d = INFINITY - elseif (c < 0) then + elseif (c < ZERO) then ! particle is inside the cylinder, thus one distance ! must be negative and one must be positive. The ! positive distance will be the one with negative sign @@ -392,7 +398,7 @@ contains ! positive sign on sqrt(quad) d = -(k + sqrt(quad))/a - if (d < 0) d = INFINITY + if (d < ZERO) d = INFINITY end if end if @@ -410,12 +416,12 @@ contains c = x**2 + y**2 + z**2 - r**2 quad = k**2 - c - if (quad < 0) then + if (quad < ZERO) then ! no intersection with sphere d = INFINITY - elseif (c < 0) then + elseif (c < ZERO) then ! particle is inside the sphere, thus one distance ! must be negative and one must be positive. The ! positive distance will be the one with negative sign @@ -430,7 +436,7 @@ contains ! positive sign on sqrt(quad) d = -(k + sqrt(quad)) - if (d < 0) d = INFINITY + if (d < ZERO) d = INFINITY end if diff --git a/src/global.f90 b/src/global.f90 index ccc63fccc..d51fc5b9b 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -56,13 +56,15 @@ module global ! Physical constants real(8), parameter :: & - & PI = 2.*acos(0.0), & ! pi - & MASS_NEUTRON = 1.0086649156, & ! mass of a neutron - & MASS_PROTON = 1.00727646677, & ! mass of a proton - & AMU = 1.66053873e-27, & ! 1 amu in kg - & N_AVOGADRO = 0.602214179, & ! Avogadro's number in 10^24/mol - & K_BOLTZMANN = 8.617342e-5, & ! Boltzmann constant in eV/K - & INFINITY = huge(0.0_8) ! positive infinity + & PI = 2.0_8*acos(0.0_8), & ! pi + & MASS_NEUTRON = 1.0086649156, & ! mass of a neutron + & MASS_PROTON = 1.00727646677, & ! mass of a proton + & AMU = 1.66053873e-27, & ! 1 amu in kg + & N_AVOGADRO = 0.602214179, & ! Avogadro's number in 10^24/mol + & K_BOLTZMANN = 8.617342e-5, & ! Boltzmann constant in eV/K + & INFINITY = huge(0.0_8), & ! positive infinity + & ZERO = 0.0_8, & + & ONE = 1.0_8 ! Boundary conditions integer, parameter :: & diff --git a/src/input_sample b/src/input_sample index 4975c01d2..be9c70d37 100644 --- a/src/input_sample +++ b/src/input_sample @@ -1,17 +1,20 @@ # test input file cell 100 40 -1 -cell 200 40 1 -2 +cell 200 41 1 -2 surface 1 sph 0 0 0 3 surface 2 sph 0 0 0 5 -material 40 -20.0 & - 3007.03c 1.0 & +material 40 -4.5 & 92238.03c 1.0 -source box -3 -3 -3 3 3 3 +material 41 -1.0 & + 1001.03c 2.0 & + 8016.03c 1.0 + +source box -1 -1 -1 1 1 1 xs_library endfb7 xs_data /opt/serpent/xsdata/endfb7 -criticality 1 1 50 +criticality 1 1 15 diff --git a/src/main.f90 b/src/main.f90 index 1a9fac72a..df74c88bb 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -20,7 +20,7 @@ program main ! Print the OpenMC title and version/date/time information call title() - verbosity = 9 + verbosity = 10 ! Initialize random number generator call RN_init_problem( 3, 0_8, 0_8, 0_8, 0 ) @@ -68,7 +68,9 @@ program main surfaces(2)%bc = BC_VACUUM call run_problem() - + ! deallocate arrays + call free_memory() + contains !===================================================================== @@ -115,8 +117,6 @@ contains ! print run time - call free_memory() - end subroutine run_problem end program main diff --git a/src/physics.f90 b/src/physics.f90 index 2146e6f96..99f1bb209 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -6,6 +6,7 @@ module physics use mcnp_random, only: rang use output, only: error, message use search, only: binary_search + use endf, only: reaction_name implicit none @@ -29,6 +30,7 @@ contains real(8) :: distance real(8) :: Sigma ! total cross-section real(8) :: f ! interpolation factor + real(8) :: tmp(3) integer :: IE ! index on energy grid ! determine what cell the particle is in @@ -91,18 +93,14 @@ contains type(Neutron), pointer :: neut - type(AceContinuous), pointer :: table + type(AceReaction), pointer :: rxn real(8) :: r1 - real(8) :: phi ! azimuthal angle - real(8) :: mu ! cosine of polar angle character(250) :: msg - integer :: cell_num integer :: i,j integer :: n_isotopes integer :: IE - - real(8) :: f, Sigma + real(8) :: f, Sigma, total real(8) :: density, density_i real(8) :: p real(8), allocatable :: Sigma_t(:) @@ -127,40 +125,149 @@ contains Sigma_t(i) = Sigma end do - ! normalize to create a discrete pdf - Sigma_t = Sigma_t / sum(Sigma_t) - ! sample nuclide r1 = rang() p = 0.0_8 + total = sum(Sigma_t) do i = 1, n_isotopes - p = p + Sigma_t(i) + p = p + Sigma_t(i) / total if (r1 < p) exit end do - table => xs_continuous(cMaterial%table(i)) - ! print *, 'sampled nuclide ', i - ! select collision type - r1 = rang() - if (r1 <= 0.5) then - ! scatter - phi = 2.*pi*rang() - mu = 2.*rang() - 1 - neut%uvw(1) = mu - neut%uvw(2) = sqrt(1. - mu**2) * cos(phi) - neut%uvw(3) = sqrt(1. - mu**2) * sin(phi) - else - neut%alive = .false. - if (verbosity >= 10) then - cell_num = cells(neut%cell)%uid - msg = " Absorbed in cell " // trim(int_to_str(cell_num)) - call message(msg, 10) - end if - return + ! Get table, total xs, interpolation factor + table => xs_continuous(cMaterial%table(i)) + Sigma = Sigma_t(i) + IE = table%grid_index(neut%IE) + f = (neut%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) + density = cMaterial%atom_percent(i)*density + + ! free memory + deallocate(Sigma_t) + + ! sample reaction channel + r1 = rang()*Sigma + p = 0.0_8 + do i = 1, table%n_reaction + rxn => table%reactions(i) + if (rxn%MT >= 200) cycle + if (IE < rxn%IE) cycle + p = p + density * (f*rxn%sigma(IE-rxn%IE+1) + (1-f)*(rxn%sigma(IE-rxn%IE+2))) + if (r1 < p) exit + end do + if (verbosity >= 10) then + msg = " " // trim(reaction_name(rxn%MT)) // " with nuclide " // & + & trim(table%name) + call message(msg, 10) end if + + ! call appropriate subroutine + select case (rxn%MT) + case (2) + call elastic_scatter(neut, table%awr) + case (102) + call n_gamma(neut) + case default + call elastic_scatter(neut, table%awr) + end select end subroutine collision +!===================================================================== +! ELASTIC_SCATTER +!===================================================================== + + subroutine elastic_scatter(neut, awr) + + type(Neutron), pointer :: neut + real(8), intent(in) :: awr + + real(8) :: phi ! azimuthal angle + real(8) :: mu ! cosine of polar angle + real(8) :: vx, vy, vz + real(8) :: vcx, vcy ,vcz + real(8) :: vel + real(8) :: u, v, w + real(8) :: E + integer :: IE + + vel = sqrt(neut%E) + + vx = vel*neut%uvw(1) + vy = vel*neut%uvw(2) + vz = vel*neut%uvw(3) + + + vcx = vx/(awr + 1.0_8) + vcy = vy/(awr + 1.0_8) + vcz = vz/(awr + 1.0_8) + + ! Transform to CM frame + vx = vx - vcx + vy = vy - vcy + vz = vz - vcz + + vel = sqrt(vx*vx + vy*vy + vz*vz) + + ! Select isotropic direcion -- this is only valid for s-wave + ! scattering + phi = 2.0_8*PI*rang() + mu = 2.0_8*rang() - 1.0_8 + u = mu + v = sqrt(1.0_8 - mu**2) * cos(phi) + w = sqrt(1.0_8 - mu**2) * sin(phi) + + vx = u*vel + vy = v*vel + vz = w*vel + + ! Transform back to LAB frame + vx = vx + vcx + vy = vy + vcy + vz = vz + vcz + + E = vx*vx + vy*vy + vz*vz + vel = sqrt(E) + + neut%E = E + neut%uvw(1) = vx/vel + neut%uvw(2) = vy/vel + neut%uvw(3) = vz/vel + + ! find energy index, interpolation factor + IE = binary_search(e_grid, n_grid, E) + neut%IE = IE + neut%interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) + + end subroutine elastic_scatter + +!===================================================================== +! LEVEL_INELASTIC +!===================================================================== + + subroutine level_inelastic + + end subroutine level_inelastic + +!===================================================================== +! N_GAMMA +!===================================================================== + + subroutine n_gamma(neut) + + type(Neutron), pointer :: neut + + integer :: cell_num + character(250) :: msg + + neut%alive = .false. + if (verbosity >= 10) then + cell_num = cells(neut%cell)%uid + msg = " Absorbed in cell " // trim(int_to_str(cell_num)) + call message(msg, 10) + end if + + end subroutine n_gamma + !===================================================================== ! MAXWELL_SPECTRUM samples an energy from the Maxwell fission ! distribution based on a rejection sampling scheme. This is described diff --git a/src/types.f90 b/src/types.f90 index 0c0d7cce6..0e4cc82cf 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -107,7 +107,7 @@ module types integer :: MT real(8) :: Q_value real(8) :: TY - integer :: energy_index + integer :: IE real(8), allocatable :: sigma(:) real(8), allocatable :: ang_cos(:,:) real(8), allocatable :: ang_pdf(:,:) @@ -145,7 +145,8 @@ module types real(8), allocatable :: nu_d_precursor_const(:,:) real(8), allocatable :: nu_d_precursor_energy(:,:) real(8), allocatable :: nu_d_precursor_prob(:,:) - type(AceReaction), pointer :: reactions(:) + integer :: n_reaction + type(AceReaction), pointer :: reactions(:) => null() end type AceContinuous @@ -214,7 +215,7 @@ module types !===================================================================== type ListKeyValueCI - type(ListKeyValueCI), pointer :: next + type(ListKeyValueCI), pointer :: next => null() type(KeyValueCI) :: data end type ListKeyValueCI @@ -224,7 +225,7 @@ module types !===================================================================== type ListKeyValueII - type(ListKeyValueII), pointer :: next + type(ListKeyValueII), pointer :: next => null() type(KeyValueII) :: data end type ListKeyValueII @@ -234,7 +235,7 @@ module types !===================================================================== type ListReal - type(ListReal), pointer :: next + type(ListReal), pointer :: next => null() real(8) :: data end type ListReal @@ -243,7 +244,7 @@ module types !===================================================================== type ListInt - type(ListInt), pointer :: next + type(ListInt), pointer :: next => null() integer :: data end type ListInt @@ -253,7 +254,7 @@ module types !===================================================================== type HashListCI - type(ListKeyValueCI), pointer :: list + type(ListKeyValueCI), pointer :: list => null() end type HashListCI !===================================================================== @@ -262,7 +263,7 @@ module types !===================================================================== type HashListII - type(ListKeyValueII), pointer :: list + type(ListKeyValueII), pointer :: list => null() end type HashListII !===================================================================== @@ -271,7 +272,7 @@ module types !===================================================================== type DictionaryCI - type(HashListCI), pointer :: table(:) + type(HashListCI), pointer :: table(:) => null() end type DictionaryCI !===================================================================== @@ -280,7 +281,7 @@ module types !===================================================================== type DictionaryII - type(HashListII), pointer :: table(:) + type(HashListII), pointer :: table(:) => null() end type DictionaryII end module types