From 66f9daba06c7f1f3989b8bfe5eac94f1880bfe57 Mon Sep 17 00:00:00 2001
From: Wiarda <wiardada@ornl.gov>
Date: Fri, 10 Sep 2021 11:03:34 -0400
Subject: [PATCH] Move file 3 reading for Maxwellian average data to C++
---
sammy/src/blk/Exploc_common.f90 | 28 +--
sammy/src/cro/mcro0.f | 30 +--
sammy/src/cro/mcro1.f | 105 +++-------
sammy/src/dat/mdat0.f90 | 17 +-
sammy/src/dat/mdat4.f90 | 16 +-
sammy/src/mlb/mmlb0.f | 16 +-
sammy/src/mlb/mmlb1.f | 8 -
sammy/src/mxw/mmxw5.f | 7 +-
sammy/src/mxw/mmxw6.f | 81 +++++---
sammy/src/rec/mrec0.f | 37 ++--
sammy/src/rec/mrec7.f | 170 ---------------
sammy/src/rec/mrec8.f90 | 52 +++--
sammy/src/sammy/CMakeLists.txt | 1 -
sammy/src/the/ZeroKCrossCorrections_M.f90 | 27 +++
sammy/src/xct/mxct0.f90 | 51 ++---
sammy/src/xct/mxct02.f90 | 242 +++++++++-------------
16 files changed, 297 insertions(+), 591 deletions(-)
delete mode 100644 sammy/src/rec/mrec7.f
diff --git a/sammy/src/blk/Exploc_common.f90 b/sammy/src/blk/Exploc_common.f90
index 751db21d4..5df315dba 100644
--- a/sammy/src/blk/Exploc_common.f90
+++ b/sammy/src/blk/Exploc_common.f90
@@ -110,12 +110,8 @@ module exploc_common_m
integer :: Nuncer ! number of uncertainties, not an array index
real(kind=8),allocatable,dimension(:)::A_Iuncs
integer,allocatable,dimension(:)::I_Ijuncs
- real(kind=8),allocatable,dimension(:)::A_Iadder
- real(kind=8),allocatable,dimension(:)::A_Iaddcr
- integer,allocatable,dimension(:)::I_Inbt
- ! old group 9
- integer,allocatable,dimension(:)::I_Iint
+ ! old group 9
real(kind=8),allocatable,dimension(:)::A_Iresol
integer,allocatable,dimension(:)::I_Iisopa
real(kind=8),allocatable,dimension(:)::A_Iccoul
@@ -505,28 +501,6 @@ module exploc_common_m
call allocate_integer_data(I_Ijuncs,want)
end subroutine make_I_Ijuncs
-
- subroutine make_I_Iint(want)
- integer::want
- call allocate_integer_data(I_Iint,want)
- end subroutine make_I_Iint
-
- subroutine make_I_Inbt(want)
- integer::want
- call allocate_integer_data(I_Inbt,want)
- end subroutine make_I_Inbt
-
-
- subroutine make_A_Iadder(want)
- integer::want
- call allocate_real_data(A_Iadder,want)
- end subroutine make_A_Iadder
-
- subroutine make_A_Iaddcr(want)
- integer::want
- call allocate_real_data(A_Iaddcr,want)
- end subroutine make_A_Iaddcr
-
subroutine make_A_Iresol(want)
integer::want
call allocate_real_data(A_Iresol,want)
diff --git a/sammy/src/cro/mcro0.f b/sammy/src/cro/mcro0.f
index f2a54ccd1..0cb6c53dc 100644
--- a/sammy/src/cro/mcro0.f
+++ b/sammy/src/cro/mcro0.f
@@ -52,19 +52,15 @@ C This is for eta-data only
C
C
C *** organize for broadening
- Icr = Kcros
IF (Nnniso.NE.1) STOP '[STOP in Samcro_0 in cro/mcro0.f]'
C
C *** find array sizes for xx
IF (resParData%getNumResonances().GT.0) CALL Findmx (Mxany)
C
c *** find array sizes for added cross section from endf/b-vi
- IF (Kaddcr.EQ.1) THEN
- CALL Read3 (Nrfil3, Npfil3)
- ELSE
- Npfil3 = 1
- Nrfil3 = 1
- END IF
+C No done in C++
+ Npfil3 = 1
+ Nrfil3 = 1
C *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-Cross section
call setAuxGridOffset(1) ! reset auxillary grid offset
call setAuxGridRowMax(0)
@@ -107,16 +103,6 @@ C
C
C *** four ***
C - - - - - - - - - - - - - - - - - - - - - - <
-C
- IF (Kaddcr.EQ.1) THEN
-C *** Andy ***
- call make_A_Iadder(Np)
- call make_A_Iaddcr(Np)
- call make_I_Inbt(Nr)
- call make_I_Iint(Nr)
- CALL Rdndf3 (A_Iadder , A_Iaddcr , I_Inbt , I_Iint , Nrfil3,
- * Npfil3)
- END IF
C
C *** five ***
call allocate_real_data(A_Ialphr, Nres)
@@ -167,12 +153,9 @@ C
C
C *** eight ***
call allocate_real_data(A_Idum, N8)
- CALL Work_Cro (I_Iflmsc , A_Iprnbk , I_Iflnbk , A_Iprbgf ,
- * I_Iflbgf , I_Indbgf , A_Ibgfmi , A_Ibgfma , A_Itexbg ,
- * A_Iteabg , A_Ith ,
+ CALL Work_Cro (I_Iflmsc ,
* A_Isigxx , A_Idasig , A_Idbsig ,
- * A_Ipiece , A_Idum , A_Iadder , A_Iaddcr ,
- * I_Inbt , I_Iint)
+ * A_Ipiece , A_Idum)
C *** Sbroutine Work_Cro generates theory and derivatives
deallocate(A_Idum)
deallocate(A_Ics)
@@ -180,9 +163,6 @@ C *** Sbroutine Work_Cro generates theory and derivatives
deallocate(A_Idphi)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
C
-C *** If making fake data, write it out, then stop
- IF (Kfake.EQ.1) CALL Newdat (A_Ith )
- IF (Kfake.EQ.1) STOP '[Normal stop in mcro0.f]'
C
Nnnsig = 1
Kiniso = Kkkiso
diff --git a/sammy/src/cro/mcro1.f b/sammy/src/cro/mcro1.f
index 2283ba4e2..33fc1c319 100644
--- a/sammy/src/cro/mcro1.f
+++ b/sammy/src/cro/mcro1.f
@@ -1,10 +1,9 @@
C
C --------------------------------------------------------------
C
- SUBROUTINE Work_Cro (Iflmsc, Parnbk, Iflnbk, Parbgf, Iflbgf,
- * Kndbgf, Bgfmin, Bgfmax, Texbgf, Teabgf, Theory,
+ SUBROUTINE Work_Cro (Iflmsc,
* Sigxxx, Dasigx, Dbsigx,
- * Pieces, Dum, Adderg, Addcro, Nbt, Int)
+ * Pieces, Dum)
C
C *** PURPOSE -- GENERATE THEORETICAL DATA "theory" AND PARTIAL
C *** DERIVATIVES "GB"
@@ -20,22 +19,20 @@ C
use EndfData_common_m
use SammyGridAccess_M
use array_sizes_common_m, only : calcData
+ use convert_to_transmission_m
+ use normalize_and_background
IMPLICIT DOUBLE PRECISION (a-h,o-z)
LOGICAL Ywhich
C
type(SammyGridAccess)::grid
- DIMENSION Iflmsc(*), Parnbk(*), Iflnbk(*), Parbgf(*), Iflbgf(*),
- * Kndbgf(*), Bgfmin(*), Bgfmax(*), Texbgf(Ntepnt,*),
- * Teabgf(Ntepnt,*),
- * Theory(Nnnsig,*), Pieces(Ngroup,*), Dum(*), Adderg(*),
- * Addcro(*), Nbt(*), Int(*)
+ logical::wantDeriv
+ DIMENSION Iflmsc(*),
+ * Pieces(Ngroup,*), Dum(*)
DIMENSION Sigxxx(Nnnsig,*),
* Dasigx(Nnnsig,*), Dbsigx(Nnnsig,Ndbxxx,*)
C
-C DIMENSION Iflmsc(Nummsc), Parnbk(Numnbk), Iflnbk(Numnbk),
-c * Parbgf(Numbgf), Iflbgf(*), Kndbgf(*), Bgfmin(*), Bgfmax(*),
-C * Theory(Numcro*Ndat), Pieces(Ngroup,Ndatb), Dum(Ndatb),
-C * Adderg(Np), Addcro(Np), Nbt(Nrfil3), Int(Nrfil3)
+C DIMENSION Iflmsc(Nummsc),
+C * Pieces(Ngroup,Ndatb), Dum(Ndatb)
C
DATA Zero /0.0d0/
C
@@ -45,6 +42,9 @@ C
numE = grid%getNumEnergies(expData)
C
Ywhich = Ydoppr.OR.Yresol
+
+ wantDeriv = Ksolve.NE.2
+ if ( covData%getPupedParam().gt.0) wantDeriv = .true.
C
IF (Kpiece.EQ.1) CALL Zero_Array (Pieces, Ngroup*numE)
Kkkkkk = 0
@@ -58,8 +58,7 @@ C *** positive energies
Jdat = J
do Jsig = 1,Nnnsig
Jcount = (Jdat-1)*Nnnsig + Jsig
- call calcData%setToZero(Jcount, Ndasig + Ndbsig)
- Theory(Jsig, Jdat) = 0.0d0
+ call calcData%setToZero(Jcount, Ndasig + Ndbsig)
end do
CALL Zero_Array (Sigxxx, Nnnsig*Nnniso)
@@ -75,11 +74,6 @@ C
Squ = dSQRT(Su)
Ipoten = 0
IF (Kpoten.NE.0 .AND. (J.EQ.1 .OR. J.EQ.numE)) THEN
-10000 FORMAT(/, ' Potential Scattering for energy',
- * 1PE13.6, /,
- * ' group l for this group for this l')
- WRITE (6,10000) Su
- WRITE (21,10000) Su
Ipoten = 1
END IF
C
@@ -92,68 +86,21 @@ C ********* Generate cross sections and derivatives
ENDIF
C
END IF
-C
-C
-C *** add endf/b file 3 to cross sections
- IF (Kaddcr.EQ.1) THEN
- Gbx = Zero
- CALL Adcros (grid%getEnergy(Jdat, expData),
- * Gbx, Adderg, Addcro,
- * Nbt, Int, Nrfil3, Npfil3, Jdat, numE)
- Sigxxx(1,1) = Sigxxx(1,1) + Gbx
- END IF
-C
-C ********* if want Maxwellian averages or energy-averages with no
-C ********* correction terms, do same as if there is broadening
-C ********* except don't convert to transmission unless needed
- IF (Maxwel.EQ.1 .OR. Knocor.EQ.1) THEN
- GO TO 90
- END IF
-C
-C ********* if there is doppler broadening, GO TO 90
- IF (Ydoppr) GO TO 90
-C
-C ********* if no Doppler, and this is transmission, make conversion
-C
-C ********* if there is resolution broadening or angular
-C distributions, do nothing else here
- IF (Yresol .OR. Yangle) GO TO 90
-C
-C ********* if ETA is being calculated and there is no broadening,
-C ********* find derivative with respect to nu
-C
-C ********* If this is fake data, store it and move on
- IF (Kfake.EQ.1) THEN
- Theory(1, J) = Sigxxx(1,1)
- cycle
- END IF
-C
-C ********* if there is normalization or background, include it
-C
-C ********* Write results onto theory if there is no broadening etc
- IF (Jjjdop.NE.1) THEN
- Kkkkkk = Kkkkkk + 1
- ! cro doesn't have isotopes, therefore the
- ! isotope is denoted as -1, to indicate that
- ! this already is the sum over all isotopes
- call calcData%addCalculatedData(Kkkkkk, Nnnsig, ndasig,
+
+
+ IF (Maxwel.EQ.1 .OR. Knocor.EQ.1.or.
+ * Ydoppr.or.Yresol .OR. Yangle.or.
+ * Kfake.eq.1) then
+ Kkkkkk = Jdat
+ else
+ Kkkkkk = Kkkkkk + 1
+ end if
+
+ call calcData%addCalculatedData(Kkkkkk, Nnnsig, ndasig,
* ndbsig, -1, Sigxxx, Dasigx, Dbsigx)
- END IF
- cycle
-C
- 90 CONTINUE
-C ********* copy results if there is broadening
- Kkkkkk = Kkkkkk + 1
- IF (Kkkmin.EQ.0) Kkkmin = Jdat
- ! cro doesn't have isotopes
- call calcData%addCalculatedData(jdat, Nnnsig, ndasig,
- * ndbsig, -1, Sigxxx, Dasigx, Dbsigx)
- end do
-C ****** end of do-loop on energies
-C
-C
C
- IF (Kpiece.EQ.1) CALL Odfpcs (Dum, Pieces)
+ end do
+ IF (Kpiece.EQ.1) CALL Odfpcs (Dum, Pieces)
C
call grid%destroy()
RETURN
diff --git a/sammy/src/dat/mdat0.f90 b/sammy/src/dat/mdat0.f90
index bc26a117c..9597460aa 100644
--- a/sammy/src/dat/mdat0.f90
+++ b/sammy/src/dat/mdat0.f90
@@ -61,12 +61,8 @@ module Samdat_0_M
!
!
! *** find array sizes for reading endf/b-vi file 3 cross sections
- IF (Kaddcr.EQ.1) THEN
- CALL Read3 (Nr, Np)
- ELSE
- Nr = 1
- Np = 1
- END IF
+ Nr = 1
+ Np = 1
!
IF (Medrpi.GT.0) THEN
N = Medrpi
@@ -279,18 +275,13 @@ module Samdat_0_M
Mnp = Np
Mnr = Nr
IF (Mnp.EQ.0) Mnp = 1
- IF (Mnr.EQ.0) Mnr = 1
- call make_A_Iadder(Mnp)
- call make_A_Iaddcr(Mnp)
- call make_I_Inbt(Mnr)
- call make_I_Iint(Mnr)
+ IF (Mnr.EQ.0) Mnr = 1
Mnr = 100*Mnr
IF (Nr.EQ.0) Mnr = 1
CALL Escale ( &
A_Ie, A_Ienerb, A_Ispken , &
- A_Iswidt , A_Isadd , A_Iadder , A_Iaddcr , &
- I_Inbt , I_Iint , Nr, Np)
+ A_Iswidt , A_Isadd)
deallocate(A_Ispken)
deallocate(A_Iswidt)
deallocate(A_Isadd)
diff --git a/sammy/src/dat/mdat4.f90 b/sammy/src/dat/mdat4.f90
index 7013137c8..dcd338baf 100644
--- a/sammy/src/dat/mdat4.f90
+++ b/sammy/src/dat/mdat4.f90
@@ -8,9 +8,7 @@ contains
!
! --------------------------------------------------------------
!
- SUBROUTINE Escale (Ee, &
- Energb, Spken, Swidth, Sadd, Adder, Addcr, Nbt, Int, &
- Nr, Np)
+ SUBROUTINE Escale (Ee, Energb, Spken, Swidth, Sadd)
!
! *** Purpose -- Generate Energb, the set of energies to be used for
! *** broadening the data (i.e., the auxiliary grid)
@@ -21,18 +19,19 @@ contains
use broad_common_m
use EndfData_common_m
use SammyGridAccess_M
+ use mthe0_M
IMPLICIT DOUBLE PRECISION (a-h,o-z)
!
type(SammyGridAccess)::grid
type(GridData)::grid_gb, grid_gm
real(kind=8),allocatable,dimension(:)::Ee(:), Energb
DIMENSION &
- Spken(*), Swidth(*), Sadd(*), Addcr(*), Adder(*), Nbt(*), Int(*)
+ Spken(*), Swidth(*), Sadd(*)
!
! DIMENSION
! * Ee(Ndatb),
! * Spken(Nres), Swidth(Nres),
-! * Sadd(*), Addcr(Np), Adder(Nr), Nbt(Nr), Int(Nr)
+! * Sadd(*)
!
DATA Zero /0.0d0/, One /1.0d0/
!
@@ -151,12 +150,11 @@ contains
!
! *** Add points for ENDF file 3 contributions, if needed
IF (Kaddcr.EQ.1) THEN
- M = Ndatb
- CALL Rdndf3 (Adder, Addcr, Nbt, Int, Nr, Np)
+ call setup_zeroK()
+ M = Ndatb
Mm = M
Eee = grid_gb%getData(Ndatb, 1)
- CALL Urrgrd (Adder, Nbt, Int, Nr, Np, Eee, &
- grid_gb, Npurr, Mm)
+ CALL Urrgrd ( Eee, grid_gb, Npurr, Mm)
Ndatb = M + Npurr
END IF
!
diff --git a/sammy/src/mlb/mmlb0.f b/sammy/src/mlb/mmlb0.f
index 6b1f42af2..624e48073 100644
--- a/sammy/src/mlb/mmlb0.f
+++ b/sammy/src/mlb/mmlb0.f
@@ -59,12 +59,8 @@ C
C *** organize for broadening
C
Icr = Kcros
- IF (Kaddcr.EQ.1) THEN
- CALL Read3 (Nrfil3, Npfil3)
- ELSE
- Npfil3 = 1
- Nrfil3 = 1
- END IF
+ Npfil3 = 1
+ Nrfil3 = 1
C
C *** guesstimate size of array needed for SAMMY-MLBW
call setAuxGridOffset(1) ! reset starting point for auxillary grid
@@ -93,14 +89,6 @@ C ### three ###
call allocate_real_data(A_Idum, N7)
C
C
- IF (Kaddcr.EQ.1) THEN
- call make_A_Iadder(Npfil3)
- call make_A_Iaddcr(Npfil3)
- call make_I_Inbt(Nrfil3)
- call make_I_Iint(Nrfil3)
- CALL Rdndf3 (A_Iadder , A_Iaddcr , I_Inbt , I_Iint , Nrfil3,
- * Npfil3)
- END IF
C
C *** Sub Work generates Theory and derivatives
CALL Work_Mlb ( I_Iflmsc , A_Iprnbk , I_Iflnbk , A_Iprbgf ,
diff --git a/sammy/src/mlb/mmlb1.f b/sammy/src/mlb/mmlb1.f
index e9d13eab7..50fe99457 100644
--- a/sammy/src/mlb/mmlb1.f
+++ b/sammy/src/mlb/mmlb1.f
@@ -105,14 +105,6 @@ C ********* If this is fake data, store it and move on
END IF
C
C
-C *** Add endf/b file 3 to cross sections
- IF (Iskip.EQ.0 .AND. Kaddcr.EQ.1) THEN
- Gbx = Zero
- ee = grid%getEnergy(Jdat, expData)
- CALL Adcros (ee, Gbx, Adderg, Addcro,
- * Nbt, Int, Nrfil3, Npfil3, Jdat, numEl)
- Sigxxx(1,1) = Sigxxx(1,1) + Gbx
- END IF
C
C ********* If want Maxwellian averages or energy-averages with no
C ********* correction terms, do same as if there is broadening
diff --git a/sammy/src/mxw/mmxw5.f b/sammy/src/mxw/mmxw5.f
index 44bb8d3c2..061ea1a68 100644
--- a/sammy/src/mxw/mmxw5.f
+++ b/sammy/src/mxw/mmxw5.f
@@ -76,8 +76,8 @@ C
C *** Findem_Mxw locates the highest energy at which there is a resonance
CALL Findem_Mxw (Caytmp, Resmax, Ndatq)
C
- Addcrx = Zero
- IF (Kaddcr.EQ.1) CALL Fil3xx (Addcrx, A_Iadder , A_Iaddcr ,Npfil3)
+ Addcrx = Zer
+ IF (Kaddcr.EQ.1) CALL Fil3xx (Addcrx)
C
Amass = Aaawww
Amass = Amass/(Amass+A_Mass_Small)
@@ -223,8 +223,7 @@ C
C ??? isotope-dependent ?
DO J=1,Ndatq
Emm = Caytmp(J)/Amass
- CALL Fil3mx (A_Iadder , A_Iaddcr , I_Inbt ,
- * I_Iint , Nrfil3, Npfil3, Emm, Sumfl3(J), Tosp)
+ CALL Fil3mx (Emm, Sumfl3(J), Tosp)
END DO
END IF
C
diff --git a/sammy/src/mxw/mmxw6.f b/sammy/src/mxw/mmxw6.f
index 809520fe7..83ce32589 100644
--- a/sammy/src/mxw/mmxw6.f
+++ b/sammy/src/mxw/mmxw6.f
@@ -128,25 +128,52 @@ C
C
C ------------------------------------------------------------------
C
- SUBROUTINE Fil3xx (Addcrx, Eendf, Cendf, Np)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION Eendf(*), Cendf(*)
+ SUBROUTINE Fil3xx (Addcrx)
+ use array_sizes_common_m, only : zeroKCalc
+ use Tab1_M
+ IMPLICIT None
+ real(kind=8)::Addcrx
+ integer::np
+ type(Tab1) :: data
+ real(kind=8)::x(1),e(1)
+ real(kind=8)::e3,x3, e4,x4, e2, x2, e1,x1
C
C *** purpose -- figure energy at which file-3 contribution becomes
C *** important (so do not stop integrating before THEN!)
C
+ call zeroKCalc%endfReader%crossAddedToStellarAverage(data)
+ Np = data%getPoints()
+
IF (Np.GT.3) THEN
- Addcrx = Eendf(3)
- IF (Cendf(3).LT.Cendf(4)) Addcrx = Eendf(4)
+ call data%getValue(3, e, x)
+ e3 = e(1)
+ x3 = x(1)
+ call data%getValue(4, e, x)
+ e4 = e(1)
+ x4 = x(1)
+ Addcrx = e3
+ IF (x3.LT.x4) Addcrx = e4
ELSE IF (Np.EQ.3) THEN
- Addcrx = Eendf(2)
- IF (Cendf(2).LT.Cendf(3)) Addcrx = Eendf(3)
+ call data%getValue(2, e, x)
+ e2 = e(1)
+ x2 = x(1)
+ call data%getValue(3, e, x)
+ e3 = e(1)
+ x3 = x(1)
+ Addcrx = e2
+ IF (x2.LT.x3) Addcrx = e3
ELSE IF (Np.EQ.2) THEN
- Addcrx = Eendf(2)
+ call data%getValue(2, e, x)
+ e2 = e(1)
+ x2 = x(1)
+ Addcrx = e2
ELSE
Write (6,10000) Np
10000 FORMAT (' Possible problem -- Np =', i3)
- Addcrx = Eendf(1)
+ call data%getValue(1, e, x)
+ e1 = e(1)
+ x1 = x(1)
+ Addcrx = e1
END IF
RETURN
END
@@ -154,22 +181,28 @@ C
C
C ------------------------------------------------------------------
C
- SUBROUTINE Fil3mx (Eendf, Cendf, Nbt, INT, Nr, Np, Emm, Sumfl3,
- * Tosp)
+ SUBROUTINE Fil3mx (Emm, Sumfl3,Tosp)
use abcerf_m
+ use array_sizes_common_m, only : zeroKCalc
+ use Tab1_M
IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION Eendf(*), Cendf(*), Nbt(*), Int(*)
+ type(Tab1) :: data
+ real(kind=8)::x(1),e(1)
DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
C
C *** purpose -- figure contribution from cross sections from endf/b file 3
C
Sum = Zero
Istart = 1
- IF (Eendf(1).EQ.Zero) Istart = 2
+ call zeroKCalc%endfReader%crossAddedToStellarAverage(data)
+ np = data%getPoints()
+ call data%getValue(1, e, x)
+ IF (e(1).EQ.Zero) Istart = 2
C *** do integral from Zero to Eendf(Istart), assuming 1/v
- Y1 = Cendf(Istart)
+ call data%getValue(Istart, e, x)
+ Y1 = x(1)
IF (Y1.ne.Zero) THEN
- X1 = Eendf(Istart)/Emm
+ X1 = e(1)/Emm
Aa = dSQRT(X1)
Aa = -X1 + Aa * Abcerf (Aa, Aa, Qa, Qb, Qc, N)
E1 = dEXP(-X1)
@@ -178,19 +211,15 @@ C *** do integral from Zero to Eendf(Istart), assuming 1/v
Istart = Istart + 1
C *** now do rest of integration, assuming linear (or constant)
DO I=Istart,Np
- X1 = Eendf(I-1)/Emm
- X2 = Eendf(I )/Emm
- Y1 = Cendf(I-1)
- Y2 = Cendf(I )
+ call data%getValue(I-1, e, x)
+ X1 = e(1)/Emm
+ Y1 = x(1)
+ call data%getValue(I, e, x)
+ X2 = e(1)/Emm
+ Y2 = x(1)
IF (X1.NE.X2) THEN
C
- DO L=1,NR
- IF (Nbt(L).GE.I) GO TO 20
- END DO
- STOP '[STOP -- I > Nbt(NR) in Fil3mx in mxw/mmxw6.f]'
-C
- 20 CONTINUE
- Intl = Int(l)
+ Intl = data%getInter(I)
C
IF (Intl.EQ.1) THEN
C *** Y is constant in X
diff --git a/sammy/src/rec/mrec0.f b/sammy/src/rec/mrec0.f
index 2b38b0f4d..2dabca5cf 100644
--- a/sammy/src/rec/mrec0.f
+++ b/sammy/src/rec/mrec0.f
@@ -20,10 +20,15 @@ C
use mrec8_m
use mdata_M
use mrec2_m
+ use mthe0_M
+ use Tab1_M
+ use array_sizes_common_m, only : zeroKCalc, zeroKCalcInit
IMPLICIT DOUBLE PRECISION (a-h,o-z)
type(GridData)::energb
real(kind=8),allocatable,dimension(:)::A_Ietab2, A_Ietx
real(kind=8),allocatable,dimension(:)::A_Isig2(:,:)
+ type(Tab1)::data
+ real(kind=8)::val
real(kind=8),allocatable,dimension(:)::A_ix1, A_Ienode, A_Iwnode
C
C
@@ -38,23 +43,26 @@ C
C
C *** find array sizes for xx
CALL Findmx (Mxany)
+ call setup_zeroK()
+
C
C *** find array sizes for endf/b-vi file 3 cross sections
IF (Kaddcr.EQ.1) THEN
- CALL Read3 (Nrfil3, Npfil3)
- ELSE
- Nrfil3 = 1
- Npfil3 = 1
- END IF
+ if (.not.zeroKCalc%readerInit) then
+ call zeroKCalc%endfReader%initialize()
+ call zeroKCalc%endfReader%setFileName(
+ * trim(zeroKCalc%endfFile))
+ zeroKCalc%readerInit = .true.
+ end if
+ call zeroKCalc%endfReader%crossAddedToStellarAverage(data)
+ End if
+ Nrfil3 = 1
+ Npfil3 = 1
C
C *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-RECONSTRUCT
CALL Estrec (N1, N2, NA, N3, NN, N5, N6, Nesave, Nesub, Nreact,
* Mxany, Nemax)
C *** zero ***
- call make_A_Iadder(Npfil3+1) ! adcros sometimes uses one extra index
- call make_A_Iaddcr(Npfil3+1)
- call make_I_Inbt(Nrfil3)
- call make_I_Iint(Nrfil3)
C
C *** one ***
C *** initialize difmax
@@ -161,9 +169,7 @@ C
C
C *** add cross sections for endf/b file 3
C - - - - - - - - - - - - - - -<
- IF (Kaddcr.EQ.1) THEN
- CALL Rdndf3 (A_Iadder , A_Iaddcr , I_Inbt , I_Iint ,
- * Nrfil3, Npfil3)
+ IF (Kaddcr.EQ.1) THEN
Maxxur = Maxnpu - Nemax
Npurr = Maxxur
I = Idimen (1, 1, '1, 1')
@@ -171,8 +177,7 @@ C - - - - - - - - - - - - - - -<
C - - - - - - - - <
C *** Urrgrd creates energy grid for unresolved resonance region
call energb%initialize()
- CALL Urrgrd (A_Iadder , I_Inbt , I_Iint , Nrfil3, Npfil3,
- * Emax, energb, Npurr, 0)
+ CALL Urrgrd (Emax, energb, Npurr, 0)
call reallocate_real_data(A_Ietab2, Nemax + Npurr)
do i = 1, npurr
A_Ietab2(Nemax + i) = Energb%getData(i, 1)
@@ -227,8 +232,8 @@ C
DO J=1,Nemax
C *** Adcros will add contribution of capture cross section from
C *** ENDF/B file 3 to Sig2(2,*)
- CALL Adcros (A_Ietab2(J),A_Isig2(2, J), A_Iadder ,
- * A_Iaddcr , I_Inbt , I_Iint , Nrfil3, Npfil3, J, Nemax)
+ val = data%interValue(A_Ietab2(J))
+ A_Isig2(2, J) = A_Isig2(2, J) + val
END DO
C - - - - - - - - - - - - - - ->
END IF
diff --git a/sammy/src/rec/mrec7.f b/sammy/src/rec/mrec7.f
deleted file mode 100644
index 4fe0d38d4..000000000
--- a/sammy/src/rec/mrec7.f
+++ /dev/null
@@ -1,170 +0,0 @@
-C
-C -----------------------------------------------------------------
-C
- SUBROUTINE Read3 (NR, NP)
- use samxxx_common_m
- use namfil_common_m
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- CHARACTER*75 Alf
-C
-C *** Purpose -- read first line of endf/b file 3 for dimensions
-C
-cx CALL Oldopn (60, Sam60x, 0)
-cx READ (60, 10000, END=100) Fxxxxx
-cx10000 FORMAT (A70)
-cx 100 CLOSE (UNIT=60)
-CX CALL Filopn (10, Fxxxxx, 0)
-C
- CALL Filopn (10, Faddcr, 0)
- READ (10,10100) Alf, Ns
-10100 FORMAT (A75, I5)
- IF (Ns.EQ.1) REWIND (UNIT=10)
- 10 CONTINUE
- READ (10, 99999, ERR=10) C1, C2, L1, L2, N1, N2,
- * Mat, Mf, Mt, Ns
-99999 FORMAT (2E11.0, 4I11, I4, I2, I3, I5)
-C
- READ (10, 99998) C1, C2, L1, L2, Nr, Np, Mat, Mf, Mt, Ns
-99998 FORMAT (2E11.0, 4I11, I4, I2, I3, I5)
- RETURN
- END
-C
-C
-C -----------------------------------------------------------------
-C
- SUBROUTINE Rdndf3 (Eendf, Cendf, Nbt, Int, Nr, Np)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION Nbt(*), Int(*), Eendf(*), Cendf(*)
-C *** Purpose -- read endf/b file 3
-C
- Nr6 = Nr/3
- Min = 1
- Max = 3
- IF (Nr6.GT.0) THEN
- DO I=1,Nr6
- READ (10,99998) (Nbt(N),Int(N),N=Min,Max), Mat, Mf, Mt, Ns
- Min = Max + 1
- Max = Max + 3
- IF (Max.GT.Nr) GO TO 20
- END DO
- 20 CONTINUE
- END IF
- IF (Nr-Min+1.EQ.2) READ (10,99994) (Nbt(N),Int(N),N=Min,Nr),
- * Mat, Mf, Mt, Ns
- IF (Nr-Min+1.EQ.1) READ (10,99992) (Nbt(N),Int(N),N=Min,Nr),
- * Mat, Mf, Mt, Ns
-99998 FORMAT (6I11, I4, I2, I3, I5)
-99994 FORMAT (4I11, 22X, I4, I2, I3, I5)
-99992 FORMAT (2I11, 44X, I4, I2, I3, I5)
- WRITE (6,66666) (Int(I),I=1,Nr)
-66666 FORMAT (' Integration type in ENDF/B file 3 =', 8I5)
-C
-C
- Np3 = Np/3
- Min = 1
- Max = 3
- IF (Max.LE.Np) THEN
- DO I=1,Np3
- READ (10,99983) (Eendf(N), Cendf(N), N=Min,Max), Mat, Mf,
- * Mt, Ns
- Min = Max + 1
- Max = Max + 3
- IF (Max.GT.Np) GO TO 40
- END DO
- 40 CONTINUE
- END IF
- IF (Min+1.EQ.Np) READ (10,99982) (Eendf(N),Cendf(N),N=Min,Np),
- * Mat, Mf, Mt, Ns
- IF (Min.EQ.Np) READ (10,99981) (Eendf(N),Cendf(N),N=Min,Np),
- * Mat, Mf, Mt, Ns
-99983 FORMAT (6E11.0, I4, I2, I3, I5)
-99982 FORMAT (4E11.0, 22X, I4, I2, I3, I5)
-99981 FORMAT (2E11.0, 44X, I4, I2, I3, I5)
-C
- CLOSE (UNIT=10)
- RETURN
- END
-C
-C
-C ------------------------------------------------------------------
-C
- SUBROUTINE Adcros (Energb, Gb, Eendf, Cendf, Nbt, INT, Nr, Np, J,
- * Ndatb)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION Eendf(*), Cendf(*), Nbt(*), Int(*)
-C
-C *** Purpose -- add cross sections from endf/b file 3 to GB
-C
- X = Energb
- DO I=1,Np
- IF (X.EQ.Eendf(I)) GO TO 15
- IF (X.LT.Eendf(I)) GO TO 20
- END DO
- WRITE (6,10000) Np, Energb, Eendf(np)
-10000 FORMAT (' Energb, Eendf(', i5, ') =', 1P2G12.4)
- STOP '[STOP -- Energb > Eendf(NP) in Adcros in rec/mrec7.f]'
-C
-C *** for Energb = Eendf(I)
- 15 CONTINUE
-C *** If discontinuous (two points with same Eendf) in endf/b file 3,then
-C *** if discontinuity occurs at beginning or end of Eendf grid use
-C *** cross section of first point, else use average cross section of
-C *** points. (This should be insignificant anyway when calculating
-C *** the Maxwellian averages)
- IF ((Eendf(I).EQ.Eendf(I+1)) .AND. (J.NE.1) .AND.
- * (J.NE.Ndatb)) THEN
- Gb = Gb + (Cendf(I)+Cendf(I+1))/2.0
- ELSE
- Gb = Gb + Cendf(I)
- END IF
- RETURN
-C
- 20 CONTINUE
-C *** for Energb < Eendf(I)
- X1 = Eendf(I-1)
- X2 = Eendf(I)
- Y1 = Cendf(I-1)
- Y2 = Cendf(I)
-C *** X1 < X < X2
-C
- DO L=1,Nr
- IF (Nbt(L).GE.I) GO TO 35
- END DO
- STOP '[STOP -- I > Nbt(NR) in Adcros in rec/mrec7.f # 2]'
-C
-C
- 35 CONTINUE
- IF (Int(L).EQ.1) THEN
-C *** y is constant in x
- Y = Y1
- Gb = Gb + Y
-C
- ELSE IF (Int(L).EQ.2) THEN
-C *** y is linear in x
- Y = Y1 + (X-X1) * (Y2-Y1)/(X2-X1)
- Gb = Gb + Y
-C
- ELSE IF (Int(L).EQ.3) THEN
-C *** y is linear in ln(x)
- Y = Y1 + (dLOG(X)-dLOG(X1)) * (Y2-Y1)/ (dLOG(X2)-dLOG(X1))
- Gb = Gb + Y
-C
- ELSE IF (Int(L).EQ.4) THEN
-C *** ln(y) is linear in x
- YL = dLOG(Y1) + (X-X1) * (dLOG(Y2)-dLOG(Y1))/(X2-X1)
- Y = dEXP(YL)
- Gb = Gb + Y
-C
- ELSE IF (Int(L).EQ.5) THEN
-C *** ln(y) is linear in ln(x)
- YL = dLOG(Y1) + (dLOG(X)-dLOG(X1)) * (dLOG(Y2)-dLOG(Y1))/
- * (dLOG(X2)-dLOG(X1))
- Y = dEXP(YL)
- Gb = Gb + Y
-C
- ELSE
- STOP '[STOP -- Int(L) > 5 in Adcros in rec/mrec7.f # 3]'
-C
- END IF
- RETURN
- END
diff --git a/sammy/src/rec/mrec8.f90 b/sammy/src/rec/mrec8.f90
index 449b8673b..f3bdbe63b 100755
--- a/sammy/src/rec/mrec8.f90
+++ b/sammy/src/rec/mrec8.f90
@@ -5,25 +5,36 @@ module mrec8_m
!
! -----------------------------------------------------------------
!
- SUBROUTINE Urrgrd (Eendf, Nbt, Int, Nr, Np, Emax, Energb, Npurr, &
- Nemax)
+ SUBROUTINE Urrgrd (Emax, Energb, Npurr,Nemax)
!
! *** Purpose -- Create energy grid for unresolved resonance region
! *** in endf/b file 3.
!
!
use AllocateFunctions_m
+ use array_sizes_common_m, only : zeroKCalc
+ use Tab1_M
IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION Eendf(*), Nbt(*), Int(*)
type(GridData)::energb
real(kind=8),allocatable,dimension(:)::X
+ type(Tab1) :: data
+ real(kind=8)::xx(1),ee(1),e1,e2,x1,x2
+ integer::inter
!
+ if (.not.zeroKCalc%readerInit) then
+ call zeroKCalc%endfReader%initialize()
+ call zeroKCalc%endfReader%setFileName(trim(zeroKCalc%endfFile))
+ zeroKCalc%readerInit = .true.
+ end if
+ call zeroKCalc%endfReader%crossAddedToStellarAverage(data)
+ np = data%getPoints()
Noemax = 0
call allocate_real_data(X, Np)
DO J=1,Np
- IF (Eendf(J).EQ.Emax) GO TO 15
- IF (Eendf(J).GT.Emax) GO TO 14
+ call data%getValue(J, ee, xx)
+ IF (ee(1).EQ.Emax) GO TO 15
+ IF (ee(1).GT.Emax) GO TO 14
END DO
STOP '[STOP -- NO Eendf .GE. Emax in Urrgrd in rec/mrec8.f]'
!
@@ -36,22 +47,23 @@ module mrec8_m
15 CONTINUE
L = 1
DO I=J,NP-1
- IF (Eendf(I).NE.Eendf(I+1)) THEN
-!
- X1 = Eendf(I)
+ call data%getValue(I, ee, xx)
+ e1 = ee(1)
+ x1 = xx(1)
+ call data%getValue(I+1, ee, xx)
+ e2 = ee(1)
+ x2 = xx(1)
+ inter = data%getInter(I+1)
+ IF (e1.NE.e2) THEN
+!
+ X1 = e1
IF (Noemax.EQ.1 .AND. I.EQ.J) X1 = Emax
- X2 = Eendf(I+1)
+ X2 = e2
Nd = 1
D = dFLOAT(Nd)
!
- DO K=1,Nr
- IF (Nbt(K).GE.I) GO TO 30
- END DO
- STOP &
- '[STOP -- I > Nbt(NR) ... error using endf energy regions]'
-!
- 30 CONTINUE
- IF ((Int(K).EQ.1) .OR. (Int(K).EQ.2)) THEN
+
+ IF ((Inter.EQ.1) .OR. (Inter.EQ.2)) THEN
! *** Int(1) => y constant in x
! *** Int(2) => y linear in x
! *** Here grid is equally spaced in energy
@@ -59,7 +71,7 @@ module mrec8_m
X(L) = X2
L = L + 1
!
- ELSE IF ((Int(K).EQ.3) .OR. (Int(K).EQ.5)) THEN
+ ELSE IF ((Inter.EQ.3) .OR. (Inter.EQ.5)) THEN
! *** Int(3) => y linear in ln(x)
! *** Int(5) => ln(y) linear in ln(x)
! *** Here grid is equally spaced in ln(energy)
@@ -67,7 +79,7 @@ module mrec8_m
X(L) = X2
L = L + 1
!
- ELSE IF (Int(K).EQ.4) THEN
+ ELSE IF (Inter.EQ.4) THEN
! *** Int(4) => ln(y) linear in x
! *** Here grid is equally spaced in exp(energy)
call reallocate_real_data(X, L, 100)
@@ -75,7 +87,7 @@ module mrec8_m
L = L + 1
!
ELSE
- WRITE (6,10000) K, Int(K)
+ WRITE (6,10000) K, Inter
10000 FORMAT(' Int(', i3, ') =', I10)
STOP &
'[STOP -- Int(K) > 5 ... endf interpolation scheme not clear]'
diff --git a/sammy/src/sammy/CMakeLists.txt b/sammy/src/sammy/CMakeLists.txt
index d4e228547..e22d82893 100644
--- a/sammy/src/sammy/CMakeLists.txt
+++ b/sammy/src/sammy/CMakeLists.txt
@@ -398,7 +398,6 @@ APPEND_SET(SAMMY_SOURCES
../rec/mrec4.f
../rec/mrec5.f
../rec/mrec6.f90
- ../rec/mrec7.f
../rec/mrec8.f90
../rpi/mrpi0.f90
diff --git a/sammy/src/the/ZeroKCrossCorrections_M.f90 b/sammy/src/the/ZeroKCrossCorrections_M.f90
index c68da6af7..035488720 100644
--- a/sammy/src/the/ZeroKCrossCorrections_M.f90
+++ b/sammy/src/the/ZeroKCrossCorrections_M.f90
@@ -348,10 +348,37 @@ contains
subroutine ZeroKCrossCorrections_addFile3(this, grid, expData)
use SammyGridAccess_M
use GridData_M
+ use Tab1_M
class(ZeroKCrossCorrections)::this
type(SammyGridAccess)::grid
type(GridDataList)::expData
+ real(kind=8)::ener
+ integer::numEl
+ integer::ipos, iel
+ type(Tab1) :: data
+ real(kind=8)::val, tmp
+
+ if (.not.this%addEndf) return
+
+ if (.not.this%readerInit) then
+ call this%endfReader%initialize()
+ call this%endfReader%setFileName(trim(this%endfFile))
+ this%readerInit = .true.
+ end if
+ call this%endfReader%crossAddedToStellarAverage(data)
+
+ numEl = grid%getNumEnergies(expData)
+ ipos = 0
+ do iel = 1, numEl
+ ener = grid%getEnergy(iel, expData)
+ if (ener.lt.0.0d0.and..not.this%wantNeg) cycle
+
+ ipos = ipos + 1
+
+ val = data%interValue(ener) + this%driver%calcData%getSharedValNs(ipos, 1, 0)
+ call this%driver%calcData%setSharedValNs(ipos, 1, 0, val)
+ end do
end subroutine
subroutine ZeroKCrossCorrections_destroy(this)
diff --git a/sammy/src/xct/mxct0.f90 b/sammy/src/xct/mxct0.f90
index a75d6bacd..ad40866cf 100644
--- a/sammy/src/xct/mxct0.f90
+++ b/sammy/src/xct/mxct0.f90
@@ -66,7 +66,7 @@ module xct_m
END IF
!
!
- Icr = Kcros
+! *** organize for broadening
Nnparx = Nnpar
Nnpar = Ndasig + Ndbsig
!
@@ -76,15 +76,12 @@ module xct_m
IF (resParData%getNumResonances().GT.0) CALL Findmx (Mxany)
!
! *** find Nrfil3, Npfil3 for added cross sections (endf/b-vi file 3)
- IF (Kaddcr.EQ.1) THEN
- CALL Read3 (Nrfil3, Npfil3)
- ELSE
- Npfil3 = 1
- Nrfil3 = 1
- END IF
+! Now done in C++
+ Npfil3 = 1
+ Nrfil3 = 1
!
! *** Count how many non-zero elements are in Xlmn
- IF (Icr.EQ.7 .OR. Icr.EQ.11 .OR. Kssmsc.GT.0) THEN
+ IF (Kcros.EQ.7 .OR. Kcros.EQ.11 .OR. Kssmsc.GT.0) THEN
ng = resParData%getNumSpinGroups()
call allocate_real_data(A_Ialj, Ntriag*ng)
CALL Kount_Clbsch (A_Izke , A_Ialj, Kslow)
@@ -103,19 +100,13 @@ module xct_m
CALL Find_If_Coulomb (A_Izeta , IfCoul, Ifdif)
CALL Estxct (Ntwo1, Nthr1, Nthr2, Nthr3, Nfour, Nfour1, Nfour2, &
Nfive1, Nfive2, Nfive3, Nfive4, Nfive1x, Nfive3x, Nsix, Neight, &
- Nnine, Nw1, Icr, Mxany, Nfprrr, K_Coul_N, numElAux)
+ Nnine, Nw1, Mxany, Nfprrr, K_Coul_N, numElAux)
!
! *** Zero ***
N = Ndasig
IF (N.ne.0) THEN
call make_I_Iisopa(N)
END IF
- IF (Npfil3.ne.0) THEN
- call make_A_Iadder(Npfil3)
- call make_A_Iaddcr(Npfil3)
- call make_I_Inbt(Nrfil3)
- call make_I_Iint(Nrfil3)
- END IF
!
IF (Ifcoul.NE.0) THEN
call make_A_Iccoul(K_Coul_N)
@@ -164,13 +155,13 @@ module xct_m
!
! *** Generate coefficients of Legendre polynomials for the various spin
! *** group pairs, if differential cross sections are needed (if
-! *** Icr=7 .OR. Icr=11 .OR. Kssmsc>0)
+! *** Kcros=7 .OR. Kcros=11 .OR. Kssmsc>0)
! *** three ***
call allocate_real_data(A_Icccll, Nthr1)
call allocate_real_data(A_Idddll, Nthr2)
call allocate_real_data(A_Ixlmn, Nthr3)
call allocate_integer_data(I_Ikxlmn, Nthr3*3)
- IF (Icr.EQ.7 .OR. Icr.EQ.11 .OR. Kssmsc.GT.0) THEN
+ IF (Kcros.EQ.7 .OR. Kcros.EQ.11 .OR. Kssmsc.GT.0) THEN
IF (Kslow.EQ.0) THEN
CALL Clbsch ( A_Ixlmn, A_Ialj , I_Ikxlmn)
ELSE
@@ -231,9 +222,6 @@ module xct_m
Nnndrc = 1
END IF
!
-! *** Andy ***
- IF (Kaddcr.EQ.1) CALL Rdndf3 (A_Iadder , A_Iaddcr , I_Inbt , &
- I_Iint , Nrfil3, Npfil3)
!
! - - - - - - - - - - - - - - - - - - - - - - <
! *** six ***
@@ -344,22 +332,14 @@ module xct_m
Lllmmm = Lllmax
IF (Lllmax.EQ.0) Lllmmm = 1
CALL Work ( calcData , calcDataSelf, &
- I_Iflmsc , A_Iprnbk , I_Iflnbk , &
- A_Iprbgf , I_Iflbgf , I_Indbgf , A_Ibgfmi , A_Ibgfma , &
- A_Itexbg , A_Iteabg , A_Ith , &
A_Isigxx , A_Idasig , A_Idbsig , A_Isigsi , A_Idasis , &
- A_Idbsis , I_Iisopa , A_Ipiece , A_Idum , A_Iadder , &
- A_Iaddcr , I_Inbt , I_Iint , A_Iedrcp , A_Icdrcp , &
+ A_Idbsis , I_Iisopa , A_Ipiece , A_Idum , &
+ A_Iedrcp , A_Icdrcp , &
A_Ixdrcp , I_Indrcp , Nnndrc , Lllmmm , Kslow)
deallocate(A_Idum)
! *** SBROUTINE Work generates theory and derivatives
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
!
-! *** if making fake data write it out, then stop
- IF (Kfake.EQ.1) THEN
- CALL Newdat (A_Ith )
- STOP '[STOP -- Normal finish for fake data]'
- END IF
!
Ks_Res = Ksolve
Ksolve = 0
@@ -369,10 +349,10 @@ module xct_m
!
Ifinal = Idimen (Ifinal, -1, ' Ifinal, -1')
Nnpar = Nnparx
- IF (Ksindi.GT.0 .AND. Kcros.EQ.8) Jwwwww = 5
Kiniso = Kkkiso
Niniso = Nnniso
IF (Iq_Val.NE.0) Niniso = Iq_Val
+
call grid%destroy()
call auxgrid%destroy()
RETURN
@@ -384,7 +364,7 @@ module xct_m
!
SUBROUTINE Estxct (Ntwo1, Nthr1, Nthr2, Nthr3, Nfour, Nfour1, &
Nfour2, Nfive1, Nfive2, Nfive3, Nfive4, Nfive1x, Nfive3x, &
- Nsix, Neight, Nnine, Nw1, Icr, Mxany, Nfprrr, K_Coul_N, &
+ Nsix, Neight, Nnine, Nw1, Mxany, Nfprrr, K_Coul_N, &
numElAux)
!
! *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-XCT
@@ -398,7 +378,6 @@ module xct_m
use rsl7_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
!
- Kcros = Icr
! *** Andy
K = 2*(Npfil3+Nrfil3)
!
@@ -413,8 +392,8 @@ module xct_m
K2 = Ntwo1
!
! *** three
- Iq_Iso = Nnniso
- IF ( (Icr.EQ.7 .OR. Icr.EQ.11) .OR. Kssmsc.GT.0) THEN
+ Iq_Iso = Nnniso
+ IF ( (Kcros.EQ.7 .OR. Kcros.EQ.11) .OR. Kssmsc.GT.0) THEN
IF (Iq_Val.GT.Iq_Iso) Iq_Iso = Iq_Val
Nthr1 = Lllmax*Iq_Iso
Nthr2 = Nthr1*Ndasig
@@ -452,7 +431,7 @@ module xct_m
Nfive2 = Nfive1*(Ndasig+Ndbsig)
IF (Nfive2.EQ.0) Nfive2 = 1
IF (Ncrssx.EQ.0) Nfive2 = 1
- IF (Icr.NE.7 .AND. Icr.NE.11 .AND. Kssmsc.EQ.0) THEN
+ IF (Kcros.NE.7 .AND. Kcros.NE.11 .AND. Kssmsc.EQ.0) THEN
Nfive3 = 1
IfCoul = 0
ELSE
diff --git a/sammy/src/xct/mxct02.f90 b/sammy/src/xct/mxct02.f90
index c22f55d6a..479c94664 100644
--- a/sammy/src/xct/mxct02.f90
+++ b/sammy/src/xct/mxct02.f90
@@ -5,10 +5,8 @@ module xct2_m
! --------------------------------------------------------------
!
SUBROUTINE Work (derivs, derivsSelf, &
- Ifl_Msc, Parnbk, Ifl_Nbk, &
- Parbgf, Ifl_Bgf, Kndbgf, Bgfmin, Bgfmax, Texbgf, Teabgf, &
- Theory, Sigxxx, Dasigx, Dbsigx, Sigsin, Dasigs, &
- Dbsigs, Isopar, Pieces, Dum, Adderg, Addcro, Nbt, Int, &
+ Sigxxx, Dasigx, Dbsigx, Sigsin, Dasigs, &
+ Dbsigs, Isopar, Pieces, Dum, &
Edrcpt, Cdrcpt, Xdrcpt, Ndrcpt, Nnndrc, Lllmmm, Kslow)
!
! *** PURPOSE -- Generate theoretical cross sections "theory" and partial
@@ -31,21 +29,20 @@ module xct2_m
use DerivativeHandler_M
use xct3_m
use xct4_m
+ use convert_to_transmission_m
! use xct6_m
use mxct27_m
IMPLICIT none
- real(8), intent(in):: Parnbk, &
- Parbgf, Bgfmin, Bgfmax, Texbgf, Teabgf, &
- Dasigx, Dbsigx, Sigsin, Dasigs, &
- Dbsigs, Pieces, Dum, Adderg, Addcro
- real(8), intent(out):: Sigxxx, Theory, Edrcpt, Cdrcpt, Xdrcpt
+ real(8), intent(in):: &
+ Dasigx, Dbsigx, Sigsin, Dasigs, &
+ Dbsigs, Pieces, Dum
+ real(8), intent(out):: Sigxxx, Edrcpt, Cdrcpt, Xdrcpt
- integer(4), intent(in):: Ifl_Msc,Ifl_Nbk,Ifl_Bgf, Kndbgf,Isopar, &
- Nbt, Int, Nnndrc, Lllmmm, Kslow
+ integer(4), intent(in):: Isopar,Nnndrc, Lllmmm, Kslow
integer(4), intent(out):: Ndrcpt
real(8):: Zero, A, Gbx, Theoryx
- integer(4):: Jdat, Idrcp, Ipoten, Iskip, Iw, Jjdat, Kkkkkk, Kkkmin, &
+ integer(4):: Jdat, Idrcp, Ipoten, Iw, irow, istart, &
Kount_Helmut, ng, numEl, TotalNdasig
integer(4) :: Iipar, iso, Jcount, Jsig,isoReal
real(8) :: val
@@ -54,22 +51,18 @@ module xct2_m
type(SammyGridAccess)::grid
type(DerivativeHandler)::derivs, derivsSelf
integer::iflagMatch
-
- DIMENSION Ifl_Msc(*), Parnbk(*), Ifl_Nbk(*), Parbgf(*), &
- Ifl_bgf(*), Kndbgf(*), Bgfmin(*), Bgfmax(*), Texbgf(*), &
- Teabgf(*), Theory(Nnnsig,*), Pieces(Ngroup,*), &
- Dum(*), Adderg(*), Addcro(*), Nbt(*), Int(*), &
+ logical::wantNegative, wantTrans, wantDeriv
+
+ DIMENSION Pieces(Ngroup,*), Dum(*), &
Edrcpt(Nnndrc,*), Cdrcpt(Nnndrc,*), Xdrcpt(*), Ndrcpt(*)
DIMENSION &
Sigxxx(Nnnsig,*), Dasigx(Nnnsig,*), Dbsigx(Nnnsig,Ndbxxx,*), &
Sigsin(*), Dasigs(*), Dbsigs(Ndbxxx,*), Isopar(*)
!
-! DIMENSION Ifl_msc(Nummsc), Parnbk(Numnbk), Ifl_nbk(Numnbk),
-! * Parbgf(Numbgf), Ifl_Bgf(), Kndbgf(), Bgfmin(), Bgfmax()
-! * Theory(Nnnsig,Ndat), Pieces(Ngroup,Ndatb),
-! * Dum(Ndatb), Adderg(Npfil3), Addcro(Npfil3), Nbt(Nrfil3),
-! * Int(Nrfil3), W...(...Ndatb)
+! DIMENSION
+! * Pieces(Ngroup,Ndatb),
+! * Dum(Ndatb), W...(...Ndatb)
!
!
DATA Zero /0.0d0/
@@ -80,7 +73,30 @@ module xct2_m
call grid%setToAuxGrid(expData)
numEl = grid%getNumEnergies(expData)
- Ywhich = Ydoppr.OR.Yresol.OR.Yangle.OR.Yssmsc.OR.Yaverg
+ ! if using Leal-Hwang doppler broadening (Kkkdop.eq.1)
+ ! don't need to extend to negative energies
+ wantNegative = .true.
+ if (Kkkdop.eq.1) wantNegative = .false.
+
+ ! do we need to convert to transmission
+ ! will only be done if no further doppler broadening is needed
+ wantTrans = .false.
+ IF (Ytrans .AND. Kcros.LT.7 .AND. Ktruet.EQ.0) wantTrans = .true.
+ if (Ydoppr) then
+ ! don't do transmission conversion
+ ! except if Maxwellian averages or energy-averages with no
+ ! correction terms
+ if (.not.(Maxwel.EQ.1 .OR. Knocor.EQ.1)) then
+ wantTrans = .false.
+ end if
+ end if
+
+
+ Ywhich = Ydoppr.OR.Yresol.OR.Yangle.OR.Yssmsc.OR.Yaverg.OR. Maxwel.EQ.1 .OR. Knocor.EQ.1
+
+ wantDeriv = Ksolve.NE.2
+ if ( covData%getPupedParam().gt.0) wantDeriv = .true.
+
!
Kount_Helmut = 0
Iw = 0
@@ -95,11 +111,12 @@ module xct2_m
iflagMatch = radFitFlags%matchFitFlag()
if(iflagMatch.gt.0) TotalNdasig = TotalNdasig + 1
- DO Iipar=1,TotalNdasig
+ DO Iipar=1,TotalNdasig
call derivs%addSharedColumn(Iipar, 1)
- end do
+ end do
end if
call derivs%setUpList(resParData, radFitFlags, Iq_Iso)
+ call derivs%reserve(numEl * Nnnsig, Ndasig + Ndbsig + 1)
IF (Iw.EQ.1.or.Ksitmp.gt.0) THEN
call derivsSelf%nullify()
@@ -110,8 +127,11 @@ module xct2_m
end do
end if
call derivsSelf%setUpList(resParData, radFitFlags, Iq_Iso)
+ call derivsSelf%reserve(numEl, Ndasig + Ndbsig + 1)
end if
+
+
CALL Zero_Integer (Isopar, Ndasig)
!
IF (Kadddc.NE.0) THEN
@@ -123,8 +143,8 @@ module xct2_m
! *** Organize for which derivative routines need to be called in Crosss
CALL Which_Derivs ()
!
- Kkkkkk = 0
- Kkkmin = 0
+ irow = 0
+ istart = 0
!
! *** (ordering is as variable # is assigned)
!
@@ -144,22 +164,36 @@ module xct2_m
! *** positive energies
DO Jdat=1,numEl
!
- do Jsig = 1,Nnnsig
- Theory(Jsig, Jdat) = 0.0d0
- end do
- Iskip = 0
+
IF ((Jdat/100000)*100000.EQ.Jdat) WRITE (6,10000) Jdat
10000 FORMAT (' *** on data point number', I10)
- Jjdat = Jdat
Su = grid%getEnergy(Jdat, expData)
- IF (Kkkdop.NE.1 .OR. Su.GT.Zero) THEN
+
+ if (Su.le.0) then
+ if (.not.wantNegative) then
+ if (Ywhich) then
+ irow = irow + 1
+ do iso = 1, Nnnsig
+ call derivs%setToZero((irow-1)*Nnnsig+iso, Ndasig+ndbsig+1)
+ if(Iw.eq.1.or.Ksitmp.gt.0) THEN
+ call derivsSelf%setToZero((irow-1)*Nnnsig+iso, Ndasig+ndbsig+1)
+ end if
+ end do
+ end if
+
+ cycle
+ end if
+ end if
+ irow = irow + 1
+
+
!
IF (Su.GT.Emax .AND. Kartgd.EQ.1) THEN
Sigxxx(1,1) = Zero
GO TO 20
END IF
IF (Su.LT.Zero) Su = - Su
-!
+
IF (Kadddc.NE.0) CALL Find_Drcpt (Edrcpt, Cdrcpt, Xdrcpt, &
Ndrcpt, Nnndrc, Idrcp, Su)
!
@@ -196,7 +230,7 @@ module xct2_m
! ************ Store Coul if needed
IF (IfCoul.GT.0) THEN
CALL Store_Coul (A_Iccoul , A_Idcoul , A_Icrssx , &
- A_Idervx , A_Icx , Jjdat)
+ A_Idervx , A_Icx , Jdat)
END IF
!
! ************ Set the particular cross sections needed for this run
@@ -204,122 +238,44 @@ module xct2_m
Dasigs, Dbsigs, Theoryx, Pieces(1,Jdat), Su, &
grid%getEnergy(Jdat, expData), Lllmmm, Kslow)
IF (Kfake.EQ.1) THEN
- Theory(1,Jdat) = Theoryx
- Iskip = 2
+ call derivs%addDataNs(Jdat, 1, 0, 1, Theoryx)
+ cycle
END IF
END IF
!
20 CONTINUE
-!
-!
-! *** Add endf/b file 3 to cross sections
- IF (Iskip.EQ.0 .AND. Kaddcr.EQ.1) THEN
- Gbx = Zero
- CALL Adcros (grid%getEnergy(Jdat, expData), &
- Gbx, Adderg, Addcro, &
- Nbt, Int, Nrfil3, Npfil3, Jjdat, numEl)
- Sigxxx(1,1) = Sigxxx(1,1) + Gbx
- END IF
-!
-!
-! ********* If want Maxwellian averages or energy-averages with no
-! ********* correction terms, do same as if there were broadening
-! ********* except don't convert to transmission unless needed
- IF (Maxwel.EQ.1 .OR. Knocor.EQ.1) THEN
- Iskip = 1
- END IF
-!
-!
-! ********* If there is Doppler-broadening, skip rest of this loop
- IF (Ydoppr) Iskip = 1
-!
-! ********* If this is angular distribution (differential), skip
- IF (Yangle) Iskip = 1
-!
-! ********* If this is capture (etc) with self-shielding and multiple-
-! ********* scattering corrections, skip the rest of this loop
- IF (Yssmsc) Iskip = 1
-!
-! ********* If no Doppler, and this is transmission, make conversion
-!
-! ********* If there is resolution, we're done here.
- IF (Yresol) Iskip = 1
-!
-!
- IF (Iskip.EQ.0) THEN
-!
-!
-!xC ********* If ETA is being calculated and there is no broadening,
-!xC ********* find derivative with respect to nu
-!x IF (Kcros.EQ.6) CALL Nnneta (Ifl_Msc, Sigxxx, Dbsigx,
-!x * Nnnsig, Nnniso)
-!
-! ********* If there is normalization or background, include it
-!
-! ********* write results onto theory if there is no broadening etc
- IF (Jjjdop.NE.1) THEN
- Kkkkkk = Kkkkkk + 1
- do Iipar = 1, Ndasig
- ! Make sure the isotope indices in derivs are consistent
- ! with the ones given in Isopar.
- ! This only matters if there are more than one isotope
- ! (Iq_Iso > 1) and if we don't already have the same value
- ! (derivs%getIsotopeForShared(Iipar).eq.Isopar(Iipar))
- ! The value for Isopar(Iipar) is not always set, i.e. is 0, in the
- ! subroutines called from this routine (derivative is zero for
- ! example). But the derivs object needs to have a value between 1 <= Iso <= Iq_Iso,
- ! for all parameters <= Ndasig, so we populated it to 1 at the beginning
- ! of this subroutine
- if (Iq_Iso.gt.1.and.Isopar(Iipar).gt.1.and. &
- derivs%getIsotopeForShared(Iipar).ne.Isopar(Iipar)) then
- call derivs%addSharedColumn(Iipar, Isopar(Iipar))
- if (Iw.eq.1.or.Ksitmp.gt.0) THEn
- call derivsSelf%addSharedColumn(Iipar, Isopar(Iipar))
- end if
- end if
- end do
- call derivs%addCalculatedData(Kkkkkk, Nnnsig, ndasig, &
- ndbsig, -1, Sigxxx(1:Nnnsig,1), Dasigx, Dbsigx(1:Nnnsig,1:Ndbxxx,1))
- END IF
-!
- ELSE IF (Iskip.EQ.1) THEN
-!
-! ********* copy results if there is broadng or if ang dist is wanted
- IF (Kkkmin.EQ.0) Kkkmin = Jdat
- do Iipar = 1, Ndasig
- if (Iq_Iso.gt.1.and.Isopar(Iipar).gt.1.and. &
- derivs%getIsotopeForShared(Iipar).ne.Isopar(Iipar)) then
- call derivs%addSharedColumn(Iipar, Isopar(Iipar))
- if (Iw.eq.1.or.Ksitmp.gt.0) THEn
- call derivsSelf%addSharedColumn(Iipar, Isopar(Iipar))
- end if
- end if
- end do
- do Iso=1,Iq_Iso
- call derivs%addCalculatedData(JJdat, Nnnsig, ndasig, &
- ndbsig, iso, Sigxxx(1:Nnnsig,Iso), Dasigx, Dbsigx(1:Nnnsig,1:Ndbxxx,Iso))
- end do
- IF (Iw.EQ.1) THEN
- do Iso=1,Iq_Iso
- call derivsSelf%addCalculatedData(JJdat, 1, ndasig, &
- ndbsig, Iso, Sigsin(Iso), Dasigs, Dbsigs(1:Ndbxxx,Iso))
- end do
- END IF
-!
- ELSE IF (Iskip.EQ.2) THEN
- do Jsig = 1,Nnnsig
- Jcount = (Jjdat-1)*Nnnsig + Jsig
- call derivs%setToZero(Jcount, Ndasig + Ndbsig)
+ do Iipar = 1, Ndasig
+ ! Make sure the isotope indices in derivs are consistent
+ ! with the ones given in Isopar.
+ ! This only matters if there are more than one isotope
+ ! (Iq_Iso > 1) and if we don't already have the same value
+ ! (derivs%getIsotopeForShared(Iipar).eq.Isopar(Iipar))
+ ! The value for Isopar(Iipar) is not always set, i.e. is 0, in the
+ ! subroutines called from this routine (derivative is zero for
+ ! example). But the derivs object needs to have a value between 1 <= Iso <= Iq_Iso,
+ ! for all parameters <= Ndasig, so we populated it to 1 at the beginning
+ ! of this subroutine
+ if (Iq_Iso.gt.1.and.Isopar(Iipar).gt.1.and. &
+ derivs%getIsotopeForShared(Iipar).ne.Isopar(Iipar)) then
+ call derivs%addSharedColumn(Iipar, Isopar(Iipar))
+ if (Iw.eq.1.or.Ksitmp.gt.0) THEn
+ call derivsSelf%addSharedColumn(Iipar, Isopar(Iipar))
+ end if
+ end if
+ end do
+ do Iso=1,Iq_Iso
+ call derivs%addCalculatedData(irow, Nnnsig, ndasig, &
+ ndbsig, iso, Sigxxx(1:Nnnsig,Iso), Dasigx, Dbsigx(1:Nnnsig,1:Ndbxxx,Iso))
+ end do
+
+ IF (Iw.EQ.1) THEN
+ do Iso=1,Iq_Iso
+ call derivsSelf%addCalculatedData(irow, 1, ndasig, &
+ ndbsig, Iso, Sigsin(Iso), Dasigs, Dbsigs(1:Ndbxxx,Iso))
end do
END IF
- else
- do Jsig = 1,Nnnsig
- Jcount = (Jjdat-1)*Nnnsig + Jsig
- call derivs%setToZero(Jcount, Ndasig + Ndbsig)
- end do
- END IF
!
END DO
--
GitLab