diff --git a/sammy/src/blk/Clm_common.f90 b/sammy/src/blk/Clm_common.f90
index ddfd9d60eeac7e41fbe4a5416dbf764423c7a7b3..0efdb4fa08cee6ddfd24438c999e16a42977fd94 100644
--- a/sammy/src/blk/Clm_common.f90
+++ b/sammy/src/blk/Clm_common.f90
@@ -4,30 +4,5 @@ module clm_common_m
! TODO: Decrypt comment below
! *** used only in clm
!
- IMPLICIT NONE
- double precision,save :: Del_Phonon
- double precision,save :: Twt
- double precision,save :: C_Trans
- double precision,save :: Tbeta
- double precision,save :: Sub
- double precision,save :: Xdop
- double precision,save :: Eps
- double precision,save :: Epsc
- double precision,save :: F0
- double precision,save :: Tev
- double precision,save :: Dwpix
- double precision,save :: Tbar
- double precision,save :: Del_S_B
- double precision,save :: Dw0
- double precision,save :: Sum_Osc_Wts
- double precision,save :: Tsave
-
- integer,save :: Mode_S_Norm
- integer,save :: Nphon
- integer,save :: N_Contin
- integer,save :: N_Osc
- integer,save :: Nbeta
- integer,save :: Nbx
- integer,save :: Nbeta_Max
-
+ IMPLICIT NONE
end module clm_common_m
diff --git a/sammy/src/clm/CrystalLatticeBroadening_M.f90 b/sammy/src/clm/CrystalLatticeBroadening_M.f90
index ed0d5757cefc97bbe63f7f319f86c21b991e3299..cb406d4a84521984d123e89bc4db426718450ea6 100644
--- a/sammy/src/clm/CrystalLatticeBroadening_M.f90
+++ b/sammy/src/clm/CrystalLatticeBroadening_M.f90
@@ -9,6 +9,49 @@ use, intrinsic :: ISO_C_BINDING
implicit none
type, extends(FreeGasDopplerBroadening) :: CrystalLatticeBroadening
+double precision :: Del_Phonon = 0.0d0
+double precision :: Twt = 0.0d0
+double precision :: C_Trans = 0.0d0
+double precision :: Tbeta = 0.0d0
+double precision :: Sub = 0.0d0
+double precision :: Xdop = 0.0d0
+double precision :: Eps = 0.0d0
+double precision :: Epsc = 0.0d0
+double precision :: F0 = 0.0d0
+double precision :: Tev = 0.0d0
+double precision :: Dwpix = 0.0d0
+double precision :: Tbar = 0.0d0
+double precision :: Del_S_B = 0.0d0
+double precision :: Dw0 = 0.0d0
+double precision :: Sum_Osc_Wts = 0.0d0
+double precision :: Tsave = 0.0d0
+integer :: Mode_S_Norm = 0
+integer :: Nphon = 0
+integer :: N_Contin = 0
+integer :: N_Osc = 0
+integer :: Nbeta = 0
+integer :: Nbx = 0
+integer :: Nbeta_Max = 25999
+
+real(kind=8),dimension(:,:),allocatable::Save
+real(kind=8),dimension(:,:),allocatable::Savex
+real(kind=8),dimension(:),allocatable::Osc_Wts
+real(kind=8),dimension(:),allocatable::Osc_Eng
+real(kind=8),dimension(:),allocatable::Phonon
+real(kind=8),dimension(:),allocatable::Ppp
+real(kind=8),dimension(:),allocatable::Beta
+real(kind=8),dimension(:),allocatable::Osc_Enx
+real(kind=8),dimension(:),allocatable::Osc_Snth
+real(kind=8),dimension(:),allocatable::Osc_Coth
+real(kind=8),dimension(:),allocatable::Bex
+real(kind=8),dimension(:),allocatable::Tlast
+real(kind=8),dimension(:),allocatable::Tnow
+real(kind=8),dimension(:),allocatable::S_Expb
+real(kind=8),dimension(:),allocatable::S_Ex
+real(kind=8),dimension(:),allocatable::Bs
+real(kind=8),dimension(:),allocatable::Ss
+real(kind=8),dimension(:),allocatable::Sab
+integer,dimension(:),allocatable::Max_T
contains
procedure, pass(this) :: initialize => CrystalLatticeBroadening_initialize
procedure, pass(this) :: broaden => CrystalLatticeBroadening_broaden
@@ -29,6 +72,25 @@ subroutine CrystalLatticeBroadening_destroy(this)
implicit none
class(CrystalLatticeBroadening) :: this
call FreeGasDopplerBroadening_destroy(this)
+ if (allocated(this%Save)) deallocate(this%Save)
+ if (allocated(this%Savex)) deallocate(this%Savex)
+ if (allocated(this%Osc_Wts)) deallocate(this%Osc_Wts)
+ if (allocated(this%Osc_Eng)) deallocate(this%Osc_Eng)
+ if (allocated(this%Phonon)) deallocate(this%Phonon)
+ if (allocated(this%Ppp)) deallocate(this%Ppp)
+ if (allocated(this%Beta)) deallocate(this%Beta)
+ if (allocated(this%Osc_Enx)) deallocate(this%Osc_Enx)
+ if (allocated(this%Osc_Snth)) deallocate(this%Osc_Snth)
+ if (allocated(this%Osc_Coth)) deallocate(this%Osc_Coth)
+ if (allocated(this%Bex)) deallocate(this%Bex)
+ if (allocated(this%Tlast)) deallocate(this%TLast)
+ if (allocated(this%TNow)) deallocate(this%TNow)
+ if (allocated(this%S_Ex)) deallocate(this%S_Ex)
+ if (allocated(this%S_Expb)) deallocate(this%S_Expb)
+ if (allocated(this%Bs)) deallocate(this%Bs)
+ if (allocated(this%Ss)) deallocate(this%Ss)
+ if (allocated(this%Sab)) deallocate(this%Sab)
+ if (allocated(this%Max_T)) deallocate(this%Max_T)
end subroutine
@@ -44,6 +106,10 @@ subroutine CrystalLatticeBroadening_broaden(this)
call this%getData(data)
call data%nullify()
call data%reserve(ndatb*data%getNnnsig(), this%numPar+1)
+ if( this%hasSelf) then
+ call this%dataSelf%nullify()
+ call this%dataSelf%reserve(ndatb*this%dataSelf%getNnnsig(), this%numPar+1)
+ end if
end subroutine
end module CrystalLatticeBroadening_M
diff --git a/sammy/src/clm/dopush.f b/sammy/src/clm/dopush.f
index bc9958d63892d5037b681aab3905ecb085c31307..3defe7efe46f3eab3677180c39c729ce5534e746 100644
--- a/sammy/src/clm/dopush.f
+++ b/sammy/src/clm/dopush.f
@@ -19,21 +19,23 @@ C
C
use fixedr_m
use broad_common_m
- use clm_common_m
use constn_common_m
use Isotop_common_m
use EndfData_common_m
use xxxb
+ use mclm2_M
+ use CrystalLatticeBroadening_M
+ use DopplerAndResolutionBroadener_M
+ use GridData_M
+ use array_sizes_common_m, only : calcData, calcDataInit
+ use AllocateFunctions_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
C
&
- DOUBLE PRECISION Osc_Eng(50), Osc_Wts(50), Phonon(1000),
- & Ar(1000), Beta(25999), Bex(25999), Sigp(4), Eref(500),
+ DOUBLE PRECISION
+ & Ar(1000), Sigp(4), Eref(500),
* Ecalc(5000),
- & Sigg(5000), Siggf(5000), Sigc(5000), E(100), Sig(100),
- & Save(25999,1000), Ppp(25999), Osc_Snth(50),
- * Osc_Enx(50),Osc_Coth(50), Tlast(25999), Tnow(25999),
- * Savex(25999,1000)
+ & Sigg(5000), Siggf(5000), Sigc(5000), E(100), Sig(100)
DOUBLE PRECISION Awr, Sigl, Dsig, Tempf,
& Rn, Toch, Kn, Sig1, Sig2, Sumgf, Arat, Rho, Ser,
& Per, EE, Sumg, Em, Sigm, Err
@@ -48,6 +50,9 @@ C & Per, EE, Sumg, Aaawww, Spi, Ap, Em, Sigm, Err
INTEGER CH, Nsyso, Nsysd, K, J, JJ,
& Ne, Nsysi, I, Nsysn, Nref
CHARACTER*40 Title
+ type(GridData)::grid
+ integer :: Mode_S_Norm
+ class(CrystalLatticeBroadening),pointer::calc
C
C
C
@@ -56,11 +61,26 @@ C
99999 FORMAT (' *** SAMMY-DOPUSH 9 Dec 05 ***')
CALL Timer (1)
CALL Timer (2)
- Nbeta_Max = 25999
C = size of Beta array (& others?)
C
C *** INITIALIZATION OF CONSTANTS
CALL Setcns (1,1)
+
+ if (.not.calcDataInit) then
+ calcDataInit = .true.
+ call calcData%initialize()
+ end if
+ if (.not.workArrayInit) then
+ call workArray%initialize()
+ workArrayInit = .true.
+ end if
+ call expData%initialize()
+ call grid%initialize()
+ call expData%addGrid(grid)
+ allocate(dopplerOption%crystalLattice)
+ calc => dopplerOption%crystalLattice
+ call calc%initialize(calcData,
+ & expData, workArray)
C
C
C
@@ -93,25 +113,29 @@ C *** NML added January 2006
C
C *** READ GLOBAL USER'S INPUT
READ (Nsysi,'(A40)') Title
- READ (Nsysi,*) CH, Mode_S_Norm, Sub
- READ (Nsysi,*) Emin, Emax, Xdop, Nphon
+ READ (Nsysi,*) CH, Mode_S_Norm, calc%Sub
+ READ (Nsysi,*) Emin, Emax, calc%Xdop, calc%Nphon
READ (Nsysi,*) Temp, Awr
- READ (Nsysi,*) Eps, Epsc, Err
+ READ (Nsysi,*) calc%Eps, calc%Epsc, Err
C
C *** READ IN CONTINUOUS DISTRIBUTION
- READ (Nsysi,*) Del_Phonon, N_Contin
- READ (Nsysi,*) (Phonon(I),I=1,N_Contin)
+ READ (Nsysi,*) calc%Del_Phonon, calc%N_Contin
+ call allocate_real_data(calc%Phonon, calc%N_Contin)
+ READ (Nsysi,*) (calc%Phonon(I),I=1,calc%N_Contin)
C
C *** READ IN DISCRETE MODES
- READ (Nsysi,*) N_Osc
- IF (N_Osc.NE.0) THEN
- READ (Nsysi,*) (Osc_Eng(I),I=1,N_Osc)
- READ (Nsysi,*) (Osc_Wts(I),I=1,N_Osc)
+ READ (Nsysi,*) calc%N_Osc
+ IF (calc%N_Osc.NE.0) THEN
+ call allocate_real_data(calc%Osc_Wts, calc%N_Osc)
+ call allocate_real_data(calc%Osc_Eng, calc%N_Osc)
+ READ (Nsysi,*) (calc%Osc_Eng(I),I=1,calc%N_Osc)
+ READ (Nsysi,*) (calc%Osc_Wts(I),I=1,calc%N_Osc)
END IF
+
C
C *** READ IN TRANSLATIONAL PART
- READ (Nsysi,*) Twt, C_Trans, Tbeta
- WRITE (Nsyso,60000) Twt, C_Trans, Tbeta
+ READ (Nsysi,*) calc%Twt, calc%C_Trans, calc%Tbeta
+ WRITE (Nsyso,60000) calc%Twt,calc%C_Trans,calc%Tbeta
60000 FORMAT (//, 'Twt, C, Tbeta=', 1p6g14.6)
C
C *** PRINTING
@@ -126,7 +150,8 @@ C *** PRINTING
& '' NORMALISATION PRECISION .................. '', F10.3, ''%'',/
& '' ERROR OF CROSS-SECTION RECONSTRuCTION .... '', F10.3, ''%'',/
& '' INTEGRATION PRECISION .................... '', F10.3, ''%'')')
- & Nphon, Temp, Awr, Emin, Emax, Xdop, Eps*100, Err*100, Epsc*100
+ & calc%Nphon, Temp, Awr, Emin, Emax,
+ & calc%Xdop, calc%Eps*100, Err*100, calc%Epsc*100
WRITE (Nsyso, '(//5X, ''RESULTS OF CONVOLUTION...'')')
WRITE (Nsyso, '(/5X, ''ENERGY'', 8X, ''SIG(NUCLEAR)'', 4X,
&''SIG(CRySTAL)'', 4X, ''SIG(fgm)'', 8X, ''SIG(fgmEFF)''/)')
@@ -138,8 +163,8 @@ c nml converted from "ratio to neutron" to "amu"
c nml so Arat should use Awr not Aaawww
Arat = Aaawww/(Aaawww+Aneutr)
Rn = Aneutr/(Aaawww+Aneutr)
- Tev = Temp*Boltzm
- Dopple = dSQRT(Tev/Awr)
+ calc%Tev = Temp*Boltzm
+ Dopple = dSQRT(calc%Tev/Awr)
C
C *** READ IN NUCLEAR CROSS-SECTION PARAMETERS
READ (Nsysn,*) Spi, Ap_Endf
@@ -149,32 +174,34 @@ C *** Endf version is in different units from sammy version
C
Zke = Twomhb*Arat
Zkte = Zke * Ap
+
+ dopplerOption%crystalLattice%Mode_S_Norm = Mode_S_Norm
c
C
C
C *** EFFECTIVE TEMPERATURE CALCULATION FOR CREATION OF Beta GRID
- Del_Beta = Del_Phonon/Tev
- CALL Start_Clm (Phonon, Ppp, Del_Beta)
+ Del_Beta = calc%Del_Phonon/calc%Tev
+ CALL Start_Clm (dopplerOption%crystalLattice,
+ & Del_Beta)
C *** INPUT -- Phonon, Tbeta, Del_Beta, N_Contin
C *** OUTPUT -- Ppp, F0, Tbar
- Tevf = Tev*Tbar
- Dwpix = F0
+ Tevf = calc%Tev*calc%Tbar
+ calc%Dwpix = calc%F0
Tevf_true = Tevf
C
C *** INITIAL Beta Mesh CALCULATION
- CALL Mesh (Beta, Tevf)
-C *** Input -- Tevf, Tev, Emax, Xdop, Del_Phonon, Sub, Nbeta_Max
+ CALL Mesh (calc, Tevf)
+C *** Input -- Tevf, Tev, Emax, Xdop, calc%Del_Phonon, Sub, Nbeta_Max
C *** Output -- Beta(Nbeta), Del_S_B, Nbeta
C
C *** Contin_X initializes things, generates Save
- CALL Contin_X (Ppp, Save, Tlast, Tnow, Savex, Del_Beta)
+ CALL Contin_X (calc, Del_Beta)
C *** Input -- Ppp(.), F0, Tbeta, Del_Beta, Sub,
C *** N_Contin, Nphon, Nbeta
C *** Output -- Save(?,Nphon)
C
C *** Discre_X initializes things for Discre
- CALL Discre_X (Osc_Eng, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- & Bex, Beta)
+ CALL Discre_X (calc)
C
C
C *** CREATION OF ARRAY OF RESONANCES ENERGIES (for starters)
@@ -207,15 +234,15 @@ C *** INVERTED STACK MODEL FOR RECONSTRuCTION OF RESONANCE LINE-SHAPE
E(J+1) = Eref(I )
E(J ) = Eref(I+1)
Tevf = Tevf_True
- CALL Sigcri (Ar, Beta, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, Save, E(J+1), Sig2, Temp, Rn, Tevf, Ch)
+ CALL Sigcri (Ar,
+ * E(J+1), Sig2, Temp, Rn, Tevf, Ch)
Sig(J+1) = Sig2
- CALL Sigcri (Ar, Beta, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, Save, E(J ), Sig1, Temp, Rn, Tevf, Ch)
+ CALL Sigcri (Ar,
+ * E(J ), Sig1, Temp, Rn, Tevf, Ch)
Sig(J) = Sig1
20 Em = (E(J+1)+E(J))/2.
- CALL Sigcri (Ar, Beta, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, Save, Em, Sigm, Temp, Rn, Tevf, Ch)
+ CALL Sigcri (Ar,
+ * Em, Sigm, Temp, Rn, Tevf, Ch)
Sigl = Sig(J+1) + (Sig(J)-Sig(J+1))*(EM-E(J+1))/(E(J)-E(J+1))
Dsig = dABS(Sigm-Sigl)
Toch = Err*Sigm
@@ -228,8 +255,8 @@ C
Sig(J+1) = Sig(J)
Sig(J ) = Sigm
Em = (E(J+1)+E(J))/2.
- CALL Sigcri (Ar, Beta, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, Save, Em, Sigm, Temp, Rn, Tevf, Ch)
+ CALL Sigcri (Ar,
+ * Em, Sigm, Temp, Rn, Tevf, Ch)
Sigl = Sig(J+1) +(Sig(J)-Sig(J+1))*(Em-E(J+1))/(E(J)-E(J+1))
Dsig = dABS(Sigm-Sigl)
Toch = Err*Sigm
@@ -247,7 +274,7 @@ C
END IF
END IF
C
- Tevf = Tbeta*Tevf + 0.5d0*Tsave
+ Tevf = calc%Tbeta*Tevf + 0.5d0*calc%Tsave
Ne = K - 1
Dopplef = dSQRT(Tevf/Awr)
dopplerOption%wantBroaden = .true.
@@ -258,13 +285,13 @@ C
DO I=1,NE
C
C *** CONVOLUTION FOR FREE GAZ
- CALL Siggaz (Ecalc(I), AR, CH, Sumg, Epsc)
+ CALL Siggaz (Ecalc(I), AR, CH, Sumg, calc%Epsc)
Sigg(I) = Sumg
C
Dummmm = Dopple
Dopple = Dopplef
C *** CONVOLUTION FOR FREE GAZ AT EFFECTIVE TEMPERATURE
- CALL Siggaz (Ecalc(I), AR, CH, Sumgf, Epsc)
+ CALL Siggaz (Ecalc(I), AR, CH, Sumgf, calc%Epsc)
Siggf(I) = Sumgf
Dopple = Dummmm
C
@@ -290,7 +317,7 @@ C
& '' ADJUSTED XDOP INTERVAL ............... '', F10.3/
& '' EFFECTIVE TEMPERATURE................. '', F10.3/
& '' DEBAY-WALLER FACTOR................... '', F10.3/)')
- & K, Nbeta, Nphon, Xdop, Tempf, Dwpix
+ & K, calc%Nbeta, calc%Nphon, calc%Xdop, Tempf, Dwpix
WRITE (6,*) 'correct run!!! '
C
CALL Timer (3)
diff --git a/sammy/src/clm/dopush3.f b/sammy/src/clm/dopush3.f
index 35e9c0586a77cca72b4caaec8f7ea800f8c76cfe..961e22abe8afbad3a812676b10921c60eee8608b 100644
--- a/sammy/src/clm/dopush3.f
+++ b/sammy/src/clm/dopush3.f
@@ -2,10 +2,15 @@ C
C
C ----------------------------------------------------------------------
C
- SUBROUTINE Sigcri (Ar, Beta, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, Save, Ei, Sigma, Temp, Rn, Tevf, Ch)
- use clm_common_m
+ SUBROUTINE Sigcri (Ar,
+ * Ei, Sigma, Temp, Rn, Tevf, Ch)
use constn_common_m
+ use mclm5_m
+ use mcml6_m
+ use mclm7_m
+ use broad_common_m
+ use CrystalLatticeBroadening_M
+ use AllocateFunctions_m
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C *** Purpose -- CRYSTAL MODEL CROSS-SECTION CALCULATION
@@ -13,64 +18,63 @@ C *** Input -- Ei,
C Rn, Recul, Tev
C *** Output -- Sigma
C
- DOUBLE PRECISION Ar, Beta, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, Save
- DIMENSION Ar(*), Beta(*), Osc_Wts(*),
- & Osc_Enx(*), Osc_Snth(50), Osc_Coth(*),
- & Bex(*), Save(Nbeta_Max,*), Max_T(Nbeta_Max)
+ DOUBLE PRECISION Ar
+ DIMENSION Ar(*)
DOUBLE PRECISION Ei, Sigma, Temp, Rn
INTEGER Ch
C
- DOUBLE PRECISION Sab, Bs, Ss, Sign, S_Expb, S_Ex
- DIMENSION Sab(25999), Bs(25999), Ss(25999), Sign(4),
- * S_Expb(25999), S_Ex(25999)
+ DOUBLE PRECISION Ss, Sign
+ DIMENSION Bs(25999), Ss(25999), Sign(4)
DOUBLE PRECISION Norm1, Sumr1, Norm, Sumr, Sintc, Sel,
& Sumc, Recul
INTEGER K, Ik, Ns, J
DOUBLE PRECISION Zero, One, Five
PARAMETER (Zero=0.0d0, One=1.0d0, Five=5.0d0)
+ class(CrystalLatticeBroadening),pointer::calc
C
C
C EFFECTIVE TEMPERATURE CALCULATION FOR CREATION OF Beta GRID
C
+ calc => dopplerOption%crystalLattice
Recul = Ei*Rn
- Alpha = Recul/Tev
+ Alpha = Recul/calc%Tev
C
C SCATTERING LaW CALCULATION
C Ei is in eV; Beta is unitless
C
C
- ix=0
+ ix=0
C *** Continuous part of distribution
- CALL Cconti (Save, Sab, Beta, Alpha, Norm, Sumr, Max_T)
+ call allocate_real_data(calc%Sab, calc%Nbeta)
+ CALL Cconti (calc, Alpha, Norm, Sumr)
C
C
Tra0 = Norm
Tra1 = Sumr
- Tbarx = Tbar
+ Tbarx = calc%Tbar
C *** Translational part, if any
- IF (Twt.NE.Zero) THEN
+ IF (calc%Twt.NE.Zero) THEN
Ndmax = 2000
- CALL Trans (Beta, Sab, Alpha, Deltab, Tra0, Tra1, Ndmax)
+ CALL Trans (calc, Alpha, Deltab, Tra0, Tra1, Ndmax)
C *** UPDATE EFfECTIVE TEMPERATURE (which is used only in Discre)
- Tevf = (Tbeta*Tevf+Twt*Tev)/(Tbeta+Twt)
- Tbarx = Tevf/Tev
+ Tevf = (calc%Tbeta*Tevf+calc%Twt*calc%Tev)/
+ & (calc%Tbeta+calc%Twt)
+ Tbarx = Tevf/calc%Tev
END IF
C
C
Osc0 = Tra0
Osc1 = Tra1
C *** Discrete oscillators, if any
- IF (N_Osc.NE.0) THEN
- CALL Discre (Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth, Sab, Beta,
- * Bex, S_Expb, S_Ex, Alpha, Tbarx, Osc0, Osc1)
+ IF (calc%N_Osc.NE.0) THEN
+ CALL Discre (calc, Alpha, Tbarx, Osc0, Osc1)
END IF
C
C
Norm1 = dABS(Norm-One)
Sumr1 = dABS(Sumr-One)
- IF (Norm1.GT.Eps .OR. Sumr1.GT.Eps) THEN
- Nphon = Nphon + 5
+ IF (Norm1.GT.calc%Eps .OR. Sumr1.GT.calc%Eps) THEN
+ calc%Nphon = calc%Nphon + 5
END IF
C
C
@@ -78,15 +82,15 @@ C Dimensioning of Alpha and Beta values
C Centering of Beta in Recul=<Beta>=Ei/(Awr+1)
C
K=0
- DO IK=Nbeta,1,-1
+ DO IK=calc%Nbeta,1,-1
K = K + 1
- Bs(K) = Recul - Beta(IK)*Tev
- Ss(K) = Sab(IK)/Tev
+ Bs(K) = Recul - calc%Beta(IK)*calc%Tev
+ Ss(K) = calc%Sab(IK)/calc%Tev
END DO
- DO IK=2,Nbeta
+ DO IK=2,calc%Nbeta
K = K + 1
- Bs(K) = Recul + Beta(IK)*Tev
- Ss(K) = Sab(IK)*dEXP(-Beta(IK))/Tev
+ Bs(K) = Recul + calc%Beta(IK)*calc%Tev
+ Ss(K) = calc%Sab(IK)*dEXP(-calc%Beta(IK))/calc%Tev
END DO
C
C
@@ -103,7 +107,7 @@ C
C Convolution for crystal between Ss and Sigman
C
DO J=1,Ns-1
- CALL Qtrapc (Ei, J, Sintc, Bs, Ss, Ns, Ar, Ch, Epsc)
+ CALL Qtrapc (Ei, J, Sintc, Bs, Ss, Ns, Ar, Ch, calc%Epsc)
Sumc = Sumc + Sintc
END DO
C
@@ -113,10 +117,10 @@ C Sigman(Ei+Recul)*EXP(-W*A)
C
CALL Csrmat (Ei+Recul, Sign, Ar)
C
- Sel = dEXP(-Dwpix*Alpha)
- IF (Mode_S_Norm.EQ.0) THEN
+ Sel = dEXP(-calc%Dwpix*Alpha)
+ IF (calc%Mode_S_Norm.EQ.0) THEN
Sigma = Sumc + Sign(Ch)*Sel
- ELSE IF (Mode_S_Norm.EQ.1) THEN
+ ELSE IF (calc%Mode_S_Norm.EQ.1) THEN
SIGMA = Sumc * (One-Sel)/(Norm-Sel) + Sign(Ch)*Sel
END IF
RETURN
diff --git a/sammy/src/clm/mclm0.f90 b/sammy/src/clm/mclm0.f90
index 47807003c77807b00cfc9847bf4fe8988537d00c..125c380deefa5f434ed3aa9753b60088279d4a3c 100644
--- a/sammy/src/clm/mclm0.f90
+++ b/sammy/src/clm/mclm0.f90
@@ -1,159 +1,72 @@
-C
-C
- Subroutine Samclm_0
-C
-C *** Purpose -- Calculate Doppler broadening using Crystal-Lattice Model
-C *** Coding based on Dmitri Naberejnev's DOPUSH code,
-C *** modified greatly by NML to conform to SAMMY conventions
-C
- use fixedi_m, only : K2reso, Kkkdex, Kkkrsl,
- * Nudwhi, Numorr, Numrpi,
- * Ndatd
- use ifwrit_m, only : Kcros, Ksolve, Kvers7, Ndatb, Ntgrlq
- use exploc_common_m
- use array_sizes_common_m
- use oopsch_common_m, only : Nowwww, Segmen
- use clm_common_m, only : N_Contin, N_Osc, Nbeta_Max, Nphon
- use lbro_common_m, only : Debug, Yresol, Yssmsc
- use AllocateFunctions_m
- use rsl7_m
- use mxct27_m
+module Samclm_m
+use CrystalLatticeBroadening_M
+implicit none
+contains
+!
+!
+ Subroutine Samclm_0(broadener)
+!
+! *** Purpose -- Calculate Doppler broadening using Crystal-Lattice Model
+! *** Coding based on Dmitri Naberejnev's DOPUSH code,
+! *** modified greatly by NML to conform to SAMMY conventions
+!
+ use fixedi_m, only : Lllmax, Ndatd, numUsedPar
+ use ifwrit_m, only : Kcros, Kdebug, Ksindi, Ksitmp, Ndatb
+ use exploc_common_m, only : I_Iflmsc
+ use lbro_common_m, only : Debug
+ use array_sizes_common_m, only : calcDataSelf
+
+ use Readclm_m
+ use mclm3_m
+
+
IMPLICIT None
- real(kind=8),allocatable,dimension(:)::A_Isave
- real(kind=8),allocatable,dimension(:)::A_Ioscwt
- real(kind=8),allocatable,dimension(:)::A_Ioscex
- real(kind=8),allocatable,dimension(:)::A_Ioscsn
- real(kind=8),allocatable,dimension(:)::A_Ioscco
- real(kind=8),allocatable,dimension(:)::A_Ibex
- real(kind=8),allocatable,dimension(:)::A_Iphono
- real(kind=8),allocatable,dimension(:)::A_Ioscen
- real(kind=8),allocatable,dimension(:)::A_Ippp
- real(kind=8),allocatable,dimension(:)::A_Itlast
- real(kind=8),allocatable,dimension(:)::A_Itnowx
- real(kind=8),allocatable,dimension(:)::A_Isavex
- real(kind=8),allocatable,dimension(:)::A_Isexpb
- real(kind=8),allocatable,dimension(:)::A_Isexxx
- real(kind=8),allocatable,dimension(:)::A_Ibs
- real(kind=8),allocatable,dimension(:)::A_Iss
- integer,allocatable,dimension(:)::I_Imaxt
- real(kind=8),allocatable,dimension(:)::A_IbetaGrid
+
integer::M, N, N_Osc_Blank, Nn, Kdatb
-C
-C
+ class(CrystalLatticeBroadening)::broadener
+!
+!
WRITE (6,99999)
99999 FORMAT (' *** SAMMY-CLM 15 Nov 07 ***')
- Segmen(1) = 'C'
- Segmen(2) = 'L'
- Segmen(3) = 'M'
- Nowwww = 0
-C
+!
CALL Initix
-C
-C
- Kdatb = Ndatd
- IF (Kdatb.EQ.0) Kdatb = Ndatb
-C
-C *** Read through the CLM file to learn array dimensions
- CALL Readclm_0 (N_Osc_Blank)
-C
-C *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMCLM
- CALL Estclm (Kdatb)
-C
- Nbeta_Max = 25999
- call allocate_real_data(A_Isave, Nphon*Nbeta_Max)
- M = IABS (N_Contin)
- N = IABS (N_Osc)
- IF (M.EQ.0) M = 1
- IF (N.EQ.0) N = 1
- call allocate_real_data(A_Ioscwt, N)
- call allocate_real_data(A_IbetaGrid, Nbeta_Max)
- call allocate_real_data(A_Ioscex , N)
- call allocate_real_data(A_Ioscsn , N)
- call allocate_real_data(A_Ioscco , N)
- call allocate_real_data(A_Ibex , Nbeta_Max)
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
- call allocate_real_data(A_Iphono , M)
- call allocate_real_data(A_Ioscen , N)
- call allocate_real_data(A_Ippp , Nbeta_Max)
- call allocate_real_data(A_Itlast , Nbeta_Max)
- call allocate_real_data(A_Itnowx , Nbeta_Max)
- call allocate_real_data(A_Isavex , Nbeta_Max*Nphon)
-C *** Read the CLM file for real
- CALL Readclm (A_Iphono , A_Ioscen , A_Ioscwt , N_Osc_Blank)
-C
-C *** Initialize the crystal-lattice calculation
- CALL Initclm (A_Isave , A_Iphono , A_Ioscen , A_Ioscwt , A_Ippp ,
- * A_IbetaGrid , A_Ioscex , A_Ioscsn , A_Ioscco , A_Ibex ,
- * A_Itlast , A_Itnowx , A_Isavex )
- deallocate(A_Iphono)
- deallocate(A_Ioscen)
- deallocate(A_Ippp)
- deallocate(A_Itlast)
- deallocate(A_Itnowx)
- deallocate(A_Isavex)
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
-C
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-C *** One ***
-C
+
+ broadener%numPar = numUsedPar
+ broadener%debugOutput = .false.
+ if (Kdebug.ne.0) broadener%debugOutput = .true.
+
+ if (Kcros.EQ.8) THEN
+ call broadener%addSelfData(calcDataSelf, Ksindi.le.0.and.Ksitmp.GT.0, Lllmax+1)
+ END IF
+ call broadener%broaden()
+
+!
+! *** Read through the CLM file to learn array dimensions
+ CALL Readclm_0 (broadener, N_Osc_Blank)
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
+! *** Read the CLM file for real
+ CALL Readclm (broadener, N_Osc_Blank)
+!
+! *** Initialize the crystal-lattice calculation
+ CALL Initclm (broadener)
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
+!
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+! *** One ***
+!
N = 1
Nn = 1
-C
-C
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
- call allocate_real_data(A_Isexpb , Nbeta_Max)
- call allocate_real_data(A_Isexxx , Nbeta_Max)
- call allocate_real_data(A_Ibs , Nbeta_Max)
- call allocate_real_data(A_Iss , Nbeta_Max)
- call allocate_integer_data(I_Imaxt , Nbeta_Max) ! now allocated as real integer
-C
-C *** Dopclm performs CLM Doppler operation
- CALL Dopclm ( I_Iflmsc ,
- * A_IbetaGrid, A_Ioscwt , A_Ioscex, A_Ioscsn, A_Ioscco , A_Ibex ,
- * A_Isexpb , A_Isexxx , A_Isave , A_Ibs , A_Iss , I_Imaxt )
+!
+!
+! *** Dopclm performs CLM Doppler operation
+ CALL Dopclm (broadener, I_Iflmsc)
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-C
-C
-C
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+!
+!
+!
CALL Write_Commons_many
- deallocate(A_IbetaGrid)
- deallocate(A_Isave)
- deallocate(A_Ioscwt)
- deallocate(A_Ioscex)
- deallocate(A_Ioscsn)
- deallocate(A_Ioscco)
- deallocate(A_Ibex)
- deallocate(A_Isexpb)
- deallocate(A_Isexxx)
- deallocate(A_Ibs)
- deallocate(A_Iss)
- deallocate(I_Imaxt)
- RETURN
-C
- END
-C
-C
-C __________________________________________________________________
-C
- SUBROUTINE Estclm (Kdatb)
-C
-C *** purpose -- estimate array size for this segment
-C
- IMPLICIT None
- integer::Kdatb, Idimen
- integer::I, K
- logical Need_isotopes
- external Idimen
-C
-C
-C *** One
-cx CALL Figure_KWs_1 (Kone)
-C
-cx K = Kone + 2*Kdatb
- K = 2*Kdatb
- K = Idimen (K, 1, 'K, 1')
- I = Idimen (K, -1, 'K, -1')
- I = Idimen (0, 1, '0, 1')
RETURN
+!
END
+end module Samclm_m
diff --git a/sammy/src/clm/mclm1.f90 b/sammy/src/clm/mclm1.f90
index 70b93ac778e7b9a273192ca7ba5c56edad95ecaa..295ab48db74c75e1929eef7af28313e07ead196f 100644
--- a/sammy/src/clm/mclm1.f90
+++ b/sammy/src/clm/mclm1.f90
@@ -1,76 +1,89 @@
-C
-C
-C -----------------------------------------------------------------
-C
- SUBROUTINE Readclm_0 (N_Osc_Blank)
-C
-C *** Purpose -- Read CLM file to learn array dimensions
-C
- use namfil_common_m
- use clm_common_m
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- CHARACTER*80 Title
- CHARACTER*5 A, Blank
- CHARACTER*1 B(80)
- CHARACTER*89 B_one
- DATA Blank /' '/
-C
+module Readclm_m
+ use CrystalLatticeBroadening_M
+ use AllocateFunctions_m, only : allocate_real_data
+
+ public Readclm_0, Readclm, Initclm
+contains
+!
+!
+! -----------------------------------------------------------------
+!
+ SUBROUTINE Readclm_0 (calc, N_Osc_Blank)
+!
+! *** Purpose -- Read CLM file to learn array dimensions
+!
+ use namfil_common_m, only : Fclmxx
+ IMPLICIT None
+ class(CrystalLatticeBroadening)::calc
+ integer::N_Osc_Blank
+ CHARACTER(80)::Title
+ CHARACTER(5):: A
+ CHARACTER(1)::B(80)
+ CHARACTER(89)::B_one
+ CHARACTER(5), parameter :: Blank = ' '
+ integer::Ixx, Kpound
+!
equivalence(B,B_one)
CALL Filopn (10, Fclmxx, 0)
-C
+!
READ (10,10200,END=400,ERR=500) Title
10200 FORMAT (A80)
- READ (10,*) Mode_S_Norm, Nphon, Sub, Xdop, Eps, Epsc
+ READ (10,*) calc%Mode_S_Norm, calc%Nphon, calc%Sub, calc%Xdop, calc%Eps, calc%Epsc
READ (10,10300) A
10300 FORMAT (A5)
IF (A.NE.Blank) STOP '[Stop in Readclm_0 in clm/mclm1.f # 4]'
-C
+!
READ (10,10300) A
CALL Convert_To_Caps (A, 5, Kpound)
-C
+!
IF (A.EQ.'CONTI') THEN
READ (10,10400) B
10400 FORMAT (80A1)
- CALL Find_Con (B, B_one, N_Contin, Ixx)
-C *** Done reading Continuous Mode; Read next card
+ CALL Find_Con (B, B_one, calc%N_Contin, Ixx)
+! *** Done reading Continuous Mode; Read next card
IF (Ixx.EQ.1) GO TO 400
READ (10,10300,End=400) A
CALL Convert_To_Caps (A, 5, Kpound)
END IF
-C
+!
IF (A.EQ.'DISCR') THEN
READ (10,10400) B
- CALL Find_Osc (B, B_one, N_Osc, N_Osc_Blank, Ixx)
+ CALL Find_Osc (B, B_one, calc%N_Osc, N_Osc_Blank, Ixx)
IF (Ixx.EQ.1) GO TO 400
-C *** Done reading Discrete Mode; Read next card
+! *** Done reading Discrete Mode; Read next card
READ (10,10300,END=400) A
CALL Convert_To_Caps (A, 5, Kpound)
END IF
-C
+!
IF (A.EQ.'TRANS') THEN
READ (10,10400,END=400) B
END IF
-C
+!
400 CONTINUE
REWIND (UNIT=10)
RETURN
-C
-C
+!
+!
500 CONTINUE
CLOSE (UNIT=10)
STOP '[Stop in Readclm_0 in clm/mclm1.f problem reading CLM file]'
END
-C
-C
-C --------------------------------------------------------------
-C
+!
+!
+! --------------------------------------------------------------
+!
SUBROUTINE Find_Con (B, Bb, Nn, Ixx)
- CHARACTER*1 B(80)
- CHARACTER*80 Bb
- CHARACTER*10 A, Blank
- DIMENSION C(8)
- DATA Blank /' '/, Zero /0.0d0/
+ IMPLICIT none
+ CHARACTER(1):: B(80)
+ CHARACTER(80):: Bb
+ CHARACTER(10):: A
+ integer::Nn, Ixx
+ real(kind=8)::C(8)
+ CHARACTER(10), parameter ::Blank =' '
+ real(kind=8),parameter::Zero=0.0d0
+ real(kind=8)::Cx
+ integer::I, L, M, N
Ixx = 0
L = 0
N = 0
@@ -95,8 +108,8 @@ C
END IF
END DO
Nn = 0
-C *** User did not count the phonons; ergo
-C *** i.e., fixed-format where SAMMY counts the number of points
+! *** User did not count the phonons; ergo
+! *** i.e., fixed-format where SAMMY counts the number of points
N = 0
30 CONTINUE
READ (10,10500,END=400) C
@@ -110,17 +123,17 @@ C *** i.e., fixed-format where SAMMY counts the number of points
END DO
GO TO 30
40 CONTINUE
-C *** Note that we have read and counted the blank line
+! *** Note that we have read and counted the blank line
RETURN
-C
+!
50 CONTINUE
-C *** User did count the phonons, so now we need to count cards
+! *** User did count the phonons, so now we need to count cards
READ (10,10000) A
10000 FORMAT (A10)
IF (A.EQ.Blank) GO TO 60
GO TO 50
60 CONTINUE
-C *** So number of cards read should now be Ncard + Ncardx; includes blank
+! *** So number of cards read should now be Ncard + Ncardx; includes blank
READ (Bb,*) Cx, Nn
RETURN
@@ -128,20 +141,25 @@ C *** So number of cards read should now be Ncard + Ncardx; includes blank
Ixx = 1
RETURN
END
-C
-C
-C --------------------------------------------------------------
-C
+!
+!
+! --------------------------------------------------------------
+!
SUBROUTINE Find_Osc (B, Bb, Nn, N_Osc_Blank, Ixx)
- CHARACTER*1 B(80)
- CHARACTER*80 Bb
- DIMENSION C(8)
- CHARACTER*10 A, Blank
- DATA Blank /' '/, Zero /0.0d0/
-C
+ IMPLICIT none
+ CHARACTER(1):: B(80)
+ CHARACTER(80)::Bb
+ integer::Nn, N_Osc_Blank, Ixx
+ real(kind=8)::C(8)
+ CHARACTER(10)::A
+ CHARACTER(10), parameter :: Blank = ' '
+ real(kind=8),parameter:: Zero=0.0d0
+ integer::I, L, N
+
+!
Ixx = 0
N_Osc_Blank = 0
-C *** Is card entirely blank?
+! *** Is card entirely blank?
DO I=1,80
IF (B(I).NE.' ') THEN
L = I
@@ -153,22 +171,22 @@ C *** Is card entirely blank?
N_Osc_Blank = 1
READ (10,10500,END=400) C
READ (10,10500,END=400) C
-C *** User did not count the phonons; nothing is on the card
+! *** User did not count the phonons; nothing is on the card
GO TO 50
-C
+!
10 CONTINUE
-C *** Does card have a decimal point?
+! *** Does card have a decimal point?
DO I=L+1,80
IF (B(I).EQ.'.') THEN
-C *** User did not count the phonons
+! *** User did not count the phonons
GO TO 40
ELSE IF (B(I).EQ.' ') THEN
GO TO 20
END IF
END DO
-C
+!
20 CONTINUE
-C *** User did count the oscillators, so now we need to count cards
+! *** User did count the oscillators, so now we need to count cards
READ (10,10000) A
10000 FORMAT (A10)
IF (A.EQ.Blank) GO TO 30
@@ -176,7 +194,7 @@ C *** User did count the oscillators, so now we need to count cards
30 CONTINUE
READ (Bb,*) Nn
RETURN
-C
+!
40 CONTINUE
N = 0
Nn = 0
@@ -184,7 +202,7 @@ C
READ (10,10500) C
10500 FORMAT (8F10.1)
50 CONTINUE
-C *** Fixed-format where SAMMY counts the number of points
+! *** Fixed-format where SAMMY counts the number of points
DO I=1,8
N = N + 1
IF (C(I).EQ.Zero) THEN
@@ -197,180 +215,184 @@ C *** Fixed-format where SAMMY counts the number of points
GO TO 50
60 CONTINUE
RETURN
-C
+!
400 CONTINUE
Nn = - N
Ixx = 1
RETURN
END
-C
-C
-C -----------------------------------------------------------------
-C
- SUBROUTINE Readclm (Phonon, Osc_Eng, Osc_Wts, N_Osc_Blank)
-C
-C *** Purpose -- Read CLM file for real
-C
- use fixedr_m
- use namfil_common_m
- use clm_common_m
- use constn_common_m
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION Phonon(*), Osc_Eng(*), Osc_Wts(*)
- CHARACTER*80 Title
- CHARACTER*5 A, Blank
- DATA Blank /' '/, Zero /0.0d0/, Half /0.5d0/, One /1.0d0/
-C
+!
+!
+! -----------------------------------------------------------------
+!
+ SUBROUTINE Readclm (calc, N_Osc_Blank)
+!
+! *** Purpose -- Read CLM file for real
+!
+ use fixedr_m, only : Emax, Emin, Temp
+ IMPLICIT None
+ class(CrystalLatticeBroadening)::calc
+ integer::N_Osc_Blank
+ CHARACTER(80):: Title
+ CHARACTER(5)::A
+ CHARACTER(5), parameter::Blank = ' '
+ real(kind=8),parameter::Zero=0.0d0, Half=0.5d0, One=1.0d0
+ real(kind=8)::Bbbb, Cccc
+ integer::I, Kpound, N
+!
READ (10,10200,END=400,ERR=400) Title
10200 FORMAT (A80)
- READ (10,*) Mode_S_Norm, Nphon, Sub, Xdop, Eps, Epsc
- IF (Sub .EQ.Zero) Sub = One
- IF (Xdop.EQ.Zero) Xdop = One
- IF (Eps .EQ.Zero) Eps = 0.080d0
- IF (Epsc.EQ.Zero) Epsc = 0.001d0
+ READ (10,*) calc%Mode_S_Norm, calc%Nphon, calc%Sub, calc%Xdop, calc%Eps, calc%Epsc
+ IF (calc%Sub .EQ.Zero) calc%Sub = One
+ IF (calc%Xdop.EQ.Zero) calc%Xdop = One
+ IF (calc%Eps .EQ.Zero) calc%Eps = 0.080d0
+ IF (calc%Epsc.EQ.Zero) calc%Epsc = 0.001d0
READ (10,10300) A
10300 FORMAT (A5)
IF (A.NE.Blank) STOP '[Stop in Readclm in clm/mclm1.f]'
-C
+!
READ (10,10300) A
CALL Convert_To_Caps (A, 5, Kpound)
-C
+!
IF (A.EQ.'CONTI') THEN
- IF (N_Contin.LE.0) THEN
- READ (10,*) Del_Phonon
-C *** Fixed-format where SAMMY counts the number of points
- N_Contin = - N_Contin
- READ (10,10500,END=400) (Phonon(I),I=1,N_Contin)
+ IF (calc%N_Contin.LE.0) THEN
+ READ (10,*) calc%Del_Phonon
+! *** Fixed-format where SAMMY counts the number of points
+ calc%N_Contin = - calc%N_Contin
+ call allocate_real_data(calc%Phonon,calc%N_Contin)
+ READ (10,10500,END=400) (calc%Phonon(I),I=1,calc%N_Contin)
10500 FORMAT (8F10.1)
ELSE
-C *** Free-format but user must count
- READ (10,*) Del_Phonon, N
- IF (N.NE.N_Contin) STOP '[Stop in Readclm in mclm1.f # 2]'
- READ (10,*) (Phonon(I),I=1,N_Contin)
+! *** Free-format but user must count
+ READ (10,*) calc%Del_Phonon, N
+ IF (N.NE.calc%N_Contin) STOP '[Stop in Readclm in mclm1.f # 2]'
+ call allocate_real_data(calc%Phonon,calc%N_Contin)
+ READ (10,*) (calc%Phonon(I),I=1,calc%N_Contin)
END IF
READ (10,10300) A
-C *** Done reading Continuous Mode; Read next card
+! *** Done reading Continuous Mode; Read next card
READ (10,10300) A
CALL Convert_To_Caps (A, 5, Kpound)
END IF
-C
+!
IF (A.EQ.'DISCR') THEN
IF (N_Osc_Blank.EQ.1) READ (10,10300) A
- IF (N_Osc.LE.0) THEN
-C *** Fixed-format where SAMMY counts the number of points
- N_Osc = - N_Osc
- READ (10,10500,END=400) (Osc_Eng(I),I=1,N_Osc)
- READ (10,10500,END=400) (Osc_Wts(I),I=1,N_Osc)
+ IF (calc%N_Osc.LE.0) THEN
+! *** Fixed-format where SAMMY counts the number of points
+ calc%N_Osc = - calc%N_Osc
+ call allocate_real_data(calc%Osc_Wts,calc%N_Osc)
+ call allocate_real_data(calc%Osc_Eng,calc%N_Osc)
+ READ (10,10500,END=400) (calc%Osc_Eng(I),I=1,calc%N_Osc)
+ READ (10,10500,END=400) (calc%Osc_Wts(I),I=1,calc%N_Osc)
ELSE
-C *** Free-format but user must count
+! *** Free-format but user must count
READ (10,*) N
- IF (N.NE.N_Osc) STOP '[Stop in Readclm in mclm1.f # 3]'
- READ (10,*) (Osc_Eng(I),I=1,N_Osc)
- READ (10,*) (Osc_Wts(I),I=1,N_Osc)
+ IF (N.NE.calc%N_Osc) STOP '[Stop in Readclm in mclm1.f # 3]'
+ call allocate_real_data(calc%Osc_Wts, calc%N_Osc)
+ call allocate_real_data(calc%Osc_Eng,calc%N_Osc)
+ READ (10,*) (calc%Osc_Eng(I),I=1,calc%N_Osc)
+ READ (10,*) (calc%Osc_Wts(I),I=1,calc%N_Osc)
END IF
READ (10,10300,END=400) A
-C *** Done reading Discrete Mode; Read next card
+! *** Done reading Discrete Mode; Read next card
READ (10,10300,END=400) A
CALL Convert_To_Caps (A, 5, Kpound)
END IF
-C
+!
IF (A.EQ.'TRANS') THEN
- READ (10,*) Twt, C_Trans, Tbeta
- WRITE (21,10600) Twt, C_Trans, Tbeta
+ READ (10,*) calc%Twt, calc%C_Trans, calc%Tbeta
+ WRITE (21,10600) calc%Twt, calc%C_Trans, calc%Tbeta
10600 FORMAT (//, 'Twt, C, Tbeta=', 1P6G14.6)
END IF
- IF (Tbeta.EQ.Zero) Tbeta = One
-C
+ IF (calc%Tbeta.EQ.Zero) calc%Tbeta = One
+!
400 CONTINUE
-C
-C *** Printing the input almost as DOPUSH does
+!
+! *** Printing the input almost as DOPUSH does
WRITE (21,10650) Title
10650 FORMAT (//, 1X, A40)
- WRITE (21,10660) Nphon, Temp, Emin, Emax, Xdop, Eps*100, Epsc*100
-10660 FORMAT (//,
- & ' REQUIRED PHONON-EXPANSION ORDER .......... ', I10, /,
- & ' TEMPERATURE............................... ', F10.3, /,
- & ' MINIMAL ENERGY ........................... ', F10.3, /,
- & ' MAXIMAL ENERGY............................ ', F10.3, /,
- & ' INITIAL XDOP INTERVAL .................... ', F10.3, /,
- & ' NORMALISATION PRECISION .................. ', F10.3, '%', /,
- & ' INTEGRATION PRECISION .................... ', F10.3, '%')
-C
+ WRITE (21,10660) calc%Nphon, Temp, Emin, Emax, calc%Xdop, calc%Eps*100, calc%Epsc*100
+10660 FORMAT (//, &
+ ' REQUIRED PHONON-EXPANSION ORDER .......... ', I10, /, &
+ ' TEMPERATURE............................... ', F10.3, /, &
+ ' MINIMAL ENERGY ........................... ', F10.3, /, &
+ ' MAXIMAL ENERGY............................ ', F10.3, /, &
+ ' INITIAL XDOP INTERVAL .................... ', F10.3, /, &
+ ' NORMALISATION PRECISION .................. ', F10.3, '%', /,&
+ ' INTEGRATION PRECISION .................... ', F10.3, '%')
+!
CLOSE (UNIT=10)
-C
- WRITE (21,10700) N_Contin
+!
+ WRITE (21,10700) calc%N_Contin
10700 FORMAT (/, ' Number of points in phonon distribution =', I5)
- Bbbb = Phonon(1)*Half
- DO I=2,N_Contin-1
- Bbbb = Bbbb + Phonon(I)
+ Bbbb = calc%Phonon(1)*Half
+ DO I=2,calc%N_Contin-1
+ Bbbb = Bbbb + calc%Phonon(I)
END DO
- Bbbb = Bbbb + Phonon(N_Contin)*Half
- Cccc = Bbbb*Del_Phonon
- WRITE (21,10800) Bbbb, Del_Phonon, Cccc
-10800 FORMAT (' Sum of points=', 1PG14.6, ' Spacing=',G14.6,
- * ' *=',G14.6)
- WRITE (21,10900) (Phonon(I),I=1,N_Contin)
+ Bbbb = Bbbb + calc%Phonon(calc%N_Contin)*Half
+ Cccc = Bbbb*calc%Del_Phonon
+ WRITE (21,10800) Bbbb, calc%Del_Phonon, Cccc
+10800 FORMAT (' Sum of points=', 1PG14.6, ' Spacing=',G14.6, &
+ ' *=',G14.6)
+ WRITE (21,10900) (calc%Phonon(I),I=1,calc%N_Contin)
10900 FORMAT (1P5G14.6)
- IF (Phonon(1).NE.Zero) THEN
- WRITE (6,11000) Phonon(1)
+ IF (calc%Phonon(1).NE.Zero) THEN
+ WRITE (6,11000) calc%Phonon(1)
11000 FORMAT ('First term must be zero, not', 1PG14.6)
STOP '[STOP in Readclm in clm/mclm1.f]'
END IF
-C
+!
RETURN
END
-C
-C
-C -----------------------------------------------------------------
-C
- SUBROUTINE Initclm (Save, Phonon, Osc_Eng, Osc_Wts, Ppp, Beta,
- * Osc_Enx, Osc_Snth, Osc_Coth, Bex, Tlast, Tnow, Savex)
-C
-C *** Purpose -- Initialize things
-C
- use ifwrit_m
- use fixedr_m
- use clm_common_m
- use constn_common_m
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION Phonon(*), Osc_Eng(*), Osc_Wts(*), Ppp(*), Beta(*),
- * Osc_Enx(*), Osc_Snth(*), Osc_Coth(*), Bex(*)
- DIMENSION Save(Nbeta_Max,*), Tlast(*), Tnow(*), Savex(Nbeta_Max,*)
-C
-C
- Tev = Boltzm*Temp
-C *** Effective temperature calculation for creation of beta grid
- Del_Beta = Del_Phonon/Tev
- IF (N_Contin.GT.0) THEN
- CALL Start_Clm (Phonon, Ppp, Del_Beta)
-C *** INPUT -- Phonon, Tbeta, Del_Beta, N_Contin
-C *** OUTPUT -- Ppp, F0, Tbar
+!
+!
+! -----------------------------------------------------------------
+!
+ SUBROUTINE Initclm (calc)
+!
+! *** Purpose -- Initialize things
+!
+ use fixedr_m, only : Temp
+ use constn_common_m, only : Boltzm
+ use mclm2_M
+ IMPLICIT None
+ class(CrystalLatticeBroadening)::calc
+ real(kind=8)::Del_Beta, Tevf, Tevf_true
+!
+!
+ calc%Tev = Boltzm*Temp
+! *** Effective temperature calculation for creation of beta grid
+ Del_Beta = calc%Del_Phonon/calc%Tev
+ IF (calc%N_Contin.GT.0) THEN
+ CALL Start_Clm (calc, Del_Beta)
+! *** INPUT -- Phonon, Tbeta, Del_Beta, N_Contin
+! *** OUTPUT -- Ppp, F0, Tbar
END IF
- Tevf = Tev*Tbar
- Dwpix = F0
+ Tevf = calc%Tev*calc%Tbar
+ calc%Dwpix = calc%F0
Tevf_true = Tevf
-C
-C *** Initial beta mesh calculation
- CALL Mesh (Beta, Tevf)
-C *** Input -- Tevf, Tev, Emax, Xdop, Del_Phonon, Sub
-C *** Output -- Beta(Nbeta), Del_S_B, Nbeta
-C
- IF (N_Contin.GT.0) THEN
-C *** Contin_X initializes things, generates array Save
- CALL Contin_X (Ppp, Save, Tlast, Tnow, Savex, Del_Beta)
-C *** Input -- Pppppp(.), Alpha, Omega, F0, Tbeta, Del_Beta, Sub,
-C *** N_Contin, Nphon, Nbeta
-C *** Output -- Save(?,Nphon)
+!
+! *** Initial beta mesh calculation
+ CALL Mesh (calc, Tevf)
+! *** Input -- Tevf, Tev, Emax, Xdop, calc%Del_Phonon, Sub
+! *** Output -- Beta(Nbeta), Del_S_B, Nbeta
+!
+ IF (calc%N_Contin.GT.0) THEN
+! *** Contin_X initializes things, generates array Save
+ CALL Contin_X (calc, Del_Beta)
+! *** Input -- Pppppp(.), Alpha, Omega, F0, Tbeta, Del_Beta, Sub,
+! *** N_Contin, Nphon, Nbeta
+! *** Output -- Save(?,Nphon)
END IF
-C
- IF (N_Osc.GT.0) THEN
-C *** Discre_X initializes things for Discre
- CALL Discre_X (Osc_Eng, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, Beta)
-C *** Input -- Osc_Eng(), Osc_Wts(), Tev, Beta(Nbeta)
-C *** Output -- Osc_Enx()
+!
+ IF (calc%N_Osc.GT.0) THEN
+! *** Discre_X initializes things for Discre
+ CALL Discre_X (calc)
+! *** Input -- Osc_Eng(), Osc_Wts(), Tev, Beta(Nbeta)
+! *** Output -- Osc_Enx()
END IF
-C
-C
+!
+!
RETURN
END
+end module Readclm_m
diff --git a/sammy/src/clm/mclm2.f90 b/sammy/src/clm/mclm2.f90
index 7f3cee7f259b8c28eeaddc8eac8492dd59394e56..b4dd8ab46c65e0f9ed6149542ca1cafe0bf88384 100644
--- a/sammy/src/clm/mclm2.f90
+++ b/sammy/src/clm/mclm2.f90
@@ -1,359 +1,361 @@
-C
-C
-C -------------------------------------------------------------------
-C
- SUBROUTINE Start_Clm (Phonon, Ppp, Del_Beta)
-C
-C *** Purpose -- Compute several integral functions of the phonon
-C *** frequency distribution.
-C *** Called by -- Initclm
-C *** Input -- Phonon(.), Tbeta, Del_Beta, N_Contin
-C *** Output -- Ppp(.), F0, Tbar
-C *** Local -- Q, Annnnn, Ubeta, V, Vv, Fs0, Fs2, W2
-C
- use clm_common_m
- IMPLICIT NONE
+module mclm2_M
+use CrystalLatticeBroadening_M
+use AllocateFunctions_m
+implicit none
+contains
+!
+!
+! -------------------------------------------------------------------
+!
+ SUBROUTINE Start_Clm (calc, Del_Beta)
+!
+! *** Purpose -- Compute several integral functions of the phonon
+! *** frequency distribution.
+! *** Called by -- Initclm
+! *** Input -- Phonon(.), Tbeta, Del_Beta, N_Contin
+! *** Output -- Ppp(.), F0, Tbar
+! *** Local -- Q, Annnnn, Ubeta, V, Vv, Fs0, Fs2, W2
+!
+ class(CrystalLatticeBroadening)::calc
INTEGER I
- DIMENSION Phonon(*), Ppp(*)
- DOUBLE PRECISION Phonon, Ppp, Del_Beta,
- * Q, Annnnn, Ubeta, V, Vv, Fs0, Fs2, W2
+ DOUBLE PRECISION Del_Beta, &
+ Q, Annnnn, Ubeta, V, Vv, Fs0, Fs2, W2
DOUBLE PRECISION Zero, Half, One
DATA Zero /0.0d0/, Half /0.5d0/, One /1.0d0/
-C
-C
-C *** ---
-C *** Copy input spectrum into Ppp array
- DO I=1,N_Contin
- Ppp(I) = Phonon(I)
+!
+!
+! *** ---
+! *** Copy input spectrum into Ppp array
+ call allocate_real_data(calc%Ppp, calc%N_Contin)
+ DO I=1,calc%N_Contin
+ calc%Ppp(I) = calc%Phonon(I)
END DO
-C
-C *** ---
-C *** Calculate normalizing constant = Annnnn
- Annnnn = Ppp(1)*Half
- DO I=2,N_Contin-1
- Annnnn = Annnnn + Ppp(I)
+!
+! *** ---
+! *** Calculate normalizing constant = Annnnn
+ Annnnn = calc%Ppp(1)*Half
+ DO I=2,calc%N_Contin-1
+ Annnnn = Annnnn + calc%Ppp(I)
END DO
- Annnnn = Annnnn + Ppp(N_Contin)*Half
-C
-C *** ---
-C *** Calculate Debye-Waller Lambda and effective temperature
-C *** at the same time as calculate Ppp = Phonon/sinth/beta
-C *** From here to next "C * ---" used to be Fsum(0...) and Fsum(2...)
-cq Ppp(1) = Ppp(2)/Del_Beta**2 Clearly this is wrong
-cq Ppp(1) = Ppp(1)/Del_Beta**2 This is also wrong
-cq Ppp(1) = Ppp(1)/Zero This is right! but ill-defined
- Ppp(1) = Ppp(1)/Del_Beta**2
+ Annnnn = Annnnn + calc%Ppp(calc%N_Contin)*Half
+!
+! *** ---
+! *** Calculate Debye-Waller Lambda and effective temperature
+! *** at the same time as calculate calc%Ppp = calc%Phonon/sinth/beta
+! *** From here to next "C * ---" used to be Fsum(0...) and Fsum(2...)
+!q calc%Ppp(1) = calc%Ppp(2)/Del_Beta**2 Clearly this is wrong
+!q calc%Ppp(1) = calc%Ppp(1)/Del_Beta**2 This is also wrong
+!q calc%Ppp(1) = calc%Ppp(1)/Zero This is right! but ill-defined
+ calc%Ppp(1) = calc%Ppp(1)/Del_Beta**2
Ubeta = Del_Beta
V = dEXP(Half*Del_Beta)
Vv = V
-cq Fs0 = Ppp(1)
+!q Fs0 = calc%Ppp(1)
Fs0 = Zero
Fs2 = Zero
- DO I=2,N_Contin-1
+ DO I=2,calc%N_Contin-1
W2 = Ubeta**2
- Ppp(I) = ( Ppp(I)/Ubeta ) / (Vv-One/Vv)
- Q = Ppp(I) * (Vv+One/Vv)
-C nml so P(I )= Phonon(I+1 ) / sinth(I Del_Beta/2) / { I Del_Beta }
-C nml P(Beta)= Phonon(Beta) / sinth( Beta/2) / { Beta }
-C nml Q = Phonon(I+1 ) * coth (I Del_Beta/2) / { I Del_Beta }
-C nml Q(Beta)= Phonon(Beta) * coth ( Beta/2) / { Beta }
+ calc%Ppp(I) = ( calc%Ppp(I)/Ubeta ) / (Vv-One/Vv)
+ Q = calc%Ppp(I) * (Vv+One/Vv)
+! nml so P(I )= calc%Phonon(I+1 ) / sinth(I Del_Beta/2) / { I Del_Beta }
+! nml P(Beta)= calc%Phonon(Beta) / sinth( Beta/2) / { Beta }
+! nml Q = calc%Phonon(I+1 ) * coth (I Del_Beta/2) / { I Del_Beta }
+! nml Q(Beta)= calc%Phonon(Beta) * coth ( Beta/2) / { Beta }
Fs0 = Fs0 + Q
Fs2 = Fs2 + Q * W2
Vv = V*Vv
Ubeta = Ubeta + Del_Beta
END DO
- Ppp(N_Contin) = ( Ppp(N_Contin)/Ubeta ) / (Vv-One/Vv)
- Q = Ppp(N_Contin) * (Vv+One/Vv)
+ calc%Ppp(calc%N_Contin) = ( calc%Ppp(calc%N_Contin)/Ubeta ) / (Vv-One/Vv)
+ Q = calc%Ppp(calc%N_Contin) * (Vv+One/Vv)
W2 = Ubeta**2
Fs0 = Fs0 + Q*Half
Fs2 = Fs2 + Q*Half*W2
-C *** ---
- F0 = ( Fs0 * Tbeta) / Annnnn
- Tbar = ( Fs2 * Half ) / Annnnn
-C
-C *** ---
-C *** Convert Ppp(Beta) into first term in phonon expansion sum
+! *** ---
+ calc%F0 = ( Fs0 * calc%Tbeta) / Annnnn
+ calc%Tbar = ( Fs2 * Half ) / Annnnn
+!
+! *** ---
+! *** Convert Ppp(Beta) into first term in phonon expansion sum
Q = Fs0*Del_Beta
- DO I=1,N_Contin
+ DO I=1,calc%N_Contin
Ubeta = Del_Beta*(I-1)
- Ppp(I) = Ppp(I)*dEXP(Half*Ubeta)/Q
+ calc%Ppp(I) = calc%Ppp(I)*dEXP(Half*Ubeta)/Q
END DO
-C
+!
RETURN
END
-C
-C
-C ----------------------------------------------------------------------
-C
- SUBROUTINE Mesh (Beta, Tevf)
-C
-C Beta GRID CALCULATION
-C
-C &&& Del_S_B - THE SOLID SPECTRUM GRID Mesh
-C Sub - NUMBER OF SubINTERVALS OF Beta GRID
-C Bmax - MAXIMAL VALUE OF Beta
-C Nbeta - NUMBER OF Beta VALUES IN THE INTERVAL (0,Bmax)
-C
-C *** Input -- Tevf, Tev, Awr, Emax, Xdop, Del_Phonon, Sub, Nbeta_Max
-C *** Output -- Beta(Nbeta), Del_S_B, Nbeta
-C
- use fixedr_m
- use clm_common_m
- use constn_common_m
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
+!
+!
+! ----------------------------------------------------------------------
+!
+ SUBROUTINE Mesh (calc, Tevf)
+!
+! Beta GRID CALCULATION
+!
+! &&& Del_S_B - THE SOLID SPECTRUM GRID Mesh
+! Sub - NUMBER OF SubINTERVALS OF Beta GRID
+! Bmax - MAXIMAL VALUE OF Beta
+! Nbeta - NUMBER OF Beta VALUES IN THE INTERVAL (0,Bmax)
+!
+! *** Input -- Tevf, Tev, Awr, Emax, Xdop, calc%Del_Phonon, Sub, Nbeta_Max
+! *** Output -- Beta(Nbeta), Del_S_B, Nbeta
+!
+ use fixedr_m, only : Aaawww, Emax
+ use constn_common_m, only : Aneutr
+ class(CrystalLatticeBroadening)::calc
INTEGER J
- DOUBLE PRECISION Beta(*), Bmax, Tevf
-C
- IF (Del_Phonon.GT.0.0d0) THEN
- Del_S_B = Del_Phonon/(Sub*Tev)
- Bmax = Xdop/Tev * dSQRT(4.0D0*Emax*Tevf*Aneutr/Aaawww)
- Nbeta = Int(Bmax/Del_S_B)
-C
- IF (Nbeta.GT.Nbeta_Max) THEN
- WRITE (6,10000) Nbeta, Nbeta_Max, Xdop
+ DOUBLE PRECISION Bmax, Tevf
+!
+ IF (calc%Del_Phonon.GT.0.0d0) THEN
+ calc%Del_S_B = calc%Del_Phonon/(calc%Sub*calc%Tev)
+ Bmax = calc%Xdop/calc%Tev * dSQRT(4.0D0*Emax*Tevf*Aneutr/Aaawww)
+ calc%Nbeta = Int(Bmax/calc%Del_S_B)
+!
+ IF (calc%Nbeta.GT.calc%Nbeta_Max) THEN
+ WRITE (6,10000) calc%Nbeta, calc%Nbeta_Max, calc%Xdop
10000 FORMAT ('Nbeta is too big', 2I10, 1P6G14.6)
END IF
-C
- Beta(1) = 0.0D0
- DO J=2,Nbeta
- Beta(J) = Beta(J-1) + Del_S_B
+!
+ call allocate_real_data(calc%Beta, calc%Nbeta)
+ calc%Beta(1) = 0.0D0
+ DO J=2,calc%Nbeta
+ calc%Beta(J) = calc%Beta(J-1) + calc%Del_S_B
END DO
ELSE
- Nbeta = 0
- END IF
+ calc%Nbeta = 0
+ END IF
RETURN
END
-C
-C
-C -------------------------------------------------------------------
-C
- SUBROUTINE Contin_X (Ppp, Save, Tlast, Tnow, Savex, Del_Beta)
-C
-C *** Purpose -- Generate array Save(Kbeta,Kphon) for use in Contin
-C *** Input -- Ppp(.), Alpha, Omega, F0, Tbeta, Del_Beta, Sub,
-C *** N_Contin, Nphon, Nbeta
-C *** Output -- Save (?,Nphon)
-C *** Dummy -- Savex(?,Nphon), Tlast(?), Tnow(?)
-C
- use fixedr_m
- use clm_common_m
- use constn_common_m
- IMPLICIT DOUBLE PRECISION (A-H,O-Z)
- DOUBLE PRECISION Ppp(*), Save(Nbeta_Max,*),
- & Tlast(*), Tnow(*), Savex(Nbeta_Max,*)
- DOUBLE PRECISION Del_Beta, Exx, Exy, Alf0, Omf0, B, C, Bt, Btp, X
- INTEGER NLast, Nnow, I, K, N
- DOUBLE PRECISION Zero, One
- PARAMETER (Zero=0.0d0, One=1.0d0)
-C
-C *** Note that Ppp(Beta) = Phonon(Beta) * dEXP(Beta/2) /
-C *** [Beta sinh(Beta/2) normalization]
-C
-C *** ------
-C *** First term in the phonon expansion sum
-C
- CALL Zero_Array (Save , Nphon*Nbeta_Max)
- CALL Zero_Array (Savex, Nphon*Nbeta_Max)
-C
- DO I=1,N_Contin
- Tlast(I) = Ppp(I)
- Savex(I,1) = Ppp(I)
+!
+!
+! -------------------------------------------------------------------
+!
+ SUBROUTINE Contin_X (calc, Del_Beta)
+!
+! *** Purpose -- Generate array Save(Kbeta,Kphon) for use in Contin
+! *** Input -- Ppp(.), Alpha, Omega, F0, Tbeta, Del_Beta, Sub,
+! *** N_Contin, Nphon, Nbeta
+! *** Output -- Save (?,Nphon)
+! *** Dummy -- Savex(?,Nphon)
+!
+ use fixedr_m, only : Aaawww, Emax, Emin
+ use constn_common_m, only : Aneutr
+ class(CrystalLatticeBroadening)::calc
+ DOUBLE PRECISION Del_Beta, Exx, Exy, Alf0, Omf0, B, C, Bt, Btp, X, Rn
+ INTEGER NLast, Nnow, I, K, N
+ real(kind=8),parameter::Zero=0.0d0, One=1.0d0
+!
+! *** Note that Ppp(Beta) = Phonon(Beta) * dEXP(Beta/2) /
+! *** [Beta sinh(Beta/2) normalization]
+!
+! *** ------
+! *** First term in the phonon expansion sum
+!
+!
+ call allocate_real_data(calc%Tlast, calc%N_Contin)
+ call reallocate_real_data_2d(calc%savex, calc%NBeta, 0, calc%Nphon, 0)
+ DO I=1,calc%N_Contin
+ calc%Tlast(I) = calc%Ppp(I)
+ calc%Savex(I,1) = calc%Ppp(I)
END DO
- N = Nphon
- NLast = N_Contin
-C
-C
-C *** ------
-C *** Add other terms in the phonon expansion sum
- IF (Tbeta.NE.Zero) THEN
+ N = calc%Nphon
+ NLast = calc%N_Contin
+!
+!
+! *** ------
+! *** Add other terms in the phonon expansion sum
+ IF (calc%Tbeta.NE.Zero) THEN
Rn = Aneutr/(Aaawww+Aneutr)
- Alf0 = Emin*Rn/Tev * F0
- Omf0 = Emax*Rn/Tev * F0
+ Alf0 = Emin*Rn/calc%Tev * calc%F0
+ Omf0 = Emax*Rn/calc%Tev * calc%F0
Exx = Alf0*dEXP(-Alf0)
Exy = Omf0*dEXP(-Omf0)
B = Exx
C = Exy
- DO N=2,Nphon
+ DO N=2,calc%Nphon
Exx = Alf0*Exx/N
Exy = Omf0*Exy/N
IF (B+Exx.EQ.B .AND. C+Exy.EQ.C) GO TO 10
B = B + Exx
C = C + Exy
- Nnow = N_Contin + NLast - 1
- CALL Convol (Ppp, Tlast, Tnow, N_Contin, NLast, Nnow,
- & Del_Beta)
+ Nnow = calc%N_Contin + NLast - 1
+ CALL Convol (calc, NLast, Nnow, Del_Beta)
+ call reallocate_real_data(calc%Tlast, Nnow, 10)
+ call reallocate_real_data(calc%Tnow, Nnow, 10)
+ call reallocate_real_data_2d(calc%savex, Nnow, 0, calc%Nphon, 0)
DO I=1,Nnow
- Tlast(I) = Tnow(I)
- Savex(I,N) = Tnow(I)
+ calc%Tlast(I) = calc%Tnow(I)
+ calc%Savex(I,N) = calc%Tnow(I)
END DO
NLast = Nnow
END DO
- IF (N.GT.Nphon) N = Nphon
+ IF (N.GT.calc%Nphon) N = calc%Nphon
END IF
-C
+!
10 CONTINUE
- IF (Nphon.EQ.N) WRITE (6,10000) Nphon
-10000 FORMAT (' ### Caution ### Number of phonons may not be enough to',
- & ' ###', /,
- & ' ### reach convergence. Try more than', I4,
- & ' ###')
- Nphon = N
-C
- DO K=1,Nbeta
- X = (K-1)/Sub
+ IF (calc%Nphon.EQ.N) WRITE (6,10000) calc%Nphon
+10000 FORMAT (' ### Caution ### Number of phonons may not be enough to', &
+ ' ###', /, &
+ ' ### reach convergence. Try more than', I4, &
+ ' ###')
+ calc%Nphon = N
+!
+ call reallocate_real_data_2d(calc%save, calc%NBeta, 0, calc%Nphon, 0)
+ calc%save = 0.0d0
+ DO K=1,calc%Nbeta
+ X = (K-1)/calc%Sub
I = X
-C *** Note that equally-spaced points in Beta are assumed throughout
+! *** Note that equally-spaced points in Beta are assumed throughout
IF (I.LT.0 .OR. I.GE.Nnow-1) THEN
-C *** Don't need to zero Save cuz is already zero
+! *** Don't need to zero Save cuz is already zero
ELSE IF (I.EQ.0) THEN
Bt = Zero
Btp = One
I = I + 1
- DO N=1,Nphon
- Save(K,N) = Savex(I,N)*(Btp-X) + X*Savex(I+1,N)
+ DO N=1,calc%Nphon
+ calc%Save(K,N) = calc%Savex(I,N)*(Btp-X) + X*calc%Savex(I+1,N)
END DO
ELSE IF (I.LT.Nnow-1) THEN
Bt = dFLOAT(I)
Btp = Bt + One
I = I + 1
- DO N=1,Nphon
- Save(K,N) = Savex(I,N)*(Btp-X) + (X-Bt)*Savex(I+1,N)
+ DO N=1,calc%Nphon
+ calc%Save(K,N) = calc%Savex(I,N)*(Btp-X) + (X-Bt)*calc%Savex(I+1,N)
END DO
ELSE
STOP '[Should never get here in Contin_X in clm/mclm3.f]'
END IF
END DO
-C
+!
RETURN
END
-C
-C
-C ----------------------------------------------------------------------
-C
- SUBROUTINE Convol (Ppp, TLast, Tnow, N_Contin, NLast, Nnow,
- & Del_Beta)
-C
-C *** PURPOSE -- Calculate the next term in the phonon expansion by
-C *** convolving P with Tlast and writing the result into
-C *** Tnow. The integral of Tnow is also checked.
-C *** METHOD -- Trapazoidal integration
-C *** CALLED BY -- Contin_X
-C
- IMPLICIT NONE
- INTEGER N_Contin, NLast, Nnow, K, J, I1, I2, Kmax
- DOUBLE PRECISION Ppp(*), Tlast(*), Tnow(*),
- & EP, Be, Del_Beta, F2, Cc, F1, Zero, Half, One, Small
- DATA Zero /0.0d0/, Half /0.5d0/, One /1.0d0/, Small /1.d-20/
-C
+!
+!
+! ----------------------------------------------------------------------
+!
+ SUBROUTINE Convol (calc, NLast, Nnow, Del_Beta)
+!
+! *** PURPOSE -- Calculate the next term in the phonon expansion by
+! *** convolving P with Tlast and writing the result into
+! *** Tnow. The integral of Tnow is also checked.
+! *** METHOD -- Trapazoidal integration
+! *** CALLED BY -- Contin_X
+!
+ class(CrystalLatticeBroadening)::calc
+ INTEGER NLast, Nnow, K, J, I1, I2, Kmax
+ DOUBLE PRECISION EP, Be, Del_Beta, F2, Cc, F1
+ real(kind=8),parameter:: Zero = 0.0d0, Half = 0.5d0, One = 1.0d0, Small = 1.d-20
+!
Kmax = 0
+ call reallocate_real_data(calc%Tnow, Nnow, 10)
DO K=1,Nnow
- Tnow(K) = Zero
+ calc%Tnow(K) = Zero
EP = One
-C
- DO J=1,N_Contin
+!
+ DO J=1,calc%N_Contin
I1 = K + J - 2
I2 = K - J
F1 = Zero
Be = (J-1)*Del_Beta
-C
- IF (Ppp(J).GT.Zero) THEN
- IF (I1+1.LE.NLast) F1 = Tlast(I1+1)*dEXP(-Be)
+!
+ IF (calc%Ppp(J).GT.Zero) THEN
+ call reallocate_real_data(calc%Tlast, NLast, 10)
+ IF (I1+1.LE.NLast) F1 = calc%Tlast(I1+1)*dEXP(-Be)
F2 = Zero
IF (I2.GE.0 .AND. I2+1.LE.NLast) THEN
- F2 = Tlast(I2+1)
+ F2 = calc%Tlast(I2+1)
ELSE IF (I2.LT.0 .AND. 1-I2.LE.NLast) THEN
Be = -I2*Del_Beta
- F2 = Tlast(1-I2)*dEXP(-Be)
+ F2 = calc%Tlast(1-I2)*dEXP(-Be)
ENDIF
- Cc = Ppp(J)*(F1+F2)
- IF (J.EQ.1 .OR. J.EQ.N_Contin) Cc = Half*Cc
- Tnow(K) = Tnow(K)+Cc
+ Cc = calc%Ppp(J)*(F1+F2)
+ IF (J.EQ.1 .OR. J.EQ.calc%N_Contin) Cc = Half*Cc
+ calc%Tnow(K) = calc%Tnow(K)+Cc
END IF
-C
+!
END DO
-C
- Tnow(K) = Tnow(K)*Del_Beta
- IF (Tnow(K).LT.Small) THEN
- Tnow(K) = Zero
+!
+ calc%Tnow(K) = calc%Tnow(K)*Del_Beta
+ IF (calc%Tnow(K).LT.Small) THEN
+ calc%Tnow(K) = Zero
ELSE
IF (K.GT.Kmax) Kmax = K
END IF
END DO
-C
+!
IF (Kmax.LT.Nnow) Nnow = Kmax
RETURN
END
-C
-C
-C ----------------------------------------------------------------------
-C
- SUBROUTINE Discre_X (Osc_Eng, Osc_Wts, Osc_Enx, Osc_Snth,
- * Osc_Coth, Bex, Beta)
-C
-C
-C *** PURPOSE -- Do energy-independent bits of Discre calculations
-C *** Called by -- Samclm_0
-C *** Uses -- no routines
-C
-C
- use fixedr_m
- use clm_common_m
- use constn_common_m
- IMPLICIT DOUBLE PRECISION (A-H,O-Z)
- DOUBLE PRECISION Osc_Eng(*), Osc_Wts(*), Osc_Enx(*), Osc_Snth(*),
- & Osc_Coth(*), Bex(*), Beta(Nbeta)
- DOUBLE PRECISION Sn, Cn, Eb
- INTEGER K, I
- DOUBLE PRECISION Zero, Half, One
- PARAMETER (Zero=0.0d0, Half=0.5d0, One=1.0d0)
-C
-C *** SET UP OSCILLATOR PARAMETERS
-C
- Sum_Osc_Wts = Zero
- DO I=1,N_Osc
- Sum_Osc_Wts = Sum_Osc_Wts + Osc_Wts(I)
+!
+!
+! ----------------------------------------------------------------------
+!
+ SUBROUTINE Discre_X (calc)
+!
+!
+! *** PURPOSE -- Do energy-independent bits of Discre calculations
+! *** Called by -- Samclm_0
+! *** Uses -- no routines
+!
+!
+ class(CrystalLatticeBroadening)::calc
+ real(kind=8)::Sn, Cn, Eb
+ INTEGER::K, I
+ real(kind=8),parameter::Zero=0.0d0, Half=0.5d0, One=1.0d0
+!
+! *** SET UP OSCILLATOR PARAMETERS
+!
+ calc%Sum_Osc_Wts = Zero
+ DO I=1,calc%N_Osc
+ calc%Sum_Osc_Wts = calc%Sum_Osc_Wts + calc%Osc_Wts(I)
END DO
-C
- DO I=1,N_Osc
- Osc_Enx(I) = Osc_Eng(I)/Tev
+!
+ call allocate_real_data(calc%Osc_Enx, calc%N_Osc)
+ DO I=1,calc%N_Osc
+ calc%Osc_Enx(I) = calc%Osc_Eng(I)/calc%Tev
END DO
-C
- Tsave = Zero
- Dw0 = Dwpix
- DO I=1,N_Osc
- Eb = dEXP(Half*Osc_Enx(I))
+!
+ calc%Tsave = Zero
+ calc%Dw0 = calc%Dwpix
+ call allocate_real_data(calc%Osc_Snth, calc%N_Osc)
+ call allocate_real_data(calc%Osc_Coth, calc%N_Osc)
+ DO I=1,calc%N_Osc
+ Eb = dEXP(Half*calc%Osc_Enx(I))
Sn = Half* (Eb-One/Eb)
Cn = Half* (Eb+One/Eb)
- Osc_Snth(I) = Osc_Wts(I)/(Sn*Osc_Enx(I))
- Osc_Coth(I) = Osc_Snth(I)*Cn
- IF (Dwpix.GT.Zero) Dwpix = Dwpix + Osc_Coth(I)
- Tsave = Tsave + Osc_Wts(I) * Osc_Eng(I) *Cn/Sn
+ calc%Osc_Snth(I) = calc%Osc_Wts(I)/(Sn*calc%Osc_Enx(I))
+ calc%Osc_Coth(I) = calc%Osc_Snth(I)*Cn
+ IF (calc%Dwpix.GT.Zero) calc%Dwpix = calc%Dwpix + calc%Osc_Coth(I)
+ calc%Tsave = calc%Tsave + calc%Osc_Wts(I) * calc%Osc_Eng(I) *Cn/Sn
END DO
-C
-C
- K = Nbeta
- DO I=1,Nbeta
- Bex(I) = - Beta(K)
+!
+!
+ K = calc%Nbeta
+ call allocate_real_data(calc%Bex, 2*calc%Nbeta)
+ DO I=1,calc%Nbeta
+ calc%Bex(I) = - calc%Beta(K)
K = K - 1
END DO
-C
- IF (Beta(1).LE.0.000000001) THEN
- Bex(Nbeta) = Zero
- K = Nbeta + 1
+!
+ IF (calc%Beta(1).LE.0.000000001) THEN
+ calc%Bex(calc%Nbeta) = Zero
+ K = calc%Nbeta + 1
ELSE
- K = Nbeta + 2
- Bex(Nbeta+1) = Beta(1)
+ K = calc%Nbeta + 2
+ calc%Bex(calc%Nbeta+1) = calc%Beta(1)
END IF
-C
- DO I=2,Nbeta
- Bex(K) = Beta(I)
+!
+ DO I=2,calc%Nbeta
+ calc%Bex(K) = calc%Beta(I)
K = K + 1
END DO
- Nbx = K - 1
- IF (Nbx.GT.Nbeta_Max) THEN
- WRITE (6,10000) Nbx, Nbeta_Max
+ calc%Nbx = K - 1
+ IF (calc%Nbx.GT.calc%Nbeta_Max) THEN
+ WRITE (6,10000) calc%Nbx, calc%Nbeta_Max
10000 FORMAT ('Nbeta_Max must be', I5, ' but is only', I5)
STOP '[Stop in Discre_X in clm/mclm2.f]'
END IF
RETURN
END
+end module mclm2_M
diff --git a/sammy/src/clm/mclm3.f90 b/sammy/src/clm/mclm3.f90
index 7dc8be21e98cc6ef92e7c39eab3fbecdb8ff011c..98cd1234a71ed4f470852a675c08feb6c6e250f3 100644
--- a/sammy/src/clm/mclm3.f90
+++ b/sammy/src/clm/mclm3.f90
@@ -1,47 +1,35 @@
-C
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Dopclm (Iflmsc,
- * Beta, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, S_Expb, S_Ex, Save, Bs, Ss, Max_T)
+module mclm3_m
+use CrystalLatticeBroadening_M
-C
-C *** PURPOSE -- FORM DOPPLER-BROADENED CROSS SECTION AND DERIVATIVES
-C
- use fixedi_m, only : Ktruet, Lllmax,
- * Numiso, numUsedPar
- use ifwrit_m, only : Kcros, Kdebug, Kfinit, Ksindi, Ksitmp,
- * Kvtemp, Kvthck, Nonu
- use fixedr_m, only : Emax, Emaxs, Emin, Emins, Temp, Thick,
- * Sitemp
+contains
+!
+!
+! --------------------------------------------------------------
+!
+ SUBROUTINE Dopclm (calc, Iflmsc)
+
+!
+! *** PURPOSE -- FORM DOPPLER-BROADENED CROSS SECTION AND DERIVATIVES
+!
+ use fixedi_m, only : Lllmax
+ use ifwrit_m, only : Kcros, Ksindi, Ksitmp,Kvtemp
+ use fixedr_m, only : Emax, Emaxs, Emin, Emins, Sitemp
use brdd_common_m, only : Iup, Kc
- use clm_common_m, only : Nbeta_Max, Tev
- use lbro_common_m, only : Yresol, Yssmsc, Ytrans
use constn_common_m, only : Boltzm
- use EndfData_common_m, only : expData, resparData,
- * radFitFlags
- use AuxGridHelper_M, only : setAuxGridOffset,
- * setAuxGridRowMax
+ use EndfData_common_m, only : expData, resparData,radFitFlags
+ use AuxGridHelper_M, only : setAuxGridOffset, setAuxGridRowMax
use SammyGridAccess_M
- use xct2_m
- use mxct27_m
- use mfgm3_M
+ use mfgm3_M, only : stetd
use Qtrap_Clm_m
+ use Getsab_m
use array_sizes_common_m, only : calcData, calcDataSelf
use convert_to_transmission_m
use DerivativeHandler_M
- use broad_common_m
- use SumIsoAndConvertToTrans_M
IMPLICIT None
- LOGICAL Need_Isotopes
- LOGICAL Another_Process_Will_Happen
type(SammyGridAccess)::grid, auxGrid
-C
+ class(CrystalLatticeBroadening)::calc
+!
integer::Iflmsc(*)
- real(kind=8):: Beta(*), Osc_Wts(*), Osc_Enx(*), Osc_Snth(*),
- * Osc_Coth(*), Bex(*), S_Expb(*), S_Ex(*), Save(Nbeta_Max,*),
- * Bs(*), Ss(*), Max_T(*)
real(kind=8)::C_Norm, eaux, em, Tevf, Tevx
integer::Isomax, Isox, Iv, Iw, J, Jdat, Kkkdat, insig
integer(C_SIZE_T)::Iso
@@ -52,73 +40,64 @@ C
call grid%setToExpGrid(expData)
call auxGrid%initialize()
call auxGrid%setToAuxGrid(expData)
- call setAuxGridOffset(1) ! reset auxillary grid starting point
- numEl = grid%getNumEnergies(expData)
- numElAux = auxGrid%getNumEnergies(expData)
+ call setAuxGridOffset(1) ! reset auxillary grid starting point
+ numElAux = calc%getNumEnergyUnbroadened()
+ numEl = calc%getNumEnergyBroadened()
call tmpCalc%initialize()
call tmpCalc%setUpList(resparData, radFitFlags, 1)
-C
-C *** Note that Sigxxx does not keep track of isotope, but calcData does
-C
-C DIMENSION Iflmsc(Nummsc), Parnbk(Numnbk), Iflnbk(Numnbk),
-C * Parbgf(Numbgf), Iflbgf(*), Kndbgf(*), Bgfmin(*), Bgfmax(*)
-C
- Tev = Boltzm*Temp
- Need_Isotopes = Yssmsc
- Another_Process_Will_Happen = Yresol.OR.Yssmsc
-C
+!
+! *** Note that Sigxxx does not keep track of isotope, but calcData does
+!
+! DIMENSION Iflmsc(Nummsc), Parnbk(Numnbk), Iflnbk(Numnbk),
+! * Parbgf(Numbgf), Iflbgf(*), Kndbgf(*), Bgfmin(*), Bgfmax(*)
+!
+ calc%Tev = Boltzm*calc%getTemperature()
+
+!
IF (Kcros.EQ.8 .AND. Ksindi.GT.0) THEN
Iv = 1
ELSE
Iv = 0
END IF
-C
- Iw = 0
- call calcData%nullify()
+!
+ Iw = 0
IF (Kcros.EQ.8 .AND. (Ksindi.GT.0 .OR. Ksitmp.GT.0)) Iw = 1
- call calcDataSelf%nullify()
- call calcData%reserve(numElAux*calcData%getNnnsig(), numUsedPar+1)
- IF ((Ksindi.GT.0.OR.Ksitmp.GT.0).AND. Kcros.EQ.8) THEN
- call calcDataSelf%reserve(numElAux, numUsedPar+1)
- end if
-C
-C
+!
+!
Now = 0
Nowx = 0
Kkkdat = 0
-C
+!
Isomax = calcData%getUsedIsotopes()
IF (Isomax.EQ.0) STOP '[STOP in Dopclm in clm/mclm3.f]'
-C
-C *** Do separately for each isotope (nuclide), since Doppler-width is
-C *** isotope-dependent
+!
+! *** Do separately for each isotope (nuclide), since Doppler-width is
+! *** isotope-dependent
DO Iso=1,Isomax
Isox = Iso
- if (Isomax.gt.1.and.
- * calcData%getRealIsotopeIndex(Iso).le.0) then
+ if (Isomax.gt.1.and. &
+ calcData%getRealIsotopeIndex(Iso).le.0) then
GO TO 80
end if
-C
+!
Kc = 1
Iup = 0
-C
+!
Kkkmin = 0
Kkkkkk = 0
Kkkdat = 0
nauxStart = 0
-C
-C *** Start of major loop over energy-points
+!
+! *** Start of major loop over energy-points
DO 70 J=1,numElAux
Jdat = J
-C
- IF ((J/1000)*1000.EQ.J .AND. Isomax.GT.1) WRITE (6,10000) J,
- * Isox
-10000 FORMAT (' *** on data point number', i10,
- * ' for nuclide number', i3)
+!
+ IF ((J/1000)*1000.EQ.J .AND. Isomax.GT.1) WRITE (6,10000) J, Isox
+10000 FORMAT (' *** on data point number', i10, ' for nuclide number', i3)
IF ((J/1000)*1000.EQ.J .AND. Isomax.EQ.1) WRITE (6,10001) J
10001 FORMAT (' *** on data point number', I10)
-C
- IF (Another_Process_Will_Happen) THEN
+!
+ IF (numElAux.eq.numEl) THEN
em = auxGrid%getEnergy(J, expData)
ELSE
IF (J.GT.numEl) THEN
@@ -127,24 +106,24 @@ C
Em = grid%getEnergy(J, expData)
END IF
END IF
-C
+!
IF (Em.LT.Emins) THEN
-C These points are the very low-energy limit,
-C of very little interest, ergo "stetd"
+! These points are the very low-energy limit,
+! of very little interest, ergo "stetd"
eaux = auxGrid%getEnergy(J + 1, expData)
IF (eaux.GE.Emin) THEN
- CALL Stetd (calcData, calcData,Kkkkkk+1,
- * Now, Kvtemp, Isox,
- * Jdat, 0)
+ CALL Stetd (calcData, calcData,Kkkkkk+1, &
+ Now, Kvtemp, Isox, &
+ Jdat, 0)
IF (Kcros.EQ.8) THEN
IF (Ksindi.GT.0) THEN
- CALL Stetd (calcDataSelf,calcDataSelf,
- * Kkkkkk+1, Nowx, Iflmsc(Ksitmp),
- * Isox, Jdat, 0)
+ CALL Stetd (calcDataSelf,calcDataSelf, &
+ Kkkkkk+1, Nowx, Iflmsc(Ksitmp), &
+ Isox, Jdat, 0)
ELSE IF (Ksitmp.GT.0) THEN
- CALL Stetd (calcData, calcDataSelf,
- * Kkkkkk+1, Nowx, Iflmsc(Ksitmp),
- * Isox, Jdat, Lllmax+1)
+ CALL Stetd (calcData, calcDataSelf, &
+ Kkkkkk+1, Nowx, Iflmsc(Ksitmp), &
+ Isox, Jdat, Lllmax+1)
END IF
END IF
GO TO 40
@@ -153,22 +132,22 @@ C of very little interest, ergo "stetd"
GO TO 70
END IF
ELSE IF (Em.GT.Emaxs) THEN
-C These points are the very high-energy limit,
-C of very little interest, ergo "stetd"
+! These points are the very high-energy limit,
+! of very little interest, ergo "stetd"
eaux = auxGrid%getEnergy(J - 1, expData)
IF (eaux.LE.Emax) THEN
- CALL Stetd (calcData, calcData,Kkkkkk+1,
- * Now, Kvtemp, Isox,
- * Jdat, 0)
+ CALL Stetd (calcData, calcData,Kkkkkk+1, &
+ Now, Kvtemp, Isox, &
+ Jdat, 0)
IF (Kcros.EQ.8) THEN
IF (Ksindi.GT.0) THEN
- CALL Stetd (calcDataSelf, calcDataSelf,
- * Kkkkkk+1, Nowx, Iflmsc(Ksitmp),
- * Isox, Jdat, 0)
+ CALL Stetd (calcDataSelf, calcDataSelf, &
+ Kkkkkk+1, Nowx, Iflmsc(Ksitmp), &
+ Isox, Jdat, 0)
ELSE IF (Ksitmp.gt.0) THEN
- Call Stetd (calcData, calcDataSelf,
- * Kkkkkk+1, Nowx, Iflmsc(Ksitmp),
- * Isox, Jdat, Lllmax+1)
+ Call Stetd (calcData, calcDataSelf, &
+ Kkkkkk+1, Nowx, Iflmsc(Ksitmp), &
+ Isox, Jdat, Lllmax+1)
END IF
END IF
GO TO 40
@@ -176,81 +155,73 @@ C of very little interest, ergo "stetd"
GO TO 80
END IF
END IF
-C
-C
-C "ELSE" but others have "GO TO" so not inside "IF" test.
-C ********* regular CLM Doppler for most cross sections
- CALL Getsab (Beta, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, S_Expb, S_Ex, Save, Bs, Ss, Max_T, Em, Tev, C_Norm,
- * Tevf, Ns)
- CALL Qtrap_Clm (tmpCalc,
- * calcData, Bs, Ss, Tev, Temp, Em,
- * C_Norm, Kvtemp, calcData%getNnnsig(),
- * Isox, 0, Ns, calcData,
- * Kkkkkk+1)
-C
-C ********* Doppler widths etc for self-indication transmission
+!
+!
+! "ELSE" but others have "GO TO" so not inside "IF" test.
+! ********* regular CLM Doppler for most cross sections
+ CALL Getsab (calc, Em, calc%Tev, C_Norm, Tevf, Ns)
+ CALL Qtrap_Clm (calc, tmpCalc, calcData, calc%Tev, calc%getTemperature(), Em, &
+ C_Norm, Kvtemp, calcData%getNnnsig(), &
+ Isox, 0, Ns, calcData, &
+ Kkkkkk+1)
+!
+! ********* Doppler widths etc for self-indication transmission
IF (Ksitmp.GT.0) THEN
Tevx = Boltzm*Sitemp
ELSE
- Tevx = Tev
+ Tevx = calc%Tev
END IF
-C
-C ********* Broaden self-indication transmission separately
+!
+! ********* Broaden self-indication transmission separately
IF (Ksindi.GT.0 .AND. Kcros.EQ.8) THEN
- CALL Getsab (Beta, Osc_Wts, Osc_Enx, Osc_Snth,
- * Osc_Coth, Bex, S_Expb, S_Ex, Save, Bs, Ss, Max_T,
- * Em, Tevx, C_Norm, Tevf, Ns)
- CALL Qtrap_Clm (tmpCalc,
- * calcDataSelf,
- * Bs, Ss, Tevx, Sitemp, Em,
- * C_Norm, Iflmsc(Ksitmp), 1, Isox, 0, Ns, calcDataSelf,
- * Kkkkkk+1)
+ CALL Getsab (calc, Em, Tevx, C_Norm, Tevf, Ns)
+ CALL Qtrap_Clm (calc, tmpCalc, &
+ calcDataSelf, &
+ Tevx, Sitemp, Em, &
+ C_Norm, Iflmsc(Ksitmp), 1, Isox, 0, Ns, calcDataSelf,&
+ Kkkkkk+1)
END IF
-C
-C ********* Broaden self-indication transm; but stored in calcData old data originally
+!
+! ********* Broaden self-indication transm; but stored in calcData old data originally
IF (Ksindi.EQ.0 .AND. Ksitmp.GT.0 .AND. Kcros.EQ.8) THEN
- CALL Getsab (Beta, Osc_Wts, Osc_Enx, Osc_Snth,
- * Osc_Coth, Bex, S_Expb, S_Ex, Save, Bs, Ss, Max_T,
- * Em, Tevx, C_Norm, Tevf, Ns)
- CALL Qtrap_Clm (tmpCalc,
- * calcData, Bs, Ss, Tevx, Sitemp, Em,
- * C_Norm, Iflmsc(Ksitmp), 1, Isox,
- * Lllmax+1, Ns, calcDataSelf, Kkkkkk+1)
+ CALL Getsab (calc, Em, Tevx, C_Norm, Tevf, Ns)
+ CALL Qtrap_Clm (calc, tmpCalc, &
+ calcData, Tevx, Sitemp, Em, &
+ C_Norm, Iflmsc(Ksitmp), 1, Isox, &
+ Lllmax+1, Ns, calcDataSelf, Kkkkkk+1)
END IF
-C
+!
40 CONTINUE
Kkkkkk = Kkkkkk + 1
-C
-C
+!
+!
Kkkdat = Kkkdat + 1
-C
-C
+!
+!
70 CONTINUE
-C *** END of do-loop on energy
-C
+! *** END of do-loop on energy
+!
80 CONTINUE
END DO
-C *** end of do-loop on isotopes (nuclides)
-C
-C
+! *** end of do-loop on isotopes (nuclides)
+!
+!
nauxStart = nauxStart + 1
call setAuxGridRowMax(Kkkdat)
IF (Now.NE.0) WRITE (21,99997) Now, Kkkdat*Isomax
- IF (Now.NE.0 .AND. Kdebug.NE.0) WRITE (06,99997) Now,Kkkdat*Isomax
-99997 FORMAT (' No Doppler broadening occured', I8,
- * ' times of a possible', I8)
-C
- call dopplerOption%crystalLattice%updateBroadenedOffset(
- * nauxStart)
- call dopplerOption%crystalLattice%setLength(Kkkkkk)
+ IF (Now.NE.0 .AND. calc%debugOutput) WRITE (06,99997) Now,Kkkdat*Isomax
+99997 FORMAT (' No Doppler broadening occured', I8, ' times of a possible', I8)
+!
+ call calc%updateBroadenedOffset( nauxStart)
+ call calc%setLength(Kkkkkk)
call grid%destroy()
call auxGrid%destroy()
call setAuxGridOffset(nauxStart)
call tmpCalc%destroy()
RETURN
-C
+!
END
+end module mclm3_m
diff --git a/sammy/src/clm/mclm4.f90 b/sammy/src/clm/mclm4.f90
index a4db678152a45cf4c2e566f20ba607b0fc06c593..0277114eb1996c3cd1af15849341bcfb4cdd44fb 100644
--- a/sammy/src/clm/mclm4.f90
+++ b/sammy/src/clm/mclm4.f90
@@ -1,101 +1,107 @@
-C
-C
-C ----------------------------------------------------------------------
-C
- SUBROUTINE Getsab (Beta, Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth,
- * Bex, S_Expb, S_Ex, Save, Bs, Ss, Max_T, Em, Tevx, C_Norm, Tevf,
- * Ns)
-C
-C *** Purpose -- CRYSTAL MODEL CROSS-SECTION CALCULATION
-C *** Input -- Em, etc.
-C *** Output -- Bs, Ss
-C
- use fixedi_m
- use ifwrit_m
- use fixedr_m
- use broad_common_m
- use clm_common_m
- use constn_common_m
- IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-C
- DIMENSION Beta(*), Osc_Wts(*), Osc_Enx(*), Osc_Snth(*),
- * Osc_Coth(*), Bex(*), S_Expb(*), S_Ex(*), Save(Nbeta_Max,*),
- * Bs(*), Ss(*), Max_T(*)
-cx DOUBLE PRECISION Em, Tevx
-C
- DIMENSION Sab(25999)
-cx DOUBLE PRECISION C_Norm1, C_Sumr1, C_Norm, C_Sumr, Recul
+module Getsab_m
+use AllocateFunctions_m
+use CrystalLatticeBroadening_M
+IMPLICIT None
+contains
+!
+!
+! ----------------------------------------------------------------------
+!
+ SUBROUTINE Getsab (calc, Em, Tevx, C_Norm, Tevf, &
+ Ns)
+!
+! *** Purpose -- CRYSTAL MODEL CROSS-SECTION CALCULATION
+! *** Input -- Em, etc.
+! *** Output -- Bs, Ss
+!
+ use fixedr_m, only : Aaawww
+ use mclm5_m
+ use mcml6_m
+ use mclm7_m
+ use constn_common_m, only : Aneutr
+!
+ class(CrystalLatticeBroadening)::calc
+ real(Kind=8)::Em, Tevx, Tevf, C_Norm
+!x DOUBLE PRECISION Em, Tevx
+!
+!x DOUBLE PRECISION C_Norm1, C_Sumr1, C_Norm, C_Sumr, Recul
INTEGER K, Ik, Ns
-cx DOUBLE PRECISION Zero, One, Five
- PARAMETER (Zero=0.0d0, One=1.0d0, Five=5.0d0)
-C
-C
-C EFFECTIVE TEMPERATURE CALCULATION FOR CREATION OF Beta GRID
-C SCATTERING LaW CALCULATION
-C Em is in eV; Beta is unitless
-C
-C
-C *** Temporarily assume (1) neutrons only, (2) one isotope only
-C *** Eventually need to fix this!
+!x DOUBLE PRECISION Zero, One, Five
+ real(kind=8),PARAMETER::Zero=0.0d0, One=1.0d0, Five=5.0d0
+ real(kind=8)::Alpha, C_Norm1, C_Sumr, C_Sumr1, Deltab
+ integer::Ndmax
+ real(kind=8)::Osc0, Osc1, Recul, Tra0, Tra1, Tbarx
+!
+!
+! EFFECTIVE TEMPERATURE CALCULATION FOR CREATION OF Beta GRID
+! SCATTERING LaW CALCULATION
+! Em is in eV; Beta is unitless
+!
+!
+! *** Temporarily assume (1) neutrons only, (2) one isotope only
+! *** Eventually need to fix this!
Recul = Em*(Aneutr/(Aaawww+Aneutr))
Alpha = Recul/Tevx
-C
-C
-C *** Continuous part of distribution
- CALL Cconti (Save, Sab, Beta, Alpha, C_Norm, C_Sumr, Max_T)
-C
-C
+ call allocate_real_data(calc%Sab, calc%Nbeta)
+!
+!
+! *** Continuous part of distribution
+ CALL Cconti (calc, Alpha, C_Norm, C_Sumr)
+!
+!
Tra0 = C_Norm
Tra1 = C_Sumr
- Tbarx = Tbar
-C *** Translational part, if any
- IF (Twt.NE.Zero) THEN
+ Tbarx = calc%Tbar
+! *** Translational part, if any
+ IF (calc%Twt.NE.Zero) THEN
Ndmax = 2000
- CALL Trans (Beta, Sab, Alpha, Deltab, Tra0, Tra1, Ndmax)
-C *** UPDATE EFfECTIVE TEMPERATURE (which is used only in Discre)
- Tevf = (Tbeta*Tevf+Twt*Tevx)/(Tbeta+Twt)
+ CALL Trans (calc, Alpha, Deltab, Tra0, Tra1, Ndmax)
+! *** UPDATE EFfECTIVE TEMPERATURE (which is used only in Discre)
+ Tevf = (calc%Tbeta*Tevf+calc%Twt*Tevx)/(calc%Tbeta+calc%Twt)
Tbarx = Tevf/Tevx
END IF
-C
-C
+!
+!
Osc0 = Tra0
Osc1 = Tra1
-C *** Discrete oscillators, if any
- IF (N_Osc.NE.0) THEN
- CALL Discre (Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth, Sab, Beta,
- * Bex, S_Expb, S_Ex, Alpha, Tbarx, Osc0, Osc1)
+! *** Discrete oscillators, if any
+ IF (calc%N_Osc.NE.0) THEN
+ CALL Discre (calc, Alpha, Tbarx, Osc0, Osc1)
END IF
-C
-C
+!
+!
C_Norm1 = dABS(C_Norm-One)
C_Sumr1 = dABS(C_Sumr-One)
- IF (C_Norm1.GT.Eps .OR. C_Sumr1.GT.Eps) THEN
- Nphon = Nphon + 5
+ IF (C_Norm1.GT.calc%Eps .OR. C_Sumr1.GT.calc%Eps) THEN
+ calc%Nphon = calc%Nphon + 5
END IF
-C
-C
-C Dimensioning of Alpha and Beta values
-C Centering of Beta in Recul=<Beta>=Em/(Awr+1)
-C
+!
+!
+! Dimensioning of Alpha and Beta values
+! Centering of Beta in Recul=<Beta>=Em/(Awr+1)
+!
K=0
- DO IK=Nbeta,1,-1
+ call allocate_real_data(calc%Bs, 2*calc%Nbeta)
+ call allocate_real_data(calc%Ss, 2*calc%Nbeta)
+ DO IK=calc%Nbeta,1,-1
K = K + 1
- Bs(K) = Recul - Beta(IK)*Tevx
- Ss(K) = Sab(IK)/Tevx
+ calc%Bs(K) = Recul - calc%Beta(IK)*Tevx
+ calc%Ss(K) = calc%Sab(IK)/Tevx
END DO
- DO IK=2,Nbeta
+ DO IK=2,calc%Nbeta
K = K + 1
- Bs(K) = Recul + Beta(IK)*Tevx
- Ss(K) = Sab(IK)*dEXP(-Beta(IK))/Tevx
+ calc%Bs(K) = Recul + calc%Beta(IK)*Tevx
+ calc%Ss(K) = calc%Sab(IK)*dEXP(-calc%Beta(IK))/Tevx
END DO
-C
-C
-C ATTENTION!!! Betas ARE NOW CENTERED AROUND Recul!!!
-C THE ELASTIC SCATTERING CONTRIBUTION
-C S(A,B) = Smult(A,B) + Delta(B)
-C IS ADDED DIRECTLY TO THE CROSS-SECTION OUT OF CONVOLUTION
-C Sigma = SigmaCONVOL + Sigma*EXP(-W*A)
-C
+!
+!
+! ATTENTION!!! Betas ARE NOW CENTERED AROUND Recul!!!
+! THE ELASTIC SCATTERING CONTRIBUTION
+! S(A,B) = Smult(A,B) + Delta(B)
+! IS ADDED DIRECTLY TO THE CROSS-SECTION OUT OF CONVOLUTION
+! Sigma = SigmaCONVOL + Sigma*EXP(-W*A)
+!
NS = K
RETURN
END
+end module Getsab_m
diff --git a/sammy/src/clm/mclm5.f90 b/sammy/src/clm/mclm5.f90
index a9ab9eb27ca3afa26e27e4651ef89a867e04d57f..bbf1ebc3c5edf051babb57ad62a2fc7bf77da9d4 100644
--- a/sammy/src/clm/mclm5.f90
+++ b/sammy/src/clm/mclm5.f90
@@ -1,49 +1,57 @@
-C
-C
-C -------------------------------------------------------------------
-C
- SUBROUTINE Cconti (Save, Sab, Beta, Alpha, C_Norm, C_Sumr, Max_T)
-C
-C
-C Phonon expansion method for scattering low calculation
-C Continuous frequency distribution case
-C Calculation of two first moments of S(Alpha,Beta)
-C
-C Main routine for calculating S(Alpha,Beta) at temperature "Tev"
-C for continuous phonon frequency distributions.
-C Called by Sigcri
-C
- use fixedr_m
- use clm_common_m
- use constn_common_m
- IMPLICIT DOUBLE PRECISION (A-H,O-Z)
- DIMENSION Save(Nbeta_Max,*), Sab(*), Beta(*), Max_T(*)
-cx DOUBLE PRECISION Alpha, C_Sumr, C_Norm,
-cx & Exx, Be, St, Add, Alw, Alp, Ssct, Ff2, Ff1, Ff1m, Ff2m,
-cx & Bem, Alf0, Sab_Max, X
+module mclm5_m
+use AllocateFunctions_m
+use constn_common_m, only : Pi
+use CrystalLatticeBroadening_M
+IMPLICIT None
+public Cconti
+contains
+!
+!
+! -------------------------------------------------------------------
+!
+ SUBROUTINE Cconti (calc, Alpha, C_Norm, C_Sumr)
+!
+!
+! Phonon expansion method for scattering low calculation
+! Continuous frequency distribution case
+! Calculation of two first moments of S(Alpha,Beta)
+!
+! Main routine for calculating S(Alpha,Beta) at temperature "Tev"
+! for continuous phonon frequency distributions.
+! Called by Sigcri
+!
+ class(CrystalLatticeBroadening)::calc
+ real(kind=8)::Alpha, C_Norm, C_Sumr
+!x DOUBLE PRECISION Alpha, C_Sumr, C_Norm,
+!x & Exx, Be, St, Add, Alw, Alp, Ssct, Ff2, Ff1, Ff1m, Ff2m,
+!x & Bem, Alf0, Sab_Max, X
INTEGER K, N, Ksmall, Kountr
-cx DOUBLE PRECISION Zero, Half, One, Four, Small, Pi
- PARAMETER (Zero=0.0d0, Half =0.5d0, One=1.0d0, Four=4.0d0,
- & Small=1.d-20)
-C
-C
+!x DOUBLE PRECISION Zero, Half, One, Four, Small, Pi
+ real(kind=8),PARAMETER::Zero=0.0d0, Half =0.5d0, One=1.0d0, Four=4.0d0, &
+ Small=1.d-20
+ real(kind=8)::Add, Alf0, Alp, Alw, Be, Bem, Exx, Ff1, Ff1m, Ff2, Ff2m
+ integer::Ntime
+ real(kind=8)::Sab_Max, Ssct, St, X
+!
+!
Sab_Max = Zero
- Alf0 = Alpha*F0
+ Alf0 = Alpha*calc%F0
Exx = Alf0*dEXP(-Alf0)
- DO K=1,Nbeta
- Add = Save(K,1)*Exx
+ DO K=1,calc%Nbeta
+ Add = calc%Save(K,1)*Exx
IF (Add.GT.Sab_Max) Sab_Max = Add
- Sab(K) = Add
+ calc%Sab(K) = Add
END DO
-C
-C
-C *** Add other terms in the phonon expansion sum
-C
- DO K=1,Nbeta
- Max_T(K) = 2
+!
+!
+! *** Add other terms in the phonon expansion sum
+!
+ call reallocate_integer_data(calc%Max_T, calc%Nbeta, 0)
+ DO K=1,calc%Nbeta
+ calc%Max_T(K) = 2
END DO
Ksmall = 1
- DO N=2,Nphon
+ DO N=2,calc%Nphon
IF (Ksmall.EQ.0) THEN
GO TO 20
END IF
@@ -53,20 +61,20 @@ C
ELSE IF (Exx.NE.Zero) THEN
Ksmall = 0
Kountr = 0
- DO K=1,Nbeta-1
+ DO K=1,calc%Nbeta-1
IF (Kountr.GT.10) GO TO 10
- St = Save(K,N)
+ St = calc%Save(K,N)
Add = St*Exx
- IF (Add+Sab(K).NE.Sab(K)) THEN
+ IF (Add+calc%Sab(K).NE.calc%Sab(K)) THEN
Ksmall = 1
- Sab(K) = Sab(K) + Add
- IF (N.EQ.Nphon .AND. Sab(K)+Sab_Max.NE.Sab_Max) THEN
- IF (Add.GT.Sab(K)*0.001d0 .AND. Max_T(K).EQ.2) THEN
-C *** IF (still adding a significant amount) THEN
-C *** switch to using asymptotic form
-C *** [NML: but why bother? we already know Sab more
-C *** accurately than this.]
- Max_T(K) = 1
+ calc%Sab(K) = calc%Sab(K) + Add
+ IF (N.EQ.calc%Nphon .AND. calc%Sab(K)+Sab_Max.NE.Sab_Max) THEN
+ IF (Add.GT.calc%Sab(K)*0.001d0 .AND. calc%Max_T(K).EQ.2) THEN
+! *** IF (still adding a significant amount) THEN
+! *** switch to using asymptotic form
+! *** [NML: but why bother? we already know Sab more
+! *** accurately than this.]
+ calc%Max_T(K) = 1
END IF
END IF
ELSE
@@ -77,53 +85,53 @@ C *** accurately than this.]
END IF
END DO
20 CONTINUE
-C
-C
-C *** ------
-C *** Start of SCT range for each Beta
-C
- DO K=2,Nbeta
- IF (Max_T(K).GT.Max_T(K-1)) Max_T(K) = Max_T(K-1)
+!
+!
+! *** ------
+! *** Start of SCT range for each Beta
+!
+ DO K=2,calc%Nbeta
+ IF (calc%Max_T(K).GT.calc%Max_T(K-1)) calc%Max_T(K) = calc%Max_T(K-1)
END DO
-C
-C
-C *** ------
-C *** Check the moments of S(Alpha,Beta); Also redefine Sab if Max_T=1
-C
+!
+!
+! *** ------
+! *** Check the moments of S(Alpha,Beta); Also redefine Sab if Max_T=1
+!
Ntime = 0
- Alw = ALpha*Tbeta
- Alp = Alw*Tbar*Four
+ Alw = ALpha*calc%Tbeta
+ Alp = Alw*calc%Tbar*Four
X = dSQRT(Pi*Alp)
K = 1
- Bem = Beta(K)
- Ff2m = Sab(K)
- Ff1m = Sab(K)*dEXP(-Bem)
+ Bem = calc%Beta(K)
+ Ff2m = calc%Sab(K)
+ Ff1m = calc%Sab(K)*dEXP(-Bem)
C_Norm = Zero
C_Sumr = Zero
- IF (Max_T(K).EQ.1) THEN
-C (asymptotic redefinition of Sab)
+ IF (calc%Max_T(K).EQ.1) THEN
+! (asymptotic redefinition of Sab)
Exx = -(Alw-Bem)**2/Alp
Ssct = Zero
IF (Exx.GT.-45.0d0) Ssct = dEXP(Exx)/X
- Sab(K) = Ssct
+ calc%Sab(K) = Ssct
END IF
- DO K=2,Nbeta
- Be = Beta(K)
- IF (Max_T(K).EQ.1) THEN
-C (asymptotic redefinition of Sab)
+ DO K=2,calc%Nbeta
+ Be = calc%Beta(K)
+ IF (calc%Max_T(K).EQ.1) THEN
+! (asymptotic redefinition of Sab)
Exx = -(Alw-Be)**2/Alp
Ssct = Zero
IF (Exx.GT.-45.0d0) Ssct = dEXP(Exx)/X
- Sab(K) = Ssct
+ calc%Sab(K) = Ssct
END IF
- IF (Sab(K).EQ.Zero) THEN
+ IF (calc%Sab(K).EQ.Zero) THEN
Ntime = Ntime + 1
IF (Ntime.GT.20) GO TO 30
ELSE
Ntime = 0
END IF
- Ff2 = Sab(K)
- Ff1 = Sab(K)*dEXP(-Be)
+ Ff2 = calc%Sab(K)
+ Ff1 = calc%Sab(K)*dEXP(-Be)
C_Norm =C_Norm +Half*(Be-Bem)*(Ff2m +Ff2 +Ff1m +Ff1 )
C_Sumr =C_Sumr +Half*(Be-Bem)*(Ff2m*Bem+Ff2*Be-Ff1m*Bem-Ff1*Be)
Ff1m = Ff1
@@ -131,16 +139,17 @@ C (asymptotic redefinition of Sab)
Bem = Be
END DO
30 CONTINUE
-C
-C
-C *** ------
-C *** Add the Del_Beta contribution to the norm
-C
- C_Norm = C_Norm + dEXP(-F0*ALpha)
- C_Sumr = C_Sumr/ALpha/Tbeta
-C
-C *** ------
-C *** Finished with continuous distribution
-C
+!
+!
+! *** ------
+! *** Add the Del_Beta contribution to the norm
+!
+ C_Norm = C_Norm + dEXP(-calc%F0*ALpha)
+ C_Sumr = C_Sumr/ALpha/calc%Tbeta
+!
+! *** ------
+! *** Finished with continuous distribution
+!
RETURN
END
+end module mclm5_m
diff --git a/sammy/src/clm/mclm6.f90 b/sammy/src/clm/mclm6.f90
index c44e2d54777a4d68dd1ffc0f9cc561e9e8c28162..f0b7eec906888ddd0248c2719477a0f32eefbca0 100644
--- a/sammy/src/clm/mclm6.f90
+++ b/sammy/src/clm/mclm6.f90
@@ -1,79 +1,82 @@
-C
-C
-C ----------------------------------------------------------------------
-C
- SUBROUTINE Trans (Beta, Sab, Alpha, Deltab, Tra0, Tra1, Ndmax)
-C
-C *** PURPOSE -- Control the addition of a translational contribution
-C *** to a continuous S(Alpha,Beta). The translational
-C *** component can be either diffusion, or a free gas.
-C *** The values of the input S(Alpha,Beta) for the
-C *** convoltion are obtained by interpolation.
-C *** CALLED BY -- Sigcri.
-C *** USES -- Stable, Sbfill, Terps.
-C *** Input -- Beta(), Sab()
-C *** -- Alpha, F0, Deltab, Twt, C_Trans, Tbeta
-C *** Output -- Sab(), Tra0, Tra1
-C *** Dummies -- Betax(), Sabx(), Sd(), Sb()
-C
- use clm_common_m
- IMPLICIT NONE
- DOUBLE PRECISION Beta(*), Sab(*)
+module mcml6_m
+use CrystalLatticeBroadening_M
+IMPLICIT NONE
+
+public Trans
+contains
+!
+!
+! ----------------------------------------------------------------------
+!
+ SUBROUTINE Trans (calc, Alpha, Deltab, Tra0, Tra1, Ndmax)
+!
+! *** PURPOSE -- Control the addition of a translational contribution
+! *** to a continuous S(Alpha,Beta). The translational
+! *** component can be either diffusion, or a free gas.
+! *** The values of the input S(Alpha,Beta) for the
+! *** convoltion are obtained by interpolation.
+! *** CALLED BY -- Sigcri.
+! *** USES -- Stable, Sbfill, Terps.
+! *** Input -- Beta(), Sab()
+! *** -- Alpha, F0, Deltab, Twt, C_Trans, Tbeta
+! *** Output -- Sab(), Tra0, Tra1
+! *** Dummies -- Betax(), Sabx(), Sd(), Sb()
+!
+ class(CrystalLatticeBroadening)::calc
DOUBLE PRECISION Alpha, Deltab, Tra0, Tra1
INTEGER Ndmax
-C
+!
DOUBLE PRECISION Betax(2000), Sabx(2000), Sd(2000), Sb(2000)
- DOUBLE PRECISION Ded, Deb, Delta, S, St, Be, Bb, Ff1, Ff2,
- & Ff1_Km1, Ff2_Km1, Be_Km1, AL, F, Ebe, Ebe_Km1, Tcal
- DOUBLE PRECISION Terps
+ DOUBLE PRECISION Ded, Deb, Delta, S, St, Be, Bb, Ff1, Ff2, &
+ Ff1_Km1, Ff2_Km1, Be_Km1, AL, F, Ebe, Ebe_Km1, Tcal
INTEGER Nsd, Iprt, Nu, I, Ibeta, Nbt
DOUBLE PRECISION Zero, Half, One, Two, Small
- PARAMETER (Zero=0.0d0, Half=0.5d0, One=1.0d0, Two=2.0d0,
- * Small=1.d-20)
-C
-C
+ PARAMETER (Zero=0.0d0, Half=0.5d0, One=1.0d0, Two=2.0d0, &
+ Small=1.d-20)
+!
+!
AL = Alpha
Iprt = 1
-C
-C
-C *** Choose Beta interval for convolution
-C
- Tcal = Twt*C_Trans*AL
- Ded = 0.4*(Tcal)/dSQRT(One+1.42*(Tcal)*C_Trans)
+!
+!
+! *** Choose Beta interval for convolution
+!
+ Tcal = calc%Twt*calc%C_Trans*AL
+ Ded = 0.4*(Tcal)/dSQRT(One+1.42*(Tcal)*calc%C_Trans)
IF (Ded.EQ.Zero ) Ded = 1000.0d0
- IF (Ded.EQ.1000.) Ded = 0.2d0*dSQRT(Twt*AL)
+ IF (Ded.EQ.1000.) Ded = 0.2d0*dSQRT(calc%Twt*AL)
Deb = 10.0d0 *AL*Deltab
Delta = dMIN1 (Deb, Ded)
Nu = 1
-C
-C
-C *** Make table of S-diffusion or S-free on this interval
-C
- CALL Stable (Sabx, Sd, Alpha, Delta, Twt, C_Trans,
- * Nu, Nsd, Ndmax)
-C
-C
+!
+!
+! *** Make table of S-diffusion or S-free on this interval
+!
+ CALL Stable (Sabx, Sd, Alpha, Delta, calc%Twt, calc%C_Trans, &
+ Nu, Nsd, Ndmax)
+!
+!
IF (Nsd.GT.1) THEN
-C
-C *** Copy original Ss(-Beta) to a temporary array
- DO I=1,Nbeta
- Betax(I) = Beta(I)
- Sabx (I) = Sab (I)
+!
+! *** Copy original Ss(-Beta) to a temporary array
+ DO I=1,calc%Nbeta
+ Betax(I) = calc%Beta(I)
+ Sabx (I) = calc%Sab (I)
END DO
-C
-C *** Loop over Beta values
- DO Ibeta=1,Nbeta
+!
+! *** Loop over Beta values
+ DO Ibeta=1,calc%Nbeta
S = Zero
Be = Betax(Ibeta)
-C
-C *** Prepare table of continuous Ss on new interval
+!
+! *** Prepare table of continuous Ss on new interval
Nbt = Nsd
- CALL Sbfill (Sb, Sabx, Betax, Delta, Be, Nbt, Nbeta, Ibeta)
-C
-C *** Convolve S-transport with S-continuous
+ CALL Sbfill (Sb, Sabx, Betax, Delta, Be, Nbt, calc%Nbeta, Ibeta)
+!
+! *** Convolve S-transport with S-continuous
DO I=1,Nbt
F = Two*(MOD(I-1,2)+1)
-C *** If I=even, F=4. If I=odd, F=2.
+! *** If I=even, F=4. If I=odd, F=2.
IF (I.EQ.1 .OR. I.EQ.Nbt) F = One
Bb = (I-1)*Delta
S = S + F*Sd(I)* ( Sb(Nbt+I-1) + Sb(Nbt-I+1)*dEXP(-Bb) )
@@ -81,27 +84,27 @@ C *** If I=even, F=4. If I=odd, F=2.
S = S*Delta/3.0D0
IF (S.LT.Small) S = Zero
St = Terps (Sd, Nbt, Delta, Be)
- IF (St.GT.Zero) S = S + dEXP(-AL*F0)*St
-C
-C *** Store results and continue Beta loop
- Sab(Ibeta) = S
+ IF (St.GT.Zero) S = S + dEXP(-AL*calc%F0)*St
+!
+! *** Store results and continue Beta loop
+ calc%Sab(Ibeta) = S
END DO
-C
-C
-C *** Check moments of calculated S(Alpha, Beta).
-C
+!
+!
+! *** Check moments of calculated S(Alpha, Beta).
+!
IF (Iprt.EQ.1) THEN
Tra0 = Zero
Tra1 = Zero
Be_Km1 = Betax(1)
Ebe_Km1 = dEXP(-Be_Km1)
- Ff1_Km1 = Sab(1) * (One+Ebe_Km1)
- Ff2_Km1 = Sab(1) * (One-Ebe_Km1) * Be_Km1
- DO Ibeta=2,Nbeta
+ Ff1_Km1 = calc%Sab(1) * (One+Ebe_Km1)
+ Ff2_Km1 = calc%Sab(1) * (One-Ebe_Km1) * Be_Km1
+ DO Ibeta=2,calc%Nbeta
Be = Betax(Ibeta)
Ebe = dEXP(-Be)
- Ff1 = Sab(Ibeta) * (One+Ebe)
- Ff2 = Sab(Ibeta) * (One-Ebe) * Be
+ Ff1 = calc%Sab(Ibeta) * (One+Ebe)
+ Ff2 = calc%Sab(Ibeta) * (One-Ebe) * Be
Tra0 = Tra0 + Half*(Be-Be_Km1) * ( Ff1_Km1 + Ff1 )
Tra1 = Tra1 + Half*(Be-Be_Km1) * ( Ff2_Km1 + Ff2 )
Ff1_Km1 = Ff1
@@ -109,131 +112,129 @@ C
Be_Km1 = Be
Ebe_Km1 = Ebe
END DO
- Tra1 = Tra1/AL/(Tbeta+Twt)
+ Tra1 = Tra1/AL/(calc%Tbeta+calc%Twt)
END IF
-C
+!
END IF
-C
+!
RETURN
END
-C
-C
-C ----------------------------------------------------------------------
-C
+!
+!
+! ----------------------------------------------------------------------
+!
SUBROUTINE Sbfill (Sb, Sabx, Betax, Delta, Be, Nbt, Nbeta, Ibeta)
-C
-C *** PURPOSE -- Generate S(Beta) on a new energy grid for convolution
-C *** with a diffusion or free-gas shape. Interpolation
-C *** is used.
-C *** For translational cases only.
-C *** CALLED BY -- Trans.
-C
- IMPLICIT NONE
+!
+! *** PURPOSE -- Generate S(Beta) on a new energy grid for convolution
+! *** with a diffusion or free-gas shape. Interpolation
+! *** is used.
+! *** For translational cases only.
+! *** CALLED BY -- Trans.
+!
INTEGER Nbt, Nbeta, Ibeta
DOUBLE PRECISION Sb(*), Sabx(*), Betax(Nbeta), Delta, Be
-C *** note that Ibeta is not used, but maybe should be ?!?
-C
+! *** note that Ibeta is not used, but maybe should be ?!?
+!
DOUBLE PRECISION Bmin, Bmax, Bet, B, St, Stm
INTEGER J, I
DOUBLE PRECISION Zero, Hundth
PARAMETER (Zero=0.0d0, Hundth=0.01d0)
-C
+!
Bmin = - Be - (Nbt-1)*Delta
Bmax = - Be + (Nbt-1)*Delta + Hundth*Delta
J = Nbeta
I = 0
Bet = Bmin
-C
+!
10 CONTINUE
-C
+!
I = I + 1
B = dABS(Bet)
-C
+!
20 CONTINUE
-C
-C *** SEARCH FOR CORRECT Beta RANGE
+!
+! *** SEARCH FOR CORRECT Beta RANGE
IF (B.GT.Betax(J)) THEN
IF (J.EQ.Nbeta) THEN
Sb(I) = Zero
GO TO 30
END IF
-C MOVE UP
+! MOVE UP
J = J + 1
GO TO 20
END IF
-C
+!
IF (B.LE.Betax(J-1)) THEN
IF (J.EQ.2) THEN
Sb(I) = Zero
GO TO 30
END IF
-C MOVE DOWN
+! MOVE DOWN
J = J - 1
GO TO 20
END IF
-C
-C
-C INTERPOLATE IN THIS RANGE
-C
+!
+!
+! INTERPOLATE IN THIS RANGE
+!
IF (Sabx(J).LE.Zero) THEN
St = -225.d0
ELSE
St = dLOG(Sabx(J))
END IF
-C
+!
IF (Sabx(J-1).LE.Zero) THEN
Stm = -225.d0
ELSE
Stm = dLOG(Sabx(J-1))
END IF
-C
+!
Sb(I) = St + (B-Betax(J))*(Stm-St) / (Betax(J-1)-Betax(J))
IF (Bet.GT.Zero) Sb(I) = Sb(I) - Bet
Sb(I) = dEXP(Sb(I))
-C
+!
30 CONTINUE
Bet = Bet + Delta
IF (Bet.LE.Bmax) GO TO 10
-C
-C *** end of loop on Bet
-C
+!
+! *** end of loop on Bet
+!
RETURN
END
-C
-C
-C ----------------------------------------------------------------------
-C
- SUBROUTINE Stable (Sabx, Sd, Alpha, Delta, Twt, C_Trans,
- * Nu, Nsd, Ndmax)
-C
-C *** PURPOSE -- Set up a table of S-Diffusion or S-Free in the array
-C *** Sd, evaluated at intervals Delta determined by Trans.
-C *** Tabulation is continued until
-C *** Sd(J) is less than 1.0d-7 * Sd(1), or Nsd is 1999
-C *** Nsd is always odd for use with Simpson's rule.
-C *** Array Sabx is used for temporary storage of Beta values.
-C *** Stable returns -Beta side of asymmetric S(Alpha,Beta).
-C *** CALLED BY- Trans.
-C *** USES -- Besk1.
-C
- IMPLICIT NONE
+!
+!
+! ----------------------------------------------------------------------
+!
+ SUBROUTINE Stable (Sabx, Sd, Alpha, Delta, Twt, C_Trans, &
+ Nu, Nsd, Ndmax)
+!
+! *** PURPOSE -- Set up a table of S-Diffusion or S-Free in the array
+! *** Sd, evaluated at intervals Delta determined by Trans.
+! *** Tabulation is continued until
+! *** Sd(J) is less than 1.0d-7 * Sd(1), or Nsd is 1999
+! *** Nsd is always odd for use with Simpson's rule.
+! *** Array Sabx is used for temporary storage of Beta values.
+! *** Stable returns -Beta side of asymmetric S(Alpha,Beta).
+! *** CALLED BY- Trans.
+! *** USES -- Besk1.
+!
INTEGER Nu, Nsd, Ndmax
DOUBLE PRECISION Sabx(Ndmax), Sd(Ndmax)
DOUBLE PRECISION Alpha, Delta, Twt, C_Trans
INTEGER Icheck, J, I
- DOUBLE PRECISION C2, C3, C4, C8, Be, C6, C7, Wal, Check0, Check1,
- & Eps, D, C5, Ex, Sfree, F, Bb, Besk1
+ DOUBLE PRECISION C2, C3, C4, C8, Be, C6, C7, Wal, Check0, Check1, &
+ Eps, D, C5, Ex, Sfree, F, Bb
DOUBLE PRECISION Pi, Zero, Quarter, Half, One, Two, Three, Four
- PARAMETER (PI=3.141592653589793238462643d0, Zero=0.0d0,
- & Quarter=0.25d0, Half=0.5d0, One=1.0d0, Two=2.0d0, Three=3.0d0,
- & Four=4.0d0)
-C
+ PARAMETER (PI=3.141592653589793238462643d0, Zero=0.0d0, &
+ Quarter=0.25d0, Half=0.5d0, One=1.0d0, Two=2.0d0, Three=3.0d0, &
+ Four=4.0d0)
+!
Icheck = 0
Eps = 1.E-7
-C
+!
IF (C_Trans.NE.Zero) THEN
-C
-C *** DIFFUSION BRANCH
+!
+! *** DIFFUSION BRANCH
D = Twt*C_Trans
C2 = dSQRT(C_Trans**2+Quarter)
C3 = Two*D*Alpha
@@ -245,33 +246,33 @@ C *** DIFFUSION BRANCH
100 CONTINUE
C6 = dSQRT(Be*Be+C4)
C7 = C6*C2
-C
+!
IF (C7.LE.One) THEN
C5 = C8*dEXP(C3+Half*Be)
ELSE
-C ELSE IF (C7.GT.One) THEN
+! ELSE IF (C7.GT.One) THEN
C5 = Zero
Ex = C3 - C7 + Half*Be
C5 = C8*dEXP(Ex)
END IF
-C
+!
Sd(J) = C5*Besk1(C7)/C6
Sabx(J) = Be
Be = Be + Delta
J = J + 1
-C
+!
IF (MOD(J,2).NE.0) THEN
GO TO 100
ELSE IF (J.LT.2000) THEN
IF (Eps*Sd(1).LT.Sd(J-1)) GO TO 100
END IF
-C
+!
J = J - 1
Nsd = J
-C
+!
ELSE
-C
-C *** FREE-GAS BRANCH
+!
+! *** FREE-GAS BRANCH
Be = Zero
J = 1
Wal = Twt*Alpha
@@ -282,22 +283,22 @@ C *** FREE-GAS BRANCH
Sabx(J) = Be
Be = Be + Delta
J = J + 1
-C
+!
IF (MOD(J,2).NE.0) THEN
GO TO 170
ELSE IF (J.LT.2000) THEN
IF (Eps*Sd(1).LT.Sd(J-1)) GO TO 170
END IF
-C
+!
J = J - 1
Nsd = J
-C
+!
END IF
-C
-C
-C
+!
+!
+!
IF (Icheck.NE.0) THEN
-C *** Check the moments of the distribution
+! *** Check the moments of the distribution
Check0 = Zero
Check1 = Zero
DO I=1,Nsd
@@ -312,93 +313,92 @@ C *** Check the moments of the distribution
Check0 = Check0*Delta/Three
Check1 = Check1*Delta/Three
Check1 = Check1/(Alpha*Twt)
-C *** Huh? Must return these values or print them or something?
-C *** Otherwise they are just ignored
+! *** Huh? Must return these values or print them or something?
+! *** Otherwise they are just ignored
END IF
-C
+!
RETURN
END
-C
-C
-C -------------------------------------------------------------------
-C
- FUNCTION Besk1 (X)
-C
-C *** PURPOSE -- Compute modified Bessel Function K1
-C *** The exponential part for X>1 is omitted (See Stable).
-C CALLED BY- Stable.
-C
- IMPLICIT NONE
- DOUBLE PRECISION X, V, U, BI1, BI3, Besk1
-C
+!
+!
+! -------------------------------------------------------------------
+!
+ real(kind=8) FUNCTION Besk1 (X)
+!
+! *** PURPOSE -- Compute modified Bessel Function K1
+! *** The exponential part for X>1 is omitted (See Stable).
+! CALLED BY- Stable.
+!
+ DOUBLE PRECISION X, V, U, BI1, BI3
+!
IF (X.LE.1.0) THEN
V = 0.125d0 * X
U = V*V
- Bi1 = ((((((((( 0.442850424d0 *U + 0.584115288d0)*U
- & + 6.070134559d0)*U + 17.864913364d0)*U
- & + 48.858995315d0)*U + 90.924600045d0)*U
- & + 113.795967431d0)*U + 85.331474517d0)*U
- & + 32.00008698 d0)*U + 3.999998802d0 )*V
- Bi3 = ((((((((( 1.304923514d0 *U + 1.47785657 d0)*U
- & + 16.402802501d0)*U + 44.732901977d0)*U
- & + 115.837493464d0)*U +198.437197312d0)*U
- & + 222.869709703d0)*U +142.216613971d0)*U
- & + 40.000262262d0)*U + 1.999996391d0)
+ Bi1 = ((((((((( 0.442850424d0 *U + 0.584115288d0)*U &
+ + 6.070134559d0)*U + 17.864913364d0)*U &
+ + 48.858995315d0)*U + 90.924600045d0)*U &
+ + 113.795967431d0)*U + 85.331474517d0)*U &
+ + 32.00008698d0)*U + 3.999998802d0 )*V
+ Bi3 = ((((((((( 1.304923514d0 *U + 1.47785657d0)*U &
+ + 16.402802501d0)*U + 44.732901977d0)*U &
+ + 115.837493464d0)*U +198.437197312d0)*U &
+ + 222.869709703d0)*U +142.216613971d0)*U &
+ + 40.000262262d0)*U + 1.999996391d0)
Besk1 = 1.0d0/X + Bi1*(dLOG(0.5d0*X)+0.5772156649d0) - V*Bi3
-C
+!
ELSE
-C
+!
U = 1.0d0/X
- Bi3 = ((((((((((((- 0.0108241775d0 *U + 0.0788000118d0)*U
- & - 0.2581303765d0)*U + 0.5050238576d0)*U
- & - 0.663229543 d0)*U + 0.6283380681d0)*U
- & - 0.4594342117d0)*U + 0.2847618149d0)*U
- & - 0.1736431637d0)*U + 0.1280426636d0)*U
- & - 0.1468582957d0)*U + 0.4699927013d0)*U
- & + 1.2533141373d0)
+ Bi3 = ((((((((((((- 0.0108241775d0 *U + 0.0788000118d0)*U &
+ - 0.2581303765d0)*U + 0.5050238576d0)*U &
+ - 0.663229543d0)*U + 0.6283380681d0)*U &
+ - 0.4594342117d0)*U + 0.2847618149d0)*U &
+ - 0.1736431637d0)*U + 0.1280426636d0)*U &
+ - 0.1468582957d0)*U + 0.4699927013d0)*U &
+ + 1.2533141373d0)
Besk1 = dSQRT(U)*Bi3
END IF
RETURN
END
-C
-C
-C ----------------------------------------------------------------------
-C
+!
+!
+! ----------------------------------------------------------------------
+!
DOUBLE PRECISION FUNCTION Terps (Sd, Nsd, Delta, Be)
-C
-C *** PURPOSE -- Interpolate for Beta=Be in table of S(Alpha,Beta).
-C *** USED IN -- Trans.
-C
- IMPLICIT NONE
+!
+! *** PURPOSE -- Interpolate for Beta=Be in table of S(Alpha,Beta).
+! *** USED IN -- Trans.
+!
INTEGER Nsd, I
DOUBLE PRECISION Sd(Nsd), Bt, Btp, ST, Delta, Stp, Stt, Be
DOUBLE PRECISION Zero
PARAMETER (Zero=0.0d0)
-C
+!
I = Be/Delta
-C
+!
IF (I.GE.Nsd-1) THEN
Terps = Zero
RETURN
END IF
-C
+!
Bt = I*Delta
Btp = Bt + Delta
I = I + 1
-C
+!
IF (Sd(I).LE.Zero) THEN
St = -225.d0
ELSE
St = dLOG(Sd(I))
END IF
-C
+!
IF (SD(I+1).LE.Zero) THEN
Stp = -225.d0
ELSE
Stp = dLOG(Sd(I+1))
END IF
-C
+!
Stt = St + (Be-Bt)*(Stp-ST)/(Btp-Bt)
Terps = dEXP(Stt)
RETURN
END
+end module mcml6_m
diff --git a/sammy/src/clm/mclm7.f90 b/sammy/src/clm/mclm7.f90
index d63b954fd39e36db48fa41cb58e9ca90b06c72dd..8b8cce6a04b906163ea956b98990e5757df1d1ec 100644
--- a/sammy/src/clm/mclm7.f90
+++ b/sammy/src/clm/mclm7.f90
@@ -1,80 +1,79 @@
-C
-C
-C ----------------------------------------------------------------------
-C
- SUBROUTINE Discre (Osc_Wts, Osc_Enx, Osc_Snth, Osc_Coth, Sab,
- & Beta, Bex, S_Expb, S_Ex, Alpha, Tbarx, Osc0, Osc1)
-C
-C
-C Scattering low calculation for discrete frequency distribution
-C
-C *** PURPOSE -- Control the convolution of discrete oscillators with
-C *** the continuous S(Alpha,Beta) computed in Contin.
-C *** Called by -- Sigcri
-C *** Uses -- Bfact, Exts, Sint
-C *** Input -- Osc_Wts(), Osc_Enx(), Osc_Snth(), Osc_Coth(),
-C *** Beta(), Bex()
-C *** -- Alpha, Dwpix, Tbeta, Tbarx
-C *** Output -- Sab()
-C *** -- Osc0, Osc1
-C *** Dummies -- Bplus(), Bminus(), S_Expb(), S_Ex(),
-C *** Bes(), Wts(), Ben(), Wtn()
-C
-C
- use clm_common_m
- use constn_common_m
- IMPLICIT DOUBLE PRECISION (A-H,O-Z)
- DOUBLE PRECISION Osc_Wts(*), Osc_Enx(*), Osc_Snth(*), Osc_Coth(*),
- * Sab(*), Beta(*), Bex(*), S_Expb(*), S_Ex(*)
- DOUBLE PRECISION Bplus(50), Bminus(50), Bes(2000), Wts(2000),
- * Ben(2000), Wtn(2000)
+module mclm7_m
+use AllocateFunctions_m
+use CrystalLatticeBroadening_M
+IMPLICIT none
+contains
+!
+!
+! ----------------------------------------------------------------------
+!
+ SUBROUTINE Discre (calc, Alpha, Tbarx, Osc0, Osc1)
+!
+!
+! Scattering low calculation for discrete frequency distribution
+!
+! *** PURPOSE -- Control the convolution of discrete oscillators with
+! *** the continuous S(Alpha,Beta) computed in Contin.
+! *** Called by -- Sigcri
+! *** Uses -- Bfact, Exts, Sint
+! *** Input -- Osc_Wts(), Osc_Enx(), Osc_Snth(), Osc_Coth(),
+! *** Beta(), Bex()
+! *** -- Alpha, Dwpix, Tbeta, Tbarx
+! *** Output -- Sab()
+! *** -- Osc0, Osc1
+! *** Dummies -- Bplus(), Bminus(), S_Expb(), S_Ex(),
+! *** Bes(), Wts(), Ben(), Wtn()
+!
+!
+ class(CrystalLatticeBroadening)::calc
+ DOUBLE PRECISION Bplus(50), Bminus(50), Bes(2000), Wts(2000), &
+ Ben(2000), Wtn(2000)
DOUBLE PRECISION Alpha, Tbarx, Osc0, Osc1
- INTEGER Maxdd
- DATA Maxdd /2000/
-C
- DOUBLE PRECISION Dwf, Wt, X, Y, Besn, Wtsn, Wsave,
- & Ff1, Ff2, Ff1m, Ff2m, Bem, Be, Qqq, AL,
- & Bzero, St, Add, Db, Sint, Enx, Sum0, Sum1, Alwt
+ INTEGER, parameter:: Maxdd = 2000
+!
+ DOUBLE PRECISION Dwf, Wt, X, Y, Besn, Wtsn, Wsave, &
+ Ff1, Ff2, Ff1m, Ff2m, Bem, Be, Qqq, AL, &
+ Bzero, St, Add, Db, Enx, Sum0, Sum1, Alwt
INTEGER I, J, Nn, N, M, K, Jj, Ibeta
- DOUBLE PRECISION Zero, Half, One
- PARAMETER (Zero=0.0d0, Half=0.5d0, One=1.0d0)
-C
-C
-C *** SET UP OSCILLATOR PARAMETERS
- CALL Exts (Sab, S_Ex, Beta, Nbeta)
-C *** initializes S_Ex as needed by Sint
-C
- Wt = Tbeta
+ real(kind=8),PARAMETER::Zero=0.0d0, Half=0.5d0, One=1.0d0
+!
+!
+! *** SET UP OSCILLATOR PARAMETERS
+ call allocate_real_data(calc%S_Ex, 2*calc%Nbeta)
+ CALL Exts (calc%Sab, calc%S_Ex, calc%Beta, calc%Nbeta)
+! *** initializes S_Ex as needed by Sint
+!
+ Wt = calc%Tbeta
Qqq = Tbarx
AL = Alpha
Dwf = Zero
- IF (AL*Dw0.LT.45.0d0) Dwf = dEXP(-AL*Dw0)
-C
- CALL Zero_Array (S_Expb,Nbeta)
-C
-C
-C *** INITIALIZE FOR DELTA FUNCTION CALCULATION
+ IF (AL*calc%Dw0.LT.45.0d0) Dwf = dEXP(-AL*calc%Dw0)
+!
+ call allocate_real_data(calc%S_Expb,calc%Nbeta)
+!
+!
+! *** INITIALIZE FOR DELTA FUNCTION CALCULATION
Ben(1) = Zero
Wtn(1) = One
Nn = 1
N = 0
-C
-C
-C *** -----
-C *** LOOP OVER ALL OSCILLATORS
-C
- DO 10 I=1,N_Osc
- X = AL*Osc_Snth(I)
- Y = AL*Osc_Coth(I)
- Enx = Osc_Enx(I)
+!
+!
+! *** -----
+! *** LOOP OVER ALL OSCILLATORS
+!
+ DO 10 I=1,calc%N_Osc
+ X = AL*calc%Osc_Snth(I)
+ Y = AL*calc%Osc_Coth(I)
+ Enx = calc%Osc_Enx(I)
CALL Bfact (Bplus, Bminus, Bzero, X, Y, Enx)
- Qqq = Qqq + Osc_Coth(I)*Enx**2
-C
-C
-C *** Do convolution for the delta functions
-C *** -----
-C *** N = 0 TERM
-C
+ Qqq = Qqq + calc%Osc_Coth(I)*Enx**2
+!
+!
+! *** Do convolution for the delta functions
+! *** -----
+! *** N = 0 TERM
+!
DO M=1,Nn
Besn = Ben(M)
Wtsn = Wtn(M)*Bzero
@@ -86,12 +85,12 @@ C
END IF
END IF
END DO
-C
-C
-C *** -----
-C *** NEGATIVE N TERMS
-C ########## DO K=1,50 ? #######
-C ### BFACT is hard-wired for max of 30, not 50
+!
+!
+! *** -----
+! *** NEGATIVE N TERMS
+! ########## DO K=1,50 ? #######
+! ### BFACT is hard-wired for max of 30, not 50
DO K=1,30
IF (Bminus(K).GT.Zero) THEN
DO M=1,Nn
@@ -105,12 +104,12 @@ C ### BFACT is hard-wired for max of 30, not 50
END DO
END IF
END DO
-C
-C
-C *** -----
-C *** POSITIVE N TERMS
-C ########## DO K=1,50 ? #######
-C ### BFACT is hard-wired for max of 30, not 50
+!
+!
+! *** -----
+! *** POSITIVE N TERMS
+! ########## DO K=1,50 ? #######
+! ### BFACT is hard-wired for max of 30, not 50
DO K=1,30
IF (Bplus(K).GT.Zero) THEN
DO M=1,Nn
@@ -124,27 +123,27 @@ C ### BFACT is hard-wired for max of 30, not 50
END DO
END IF
END DO
-C
-C
-C *** -----
-C *** CONTINUE OSCILLATOR LOOP
-C
+!
+!
+! *** -----
+! *** CONTINUE OSCILLATOR LOOP
+!
Nn = N
DO M=1,Nn
Ben(M) = Bes(M)
Wtn(M) = Wts(M)
END DO
N = 0
- Wt = Wt + Osc_Wts(I)
-C
+ Wt = Wt + calc%Osc_Wts(I)
+!
10 CONTINUE
-C *** END OF LOOP OVER ALL OSCILLATORS
+! *** END OF LOOP OVER ALL OSCILLATORS
N = Nn
-C
-C
-C
-C *** -----
-C *** Sort the discrete lines and throw out the smallest ones
+!
+!
+!
+! *** -----
+! *** Sort the discrete lines and throw out the smallest ones
Nn = N - 1
DO I=1,N-1
K = 0
@@ -162,128 +161,125 @@ C *** Sort the discrete lines and throw out the smallest ones
IF (K.EQ.0) GO TO 20
END DO
20 CONTINUE
-C
+!
DO I=1,Nn
N = I
IF (Wts(I).LT.1.E-5) GO TO 30
END DO
30 CONTINUE
-C
-C
-C
-C *** -----
-C *** Add the continuum part to the scattering law
+!
+!
+!
+! *** -----
+! *** Add the continuum part to the scattering law
Alwt = Alpha * Wt
Qqq = Qqq * Alwt * 2.0d0
DO M=1,N
IF (Wts(M).NE.Zero) THEN
- DO J=1,Nbeta
- Be = - Beta(J) - Bes(M)
- St = Sint (Beta, Bex, S_Ex, Be, Alwt, Qqq, Del_S_B, Nbx,
- * Nbeta)
+ DO J=1,calc%Nbeta
+ Be = - calc%Beta(J) - Bes(M)
+ St = Sint (calc%Beta, calc%Bex, calc%S_Ex, Be, Alwt, Qqq, calc%Del_S_B, calc%Nbx, calc%Nbeta)
IF (St.NE.Zero) THEN
Add = Wts(M)*St
IF (Add.GE.1.E-20) THEN
- S_Expb(J) = S_Expb(J) + Add
+ calc%S_Expb(J) = calc%S_Expb(J) + Add
END IF
END IF
END DO
END IF
END DO
-C
-C
-C *** -----
-C *** Add the delta functions to the scattering law
-C *** Delta(0.) is saved fro [from? for?] the incoherent elastic
-C
+!
+!
+! *** -----
+! *** Add the delta functions to the scattering law
+! *** Delta(0.) is saved fro [from? for?] the incoherent elastic
+!
JJ = 0
IF (Dwf.GE.1.E-10) THEN
- IF (Twt.LE.Zero) THEN
+ IF (calc%Twt.LE.Zero) THEN
DO M=1,N
IF (Bes(M).LT.Zero) THEN
Be = -Bes(M)
-C
- IF (Be.LE.Beta(Nbeta-1)) THEN
+!
+ IF (Be.LE.calc%Beta(calc%Nbeta-1)) THEN
Db = 1000.d0
- DO J=1,Nbeta
+ DO J=1,calc%Nbeta
JJ = J
- IF (dABS(Be-Beta(J)).GT.Db) GO TO 40
- Db = dABS(Be-Beta(J))
+ IF (dABS(Be-calc%Beta(J)).GT.Db) GO TO 40
+ Db = dABS(Be-calc%Beta(J))
END DO
40 CONTINUE
-C
+!
IF (JJ.LE.2) THEN
- Add = Wts(M)/Beta(JJ)
+ Add = Wts(M)/calc%Beta(JJ)
ELSE
- Add = 2*Wts(M)/(Beta(JJ)-Beta(JJ-2))
+ Add = 2*Wts(M)/(calc%Beta(JJ)-calc%Beta(JJ-2))
END IF
-C
+!
Add = Add*Dwf
IF (Add.GE.1.d-20) THEN
- S_Expb(JJ-1) = S_Expb(JJ-1) + Add
+ calc%S_Expb(JJ-1) = calc%S_Expb(JJ-1) + Add
END IF
-C
+!
END IF
END IF
END DO
END IF
END IF
-C
-C
-C *** -----
-C *** Record the results
- DO J=1,Nbeta
- Sab(J) = S_Expb(J)
+!
+!
+! *** -----
+! *** Record the results
+ DO J=1,calc%Nbeta
+ calc%Sab(J) = calc%S_Expb(J)
END DO
-C
-C
-C *** -----
-C *** Check moments of calculated S(Alpha,Beta).
-C Ibeta = 1
- Bem = Beta(1)
- Ff1m = Sab(1)
- Ff2m = Sab(1)*dEXP(-Bem)
+!
+!
+! *** -----
+! *** Check moments of calculated S(Alpha,Beta).
+! Ibeta = 1
+ Bem = calc%Beta(1)
+ Ff1m = calc%Sab(1)
+ Ff2m = calc%Sab(1)*dEXP(-Bem)
Sum0 = Zero
Sum1 = Zero
- DO Ibeta=2,Nbeta
- Be = Beta(Ibeta)
- Ff2 = Sab(Ibeta)
- Ff1 = Sab(Ibeta)*dEXP(-Be)
+ DO Ibeta=2,calc%Nbeta
+ Be = calc%Beta(Ibeta)
+ Ff2 = calc%Sab(Ibeta)
+ Ff1 = calc%Sab(Ibeta)*dEXP(-Be)
Sum0 = Sum0 + Half*(Be-Bem)*(Ff1m+Ff2m+Ff1+Ff2)
- Sum1 = Sum1 + Half*(Be-Bem)
- & *(Ff2m*Bem+Ff2*Be-Ff1m*Bem-Ff1*Be)
+ Sum1 = Sum1 + Half*(Be-Bem)*(Ff2m*Bem+Ff2*Be-Ff1m*Bem-Ff1*Be)
Ff1m = Ff1
Ff2m = Ff2
Bem = Be
END DO
- IF (Twt.EQ.Zero) Sum0 = Sum0 + dEXP(-AL*Dwpix)
+ IF (calc%Twt.EQ.Zero) Sum0 = Sum0 + dEXP(-AL*calc%Dwpix)
Sum1 = Sum1/AL
Osc0 = Sum0
Osc1 = Sum1
RETURN
END
-C
-C
-C ----------------------------------------------------------------------
-C
+!
+!
+! ----------------------------------------------------------------------
+!
SUBROUTINE Exts (Sab, S_Ex, Beta, Nbeta)
-C
-C *** PURPOSE -- Set up the Array S_Ex for Sint.
-C *** Sab contains the asymmetric S(alpha,beta) for negative
-C *** Beta, and S_Ex contains the asymmetrix S(alpha,beta)
-C *** extended to plus and minus Beta.
-C *** CALLED BY -- Discre.
-C
- IMPLICIT NONE
+!
+! *** PURPOSE -- Set up the Array S_Ex for Sint.
+! *** Sab contains the asymmetric S(alpha,beta) for negative
+! *** Beta, and S_Ex contains the asymmetrix S(alpha,beta)
+! *** extended to plus and minus Beta.
+! *** CALLED BY -- Discre.
+!
INTEGER Nbeta, K, I
- DOUBLE PRECISION Sab(*), S_Ex(*), Beta(Nbeta)
-C
+ real(kind=8)::Beta(:), S_Ex(:), Sab(:)
+!
K = Nbeta
DO I=1,Nbeta
S_Ex(I) = Sab(K)
K = K - 1
END DO
-C
+!
IF (Beta(1).LE.0.000000001) THEN
S_Ex(Nbeta) = Sab(1)
K = Nbeta + 1
@@ -291,40 +287,39 @@ C
K = Nbeta + 2
S_Ex(Nbeta+1) = Sab(1)
END IF
-C
+!
DO I=2,Nbeta
S_Ex(K) = Sab(I)*dEXP(-Beta(I))
K = K + 1
END DO
-C
+!
RETURN
END
-C
-C
-C ----------------------------------------------------------------------
-C
- DOUBLE PRECISION FUNCTION Sint (Beta, Bex, S_Ex, Be, Alwt,
- * Qqq, Del_S_B, Nbx, Nbeta)
-C
-C INTERPOLATES IN SCATTERING FUNCTION, OR USES SCT APPROX
-C TO EXTRAPOLATE OUTSIDE THE RANGE IN MEMORY.
-C CALLED BY DISCRE.
-C
- IMPLICIT NONE
- DOUBLE PRECISION Beta(*), Bex(*), S_Ex(*)
+!
+!
+! ----------------------------------------------------------------------
+!
+ DOUBLE PRECISION FUNCTION Sint (Beta, Bex, S_Ex, Be, Alwt, &
+ Qqq, Del_S_B, Nbx, Nbeta)
+!
+! INTERPOLATES IN SCATTERING FUNCTION, OR USES SCT APPROX
+! TO EXTRAPOLATE OUTSIDE THE RANGE IN MEMORY.
+! CALLED BY DISCRE.
+!
+ real(kind=8)::Beta(:), Bex(:), S_Ex(:)
DOUBLE PRECISION Be, Alwt, Qqq, Del_S_B
INTEGER Nbx, Nbeta
-C
+!
DOUBLE PRECISION Sv, Ex, Ss1, Ss3
INTEGER K1, K2, K3
DOUBLE PRECISION Zero, Pi
PARAMETER (Zero=0.0d0, Pi=3.141592653589793238462643d0)
-C
-C
+!
+!
IF (dABS(Be).GT.Beta(Nbeta)) THEN
-C
-C *** -----
-C *** SCT approximation [extrapolation]
+!
+! *** -----
+! *** SCT approximation [extrapolation]
IF (Alwt.LE.Zero) THEN
Sv = Zero
ELSE
@@ -332,16 +327,16 @@ C *** SCT approximation [extrapolation]
IF (Be.GT.Zero) Ex = Ex - Be
Sv = dEXP(Ex)/(Pi*Qqq)
END IF
-C
+!
ELSE
-C
-C *** -----
-C *** Interpolation
+!
+! *** -----
+! *** Interpolation
K1 = 1
K2 = Nbeta
K3 = Nbx
-C
-C ** Bisect for Be
+!
+! ** Bisect for Be
130 CONTINUE
IF (Be.EQ.Bex(K2)) THEN
Sv = S_Ex(K2)
@@ -369,61 +364,60 @@ C ** Bisect for Be
END IF
END IF
END IF
-C
+!
Sint = Sv
-C
+!
RETURN
END
-C
-C
-C -------------------------------------------------------------------
-C
+!
+!
+! -------------------------------------------------------------------
+!
SUBROUTINE Bfact (Bplus, Bminus, Bzero, X, Y, Enx)
-C
-C CALCULATES THE Bessel-FUNCTION TERMS FOR DISCRETE OSCILLATORS.
-C CALLED BY DISCRE.
-C
- IMPLICIT NONE
+!
+! CALCULATES THE Bessel-FUNCTION TERMS FOR DISCRETE OSCILLATORS.
+! CALLED BY DISCRE.
+!
DOUBLE PRECISION Bplus(*), Bminus(*), Bzero, X, Y, Enx
- DOUBLE PRECISION Xx, Yy, Xy, Xt, U, V,
- & Bessio, Bessi1, Zero, One, Half
+ DOUBLE PRECISION Xx, Yy, Xy, Xt, U, V, &
+ Bessio, Bessi1, Zero, One, Half
INTEGER Imax, I, J, K
PARAMETER (Zero=0.0d0, Half=0.5d0, One=1.0D0)
-C
-C
+!
+!
Xx = X
Yy = Y
Xy = Enx
Xt = Xx/3.75d0
-C
-C
-C *** Compute Bessio and Bessi1
-C
+!
+!
+! *** Compute Bessio and Bessi1
+!
IF (Xt.LE.One) THEN
U = Xt*Xt
- Bessio = One
- & + U * (3.5156229 + U * (3.0899424 + U * (1.2067492
- & + U * (0.2659732 + U * (0.0360768 + U * 0.0045813 )))))
- Bessi1 = (0.5
- & + U * (0.87890594 + U * (0.51498869 + U *(0.15084934
- & + U * (0.02658733 + U * (0.00301532 + U * 0.00032411))))))*Xx
+ Bessio = One &
+ + U * (3.5156229 + U * (3.0899424 + U * (1.2067492 &
+ + U * (0.2659732 + U * (0.0360768 + U * 0.0045813 )))))
+ Bessi1 = (0.5 &
+ + U * (0.87890594 + U * (0.51498869 + U *(0.15084934 &
+ + U * (0.02658733 + U * (0.00301532 + U * 0.00032411))))))*Xx
ELSE
U = One/Xt
- Bessio = (
- & 0.39894228 + U * (0.01328592 + U * ( 0.00225319
- & + U * (-0.00157565 + U * (0.00916281 + U * (-0.02057706
- & + U * ( 0.02635537 + U * (-.01647633 + U * 0.00392377
- & )))))))) / SQRT(Xx)
- Bessi1 = 0.02282967 + U * (-0.02895312 + U * ( 0.01787654
- & -U * 0.00420059 ))
- Bessi1 = 0.39894228 + U * (-0.03988024 + U * (-0.00362018
- & + U * ( 0.00163801 + U * (-0.01031555 + U * Bessi1 ))))
+ Bessio = ( &
+ 0.39894228 + U * (0.01328592 + U * ( 0.00225319 &
+ + U * (-0.00157565 + U * (0.00916281 + U * (-0.02057706 &
+ + U * ( 0.02635537 + U * (-.01647633 + U * 0.00392377 &
+ )))))))) / SQRT(Xx)
+ Bessi1 = 0.02282967 + U * (-0.02895312 + U * ( 0.01787654 &
+ -U * 0.00420059 ))
+ Bessi1 = 0.39894228 + U * (-0.03988024 + U * (-0.00362018 &
+ + U * ( 0.00163801 + U * (-0.01031555 + U * Bessi1 ))))
Bessi1 = Bessi1/SQRT(Xx)
END IF
-C
-C
-C *** Generate higher-order Bessel functions by reverse recursion
-C
+!
+!
+! *** Generate higher-order Bessel functions by reverse recursion
+!
Imax = 30
Bplus(Imax ) = Zero
Bplus(Imax-1) = One
@@ -436,17 +430,17 @@ C
Bplus(J) = Bplus(J)*1.0d-10
END DO
END IF
-C end of loop on I
+! end of loop on I
END DO
-C
+!
U = Bessi1/Bplus(1)
DO I=1,Imax-1
Bplus(I) = Bplus(I)*U
END DO
-C
-C
-C *** Apply exponential terms to Bessel Functiions
-C
+!
+!
+! *** Apply exponential terms to Bessel Functiions
+!
U = Xy*Half
IF (Xt.LT.One) THEN
V = - Yy
@@ -462,3 +456,4 @@ C
END DO
RETURN
END
+end module mclm7_m
diff --git a/sammy/src/clm/mclm8.f90 b/sammy/src/clm/mclm8.f90
index b8816fe2a46fc465fb9052dcdee146ea2a7e9dc9..9d10f02a4c5df27d6ef0c207779892b4c6aba0ce 100644
--- a/sammy/src/clm/mclm8.f90
+++ b/sammy/src/clm/mclm8.f90
@@ -1,22 +1,22 @@
-C
-C
-C -----------------------------------------------------------------
-C
- SUBROUTINE Find_L (grid, Xx, L1, L2, Lstart, L1_Error,
- * L2_Error)
-C *** Purpose -- find L1 & L2 such that Energy(L1) < Xx < Energy(L2)
-C
- use SammyGridAccess_M
- use EndfData_common_m
- IMPLICIT NONE
+module mcml8_m
+use CrystalLatticeBroadening_M
+implicit none
+contains
+!
+!
+! -----------------------------------------------------------------
+!
+ SUBROUTINE Find_L (calc, Xx, L1, L2, Lstart, L1_Error, L2_Error)
+! *** Purpose -- find L1 & L2 such that Energy(L1) < Xx < Energy(L2)
+!
+ class(CrystalLatticeBroadening)::calc
DOUBLE PRECISION Xx
- type(SammyGridAccess)::grid ! can pass even in f77 as in same compile unit
INTEGER L1, L2, Lstart, Lk, L1_Error, L2_Error
real(kind=8)::eLstart, eL1, eL2, eLk
integer::numEl
-C
- numEl = grid%getNumEnergies(expData)
- eLstart = grid%getEnergy(Lstart, expData)
+!
+ numEl = calc%getNumEnergyUnbroadened()
+ eLstart = calc%getEnergyUnbroadened(Lstart)
IF (Xx.EQ.eLstart) THEN
L1 = Lstart
L2 = Lstart
@@ -25,24 +25,24 @@ C
L1_Error = L1_Error + 1
L1 = 1
L2 = 2
-CX STOP '[Stop in Find_L in clm/mclm8.f # 1]'
+!X STOP '[Stop in Find_L in clm/mclm8.f # 1]'
ELSE
DO Lk=Lstart-1,1,-1
- eLk = grid%getEnergy(Lk, expData)
+ eLk = calc%getEnergyUnbroadened(Lk)
IF (Xx.GE.eLk) THEN
L1 = Lk
L2 = Lk + 1
- eL1 = grid%getEnergy(L1, expData)
- eL2 = grid%getEnergy(L2, expData)
+ eL1 = calc%getEnergyUnbroadened(L1)
+ eL2 = calc%getEnergyUnbroadened(L2)
IF (Xx.LT.eL2) RETURN
WRITE (6,10000) L1, eL1, Xx, L2, eL2
-10000 FORMAT (' Eb(', I5, ')=',1PG14.6,' Xx=', 1PG14.6,
- * ' Eb(', I5, ')=', 1PG14.6)
+10000 FORMAT (' Eb(', I5, ')=',1PG14.6,' Xx=', 1PG14.6, &
+ ' Eb(', I5, ')=', 1PG14.6)
STOP '[Stop in Find_L in clm/mclm8.f # 2]'
END IF
END DO
Lk = 1
- eLk = grid%getEnergy(Lk, expData)
+ eLk = calc%getEnergyUnbroadened(Lk)
WRITE (6,10000) Lstart, eLstart, Xx, Lk, eLk
STOP '[Stop in Find_L in clm/mclm8.f # 3]'
END IF
@@ -53,24 +53,25 @@ CX STOP '[Stop in Find_L in clm/mclm8.f # 1]'
L2 = numEl
ELSE
DO Lk=Lstart+1,numEl
- eLk = grid%getEnergy(Lk, expData)
+ eLk = calc%getEnergyUnbroadened(Lk)
IF (Xx.LE.eLk) THEN
L2 = Lk
L1 = Lk - 1
- eL1 = grid%getEnergy(L1, expData)
- eL2 = grid%getEnergy(L2, expData)
+ eL1 = calc%getEnergyUnbroadened(L1)
+ eL2 = calc%getEnergyUnbroadened(L2)
IF (Xx.GE.eL1) RETURN
WRITE (6,10000) L1, eL1, Xx, L2, eL2
STOP '[Stop in Find_L in clm/mclm8.f # 5]'
END IF
END DO
- eL2 = grid%getEnergy(numEl, expData)
+ eL2 = calc%getEnergyUnbroadened(numEl)
WRITE (6,10100) numEl, eL2, Xx
10100 FORMAT (' Eb(', I5, ')=',1PG14.6,' Xx=', 1PG14.6)
L1 = numEl - 1
L2 = numEl
-cx STOP '[Stop in Find_L in clm/mclm8.f # 6]'
+!x STOP '[Stop in Find_L in clm/mclm8.f # 6]'
END IF
END IF
RETURN
END
+end module mcml8_m
diff --git a/sammy/src/clm/mclm8a.f90 b/sammy/src/clm/mclm8a.f90
index 663426350013c5c8761ad222a414a6abfe1640af..af60f35106e3108ce5c57fa804769a194a05f750 100644
--- a/sammy/src/clm/mclm8a.f90
+++ b/sammy/src/clm/mclm8a.f90
@@ -1,4 +1,12 @@
module Qtrap_Clm_m
+ use CrystalLatticeBroadening_M
+ use fixedr_m, only : Aaawww
+ use constn_common_m, only : Aneutr
+ use DerivativeHandler_M
+ use mcml8_m
+ use, intrinsic :: ISO_C_BINDING
+ IMPLICIT none
+
private
public Qtrap_Clm
contains
@@ -6,8 +14,8 @@ module Qtrap_Clm_m
!
! -----------------------------------------------------------------
!
- SUBROUTINE Qtrap_Clm (tmpCalc, derivs, &
- Bs, Ss, Tevx, Tempx, Em, C_Norm, Ktempx, Na, &
+ SUBROUTINE Qtrap_Clm (calc, tmpCalc, derivs, &
+ Tevx, Tempx, Em, C_Norm, Ktempx, Na, &
Isox, Locate, Ns, derivsNew, irow)
!
! *** INTEGRATION BY TRAPEZOIDAL METHOD
@@ -19,36 +27,23 @@ module Qtrap_Clm_m
! *** Bs - Beta array (variable of integration)
! *** Ss - Weighting function (S_alpha_beta)
!
- use fixedi_m, only : numUsedPar
- use fixedr_m, only : Aaawww
- use clm_common_m, only : Dwpix, Eps, Mode_S_Norm
- use constn_common_m, only : Aneutr
- use SammyGridAccess_M
- use EndfData_common_m
- use DerivativeHandler_M
- use, intrinsic :: ISO_C_BINDING
- IMPLICIT None
+ class(CrystalLatticeBroadening)::calc
INTEGER Ktempx, Na, Isox, Locate, Ns
- DOUBLE PRECISION Bs(*), Ss(*)
DOUBLE PRECISION Tevx, Tempx, Em, C_Norm
DOUBLE PRECISION S, Olds, Xxa, Xxb, Del, Sel, E_Prime
DOUBLE PRECISION Aaa, Xx, Yy, val
INTEGER Jmax, Ibet, Lstart, L1, L2, irow, Ipar, N
- type(SammyGridAccess)::grid
type(DerivativeHandler)::derivs, derivsNew, tmpCalc
real(kind=8),parameter:: Zero=0.0d0, Half=0.5d0, One=1.0d0
real(kind=8)::Alpha, Recul
integer::It, J, Jj, L1er, L2er, Nax, numEl
logical(C_BOOL)::accu
-!
- call grid%initialize()
- call grid%setToAuxGrid(expData)
!
L1er = 0
L2er = 0
Nax = derivsNew%getNnnsig()
call tmpCalc%setNnsig(Nax)
- call tmpCalc%reserve(Nax, numUsedPar+1)
+ call tmpCalc%reserve(Nax, calc%numPar+1)
IF (Locate.NE.0) Nax = 1
!
Lstart = 1
@@ -58,20 +53,20 @@ module Qtrap_Clm_m
!
call tmpCalc%nullify()
!
- Xxa = Bs(Ibet)
- Xxb = Bs(Ibet+1)
+ Xxa = calc%Bs(Ibet)
+ Xxb = calc%Bs(Ibet+1)
Del = Xxb - Xxa
Aaa = Del*Half
!
! *** First Point (Xxa)
Xx = Em + Xxa
- Call Find_L (grid, Xx, L1, L2, Lstart, L1er, L2er)
+ Call Find_L (calc, Xx, L1, L2, Lstart, L1er, L2er)
Lstart = L1
- Yy = Ss(Ibet)
+ Yy = calc%Ss(Ibet)
IF (Yy.EQ.Zero) THEN
ELSE
Yy = Yy*Aaa
- Call Sig_Int (grid, tmpCalc, derivs, &
+ Call Sig_Int (calc, tmpCalc, derivs, &
Tempx, Xx, Yy, Ktempx, Na, Isox, &
Locate, L1, L2, 1, 1)
END IF
@@ -79,12 +74,12 @@ module Qtrap_Clm_m
!
! *** Last Point (Xxb)
Xx = Em + Xxb
- Call Find_L (grid, Xx, L1, L2, Lstart, L1er, L2er)
- Yy = Ss(Ibet+1)
+ Call Find_L (calc, Xx, L1, L2, Lstart, L1er, L2er)
+ Yy = calc%Ss(Ibet+1)
IF (Yy.EQ.Zero) THEN
ELSE
Yy = Yy*Aaa
- Call Sig_Int (grid, tmpCalc, derivs, &
+ Call Sig_Int (calc, tmpCalc, derivs, &
Tempx, Xx, Yy, Ktempx, Na, Isox, &
Locate, L1, L2, 1, 1)
END IF
@@ -97,7 +92,7 @@ module Qtrap_Clm_m
E_Prime = Xxa + Half*Del
!
! *** Take half of what's already there
- do ipar = 0, numUsedPar
+ do ipar = 0, calc%numPar
do N = 1, Nax
val = tmpCalc%getDataNs(1, N, ipar, 1)*Half
if (val.eq.0.0d0) cycle
@@ -108,14 +103,14 @@ module Qtrap_Clm_m
DO Jj=1,It
!
! *** Find multiplier Yy
- IF (Ss(Ibet).EQ.Zero) THEN
+ IF (calc%Ss(Ibet).EQ.Zero) THEN
Yy = Zero
ELSE
- Yy = Ss(Ibet)*(Xxb-E_Prime)
+ Yy = calc%Ss(Ibet)*(Xxb-E_Prime)
END IF
- IF (Ss(Ibet+1).EQ.Zero) THEN
+ IF (calc%Ss(Ibet+1).EQ.Zero) THEN
ELSE
- Yy = Yy + Ss(Ibet+1)*(E_Prime-Xxa)
+ Yy = Yy + calc%Ss(Ibet+1)*(E_Prime-Xxa)
END IF
!
! *** Add new piece
@@ -123,9 +118,9 @@ module Qtrap_Clm_m
ELSE
Yy = Yy*Aaa
Xx = Em + E_Prime
- Call Find_L (grid, Xx, L1, L2, Lstart, L1er, L2er)
+ Call Find_L (calc, Xx, L1, L2, Lstart, L1er, L2er)
- Call Sig_Int (grid, tmpCalc, derivs, &
+ Call Sig_Int (calc, tmpCalc, derivs, &
Tempx, Xx, Yy, Ktempx, Nax, Isox, &
Locate, L1, L2, 1, 1)
END IF
@@ -134,7 +129,7 @@ module Qtrap_Clm_m
END DO
!
S = tmpCalc%getDataNs(1,1,0,1)
- IF (ABS(S-Olds).LT.Eps*ABS(Olds)) GO TO 10
+ IF (ABS(S-Olds).LT.calc%Eps*ABS(Olds)) GO TO 10
It = It*2
Aaa = Aaa*Half
Del = Del*Half
@@ -147,7 +142,7 @@ module Qtrap_Clm_m
accu = .true.
call derivsNew%setAccumulate(accu)
- do ipar = 0, numUsedPar
+ do ipar = 0, calc%numPar
do N = 1, Nax
val = tmpCalc%getDataNs(1, N, ipar, 1)
if (val.eq.0.0d0) cycle
@@ -166,11 +161,11 @@ module Qtrap_Clm_m
Recul = Em*(Aneutr/(Aaawww+Aneutr))
Alpha = Recul/Tevx
!
- Sel = dEXP(-Dwpix*Alpha)
- IF (Mode_S_Norm.EQ.0) THEN
- ELSE IF (Mode_S_Norm.EQ.1) THEN
+ Sel = dEXP(-calc%Dwpix*Alpha)
+ IF (calc%Mode_S_Norm.EQ.0) THEN
+ ELSE IF (calc%Mode_S_Norm.EQ.1) THEN
Yy = (One-Sel)/(C_Norm-Sel)
- do ipar = 0, numUsedPar
+ do ipar = 0, calc%numPar
do N = 1, Nax
val = Yy * derivsNew%getDataNs( irow, N, Ipar, Isox)
if (val.eq.0.0d0) cycle
@@ -183,25 +178,23 @@ module Qtrap_Clm_m
! Sigman(Em+Recul)*EXP(-W*A)
Xx = Em + Recul
Yy = Sel
- Call Find_L (grid, Xx, L1, L2, Lstart, L1er, L2er)
- Call Sig_Int (grid, derivsNew, derivs, &
+ Call Find_L (calc, Xx, L1, L2, Lstart, L1er, L2er)
+ Call Sig_Int (calc, derivsNew, derivs, &
Tempx, Xx, Yy, Ktempx, Na, Isox, &
Locate, L1, L2, irow, Isox)
!
IF (L1er.GT.0) THEN
- WRITE (6,10100) L1er, Em, grid%getEnergy(1, expData)
+ WRITE (6,10100) L1er, Em, calc%getEnergyUnbroadened(1)
10100 FORMAT ('###', I5, '= number of omitted terms', 1X, &
'for E=', 1PG13.5, ' at Eb(1)=', 1PG13.5)
END IF
IF (L2er.GT.0) THEN
- numEl = grid%getNumEnergies(expData)
+ numEl = calc%getNumEnergyUnbroadened()
WRITE (6,10200) L2er, numEl, Em, &
- grid%getEnergy(numEl, expData)
+ calc%getEnergyUnbroadened(numEl)
10200 FORMAT ('###', I5, '= omitted terms', 1X, &
'for E=', 1PG13.5, ' at Eb(', I5, ')=', 1PG13.5)
END IF
-!
- call grid%destroy()
RETURN
END SUBROUTINE Qtrap_Clm
@@ -210,20 +203,14 @@ module Qtrap_Clm_m
!
! -----------------------------------------------------------------
!
- SUBROUTINE Sig_Int_L (grid, tmpCalc, derivs, &
+ SUBROUTINE Sig_Int_L (calc, tmpCalc, derivs, &
Tempx, Xx, Yy, Ktempx, Na, Isox, &
Locate, L1, L2, irow, ourIso)
!
! *** Purpose -- Linear interpolate to get CrossSection(Xx) using
! *** CrossSection(Energb(L1)) and CrossSection(Energb(L2))
!
- use fixedi_m, only : numUsedPar
- use constn_common_m
- use SammyGridAccess_M
- use EndfData_common_m
- use DerivativeHandler_M
- IMPLICIT none
- type(SammyGridAccess)::grid
+ class(CrystalLatticeBroadening)::calc
INTEGER Ktempx, Na, Isox, Locate, L1, L2
type(DerivativeHandler)::derivs, tmpCalc
real(kind=8)::Tempx, Xx, Yy
@@ -237,8 +224,8 @@ module Qtrap_Clm_m
Aa(1) = Yy
Aa(2) = Zero
ELSE
- E1 = grid%getEnergy(L1, expData)
- E2 = grid%getEnergy(L2, expData)
+ E1 = calc%getEnergyUnbroadened(L1)
+ E2 = calc%getEnergyUnbroadened(L2)
D = E2 - E1
Aa(1) = (E2-Xx)/D * Yy
Aa(2) = (Xx-E1)/D * Yy
@@ -249,7 +236,7 @@ module Qtrap_Clm_m
! ***
! ---
! update cross section (ipar=0) as well as derivatives
- do Ipar = 0, numUsedPar
+ do Ipar = 0, calc%numPar
DO N=1,Na
L = 0
val = 0.0d0
@@ -267,7 +254,7 @@ module Qtrap_Clm_m
! ***
! ---
! update cross section (ipar=0) as well as derivatives
- do Ipar = 0, numUsedPar
+ do Ipar = 0, calc%numPar
L = 0
val = 0.0d0
DO LL=L1,L2
@@ -287,21 +274,14 @@ module Qtrap_Clm_m
!
! -----------------------------------------------------------------
!
- SUBROUTINE Sig_Int (grid, tmpCalc, derivs, &
+ SUBROUTINE Sig_Int (calc, tmpCalc, derivs, &
Tempx, Xx, Yy, Ktempx, Na, Isox, &
Locate, L1, L2, irow, ourIso)
!
! *** Purpose -- Quadratic interpolate to get CrossSection(Xx) using
! *** CrossSection(Energb(L1-1)) to CrossSection(Energb(L2+1))
!
- use fixedi_m, only : numUsedPar
- use constn_common_m
- use SammyGridAccess_M
- use EndfData_common_m
- use DerivativeHandler_M
- use, intrinsic :: ISO_C_BINDING
- IMPLICIT none
- type(SammyGridAccess)::grid
+ class(CrystalLatticeBroadening)::calc
INTEGER Ktempx, Na, Isox, Locate, L1, L2
real(kind=8)::Tempx, Xx, Yy
type(DerivativeHandler)::derivs, tmpCalc
@@ -313,19 +293,19 @@ module Qtrap_Clm_m
logical(C_BOOL)::accu
!
- Ndatb = grid%getNumEnergies(expData)
+ Ndatb = calc%getNumEnergyUnbroadened()
IF (L1.EQ.1 .OR. L2.EQ.Ndatb .OR. L1.EQ.L2) THEN
- CALL Sig_Int_L (grid, tmpCalc, derivs, &
+ CALL Sig_Int_L (calc, tmpCalc, derivs, &
Tempx, Xx, Yy, Ktempx, Na, Isox, &
Locate, L1, L2, irow, ourIso)
RETURN
END IF
!
Iso = Isox
- E1 = grid%getEnergy(L1-1, expData)
- E2 = grid%getEnergy(L1 , expData)
- E3 = grid%getEnergy(L2 , expData)
- E4 = grid%getEnergy(L2+1, expData)
+ E1 = calc%getEnergyUnbroadened(L1-1)
+ E2 = calc%getEnergyUnbroadened(L1)
+ E3 = calc%getEnergyUnbroadened(L2)
+ E4 = calc%getEnergyUnbroadened(L2+1)
!
D21 = E2 - E1
D32 = E3 - E2
@@ -353,7 +333,7 @@ module Qtrap_Clm_m
! ***
! ---
! update cross section (ipar=0) as well as derivatives
- Do Ipar = 0, numUsedPar
+ Do Ipar = 0, calc%numPar
DO N=1,Na
L = 0
val = 0.0d0
@@ -371,7 +351,7 @@ module Qtrap_Clm_m
! IF (Locate.NE.0) THEN
! ***
! update cross section (ipar=0) as well as derivatives
- Do Ipar = 0, numUsedPar
+ Do Ipar = 0, calc%numPar
L = 0
val = 0.0d0
DO LL=L1-1,L2+1
diff --git a/sammy/src/cro/mcro5.f90 b/sammy/src/cro/mcro5.f90
index b96d0626a33445d2f05f011b6fdb9d0a0a8aa3b0..c708f3fcf0bcd289bdd0ec783e150c5e9ba681a3 100644
--- a/sammy/src/cro/mcro5.f90
+++ b/sammy/src/cro/mcro5.f90
@@ -163,6 +163,7 @@ end module sumAndConvertAfterBroadening_m
use broad_common_m
use samdbd_m
use Samfgm_0_m
+ use Samclm_m
use FreeGasDopplerBroadening_M
use DopplerAndResolutionBroadener_M
use GridData_M
@@ -257,7 +258,7 @@ end module sumAndConvertAfterBroadening_m
else IF ( dopplerOption%bType.EQ.2) then
call Samfgm_0(dopplerOption%freeGas)
else IF ( dopplerOption%bType.EQ.3) then
- call Samclm_0
+ call Samclm_0(dopplerOption%crystalLattice)
end if
call sumAndConvertAfterBroadening(moreBroadening, Yssmsc, dopplerOption%broadener)
diff --git a/sammy/src/dbd/HighEnergyFreeGas_m.f90 b/sammy/src/dbd/HighEnergyFreeGas_m.f90
index 293fbcd84ea30bc89b769bf6655229c31931418a..65d4210abd11ebc35e175a407935887859d862be 100644
--- a/sammy/src/dbd/HighEnergyFreeGas_m.f90
+++ b/sammy/src/dbd/HighEnergyFreeGas_m.f90
@@ -41,7 +41,6 @@ subroutine HighEnergyFreeGas_broaden(this)
! reserve the data
ndatb = this%getNumEnergyBroadened()
- write(0,*)" I am called !!!!"
call this%getData(data)
call data%nullify()
call data%reserve(ndatb*data%getNnnsig(), this%numPar+1)
diff --git a/sammy/src/fgm/FreeGasDopplerBroadening_M.f90 b/sammy/src/fgm/FreeGasDopplerBroadening_M.f90
index ef6b5e421b7c27397e1952aa0a91d4afcfa65124..bee13b3510fca9e246f6a7fb5d6cbc7388ef1707 100644
--- a/sammy/src/fgm/FreeGasDopplerBroadening_M.f90
+++ b/sammy/src/fgm/FreeGasDopplerBroadening_M.f90
@@ -36,6 +36,7 @@ subroutine FreeGasDopplerBroadening_initialize(this, hand, list, work)
call DopplerBroadening_initialize(this, hand, list, work)
this%debugOutput = .false.
this%Elowbr = 0.0d0
+ this%numPar = 0
end subroutine
subroutine FreeGasDopplerBroadening_destroy(this)