diff --git a/sammy/src/amr/mamr6.f90 b/sammy/src/amr/mamr6.f90
index d3c3075d15ab06b47226b75a997352f45f7cc376..b08cb957665cf70897a931892c2675bf8ba5b6dc 100755
--- a/sammy/src/amr/mamr6.f90
+++ b/sammy/src/amr/mamr6.f90
@@ -552,7 +552,7 @@ module amr6
use fixedi_m
use fxxxdx_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION Parpmc(*), Jflpmc(*), Delpmc(*), Dum(*), Idum(*)
+ DIMENSION Parpmc(4,*), Jflpmc(4,*), Delpmc(*), Dum(*), Idum(*)
WRITE (6,10000)
10000 FORMAT (' ### Oops... SAMAMR cannot accomodate PMC yet')
STOP '[STOP in Pmcfix in amr/mamr6.f]'
diff --git a/sammy/src/blk/Exploc_common.f90 b/sammy/src/blk/Exploc_common.f90
index 3c9fa46630d4019794341b02ad3e3dd3bb88e962..b41b922da0e0c2029fe358807d8e198786e47923 100644
--- a/sammy/src/blk/Exploc_common.f90
+++ b/sammy/src/blk/Exploc_common.f90
@@ -46,8 +46,8 @@ module exploc_common_m
real(kind=8),allocatable,dimension(:)::A_Idemsc
integer,allocatable,dimension(:)::I_Irdmsc
integer,allocatable,dimension(:)::I_Ijkmsc
- real(kind=8),allocatable,dimension(:)::A_Iprpmc
- integer,allocatable,dimension(:)::I_Iflpmc
+ real(kind=8),allocatable,dimension(:,:)::A_Iprpmc
+ integer,allocatable,dimension(:,:)::I_Iflpmc
! restart group counting
! old group 1
@@ -284,12 +284,22 @@ module exploc_common_m
subroutine make_A_Iprpmc(want)
integer::want
- call allocate_real_data(A_Iprpmc,want)
+ call reallocate_real_data_2d(A_Iprpmc, 4, 0, want, 0)
end subroutine make_A_Iprpmc
subroutine make_I_Iflpmc(want)
integer::want
- call allocate_integer_data(I_Iflpmc,want)
+
+ if (want.eq.0) return
+ if (allocated(I_Iflpmc)) then
+ if (size(I_Iflpmc,dim=2).lt.want) then
+ deallocate(I_Iflpmc)
+ end if
+ end if
+ if (.not.allocated(I_Iflpmc)) then
+ allocate(I_Iflpmc(4,want))
+ end if
+ I_Iflpmc = 0
end subroutine make_I_Iflpmc
diff --git a/sammy/src/cro/mcro0.f90 b/sammy/src/cro/mcro0.f90
index 6022d3187afbae205901e0e5dbae4fe5e892e643..0435f407c0403947014b29de3bcf58b6f0d137cb 100644
--- a/sammy/src/cro/mcro0.f90
+++ b/sammy/src/cro/mcro0.f90
@@ -3,7 +3,7 @@
!
SUBROUTINE Samcro_0
!
- use fixedi_m, only : Iq_Iso, Jcros, Kiniso, Kkkiso, Niniso, Nnniso, Nnnsig, Npfil3, Nres, Nrext, Nrfil3, Ntotc, needResDerivs
+ use fixedi_m, only : Iq_Iso, Jcros, Kiniso, Kkkiso, Niniso, Nnniso, Nnnsig, Npfil3, Nres, Nrext, Nrfil3, Ntotc
use ifwrit_m, only : Ks_Res, Ksolve
use exploc_common_m, only : A_Ibound, A_Iechan, A_Ipolar, A_Izke
use oopsch_common_m, only : Nowwww, Segmen
diff --git a/sammy/src/cro/mcro2.f90 b/sammy/src/cro/mcro2.f90
index 97d66ef35e425564be4ec9cefb4762a2ac96d8a2..9719f6c6c4cb763d8061ed324125a0876fd009aa 100644
--- a/sammy/src/cro/mcro2.f90
+++ b/sammy/src/cro/mcro2.f90
@@ -116,7 +116,7 @@ contains
!
SUBROUTINE Setr_Cro (calc, Ntot, Bound, Echan, &
Min, igr, &
- Z, Rmat, Sphr, Sphi, Phr, Phi, Zke, Krext, Lrmat)
+ Z, Rmat, Sphr, Sphi, Phr, Phi, Zke, Lrmat)
!
! *** PURPOSE -- GENERATE Rmat = 1/(S-B+IP)
! *** - Sum Beta*Beta/((DEL E)**2-(Gamgam/2)**2)
diff --git a/sammy/src/cro/mcro2a.f90 b/sammy/src/cro/mcro2a.f90
index 09b7af13e808593d8987c6161f91f6e0d68848af..dee459d3b774e66a4bce5d8e34c2ffc4e52a94ca 100644
--- a/sammy/src/cro/mcro2a.f90
+++ b/sammy/src/cro/mcro2a.f90
@@ -18,8 +18,7 @@ contains
class(CroCrossCalc)::calc
!
!
- CALL Abpart_Cro ( calc, A_Ixden , &
- A_Idifma , I_Inotu)
+ CALL Abpart_Cro ( calc, A_Ixden, A_Idifma)
!
CALL Parsh ( calc, &
A_Izke , &
@@ -30,7 +29,7 @@ contains
!
! --------------------------------------------------------------
!
- SUBROUTINE Abpart_Cro (calc, Xden, Difmax, Notu)
+ SUBROUTINE Abpart_Cro (calc, Xden, Difmax)
!
! *** PURPOSE -- GENERATE Upr AND Upi = ENERGY-DEPENDENT Pieces OF
! *** PR AND PI = PARTIAL OF R WRT U-PARAMETERS
@@ -53,7 +52,6 @@ contains
real(kind=8):: &
Xden(*), &
Difmax(*)
- integer:: Notu(*)
real(kind=8)::Zero, Two
real(kind=8)::Aa, G2, G3
integer::I, Igam, igr, Ij, Ipar, J, K, M, N2, N, Iflr, ires
@@ -99,7 +97,7 @@ contains
! *** PARTIAL DERIVATIVES
!
!
- IF (Ksolve.NE.2 .AND. needResDerivs) THEN
+ IF (calc%wantDerivs) THEN
calc%pr = 0.0d0
calc%pi = 0.0d0
@@ -194,31 +192,28 @@ contains
type(SammySpinGroupInfo)::spinInfo
type(RMatResonance)::resonance
type(SammyRExternalInfo)::rextInfo
- logical::ifcap
!
! DIMENSION
! * Zke(Ntotc,Ngroup)
!
!
!
- Krext = Nrext
- IF (Krext.EQ.0) Krext = 1
!
Nnf1 = 0
Nn2 = 0
- Kstart = 0
Jstart = Nfpres
! *** DO LOOP OVER GROUPS (IE SPIN-PARITY GROUPS) -
! *** GOES TO END OF SUBROUTINE
!
!
istart = 0
+ ipar = 0
DO N=1,resParData%getNumSpinGroups()
call resParData%getSpinGroupInfo(spinInfo, N)
iflAbund = spinInfo%getAbundanceFitFlag()
min = istart + 1
- call getParamPerSpinGroup(istart, N, Npr, resDeriv, &
- Kstart, ifcap)
+ call calc%getParamPerSpinGroup(istart, N)
+
IF (spinInfo%getIncludeInCalc()) THEN
IF (Numiso.GT.0) THEN
@@ -268,7 +263,7 @@ contains
CALL Setr_Cro (calc, Ntotnn, &
A_Ibound , A_Iechan , Min, n, &
A_Iz , A_Irmat , A_Isphr , A_Isphi , A_Iphr , A_Iphi , &
- Zke(1,N), Krext, Lrmat)
+ Zke(1,N), Lrmat)
!
IF (Lrmat.EQ.1) THEN
CALL Zeror_Cro (A_Iwr , A_Iwi , A_Ipwrr , A_Ipwri , &
@@ -292,43 +287,39 @@ contains
! *** TOTAL CROSS SECTIONS
nent = spinInfo%getNumEntryChannels()
next = spinInfo%getNumExitChannels()
- IF (Kcros.EQ.1) CALL Total (calc, Agoj, nent, Ntotnn, &
+ IF (Kcros.EQ.1) CALL Total (calc, spinInfo, &
A_Ics, A_Isi, &
A_Idphi , A_Iwr , A_Iwi , A_Ipwrr , A_Ipwri , A_Itr , &
- A_Iti , A_Iqr , A_Iqi , I_Inotu , Krext, Lrmat, &
- min , N, Zke(1,N), &
- iflAbund, &
- Iso)
+ A_Iti , A_Iqr , A_Iqi , Lrmat, &
+ N, Zke(1,N), &
+ iflAbund, ipar)
!
! *** SCATTERING (ELASTIC) CROSS SECTION
- IF (Kcros.EQ.2) CALL Elastc (calc, Agoj, Nent, Ntotnn, &
+ IF (Kcros.EQ.2) CALL Elastc (calc, spinInfo, &
A_Ics, &
A_Isi, A_Idphi , A_Iwr , A_Iwi, A_Ipwrr , A_Ipwri , &
- A_Itr , A_Iti , A_Iqr , A_Iqi , I_Inotu, Krext, Lrmat, &
- Min , N, Zke(1,N), &
- iflAbund, &
- Iso)
+ A_Itr , A_Iti , A_Iqr , A_Iqi , Lrmat, &
+ N, Zke(1,N), &
+ iflAbund, ipar)
!
! *** REACTION (FISSION, INELASTIC SCATTERING, ETC.) CROSS SECTIONS
IF (Kcros.EQ.3 .OR. Kcros.EQ.5 .OR. Kcros.EQ.6) CALL Reactn &
- ( calc, Agoj, Nent, Next, Ntotnn, &
+ ( calc, spinInfo, &
A_Iwr, A_Iwi, A_Ipwrr , &
- A_Ipwri , A_Itr , A_Iti , A_Iqr , A_Iqi , I_Inotu, &
- Krext, Lrmat, Min , N, &
- Zke(1,N), iflAbund, &
- Iso)
+ A_Ipwri , A_Itr , A_Iti , A_Iqr , A_Iqi , &
+ Lrmat, N, &
+ Zke(1,N), iflAbund, ipar)
!
! *** CAPTURE CROSS SECTION
IF (Kcros.EQ.4 .OR. Kcros.EQ.5 .OR. Kcros.EQ.6) CALL Captur &
- ( calc, Agoj, Nent, Next, Ntotnn, &
+ ( calc, spinInfo, &
A_Iwr, A_Iwi, A_Ipwrr , &
- A_Ipwri , A_Itr , A_Iti , A_Iqr , A_Iqi, I_Inotu, &
- Krext, Lrmat,Min, N, &
- Zke(1,N), iflAbund, &
- Iso)
+ A_Ipwri , A_Itr , A_Iti , A_Iqr , A_Iqi, &
+ Lrmat, N, &
+ Zke(1,N), iflAbund, ipar)
!
END IF
- Kstart = Kstart + Npr
+ ipar = ipar + calc%inumSize
END DO
!
RETURN
diff --git a/sammy/src/cro/mcro4.f90 b/sammy/src/cro/mcro4.f90
index a8421e188c19bd31d1398b919fc80907a2ea160b..9d61d79ab316ef821972f42f571e0eba4841cce8 100644
--- a/sammy/src/cro/mcro4.f90
+++ b/sammy/src/cro/mcro4.f90
@@ -90,7 +90,7 @@ contains
call needRadiusDeriv(Ksolve, igr, iflApe, iflApt)
!
- IF (resDeriv .OR. Npx.NE.0) THEN
+ IF (calc%wantDerivs) THEN
Kl = 0
DO K=1,Ntot
DO L=1,K
@@ -179,7 +179,7 @@ contains
END
subroutine needRadiusDeriv(Ksolve, igr, iflApe, iflApt)
- use EndfData_common_m, only : resparData, radFitFlags
+ use EndfData_common_m, only : resParData, radFitFlags
use SammySpinGroupInfo_M
IMPLICIT none
integer::Ksolve, igr
@@ -206,10 +206,9 @@ contains
!
! --------------------------------------------------------------
!
- SUBROUTINE Total (calc, Agoj, Nent, Ntot, &
+ SUBROUTINE Total (calc, spinInfo, &
Cs, Si, Dphi, Wr, Wi, Pwrhor, Pwrhoi, TR, TI, Qr, Qi, &
- Notu, Krext, Lrmat, Minres, igr, Zke, If_Zke, &
- iso)
+ Lrmat, igr, Zke, If_Zke, ipar)
!
! *** PURPOSE -- GENERATE CROSS SECTION and
! *** PARTIAL DERIVATIVE OF TOTAL
@@ -219,19 +218,20 @@ contains
use ifwrit_m
use varyr_common_m
use constn_common_m
- use EndfData_common_m, only : resparData, radFitFlags
+ use SammySpinGroupInfo_M
+ use EndfData_common_m, only : radFitFlags
IMPLICIT none
class(CroCrossCalc)::calc
real(kind=8)::agoj
- integer :: nent, ntot, Krext, Lrmat, Minres, igr, If_Zke, iso
+ integer :: nent, ntot, Lrmat, igr, If_Zke, ipar
!
!
+ type(SammySpinGroupInfo)::spinInfo
real(kind=8) :: &
Cs(*), Si(*), Dphi(*), &
Wr(*), Wi(*), Pwrhor(*), Pwrhoi(*), Tr(*), Ti(*), Qr(NN,*), &
Qi(NN,*), Zke(*)
- integer :: Notu(*)
real(kind=8)::val
!
@@ -246,66 +246,42 @@ contains
integer :: I, Ij, J, K, Kl
logical::iflApe, iflApt
- call needRadiusDeriv(Ksolve, igr, iflApe, iflApt)
- IF (iflApe.or.iflApt) THEN
-!
- Sum = Zero
- Sumc = Zero
- Kiso = If_Zke
- Ij = 0
- B = Two*Agoj*Pi100/Squ
- DO I=1,Nent
- Ij = Ij + I
- A = 1.0D0/Zke(I)**2
- Sumc = Sumc + A
- Sum = Sum + ( Cs(I)*Wr(Ij) + Si(I)*Wi(Ij) )*A
- Ktru = radFitFlags%getTrueFitFlag(Igr, I)
- IF (Ktru.GT.0) THEN
- val = - B* ( Cs(I)*Pwrhor(Ij) + Si(I)*Pwrhoi(Ij) )/Zke(I)
- call calc%crossData%setSharedValNs(calc%row, 1, Ktru, val)
- END IF
- Keff = radFitFlags%getEffFitFlag(Igr, I)
- IF (Keff.GT.0) THEN
- val = - Two*B* ( Cs(I)*Wi(Ij)-Si(I)*Wr(Ij) )*Dphi(I)/Zke(I)
- call calc%crossData%setSharedValNs(calc%row, 1, Keff, val)
- END IF
- END DO
- A = Two*Agoj*Pi100*(Sumc-Sum)/Su
- call calc%crossData%setSharedValNs(calc%row, 1, 0, A)
- IF (Kiso.GT.0) THEN
- val = A/VarAbn
- call calc%crossData%setSharedValNs(calc%row, 1, Kiso, val)
- END IF
-!
- ELSE
-!
- Sum = Zero
- Sumc = Zero
- Ij = 0
- DO I=1,Nent
- Ij = Ij + I
- Sum = Sum + ( Cs(I)*Wr(Ij) + Si(I)*Wi(Ij) )/Zke(I)**2
- Sumc = Sumc + One/Zke(I)**2
- END DO
- A = Two*Pi100*Agoj*(Sumc-Sum)/Su
- call calc%crossData%setSharedValNs(calc%row, 1, 0, A)
-!
- IF (Ksolve.NE.2) THEN
- Kiso = If_Zke
- IF (Kiso.GT.0) THEN
- val = A/VarAbn
- call calc%crossData%setSharedValNs(calc%row, 1, Kiso, val)
- END IF
+ Agoj = calc%getAbundance(igr)*spinInfo%getGFactor()
+ nent = spinInfo%getNumEntryChannels()
+ Ntot = spinInfo%getNumChannels()
+ Sum = Zero
+ Sumc = Zero
+ B = Two*Agoj*Pi100/Squ
+ Ij = 0
+ dO I=1,Nent
+ Ij = Ij + I
+ A = 1.0D0/Zke(I)**2
+ Sum = Sum + ( Cs(I)*Wr(Ij) + Si(I)*Wi(Ij) ) * A
+ Sumc = Sumc + A
+ Ktru = radFitFlags%getTrueFitFlag(Igr, I)
+ IF (Ktru.GT.0.and.calc%wantDerivs) THEN
+ val = - B* ( Cs(I)*Pwrhor(Ij) + Si(I)*Pwrhoi(Ij) )/Zke(I)
+ call calc%crossData%setSharedValNs(calc%row, 1, Ktru, val)
END IF
-!
-!
+ Keff = radFitFlags%getEffFitFlag(Igr, I)
+ IF (Keff.GT.0.and.calc%wantDerivs) THEN
+ val = - Two*B* ( Cs(I)*Wi(Ij)-Si(I)*Wr(Ij) )*Dphi(I)/Zke(I)
+ call calc%crossData%setSharedValNs(calc%row, 1, Keff, val)
+ END IF
+ END DO
+ A = Two*Pi100*Agoj*(Sumc-Sum)/Su
+ call calc%crossData%setSharedValNs(calc%row, 1, 0, A)
+ Kiso = If_Zke
+ IF (Kiso.GT.0.and.calc%wantDerivs) THEN
+ val = A/VarAbn
+ call calc%crossData%setSharedValNs(calc%row, 1, Kiso, val)
END IF
!
!
!
IF (Lrmat.EQ.1) RETURN
IF (Ksolve.EQ.2) RETURN
- IF (.not.resDeriv .AND. Npx.EQ.0) RETURN
+ IF (.not.calc%wantDerivs) RETURN
!
CALL Zero_Array (Tr, Ntriag)
CALL Zero_Array (Ti, Ntriag)
@@ -329,11 +305,9 @@ contains
END DO
END DO
!
- IF (resDeriv) CALL Derres_Cro (calc, Agoj, Ntot, &
- Tr, Ti, Notu, Minres, igr, Zke, Iso)
+ CALL Derres_Cro (calc, spinInfo, Igr, Tr, Ti, ipar)
!
- IF (Npx.NE.0) CALL Derext_Cro(calc, Agoj, Ntot, TR, &
- Zke, Iso)
+ CALL Derext_Cro(calc, spinInfo, Igr, TR)
!
RETURN
END
@@ -341,9 +315,9 @@ contains
!
! --------------------------------------------------------------
!
- SUBROUTINE Elastc (calc, Agoj, Nent, Ntot, &
- Cs, Si, Dphi, Wr, Wi, Pwrhor, Pwrhoi, Tr, Ti, Qr, Qi, Notu, &
- Krext, Lrmat, Minres, igr, Zke, If_Zke, Iso)
+ SUBROUTINE Elastc (calc, spinInfo, &
+ Cs, Si, Dphi, Wr, Wi, Pwrhor, Pwrhoi, Tr, Ti, Qr, Qi, &
+ Lrmat, igr, Zke, If_Zke, ipar)
!
! *** PURPOSE -- GENERATE CROSS SECTION
! *** AND PARTIAL DERIVATIVE OF SCATTERING
@@ -353,18 +327,19 @@ contains
use ifwrit_m
use varyr_common_m
use constn_common_m
- use EndfData_common_m, only : resparData, radFitFlags
+ use SammySpinGroupInfo_M
+ use EndfData_common_m, only : radFitFlags
IMPLICIT NONE
!
class(CroCrossCalc)::calc
+ type(SammySpinGroupInfo)::spinInfo
real(kind=8) :: agoj
- integer :: Nent, Ntot, Krext, Lrmat, Minres, igr, If_Zke,Iso
+ integer :: Nent, Ntot, Lrmat, igr, If_Zke, ipar
real(kind=8) :: &
Cs(*), Si(*), Dphi(*), &
Wr(*), Wi(*), Pwrhor(*), Pwrhoi(*), Tr(*), Ti(*), Qr(NN,*), &
Qi(NN,*), Zke(*)
- integer :: Notu(*)
real(kind=8) :: val
logical(C_BOOL)::accu
@@ -376,94 +351,77 @@ contains
!
real(kind=8),parameter :: Zero = 0.0d0, One = 1.0d0, Two = 2.0d0
!
- real(kind=8) :: A, B, C, Csk, Sik,Sum, Sum1, Sumc
+ real(kind=8) :: A, B, Bb, C, Csk, Sik,Sum, Sum1, Sumc
integer :: I, Ii, Ij, J, K, Kl, L
- logical::iflApe, iflApt
+ logical::iflAp
Sum = Zero
Sum1 = Zero
Sumc = Zero
- call needRadiusDeriv(Ksolve, igr, iflApe, iflApt)
- IF (iflApe.or.iflApt) THEN
-!
- C = Two*Agoj*Pi100/Squ
- Ij = 0
- Ii = 0
- DO I=1,Nent
- Ii = Ii + I
- B = One/Zke(I)**2
- Sumc = Sumc + B
- Sum = Sum + ( Cs(I)*Wr(Ii)+Si(I)*Wi(Ii) )*B
- B = C/Zke(I)
+ Agoj = calc%getAbundance(igr)*spinInfo%getGFactor()
+ nent = spinInfo%getNumEntryChannels()
+ Ntot = spinInfo%getNumChannels()
+
+
+ Ij = 0
+ C = Two*Agoj*Pi100/Squ
+ Ii = 0
+ DO I=1,Nent
+ B = One/Zke(I)**2
+ Sumc = Sumc + B
+
+ Ii = Ii + I
+ iflAp = .false.
+
+ if (calc%wantDerivs) then
Ktru = radFitFlags%getTrueFitFlag(Igr, I)
- IF (Ktru.GT.0) THEN
- val = - ( Cs(I)*Pwrhor(Ii)+Si(I)*Pwrhoi(Ii) )*B
- call calc%crossData%setSharedValNs(calc%row, 1, Ktru, val)
- END IF
Keff = radFitFlags%getEffFitFlag(Igr, I)
- IF (Keff.GT.0) THEN
- val = Two*( Cs(I)*Wi(Ii)-Si(I)*Wr(Ii) )*Dphi(I)*B
- call calc%crossData%setSharedValNs(calc%row, 1, Keff, val)
- END IF
- DO J=1,I
- Ij = Ij + 1
- A = Wr(Ij)**2 + Wi(Ij)**2
- IF (I.NE.J) A = A + A
- Sum1 = Sum1 + A*B
- IF (Ktru.GT.0) THEN
+ iflAp = Ktru.ne.0.or.Keff.ne.0
+ end if
+ if (iflAp) then
+ Bb = C/Zke(I)
+ IF (Ktru.GT.0) THEN
+ val = - ( Cs(I)*Pwrhor(Ii)+Si(I)*Pwrhoi(Ii) )*B
+ call calc%crossData%setSharedValNs(calc%row, 1, Ktru, val)
+ END IF
+ IF (Keff.GT.0) THEN
+ val = Two*( Cs(I)*Wi(Ii)-Si(I)*Wr(Ii) )*Dphi(I)*B
+ call calc%crossData%setSharedValNs(calc%row, 1, Keff, val)
+ END IF
+ else
+ Ktru = 0
+ end if
+
+ DO J=1,I
+ Ij = Ij + 1
+ A = Wr(Ij)**2 + Wi(Ij)**2
+ IF (I.NE.J) A = A + A
+ Sum1 = Sum1 + A*B
+ IF (I.EQ.J) THEN
+ Sum = Sum + ( Cs(I)*Wr(Ij) + Si(I)*Wi(Ij) )*B
+ IF (Ktru.GT.0) THEN
A = Wr(Ij)*Pwrhor(Ij) + Wi(Ij)*Pwrhoi(Ij)
IF (I.NE.J) A = A + A
- val = A*B
+ val = A*Bb
call calc%crossData%setSharedValNs(calc%row, 1, Ktru, val)
END IF
- END DO
- END DO
- A = Two*Agoj*Pi100*(Sumc-Two*Sum+Sum1)/Su
- call calc%crossData%setSharedValNs(calc%row, 1, 0, A)
- Kiso = If_Zke
- IF (Kiso.GT.0) THEN
- val = calc%crossData%getSharedValNs(calc%row, 1, Kiso)
- val = val/VarAbn
- accu = .false.
- call calc%crossData%setAccumulate(accu)
- call calc%crossData%setSharedValNs(calc%row, 1, Kiso, val)
- accu = .true.
- call calc%crossData%setAccumulate(accu)
- END IF
-!
- ELSE
-!
- Ij = 0
- DO I=1,Nent
- B = One/Zke(I)**2
- Sumc = Sumc + B
- DO J=1,I
- Ij = Ij + 1
- A = Wr(Ij)**2 + Wi(Ij)**2
- IF (I.NE.J) A = A + A
- Sum1 = Sum1 + A*B
- IF (I.EQ.J) THEN
- Sum = Sum + ( Cs(I)*Wr(Ij) + Si(I)*Wi(Ij) )*B
- END IF
- END DO
+ END IF
END DO
- B = Pi100*Agoj*(Sumc-Two*Sum+Sum1)/Su
- call calc%crossData%setSharedValNs(calc%row, 1, 0, B)
+ END DO
+ B = Pi100*Agoj*(Sumc-Two*Sum+Sum1)/Su
+ call calc%crossData%setSharedValNs(calc%row, 1, 0, B)
!
- IF (Ksolve.NE.2) THEN
- Kiso = If_Zke
- IF (Kiso.GT.0) THEN
- val = B/VarAbn
- call calc%crossData%setSharedValNs(calc%row, 1, Kiso, val)
- END IF
- END IF
+ Kiso = If_Zke
+ IF (Kiso.GT.0.and.calc%wantDerivs) THEN
+ val = B/VarAbn
+ call calc%crossData%setSharedValNs(calc%row, 1, Kiso, val)
END IF
!
!
!
IF (Lrmat.EQ.1) RETURN
IF (Ksolve.EQ.2) RETURN
- IF (.not.resDeriv .AND. Npx.EQ.0) RETURN
+ IF (.not.calc%wantDerivs) RETURN
!
CALL Zero_Array (TR, Ntriag)
CALL Zero_Array (TI, Ntriag)
@@ -492,11 +450,9 @@ contains
END DO
END DO
!
- IF (resDeriv) CALL Derres_Cro (calc, Agoj, Ntot, &
- TR, TI, Notu, Minres, igr, Zke, Iso)
+ CALL Derres_Cro (calc, spinInfo, Igr, Tr, Ti, ipar)
!
- IF (Npx.NE.0) CALL Derext_Cro (calc, Agoj, Ntot, TR, &
- Zke, Iso)
+ CALL Derext_Cro (calc, spinInfo, Igr, TR)
!
RETURN
END
@@ -504,9 +460,9 @@ contains
!
! --------------------------------------------------------------
!
- SUBROUTINE Reactn (calc, Agoj, Nent, Next, Ntot, &
- Wr, Wi, Pwrhor, Pwrhoi, Tr, Ti, Qr, Qi, Notu, Krext, Lrmat, &
- Minres, igr, Zke, If_Zke, Iso)
+ SUBROUTINE Reactn (calc, spinInfo, &
+ Wr, Wi, Pwrhor, Pwrhoi, Tr, Ti, Qr, Qi, Lrmat, &
+ igr, Zke, If_Zke, ipar)
!
! *** PURPOSE -- GENERATE CROSS SECTION
! *** AND PARTIAL DERIVATIVE OF REACTION
@@ -516,19 +472,20 @@ contains
use ifwrit_m
use varyr_common_m
use constn_common_m
- use EndfData_common_m, only : resparData, radFitFlags
+ use SammySpinGroupInfo_M
+ use EndfData_common_m, only : radFitFlags
use SammySpinGroupInfo_M
IMPLICIT none
!
class(CroCrossCalc)::calc
+ type(SammySpinGroupInfo)::spinInfo
real(kind=8) :: Agoj
- integer :: Nent, Next, Ntot, Krext, Lrmat, Minres, igr
- integer :: If_Zke, Iso
+ integer :: Nent, Next, Ntot, Lrmat, igr, ipar
+ integer :: If_Zke
real(kind=8) :: &
Wr(*), Wi(*), Pwrhor(*), &
Pwrhoi(*), Tr(*), Ti(*), Qr(NN,*), Qi(NN,*), Zke(*)
- integer :: Notu(*)
real(kind=8) :: val
!
! DIMENSION
@@ -541,9 +498,13 @@ contains
integer :: I, Ij, J, Jj, K, Kl, L, Ll
logical::iflApe, iflApt
!
- F1sum = Zero
- F2sum = Zero
+ F1sum = Zero
+ F2sum = Zero
call needRadiusDeriv(Ksolve, igr, iflApe, iflApt)
+ Agoj = calc%getAbundance(igr)*spinInfo%getGFactor()
+ nent = spinInfo%getNumEntryChannels()
+ next = spinInfo%getNumExitChannels()
+ Ntot = spinInfo%getNumChannels()
IF (iflApt) THEN
!
Sum = Zero
@@ -623,7 +584,7 @@ contains
!
IF (Lrmat.EQ.1) RETURN
IF (Ksolve.EQ.2) RETURN
- IF (.not.resDeriv .AND. Npx.EQ.0) RETURN
+ IF (.not.calc%wantDerivs) RETURN
!
CALL Zero_Array (Tr, Ntriag)
CALL Zero_Array (Ti, Ntriag)
@@ -651,11 +612,9 @@ contains
END DO
END IF
!
- IF (resDeriv) CALL Derres_Cro (calc, Agoj, Ntot, &
- Tr, Ti, Notu, Minres, igr, Zke, Iso)
+ CALL Derres_Cro (calc, spinInfo, Igr, Tr, Ti, ipar)
!
- IF (Npx.NE.0) CALL Derext_Cro (calc, Agoj, Ntot, Tr, &
- Zke, Iso)
+ CALL Derext_Cro (calc, spinInfo, igr, Tr)
!
RETURN
END
@@ -663,9 +622,9 @@ contains
!
! --------------------------------------------------------------
!
- SUBROUTINE Captur (calc, Agoj, Nent, Next, Ntot, &
- Wr, Wi, Pwrhor, Pwrhoi, Tr, Ti, Qr, Qi, Notu, Krext, &
- Lrmat, Minres, igr, Zke, If_Zke, Iso)
+ SUBROUTINE Captur (calc, spinInfo, &
+ Wr, Wi, Pwrhor, Pwrhoi, Tr, Ti, Qr, Qi, &
+ Lrmat, igr, Zke, If_Zke, ipar)
!
! *** PURPOSE -- GENERATE CROSS SECTION
! *** AND PARTIAL DERIVATIVE OF CAPTURE
@@ -675,21 +634,22 @@ contains
use ifwrit_m
use varyr_common_m
use constn_common_m
- use EndfData_common_m, only : resparData, radFitFlags
+ use SammySpinGroupInfo_M
+ use EndfData_common_m, only : radFitFlags
use SammySpinGroupInfo_M
IMPLICIT None
!
!
class(CroCrossCalc)::calc
+ type(SammySpinGroupInfo)::spinInfo
real(kind=8) :: agoj
- integer :: Nent, Next, Ntot, Krext, Lrmat, &
- igr, If_zke, iso, minres
+ integer :: Nent, Next, Ntot, Lrmat, &
+ igr, If_zke, ipar
real(kind=8):: &
Wr(*), Wi(*), &
Pwrhor(*), Pwrhoi(*), Tr(*), Ti(*), Qr(NN,*), &
Qi(NN,*), Zke(*)
- integer:: Notu(*)
real(kind=8)::val
!
! DIMENSION
@@ -702,72 +662,48 @@ contains
integer :: I, Ij, J, Jj, K, Kl, L, Ll
logical::iflApe, iflApt
!
- call needRadiusDeriv(Ksolve, igr, iflApe, iflApt)
-!
- IF (iflApt) THEN
-!
- Sum = Zero
- Sumc = Zero
- C = Two*Agoj*Pi100/Squ
- DO J=1,Nent
- Jj = (J*(J-1))/2
- B = One/Zke(J)**2
- Ktru = radFitFlags%getTrueFitFlag(Igr, J)
- A = C/Zke(J)
- DO I=1,Ntot
- IF (I.LE.J) THEN
- Ij = Jj + I
- ELSE
- Ij = (I*(I-1))/2 + J
- END IF
- Sum = Sum + ( Wr(Ij)**2 + Wi(Ij)**2 )*B
- IF (Ktru.GT.0) THEN
- val = - Two* &
- ( Wr(Ij)*Pwrhor(Ij)+Wi(Ij)*Pwrhoi(Ij) )*A
+ Agoj = calc%getAbundance(igr)*spinInfo%getGFactor()
+ nent = spinInfo%getNumEntryChannels()
+ next = spinInfo%getNumExitChannels()
+ Ntot = spinInfo%getNumChannels()
+
+ Sum = Zero
+ Sumc = Zero
+ C = Two*Agoj*Pi100/Squ
+ DO J=1,Nent
+ Jj = (J*(J-1))/2
+ B = One/Zke(J)**2
+ Sumc = Sumc + B
+ A = C/Zke(J)
+ Ktru = radFitFlags%getTrueFitFlag(Igr, J)
+ if (.not.calc%wantDerivs) Ktru = 0
+ DO I=1,Ntot
+ IF (I.LE.J) THEN
+ Ij = Jj + I
+ ELSE
+ Ij = (I*(I-1))/2 + J
+ END IF
+ Sum = Sum + ( Wr(Ij)**2 + Wi(Ij)**2 )*B
+ IF (Ktru.GT.0) THEN
+ val = - Two* &
+ ( Wr(Ij)*Pwrhor(Ij)+Wi(Ij)*Pwrhoi(Ij) )*A
! ??????????? is this right ???????????? methinks not
- call calc%crossData%setSharedValNs(calc%row, 1, Ktru, val)
- END IF
- END DO
- END DO
- B = Pi100*Agoj*(Sumc-Sum)/Su
- call calc%crossData%setSharedValNs(calc%row, 1, 0, B)
- Kiso = If_Zke
- IF (Kiso.GT.0) THEN
- val = B/VarAbn
- call calc%crossData%setSharedValNs(calc%row, 1, Kiso, val)
- END IF
-!
- ELSE
- Sum = Zero
- Sumc = Zero
- DO J=1,Nent
- Jj = (J*(J-1))/2
- B = One/Zke(J)**2
- Sumc = Sumc + B
- DO I=1,Ntot
- IF (I.LE.J) THEN
- Ij = Jj + I
- ELSE
- Ij = (I*(I-1))/2 + J
- END IF
- Sum = Sum + ( Wr(Ij)**2 + Wi(Ij)**2 )*B
- END DO
- END DO
- B = Pi100*Agoj*(Sumc-Sum)/Su
- call calc%crossData%setSharedValNs(calc%row, 1, 0, B)
- IF (Ksolve.NE.2) THEN
- Kiso = If_Zke
- IF (Kiso.GT.0) THEN
- val = B/VarAbn
- call calc%crossData%setSharedValNs(calc%row, 1, Kiso, val)
+ call calc%crossData%setSharedValNs(calc%row, 1, Ktru, val)
END IF
- END IF
+ END DO
+ END DO
+ B = Pi100*Agoj*(Sumc-Sum)/Su
+ call calc%crossData%setSharedValNs(calc%row, 1, 0, B)
+ Kiso = If_Zke
+ IF (Kiso.GT.0.and.calc%wantDerivs) THEN
+ val = B/VarAbn
+ call calc%crossData%setSharedValNs(calc%row, 1, Kiso, val)
END IF
!
!
IF (Lrmat.EQ.1) RETURN
IF (Ksolve.EQ.2) RETURN
- IF (.not.resDeriv .AND. Npx.EQ.0) RETURN
+ IF (.not.calc%wantDerivs) RETURN
!
CALL Zero_Array (Tr, Ntriag)
CALL Zero_Array (Ti, Ntriag)
@@ -814,55 +750,52 @@ contains
END DO
END IF
!
- IF (resDeriv) CALL Derres_Cro (calc, Agoj, Ntot, &
- Tr, Ti, Notu, Minres, igr, Zke, Iso)
+ CALL Derres_Cro (calc, spinInfo, Igr, Tr, Ti, ipar)
!
- IF (Npx.NE.0) CALL Derext_Cro (calc, Agoj, Ntot, Tr, &
- Zke, Iso)
+ CALL Derext_Cro (calc, spinInfo, Igr, Tr)
RETURN
END
!
!
! --------------------------------------------------------------
!
- SUBROUTINE Derres_Cro (calc, Agoj, Ntot, Tr, Ti, Notu, &
- Minres, igr, Zke, Iso)
+ SUBROUTINE Derres_Cro (calc, spinInfo, Igr, Tr, Ti, ipar)
!
- use fixedi_m
- use ifwrit_m
- use varyr_common_m
- use constn_common_m
- use EndfData_common_m
- use SammyResonanceInfo_M
+ use SammySpinGroupInfo_M
+ use constn_common_m, only : Fourpi
IMPLICIT NONE
!
class(CroCrossCalc)::calc
+ type(SammySpinGroupInfo)::spinInfo
real(kind=8) :: agoj
- integer :: Ntot, Minres, igr, Iso
+ integer :: Ntot, igr, ipar
- real(kind=8) :: Tr(*), Ti(*), Zke(*)
- integer :: Notu(*)
- type(SammyResonanceInfo)::resInfo
+ real(kind=8) :: Tr(*), Ti(*)
real(kind=8)::val
-!
-! DIMENSION
-! * Tr(NN), Ti(NN)
-!
- real(kind=8) :: Zero = 0.0d0, One = 1.0d0
-!,
- real(kind=8) :: A, S
+
+ real(kind=8),parameter :: Zero = 0.0d0, One = 1.0d0
+
+ real(kind=8) :: A, S, Su
integer :: I, Ifl, Ij, J, M, Mm
+
+ if (.not.calc%wantDerivs) return ! don't need derivatives
+ if(calc%inumSize.eq.0) return ! no resonance parameters are varied in this spin group
+
+ ntot = spinInfo%getNumChannels()
+ Agoj = calc%getAbundance(igr)*spinInfo%getGFactor()
+ Su = calc%ener
- IF (resDeriv) THEN
- DO Mm=1,Npr
- M = Kstart + Mm
- IF (Notu(M).NE.0) THEN
+
+ DO Mm=1,calc%inumSize
+ M = ipar + Mm
+ Ifl = calc%Inum(Mm,1)
+ if (.not.calc%covariance%contributes(Ifl)) cycle
+
S = Zero
Ij = 0
- Ifl = Notu(M)
- if (ifl.lt.0) ifl = -1 * ifl
+
DO I=1,Ntot
- A = 1.0D0/Zke(I)**2
+ A = 1.0D0/calc%Zke(I, Igr)**2
DO J=1,I
Ij = Ij + 1
IF (calc%Pi(Ij,M).NE.Zero) THEN
@@ -875,9 +808,9 @@ contains
END DO
val = Fourpi*Agoj*S/Su
call calc%crossData%setSharedValNs(calc%row, 1, Ifl, val)
- END IF
+
END DO
- END IF
+
RETURN
END
@@ -885,41 +818,45 @@ contains
!
! --------------------------------------------------------------
!
- SUBROUTINE Derext_Cro (calc, Agoj, Ntot, Tr, Zke, Iso)
+ SUBROUTINE Derext_Cro (calc, spinInfo, igr, Tr)
!
- use fixedi_m
- use ifwrit_m
- use varyr_common_m
- use constn_common_m
+ use constn_common_m, only : Fourpi
use SammyRExternalInfo_M
use RMatResonanceParam_M
- use EndfData_common_m, only : resParData
+ use SammyRExternalInfo_M
IMPLICIT none
!
class(CroCrossCalc)::calc
+ type(SammySpinGroupInfo)::spinInfo
real(kind=8) :: Agoj
- integer :: Ntot, Iso
+ integer :: Ntot, igr
- real(kind=8) :: Tr(*),Zke(*)
+ real(kind=8) :: Tr(*)
type(SammyRExternalInfo)::rextInfo
type(RExternalFunction)::rext
real(kind=8)::Parext(7)
- real(kind=8)::val
-! DIMENSION Tr(NN)
-!
- real(kind=8) :: Two = 2.0d0
+ real(kind=8)::val, Su
+ real(kind=8), parameter:: Two = 2.0d0
!
real(kind=8) :: B, C
- integer :: I, Ifl, Ij, J
+ integer :: I, Ifl, Ij, J, Nrext
+ if (.not.calc%wantDerivs) return ! don't need derivatives
+ if (calc%resData%getNumRext().eq.0) return ! no external R-Matrix
+
+ ntot = spinInfo%getNumChannels()
+ Agoj = calc%getAbundance(igr)*spinInfo%getGFactor()
+ Su = calc%ener
C = Fourpi*Agoj/Su
Ij = 0
DO I=1,Ntot
- B = C/Zke(I)**2
+ B = C/calc%Zke(I, igr)**2
Ij = Ij + I
- IF (resparData%hasRexInfo(Nnnn, I)) THEN
- call resparData%getRextInfoByGroup(rextInfo, Nnnn, I)
- call resparData%getRext(rext, rextInfo)
+ IF (calc%resData%hasRexInfo(igr, I)) THEN
+ call calc%resData%getRextInfoByGroup(rextInfo, igr, I)
+ call calc%resData%getRext(rext, rextInfo)
+ Nrext = rextInfo%getNrext()
+
Parext = 0.0d0
DO J = 1, rextInfo%getNrext()
Parext(J) = rext%getSammyValue(J)
diff --git a/sammy/src/inp/minp01.f b/sammy/src/inp/minp01.f
index d6d0cd32f91c94ab7bfc65d031e3362eb896053b..0f17873b193987bc9e6a2b8f7b74605d5aa9f5d3 100644
--- a/sammy/src/inp/minp01.f
+++ b/sammy/src/inp/minp01.f
@@ -407,8 +407,8 @@ C
C *** three PAR
IF (Numpmc.NE.0) THEN
M = Numpmc
- call make_A_Iprpmc(4*M)
- call make_I_Iflpmc(4*M)
+ call make_A_Iprpmc(M) ! 4 is added automatically
+ call make_I_Iflpmc(M) ! 4 is added automatically
call make_A_Idepmc(4*M)
call make_I_Isopmc(M)
END IF
diff --git a/sammy/src/ref/mrfs2.f b/sammy/src/ref/mrfs2.f
index 9e17e7d13804890989a9635a311f0e33ac095a38..c14b023a2c7e853748cf2b78d750370f93653ae8 100644
--- a/sammy/src/ref/mrfs2.f
+++ b/sammy/src/ref/mrfs2.f
@@ -177,11 +177,11 @@ C *** two
call make_I_Irdmsc(Ngroup)
C
C *** three
- M = Numpmc*4
+ M = Numpmc ! 4 is added automatically
IF (M.EQ.0) M = 1
call make_A_Iprpmc(M)
call make_I_Iflpmc(M)
- call make_A_Idepmc(M)
+ call make_A_Idepmc(M*4)
M = Numpmc
IF (M.EQ.0) M = 1
call make_I_isopmc(M)
diff --git a/sammy/src/the/ZeroKCrossCorrections_M.f90 b/sammy/src/the/ZeroKCrossCorrections_M.f90
index 89f4898a08b44df7b8404fde005edb67d63bab21..575c19ea88c4680d10918b5b509c39baecd5af12 100644
--- a/sammy/src/the/ZeroKCrossCorrections_M.f90
+++ b/sammy/src/the/ZeroKCrossCorrections_M.f90
@@ -34,6 +34,7 @@ module ZeroKCrossCorrections_M
procedure, pass(this) :: applyNorm => ZeroKCrossCorrections_applyNorm ! Normalize and apply background correction to 0K data. This function uses SAMMY global parameters
procedure, pass(this) :: convertToTrans => ZeroKCrossCorrections_convertToTrans ! Convert to transmission. This function uses SAMMY global parameters
procedure, pass(this) :: Fix_Eta => ZeroKCrossCorrections_Fix_Eta ! Calculate Eta, if desired. This function uses SAMMY global parameters
+ procedure, pass(this) :: AddParam => ZeroKCrossCorrections_AddParam ! add the paramagnetic cross section. This function uses SAMMY global parameters
end type
contains
@@ -387,6 +388,92 @@ contains
end do
end subroutine
+
+ subroutine ZeroKCrossCorrections_AddParam(this, grid, expData, covariance)
+ use SammyGridAccess_M
+ use paramagnetic_cross_m
+ use GridData_M
+ use fixedi_m, only : Numpmc
+ use ifwrit_m, only : Ksindi,Ksolve
+ use exploc_common_m, only : A_Iprpmc , I_Iflpmc , I_Isopmc
+
+ class(ZeroKCrossCorrections)::this
+ type(SammyGridAccess)::grid
+ type(GridDataList)::expData
+ type(CovarianceData)::covariance
+
+ integer::iel, numEl, ipos, iso, numIso, ns
+ real(kind=8)::ener, abund, pmc
+ integer::ourIso, ii, Ifl, i
+ logical::wantPara, wantDeriv
+ logical(C_BOOL)::accu
+
+
+ if (Numpmc.eq.0) return ! no paramagnetic cross section
+ wantPara = .false.
+ if (this%driver%calculator%reactType.eq.1) wantPara = .true.
+ if (Ksindi.gt.0.and.this%driver%calculator%reactType.eq.8) wantPara = .true.
+ if (.not.wantPara) return
+ numEl = grid%getNumEnergies(expData)
+
+ ipos = 0
+ numIso = this%driver%calculator%resData%getNumIso()
+ ns = this%driver%calcData%getNnnsig()
+ ourIso = this%driver%calcData%getNumberIsotopes()
+ accu = .true.
+ call this%driver%calcData%setAccumulate(accu)
+ call this%driver%calcDataSelf%setAccumulate(accu)
+
+ wantDeriv = this%driver%getWantDerivs()
+
+ if (.not.wantDeriv) then ! check for pup'ed
+ do ii = 1,Numpmc
+ do i = 1, 4
+ Ifl = I_Iflpmc(i, ii)
+ if (Ifl.le.0) cycle
+ if (ksolve.ne.2) then
+ wantDeriv = .true.
+ else
+ if ( covariance%isPupedParameter(Ifl)) then
+ wantDeriv = .true.
+ end if
+ end if
+ end do
+ end do
+ end if
+
+ do iel = 1, numEl
+ ener = grid%getEnergy(iel, expData)
+ if (ener.lt.0.0d0.and..not.this%wantNeg) cycle
+ ipos = ipos + 1
+
+ do iso = 1, numIso
+ abund = this%driver%calculator%resData%getAbundanceByIsotope(iso)
+ pmc = Addpmc (A_Iprpmc, I_Isopmc, Numpmc, iso, abund, ener)
+ if (pmc.eq.0.0d0) cycle
+ ii = iso
+ if (ourIso.eq.1) ii = 1 ! can already sum
+ if (this%driver%calculator%reactType.eq.1) then
+ call this%driver%calcData%addDataNs(ipos, Ns, 0, ii, pmc)
+ else if (this%driver%calculator%reactType.eq.8) then
+ call this%driver%calcDataSelf%addDataNs(ipos, 1, 0, ii, pmc)
+ end if
+
+ if (wantDeriv) then
+ if (this%driver%calculator%reactType.eq.1) then
+ call Dddpmc (A_Iprpmc, I_Iflpmc, I_Isopmc, Numpmc, iso, ii, abund, ener, this%driver%calcData, ipos, Ns)
+ else if (this%driver%calculator%reactType.eq.8) then
+ call Dddpmc (A_Iprpmc, I_Iflpmc, I_Isopmc, Numpmc, iso, ii, abund, ener, this%driver%calcData, ipos, 1)
+ end if
+ end if
+ end do
+ end do
+ accu = .false.
+ call this%driver%calcData%setAccumulate(accu)
+ call this%driver%calcDataSelf%setAccumulate(accu)
+
+ end subroutine
+
!!
!! Add File 3 data if desired
!! - grid energy grid on which to calculate the cross section
diff --git a/sammy/src/the/mthe0.f90 b/sammy/src/the/mthe0.f90
index 9c90f37adb41a0642ec242b6e342b58b7aa9d0ab..7a05f90cedd8ab7d0273f75bd777faf99c9a4d94 100644
--- a/sammy/src/the/mthe0.f90
+++ b/sammy/src/the/mthe0.f90
@@ -134,7 +134,8 @@ module mthe0_M
end if
if (Kfake.ne.1) then ! normal additional correction as needed
- call zeroKCalc%addFile3(grid, expData) ! add file 3 data if needed
+ call zeroKCalc%AddParam(grid, expData, covData) ! add paramagnetic cross section if desired
+ call zeroKCalc%addFile3(grid, expData) ! add file 3 data if needed
call zeroKCalc%convertToTrans(grid, expData, covData) ! convert to transmission if needed
call zeroKCalc%applyNorm(grid, expData, covData) ! apply normalization and background if needed
end if
diff --git a/sammy/src/xct/XctCrossCalc_M.f90 b/sammy/src/xct/XctCrossCalc_M.f90
index 60c7e4f35545a3321d7d9cca4e60a00137d627ea..c18f7aaabc884eb93d4fd1703ef94f9d1a40b6e7 100644
--- a/sammy/src/xct/XctCrossCalc_M.f90
+++ b/sammy/src/xct/XctCrossCalc_M.f90
@@ -202,9 +202,6 @@ subroutine XctCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itze
allocate(this%Ifcros(ntot))
this%Ifcros = .true. ! assume cross sections are to be calculated for all channels
- call info%destroy()
-
- call info%initialize(pars, cov, rad)
this%ntotc = info%getMaxChannel() -1 ! for most of the internal arrays the gamma channel is not used in the internal arrays
this%ntriag = (this%ntotc*(this%ntotc+1))/2
countCombined = .true.
diff --git a/sammy/src/xct/mxct18.f90 b/sammy/src/xct/mxct18.f90
index b4c5a64b1b935cae76eeb39826f763ef498f26c6..b8f933a19c89c7ba353f09be1cd373b14fcb152c 100644
--- a/sammy/src/xct/mxct18.f90
+++ b/sammy/src/xct/mxct18.f90
@@ -47,7 +47,7 @@ contains
STOP '[STOP in Zwhich in mxct18.f A_Iprmsc is too small ]'
end if
CALL Prtclr ( I_Ixciso , &
- I_Iflmsc , A_Iprpmc , I_Iflpmc , I_Isopmc , A_Ietax , &
+ I_Iflmsc , A_Ietax , &
A_Ietaee , Theory, Sigxxx, Dasigx, Dbsigx, &
A_Icrss , A_Ideriv , A_Itermf, I_Iisopa , Su, Eb)
!
@@ -57,8 +57,7 @@ contains
! *** type of cross section
CALL Diffel ( calc, A_Isiabn , &
I_Ifexcl , &
- A_Iprdet , I_Ifldet , I_Iigrde , I_Iflmsc , A_Iprpmc , &
- I_Iflpmc , I_Isopmc , &
+ A_Iprdet , I_Ifldet , I_Iigrde , I_Iflmsc , &
Sigxxx, Dasigx, Dbsigx, Sigsin, Dasigs, Dbsigs, &
I_Iisopa , A_Icccll, A_Idddll, &
A_Icrss , A_Ideriv , A_Icrssx , A_Idervx , A_Itermf, &
diff --git a/sammy/src/xct/mxct19.f90 b/sammy/src/xct/mxct19.f90
index d5427b720b6aee49f3c1e9a30979645b11a741ba..bc43961c5f6292ec793b3715766bccaaaa3316d6 100644
--- a/sammy/src/xct/mxct19.f90
+++ b/sammy/src/xct/mxct19.f90
@@ -2,8 +2,8 @@
!
! ______________________________________________________________________
!
- SUBROUTINE Prtclr (Ixciso, Iflmsc, Parpmc, &
- Iflpmc, Isopmc, Etanux, Etaeee, Theory, Sigxxx, Dasigx, Dbsigx, &
+ SUBROUTINE Prtclr (Ixciso, Iflmsc, &
+ Etanux, Etaeee, Theory, Sigxxx, Dasigx, Dbsigx, &
Crss, Deriv, Termf, Isopar, Su, Eb)
!
! *** Purpose -- Set Sigxxx(...) = the particular cross sections needed
@@ -23,13 +23,12 @@
IMPLICIT DOUBLE PRECISION (a-h,o-z)
!
DIMENSION Ixciso(*), Iflmsc(*), &
- Parpmc(4,*), Iflpmc(4,*), Isopmc(*), Etanux(*), Etaeee(*), &
+ Etanux(*), Etaeee(*), &
Sigxxx(Nnnsig,*), Dasigx(Nnnsig,*), Dbsigx(Nnnsig,Ndbxxx,*), &
Isopar(*), Crss(Ncrsss,*), Deriv(Ncrsss,Nnpar,*), &
Termf(*)
! DIMENSION
-! * Ixciso(Numiso), Iflmsc(Nummsc), Parpmc(4,Numpmc),
-! * Iflpmc(4,Numpmc), Isopmc(Numpmc), Etanux(Mjetan),
+! * Ixciso(Numiso), Iflmsc(Nummsc), Etanux(Mjetan),
! * Etaeee(Mjetan), Sigxxx(Nnnsig,Nnniso),
! * Dasigx(Nnnsig,Ndasig), Dbsigx(Nnnsig,Ndbsig,Nnniso),
! * Isopar(Nfpall), Crss(Ncrsss,Ngroup), Deriv(Ncrsss,Nnpar,Ngroup)
@@ -151,9 +150,7 @@
!
15 CONTINUE
Sigxxx(1,Iso) = Answer
- IF (Kfake.EQ.1) Theory = Answer + Theory
- IF (Kcros.EQ.1 .AND. Numpmc.GT.0) CALL Addpmc (Parpmc, &
- Isopmc, Sigxxx, Su, Eb, 1, Iso, Nnnsig, 0)
+ IF (Kfake.EQ.1) Theory = Answer + Theory
!
!
! *** Now, set the derivatives (if needed)
@@ -324,11 +321,6 @@
END IF
END IF
END IF
-!
- Nnn = 1
- IF (Kcros.EQ.1 .AND. Nfppmc.GT.0) CALL Dddpmc ( &
- Parpmc, Iflpmc, Isopmc, Dbsigx, Su, Eb, Nnnsig, Nnn, &
- Iso, 0)
!
END IF
60 CONTINUE
diff --git a/sammy/src/xct/mxct20.f90 b/sammy/src/xct/mxct20.f90
index f6e61a75cf362922aad820cf6096bfee9435e3da..eae5434c24981ae2bb2a6cc02441f24e0e56afbb 100644
--- a/sammy/src/xct/mxct20.f90
+++ b/sammy/src/xct/mxct20.f90
@@ -6,8 +6,8 @@ contains
! ______________________________________________________________________
!
SUBROUTINE Diffel (calc, Siabnd, Jfexcl, &
- Pardet, Ifldet, Igrdet, Iflmsc, Parpmc, &
- Iflpmc, Isopmc, Sigxxx, Dasigx, Dbsigx, Sigsin, Dasigs, Dbsigs, &
+ Pardet, Ifldet, Igrdet, Iflmsc, &
+ Sigxxx, Dasigx, Dbsigx, Sigsin, Dasigs, Dbsigs, &
Isopar, Ccclll, Dddlll, Crss , Deriv , Crssx , &
Derivx, Termf , Termfx, Echan , Cmlab , Iso_Qv, Lllmmm, Su, Eb, &
Kslow)
@@ -29,7 +29,6 @@ contains
class(XctCrossCalc)::calc
DIMENSION Siabnd(*), Jfexcl(Ntotc,*), &
Pardet(*), Ifldet(*), Igrdet(*), Iflmsc(*), &
- Parpmc(*), Iflpmc(*), Isopmc(*), &
Sigxxx(Nnnsig,*), Dasigx(Nnnsig,*), Dbsigx(Nnnsig,Ndbxxx,*), &
Sigsin(*), Dasigs(*), Dbsigs(Ndbxxx,*), Isopar(*), &
Termf(*), Termfx(*), Echan(*), Cmlab(3,*), Iso_Qv(*)
@@ -72,12 +71,12 @@ contains
! ****** set the cross sections (non-angle-differential):
! ****** beginning with self-indication transmission
IF (Ksindi.GT.0 .AND. Kcros.EQ.8) CALL Setsel (calc, Siabnd, &
- Parpmc, Isopmc, Sigsin, Crss, Su, Eb)
+ Sigsin, Crss, Su, Eb)
!
!
! ****** now do other terms
CALL Setoth (calc, Pardet, Igrdet, &
- Parpmc, Isopmc, Sigxxx, Crss, Termf, Termfx, Su, Eb)
+ Sigxxx, Crss, Termf, Termfx, Su, Eb)
END IF
!
20 CONTINUE
@@ -110,12 +109,11 @@ contains
! ****** Now, the derivatives of the angle-integrated cross sections
! ****** beginning with self-indication transmission
IF (Ksindi.GT.0 .AND. Kcros.EQ.8) call Dersel (Siabnd, &
- Iflmsc, Parpmc, Iflpmc, Isopmc, &
- Dasigs, Dbsigs, Crss, Deriv, Isopar, Su, Eb)
+ Iflmsc, Dasigs, Dbsigs, Crss, Deriv, Isopar, Su, Eb)
!
! ****** Now derivatives of the other cross sections
CALL Deroth (Pardet, Ifldet, &
- Igrdet, Iflmsc, Parpmc, Iflpmc, Isopmc, Dasigx, Dbsigx, &
+ Igrdet, Iflmsc, Dasigx, Dbsigx, &
Crss, Deriv, Isopar, Termf, Termfx, Su, Eb)
!
END IF
diff --git a/sammy/src/xct/mxct21.f90 b/sammy/src/xct/mxct21.f90
index 65d78ac5fbde79616347c6aa505778d5e8e9763f..10cfa6512f8b06844a695ac22c58222545d10ab9 100644
--- a/sammy/src/xct/mxct21.f90
+++ b/sammy/src/xct/mxct21.f90
@@ -131,20 +131,19 @@ contains
!
! --------------------------------------------------------------
!
- SUBROUTINE Setsel (calc, Siabnd, Parpmc, Isopmc, &
- Sigsin, Crss, Su, Eb)
+ SUBROUTINE Setsel (calc, Siabnd, Sigsin, Crss, Su, Eb)
!
! *** Purpose -- set self-indication transmission if needed
!
use fixedi_m, only : Ncrsss, Numiso, Numpmc
use constn_common_m, only : Fourpi
use SammySpinGroupInfo_M
+ use paramagnetic_cross_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
!
class(XctCrossCalc)::calc
real(kind=8)::Siabnd(*)
- real(kind=8)::Parpmc(*), Sigsin(*)
- integer::Isopmc(*)
+ real(kind=8)::Sigsin(*)
real(kind=8)::Crss(Ncrsss,*)
real(kind=8)::Su, Eb
type(SammySpinGroupInfo)::spinInfo
@@ -166,9 +165,7 @@ contains
END IF
END IF
END DO
- Sigsin(Iso) = Sigsin(Iso) + Sitota*Fourpi/Eb
- IF (Numpmc.NE.0) call Addpmc (Parpmc, Isopmc, &
- Sigsin, Su, Eb, 1, Iso, 1, 0)
+ Sigsin(Iso) = Sigsin(Iso) + Sitota*Fourpi/Eb
END DO
RETURN
END
@@ -177,7 +174,7 @@ contains
! ______________________________________________________________________
!
SUBROUTINE Setoth ( calc, Pardet, &
- Igrdet, Parpmc, Isopmc, Sigxxx, Crss, Termf, Termfx, Su, Eb)
+ Igrdet, Sigxxx, Crss, Termf, Termfx, Su, Eb)
!
! *** Purpose -- Set "other" cross sections as needed
!
@@ -186,11 +183,12 @@ contains
use lbro_common_m, only : Yangle
use constn_common_m, only : Fourpi
use SammySpinGroupInfo_M
+ use paramagnetic_cross_m
!
class(XctCrossCalc)::calc
- real(kind=8):: Pardet(*), Parpmc(*), &
+ real(kind=8):: Pardet(*), &
Sigxxx(Nnnsig,*), Crss(Ncrsss,*), Termf(*), Termfx(*)
- integer:: Igrdet(*), Isopmc(*)
+ integer:: Igrdet(*)
type(SammySpinGroupInfo)::spinInfo
real(kind=8),parameter::Zero = 0.0d0
real(kind=8)::Su, Eb, Val
@@ -254,8 +252,6 @@ contains
! *** total cross section
Sigxxx(Nnn,Iso) = (Termn+Terma)*Fourpi/Eb
Isox = Iso
- IF (Numpmc.GT.0) CALL Addpmc (Parpmc, Isopmc, &
- Sigxxx, Su, Eb, Nnn, Isox, Nnnsig, 0)
!
! *** elastic scattering cross section
IF (Kcros.EQ.2) Sigxxx(Nnnsig,Iso) = Termn*Fourpi/Eb
diff --git a/sammy/src/xct/mxct23.f90 b/sammy/src/xct/mxct23.f90
index e26ebbe25639f1cdabdb6e4dfceb82f7e74dd515..6dfca2191183bf45693714609865e674f8e9f8d6 100644
--- a/sammy/src/xct/mxct23.f90
+++ b/sammy/src/xct/mxct23.f90
@@ -3,7 +3,7 @@
! ______________________________________________________________________
!
SUBROUTINE Dersel (Siabnd,Iflmsc, &
- Parpmc, Iflpmc, Isopmc, Dasigs, Dbsigs, Crss, Deriv, Isopar, &
+ Dasigs, Dbsigs, Crss, Deriv, Isopar, &
Su, Eb)
!
! *** purpose -- set derivatives of transmission sample total cross
@@ -16,10 +16,11 @@
use EndfData_common_m
use SammySpinGroupInfo_M
use SammyIsoInfo_M
+ use paramagnetic_cross_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
!
DIMENSION Siabnd(*), &
- Iflmsc(*), Parpmc(*), Iflpmc(*), Isopmc(*), Dasigs(*), &
+ Iflmsc(*), Dasigs(*), &
Dbsigs(Ndbxxx,*), Crss(Ncrsss,*), Deriv(Ncrsss,Nnpar,*), &
Isopar(*)
type(SammySpinGroupInfo)::spinInfo
@@ -102,14 +103,6 @@
END DO
END IF
END DO
-!
-! *** si3 now derivatives wrt paramagnetic cross section parameters
- IF (Nfppmc.NE.0) THEN
- DO Iso=1,Nnniso
- CALL Dddpmc (Parpmc, Iflpmc, Isopmc, Dbsigs, &
- Su, Eb, 1, 1, Iso, 1)
- END DO
- END IF
!
RETURN
END
@@ -118,7 +111,7 @@
! ______________________________________________________________________
!
SUBROUTINE Deroth ( Pardet, Ifldet, &
- Igrdet, Iflmsc, Parpmc, Iflpmc, Isopmc, Dasigx, Dbsigx, Crss, &
+ Igrdet, Iflmsc, Dasigx, Dbsigx, Crss, &
Deriv, Isopar, Termf, Termfx, Su, Eb)
!
! *** Purpose -- Set derivatives of "other" cross sections
@@ -130,10 +123,10 @@
use constn_common_m
use EndfData_common_m
use SammySpinGroupInfo_M
+ use paramagnetic_cross_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
!
- DIMENSION Iflmsc(*), Parpmc(*), &
- Iflpmc(*), Isopmc(*), Pardet(*), Ifldet(*), &
+ DIMENSION Iflmsc(*), Pardet(*), Ifldet(*), &
Igrdet(*), Dasigx(Nnnsig,*), Dbsigx(Nnnsig,Nnpar,*), &
Crss(Ncrsss,*), Deriv(Ncrsss,Nnpar,*), Isopar(*), Termf(*), &
Termfx(*)
@@ -237,10 +230,7 @@
END IF
ELSE
Dbsigx(Nnn,Iiparn,Iso) = (Dtermn+Dterma)*Fourpi/Eb
- END IF
- IF (Nfppmc.GT.0 .AND. Ix.EQ.0) CALL Dddpmc ( &
- Parpmc, Iflpmc, Isopmc, Dbsigx, Su, Eb, Nnnsig, &
- Nnn, Iso, 0)
+ END IF
Ix = Ix + 1
END IF
!
diff --git a/sammy/src/xct/mxct24.f90 b/sammy/src/xct/mxct24.f90
index 59d1bf058d6638ed161f0b13bcba5baffe42ea17..98fae5bbdb0fd84f50548d34bce06d292396bc12 100644
--- a/sammy/src/xct/mxct24.f90
+++ b/sammy/src/xct/mxct24.f90
@@ -1,30 +1,32 @@
+module paramagnetic_cross_m
+! test case tr055
+IMPLICIT none
+contains
!
!
! --------------------------------------------------------------
!
- SUBROUTINE Addpmc (Parpmc, Isopmc, Sigxxx, Su, Eb, Nnn, &
- Isox, Kx, Itest)
+ real(kind=8) function Addpmc (Parpmc, Isopmc, Numpmc, Isox, Ab, Eb) result(Answer)
!
! *** purpose -- Add paramgnetic cross section thereof
!
- use fixedi_m
- use ifwrit_m
- use fixedr_m
- use EndfData_common_m
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
+ !use fixedi_m, only : Nnniso, Numiso, Numpmc
+ !use EndfData_common_m, only : resParData
+
!
- DIMENSION Parpmc(4,*), Isopmc(*), Sigxxx(kx,*)
- DATA Zero /0.0d0/
+ real(kind=8)::Parpmc(:,:)
+ integer::Isopmc(:)
+ real(kind=8)::Eb
+ integer::Isox, Numpmc
+ real(kind=8),parameter:: Zero=0.0d0
+ real(kind=8)::Ab, Bep, F, X, Su
+ integer::Ipmcs
+
!
Answer = Zero
+ Su = dabs(Eb)
DO Ipmcs=1,Numpmc
- IF (Itest.EQ.1 .OR. Nnniso.NE.Numiso .OR. &
- Isox.EQ.Isopmc(Ipmcs) ) THEN
- IF (Numiso.GT.0) THEN
- Ab = resParData%getAbundanceByIsotope(Isopmc(Ipmcs))
- ELSE
- Ab = 1.0d0
- END IF
+ IF ( Isox.EQ.Isopmc(Ipmcs) ) THEN
Bep = Parpmc(2,Ipmcs)*Su**Parpmc(3,Ipmcs)
F = Ab*Parpmc(4,Ipmcs)
X = Datan(Bep)/(Bep*(1.0d0+Bep**2)**4)
@@ -33,113 +35,72 @@
END IF
END DO
!
- IF (Eb.LT.Zero) Answer = - Answer
- Sigxxx(Nnn,Isox) = Sigxxx(Nnn,Isox) + Answer
+ IF (Eb.LT.Zero) Answer = - Answer
RETURN
- END
+ END function
!
!
! ______________________________________________________________________
!
- SUBROUTINE Dddpmc (Parpmc, Iflpmc, Isopmc, gx, su, Eb, &
- na, Nab, Iso, itest)
+ SUBROUTINE Dddpmc (Parpmc, Iflpmc, Isopmc, Numpmc, Isox, Iso, Ab, Eb, derivs, ipos, Ns)
!
! *** purpose -- Add derivatives for paramagnetic cross section
!
- use fixedi_m
- use ifwrit_m
- use fixedr_m
- use EndfData_common_m
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
+ use DerivativeHandler_M
!
- DIMENSION Parpmc(4,*), Iflpmc(4,*), Isopmc(*)
-! DIMENSION Parpmc(4,Numpmc), Iflpmc(4,Numpmc), Isopmc(Numpmc)
- DIMENSION Gx(Na,Ndbxxx,*)
-! DIMENSION Dbsigx(Nnnsig,Ndbsig,Nnniso) or
-! Dbsigs( 1,Ndbsig,Nnniso)
- DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
+ class(DerivativeHandler)::derivs
+ real(kind=8)::Parpmc(:,:)
+ integer::Iflpmc(:,:), Isopmc(:)
+ integer::Isox, Iso, Ipos, ns, Numpmc
+ real(kind=8)::Ab, Eb
+ real(kind=8):: Su
+ real(kind=8),parameter:: Zero = 0.0d0, One = 1.0d0, Two = 2.0d0
+
+ real(kind=8)::Abx, Bep, Df, F, X
+ integer::Ifl, Ipmcs
!
!
+ Su = dabs(Eb)
DO Ipmcs=1,Numpmc
- IF (Nnniso.NE.Numiso .OR. Iso.EQ.Isopmc(Ipmcs)) THEN
- Bep = Parpmc(2,Ipmcs)*Su**Parpmc(3,Ipmcs)
- IF (Numiso.GT.0) THEN
- Ab = resParData%getAbundanceByIsotope(Isopmc(Ipmcs))
- ELSE
- Ab = One
- END IF
+ IF (Isox.EQ.Isopmc(Ipmcs)) THEN
+ Bep = Parpmc(2,Ipmcs)*Su**Parpmc(3,Ipmcs)
Abx = Ab*Parpmc(4,Ipmcs)
F = dATAN(Bep)/(Bep*(One+Bep**2)**4)
- IF (Iflpmc(1,Ipmcs).GT.0) THEN
+ Ifl = Iflpmc(1,Ipmcs)
+ IF (Ifl.GT.0) THEN
X = Abx*Parpmc(1,Ipmcs)*F**2*Two
IF (Eb.LT.Zero) X = -X
- Ifl = Iflpmc(1,Ipmcs)
- IF (Ifl.LE.Ndasig) THEN
- STOP '[STOP Ifl.LE.Ndasig in Dddpmc in xct/mxct24.f]'
- ELSE
- Ifl = Ifl - Ndasig
- END IF
- IF (Gx(Nab,Ifl,Iso).EQ.Zero) THEN
- Gx(Nab,Ifl,Iso) = X
- ELSE
- Gx(Nab,Ifl,Iso) = Gx(Nab,Ifl,Iso) + X
- END IF
+ call derivs%addDataNs(ipos, Ns, Ifl, Iso, X)
END IF
- IF (Iflpmc(2,Ipmcs).GT.0) THEN
+ Ifl = Iflpmc(2,Ipmcs)
+ IF (Ifl.GT.0) THEN
X = Abx*(Parpmc(1,Ipmcs))**2
Df = One - (One/Bep+9.0d0*Bep)*dATAN(Bep)
Df = Df / (One+Bep**2)**5
Df = X*Two*F*Df
Df = Df/Parpmc(2,Ipmcs)
IF (Eb.LT.Zero) Df = -Df
- Ifl = Iflpmc(2,Ipmcs)
- IF (Ifl.LE.Ndasig) THEN
- STOP '[Ifl.LE.Ndasig in Dddpmc in xct/mxct24.f # 2]'
- ELSE
- Ifl = Ifl - Ndasig
- END IF
- IF (Gx(Nab,Ifl,Iso).EQ.Zero) THEN
- Gx(Nab,Ifl,Iso) = Df
- ELSE
- Gx(Nab,Ifl,Iso) = Gx(Nab,Ifl,Iso) + Df
- END IF
+ call derivs%addDataNs(ipos, Ns, Ifl, Iso, Df)
END IF
- IF (Iflpmc(3,Ipmcs).GT.0) THEN
+ Ifl = Iflpmc(3,Ipmcs)
+ IF (Ifl.GT.0) THEN
X = Abx*(Parpmc(1,Ipmcs))**2
Df = One - (One/Bep+9.0d0*Bep)*dATAN(Bep)
Df = Df / (One+Bep**2)**5
Df = X*Two*F*Df
Df = Df*dLOG(Su)
- IF (Eb.LT.Zero) Df = -Df
- Ifl = Iflpmc(3,Ipmcs)
- IF (Ifl.LE.Ndasig) THEN
- STOP '[Ifl.LE.Ndasig in Dddpmc in xct/mxct24.f # 3]'
- ELSE
- Ifl = Ifl - Ndasig
- END IF
- IF (Gx(Nab,Ifl,Iso).EQ.Zero) THEN
- Gx(Nab,Ifl,Iso) = Df
- ELSE
- Gx(Nab,Ifl,Iso) = Gx(Nab,Ifl,Iso) + Df
- END IF
+ IF (Eb.LT.Zero) Df = -Df
+ call derivs%addDataNs(ipos, Ns, Ifl, Iso, Df)
END IF
- IF (Iflpmc(4,Ipmcs).GT.0) THEN
+ Ifl = Iflpmc(4,Ipmcs)
+ IF (Ifl.GT.0) THEN
X = Abx*(Parpmc(1,Ipmcs)*F)**2
Df = X/Parpmc(4,Ipmcs)
IF (Eb.LT.Zero) Df = -Df
- Ifl = Iflpmc(4,Ipmcs)
- IF (Ifl.LE.Ndasig) THEN
- STOP '[Ifl.LE.Ndasig in Dddpmc in xct/mxct24.f # 4]'
- ELSE
- Ifl = Ifl - Ndasig
- END IF
- IF (Gx(Nab,Ifl,Iso).EQ.Zero) THEN
- Gx(Nab,Ifl,Iso) = Df
- ELSE
- Gx(Nab,Ifl,Iso) = Gx(Nab,Ifl,Iso) + Df
- END IF
+ call derivs%addDataNs(ipos, Ns, Ifl, Iso, Df)
END IF
END IF
END DO
RETURN
- END
+ END SUBROUTINE
+end module paramagnetic_cross_m