diff --git a/sammy/src/amr/mamr2.f90 b/sammy/src/amr/mamr2.f90
index 77c83c95881591efa62f4c677d9ed6c6bedcb83e..3bac4e0299525afa3dc23a4f8ed66f9fa193c49e 100644
--- a/sammy/src/amr/mamr2.f90
+++ b/sammy/src/amr/mamr2.f90
@@ -234,6 +234,11 @@ module amr2
       Co2    = Dum(9)
       Do2    = Dum(10)
       Dopple = Dum(11)
+      if (Dopple.ne.0.0d0) then
+        dopplerOption%wantBroaden = .true.
+      else
+         dopplerOption%wantBroaden = .false.
+      end if
       call radFitFlags%setMatchingK(Dum(12))
       Dist   = Dum(13)
       Anorm  = Dum(14)
diff --git a/sammy/src/blk/Abro_common.f90 b/sammy/src/blk/Abro_common.f90
index 46317cb405bef50234f972dba0fd2213c789bc23..9c70f35d297d7f53a680a98b7c1d7d374fccf47e 100644
--- a/sammy/src/blk/Abro_common.f90
+++ b/sammy/src/blk/Abro_common.f90
@@ -7,7 +7,6 @@ module abro_common_m
       double precision, pointer :: Thck     => Aabro(1)
       double precision, pointer :: Tempe    => Aabro(2)
       double precision, pointer :: Delv     => Aabro(3)
-      double precision, pointer :: Delt     => Aabro(4)
       double precision, pointer :: Dddelt   => Aabro(5)
       double precision, pointer :: Odffff   => Aabro(6)
       double precision, pointer :: Abro1    => Aabro(7)
diff --git a/sammy/src/blk/Broad_common.f90 b/sammy/src/blk/Broad_common.f90
index e45b0f84dede6ebd364c88de71bb82678bd822f6..5ecdd9d97aa3f76731acdf5682ee25c377b72c55 100644
--- a/sammy/src/blk/Broad_common.f90
+++ b/sammy/src/blk/Broad_common.f90
@@ -2,6 +2,8 @@
 module broad_common_m
      use FreeGasDopplerBroadening_M
      use HighEnergyFreeGas_m
+     use LealHwangBroadening_M
+     use CrystalLatticeBroadening_M
      use DopplerAndResolutionBroadener_M
      use DopplerBroadening_M
      use GridData_M
@@ -37,11 +39,12 @@ module broad_common_m
      type DopplerBroadenInfo
        type(FreeGasDopplerBroadening),pointer::freeGas
        type(HighEnergyFreeGas),pointer::highEnergyFreeGas
-       type(DopplerBroadening),pointer::lealHwang
-       type(DopplerBroadening),pointer::crystalLattice
+       type(LealHwangBroadening),pointer::lealHwang
+       type(CrystalLatticeBroadening),pointer::crystalLattice
 
        class(DopplerBroadening),pointer::broadener
        integer::bType=2  ! default to freeGas
+       logical::wantBroaden=.false.
      end type
      type(DopplerBroadenInfo)::dopplerOption
 
diff --git a/sammy/src/blk/Fixedi_common.f90 b/sammy/src/blk/Fixedi_common.f90
index 0e3d992cd42b72e50b6f8044ab5e0c10e394c7b8..9045efbb59c9073dcf93cf047c8f19d4d00d399a 100644
--- a/sammy/src/blk/Fixedi_common.f90
+++ b/sammy/src/blk/Fixedi_common.f90
@@ -122,8 +122,6 @@ module fixedi_m
      integer,pointer :: Nres    => lfdim(75)
      integer,pointer :: Nrext   => lfdim(76)
      integer,pointer :: Nlfdimsiz  => lfdim(77)
-     integer,pointer :: Kjjjjj  => lfdim(78)
-     integer,pointer :: K2pls1  => lfdim(79)
      integer,pointer :: Ngtvv   => lfdim(80)
      integer,pointer :: Nxtra   => lfdim(81)
      integer,pointer :: Nntype  => lfdim(83)
diff --git a/sammy/src/blk/Fixedr_common.f90 b/sammy/src/blk/Fixedr_common.f90
index dfb068d0c5ef41ad64e0abfd72b218407610b05f..1a77c5e56237e553d97761f3f43308a2a3208702 100644
--- a/sammy/src/blk/Fixedr_common.f90
+++ b/sammy/src/blk/Fixedr_common.f90
@@ -29,7 +29,6 @@ module fixedr_m
      ! old group 3
      double precision, pointer :: Backd   => Ff(21)
      double precision, pointer :: Backf   => Ff(22)
-     double precision, pointer :: Delttt  => Ff(23)
      double precision, pointer :: Delttx  => Ff(24)
      double precision, pointer :: Delvvv  => Ff(25)
      double precision, pointer :: Delvvx  => Ff(26)
diff --git a/sammy/src/blk/Ifwrit_common.f90 b/sammy/src/blk/Ifwrit_common.f90
index df2c07138e34ef5a557489b2d460f28c08a3838c..fd1bebafe5d7b2a985ab2f71e73ed6ba5fd69196 100644
--- a/sammy/src/blk/Ifwrit_common.f90
+++ b/sammy/src/blk/Ifwrit_common.f90
@@ -83,7 +83,6 @@ module ifwrit_m
      
      ! old group 8
      integer,pointer :: Kedxfw => Lwrit(63)
-     integer,pointer :: Jjjdop => Lwrit(64)
      integer,pointer :: Inmdts => Lwrit(65)
      integer,pointer :: Kexpsh => Lwrit(66)
      integer,pointer :: Kendf  => Lwrit(67)
diff --git a/sammy/src/clm/dopush.f b/sammy/src/clm/dopush.f
index 1de204dccbbf71ddebd9346540b86836fe550a81..bc9958d63892d5037b681aab3905ecb085c31307 100644
--- a/sammy/src/clm/dopush.f
+++ b/sammy/src/clm/dopush.f
@@ -250,6 +250,7 @@ C
       Tevf = Tbeta*Tevf + 0.5d0*Tsave
       Ne = K - 1
       Dopplef = dSQRT(Tevf/Awr)
+      dopplerOption%wantBroaden = .true.
 C
 C *** FREE GAZ, FREE GAZ AT EFFECTIVE TEMPERATURE AND
 C     NUCLEAR CROSS-SECTION CALCULATIONS ON ADJUSTED ENERGY GRID
diff --git a/sammy/src/cro/mcro5.f90 b/sammy/src/cro/mcro5.f90
index afa5c3abac0ae59c79f5d52f41d7b0a985943dc4..b96d0626a33445d2f05f011b6fdb9d0a0a8aa3b0 100644
--- a/sammy/src/cro/mcro5.f90
+++ b/sammy/src/cro/mcro5.f90
@@ -88,13 +88,7 @@ end module sumAndConvertAfterBroadening_m
       IMPLICIT none
       integer::icr, Ntotal
       logical::haveDerivs, doAngle
-!
-!     deltdp => delttt from fixedr_m (Ff(23) assigned to delt (from common block b39
-!     dpdelv => delvvv  from fixedr_m (Ff(25) assigned to  delv (from common block b39
-!     dpdelt =>  Gammmm from fixedr_m (Ff(27) assigned to  Dddelt (from common block b39
-!
-!      EQUIVALENCE (Deltdp,Delttt), (Dpdelv,Delvvv), (Dpdelt,Gammmm)
-!
+
       Jcros  = Kcros
       Icr    = Kcros
       doAngle = Kcros.EQ.7 .OR. Kcros.eq.11
@@ -103,7 +97,6 @@ end module sumAndConvertAfterBroadening_m
       Thck   = Thick
       Tempe  = Temp
       Delv   = Delvvv  ! Dpdelv
-      Delt   = Delttt  ! Deltdp
       Dddelt = Gammmm  ! Dpdelt
       Odffff = Odfmul
 !
@@ -144,7 +137,6 @@ end module sumAndConvertAfterBroadening_m
       Filein = 'SAM51.DAT'
       Filout = 'SAM54.DAT'
 !
-      Jjjdop = dopplerOption%bType
 
       RETURN
       END
@@ -182,11 +174,17 @@ end module sumAndConvertAfterBroadening_m
       logical::debugOut
       logical(C_BOOL)::moreBroadening, needIsos
       type(FreeGasDopplerBroadening)::broadener
-      integer::Jwwwww
+      integer::Jwwwww, Jjjdop, btypeSave
 
       if (Kresol.EQ.1) Kplotu = 0
       debugOut = debug.or.Kplotu.ne.0
 
+      ! jjjdop = 1 is for special print-out in
+      ! Samint_0 right after Leal Hwang broadening
+      ! Samint_0  will reset to 0, for normal
+      ! intermediate printout
+      Jjjdop = dopplerOption%bType
+
       ! Jwwwww sets the title for intermediate or final output
       Jwwwww = 1
       IF (Ksindi.GT.0 .and. Kcros.EQ.8) jwwwww = 5
@@ -194,10 +192,10 @@ end module sumAndConvertAfterBroadening_m
       ! do maxwellian if desired
       if (Maxwel.eq.1) then
          ! if we want intermediate output, do so
-         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww)
-         if (Ksolve.ne.2) then
+         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww, Jjjdop)
+         if (Ksolve.ne.2) then           
             call Samsqu_0
-         end if
+         end if                  
          call Sammxw_0
          call Samend_0
          return
@@ -210,7 +208,7 @@ end module sumAndConvertAfterBroadening_m
          call calcDataSelf%sumOverIsotopes(numUsedPar+1)
 
          ! if we want intermediate output, do so
-         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww)
+         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww, Jjjdop)
 
          ! if fit: set's fitoption and 'genfit'
          call Samsqu_0
@@ -220,7 +218,7 @@ end module sumAndConvertAfterBroadening_m
       ! any angles to calculate
       IF (Kcros.EQ.7 .OR. Kcros.EQ.11) then
           ! if we want intermediate output, do so
-          if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww)
+          if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww, Jjjdop)
           call Samang_0
           Jwwwww = 9
 
@@ -244,7 +242,7 @@ end module sumAndConvertAfterBroadening_m
          end if
 
          ! if we want intermediate output, do so
-         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww)
+         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww, Jjjdop)
 
          ! switch handler data grid
          call Set_Kws
@@ -254,7 +252,8 @@ end module sumAndConvertAfterBroadening_m
          IF (dopplerOption%bType.EQ.0) then
             call Samdbd_0(dopplerOption%highEnergyFreeGas)
          else IF (dopplerOption%bType.EQ.1)  then
-            call Samdop_0            
+            call Samdop_0(dopplerOption%lealHwang)
+            Jjjdop = 1
          else IF ( dopplerOption%bType.EQ.2)  then
             call Samfgm_0(dopplerOption%freeGas)            
          else IF ( dopplerOption%bType.EQ.3) then
@@ -272,19 +271,20 @@ end module sumAndConvertAfterBroadening_m
       if (Ntgrlq.EQ.1) then  
           if (Ksolve.NE.2) then
               ! if we want intermediate output, do so              
-              if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww)
+              if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww, Jjjdop)
               call Samsqu_0
           else if (Ksolve.EQ.2) then
               call Samntg_0
               Kplotu = 0
-          end if
+          end if          
           return  ! can't have multiple scattering or resolution broadening
       end if
 
       ! multiple scattering
       IF (Yssmsc) THEN
          ! if we want intermediate output, do so
-         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww)         
+         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww, Jjjdop)
+         Jjjdop = 0
          IF (Kvers7.EQ.0) then
              call Samssm_0
          else IF (Kvers7.EQ.1) then
@@ -300,7 +300,8 @@ end module sumAndConvertAfterBroadening_m
       ! resolution brodening
       IF (Yresol) THEN
          ! if we want intermediate output, do so
-         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww)
+         if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww, Jjjdop)
+         Jjjdop = 0
          Kplotu = 0
 
          ! allowed is:
@@ -429,10 +430,17 @@ end module sumAndConvertAfterBroadening_m
              end if
              call normAndBack(moreBroadening, resolutionOption%resBroad1)
 
+             ! if Leal Hwang we use a different auxillary grid
+             ! after first resolution function this is the
+             ! normal grid, so signal that to the
+             ! second resolution function
+             btypeSave = dopplerOption%bType
+             dopplerOption%bType = 2
+
              Jwwwww = 4
 
              ! if we want output after first resolution function
-             if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww)
+             if (debugOut)  call Samint_0(debugOut, Yssmsc, Jwwwww, Jjjdop)
 
 
              moreBroadening = .false.
@@ -452,6 +460,9 @@ end module sumAndConvertAfterBroadening_m
                call Samrsl_0               
             end if
 
+            ! ensure that doppler broadening is reset to the user desired type
+            dopplerOption%bType =  btypeSave
+
             Jwwwww = 4
 
             call normAndBack(moreBroadening, resolutionOption%resBroad2)
@@ -460,7 +471,9 @@ end module sumAndConvertAfterBroadening_m
 
 
       ! final output
-      call Samint_0 (.false., Yssmsc, Jwwwww) ! which will set the correct values for where to next
+      ! TODO: check whether we should set Jjjdop to 0 to make sure  we print
+      !       final output. Tr058 will fail if that is done
+      call Samint_0 (.false., Yssmsc, Jwwwww, Jjjdop) ! which will set the correct values for where to next
    end subroutine
 
 
diff --git a/sammy/src/dat/mdat0.f90 b/sammy/src/dat/mdat0.f90
index 375fe6703f117733ddde68267e11631e20c8e9f6..d4e959f774131e9741e98cb1b0afd491eea2a2a4 100644
--- a/sammy/src/dat/mdat0.f90
+++ b/sammy/src/dat/mdat0.f90
@@ -521,7 +521,7 @@ contains
     !
     use oops_common_m, only : Msize
     use fixedi_m, only : Krefit, Nres, Ntotc, Numdtp
-    use ifwrit_m, only : Iptdop, Iptwid, Jjjdop, Kartgd, Kdatv, Ktzero, Ndat, Ndatb
+    use ifwrit_m, only : Iptdop, Iptwid, Kartgd, Kdatv, Ktzero, Ndat, Ndatb
     use fixedr_m
     use broad_common_m, only :  dopplerOption
     IMPLICIT none
@@ -602,7 +602,6 @@ contains
           END IF
           ! ###       four ###
           K4 = Ndatbm
-          Jjjdop = 1
           !
        END IF
        I = Idimen (I, -1, 'I, -1')
diff --git a/sammy/src/dat/mdat7.f90 b/sammy/src/dat/mdat7.f90
index 84af315fdc9c3782c4faa27b31d9dcdd2bbc5442..2a66261933f0b45c1594d80210bacc9fdb99685e 100644
--- a/sammy/src/dat/mdat7.f90
+++ b/sammy/src/dat/mdat7.f90
@@ -22,21 +22,16 @@ module mdat7_M
       use dop2_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 !
-     type(GridData)::coefGrid
      real(kind=8),allocatable,dimension(:)::Energb
-!
-!   delt => delttt from fixedr_m (Ff(23) assigned to delt 
-!   dpdelt => gammmm from fixedr_m (Ff(27)
-!   dpdelv => delvvv  from fixedr_m (Ff(25) assigned t   
-!      
-!      EQUIVALENCE (Delttt,Delt), (Gammmm,Dpdelt), (Delvvv,Dpdelv)
+
 !   Todo: Retrieve Boltzman constant from the SAMMY constant module.
       DATA Xkb /8.617065D-05/, Zero /0.0d0/
 !
-      double precision, pointer::delt => Ff(23)
+
       double precision, pointer::Dpdelt => Ff(27)
       double precision, pointer::Dpdelv => Ff(25)
       
+      if (.not.associated(dopplerOption%lealHwang)) return
       
       Emid = (Emin+Emax)*0.5d0
       Vvmin = DSQRT(Eminr)
@@ -44,15 +39,17 @@ module mdat7_M
       Xvmin = DSQRT(Emin)
       Xvmax = DSQRT(Emax)
 !
-      Delt = Delttt
-      IF (Delt.EQ.0.0D0) Delt = 5.0D0
+      delT = dopplerOption%lealHwang%Delt
+      if (delT.eq.0.0d0) Delt = 5.0d0
+!
 ! *** a reasonable number ?  probably will want to change eventually
 !
       Itime = 0
-      AJ = Temp/Delt
+      AJ = Temp/delT
       J = AJ
       IF (J.LE.10) J = 10
       Delt = Temp/Dfloat(J)
+      dopplerOption%lealHwang%Delt = Delt
 ! *** adjust delt slightly so that "j" steps gets us exactly to Temp
 !
       Dpdelt = (Delt*Xkb)/(4.0D0*Aaawww)
@@ -70,22 +67,22 @@ module mdat7_M
               /,'           xvmin, xvmax =', 23X, 3(1pg14.5))
 !
 ! 
-      Kjjjjj = J
-      K2pls1 = J + J + 1
+      dopplerOption%lealHwang%Kjjjjj = J
+      dopplerOption%wantBroaden = .true.
+      dopplerOption%lealHwang%K2pls1 = J + J + 1
 !
 !
 ! *** Coefgn generates coefficients Coef (which we put into Energb
 ! ***   for now) for solving heat equation for performing Doppler
-! ***   broadening
-      call coefGrid%initialize()
-      CALL Coefgn (coefGrid, Kjjjjj, K2pls1, Mkkkkk, M2pls1, 0)
-      call reallocate_real_data(Energb, coefGrid%getLength())
+! ***   broadening     
+      CALL Coefgn (dopplerOption%lealHwang%coefGrid, dopplerOption%lealHwang%Kjjjjj,  Mkkkkk, M2pls1)
+      call reallocate_real_data(Energb, dopplerOption%lealHwang%coefGrid%getLength())
 
 ! *** Mkkkkk is the number actually used, since the rest are
 ! ***      effectively zero
 !
-      WRITE (21,10002) Mkkkkk, Kjjjjj
-      WRITE ( 6,10002) Mkkkkk, Kjjjjj
+      WRITE (21,10002) Mkkkkk, dopplerOption%lealHwang%Kjjjjj
+      WRITE ( 6,10002) Mkkkkk, dopplerOption%lealHwang%Kjjjjj
 10002 FORMAT (' New and actual values of J = ', 2I8)
 !
       N = (Vmax-Vvmin)/Dpdelv
@@ -96,7 +93,6 @@ module mdat7_M
       IF (Ndatb.GT.Kdatbm) STOP '[STOP in Llgrid in dat/mdat7.f]'
 !
 !      Delvvv = Dpdelv
-!      Delttt = Dpdelt
 !      Gammmm = Dpdelt
 !
 !
@@ -116,28 +112,25 @@ module mdat7_M
 !
          IF (Ngtvv.NE.0) THEN
             DO I=1,Ngtvv
-               call coefGrid%addData(I, 1, -(Vvmin + Dpdelv*Dfloat(I)) **2)
+               call dopplerOption%lealHwang%coefGrid%addData(I, 1, -(Vvmin + Dpdelv*Dfloat(I)) **2)
             END DO
-            call coefGrid%addData(Ngtvv+1, 1, Zero)
+            call dopplerOption%lealHwang%coefGrid%addData(Ngtvv+1, 1, Zero)
             Ngtvv = Ngtvv + 1
             Nnn = Ngtvv + 1
          END IF
       END IF
 !
       DO I=Nnn,Ndatb
-          call coefGrid%addData(I, 1, (Vvmin + Dpdelv*Dfloat(I)) **2)
+          call dopplerOption%lealHwang%coefGrid%addData(I, 1, (Vvmin + Dpdelv*Dfloat(I)) **2)
       END DO
 
 
-      call reallocate_real_data(Energb, coefGrid%getLength())
-      do i = 1, coefGrid%getLength()
-           Energb(i) = coefGrid%getData(i,1)
+      call reallocate_real_data(Energb, dopplerOption%lealHwang%coefGrid%getLength())
+      do i = 1, dopplerOption%lealHwang%coefGrid%getLength()
+           Energb(i) = dopplerOption%lealHwang%coefGrid%getData(i,1)
       end do
-      call coefGrid%destroy()
-
-
 !
-      Nwd = K2pls1
+      Nwd = dopplerOption%lealHwang%K2pls1
       Awid = Dpdelv*(2.0D0*DSQRT(Emid+Dpdelv))
       RETURN
 !
diff --git a/sammy/src/dex/mdex0.f b/sammy/src/dex/mdex0.f
index e1c699b79a05c7ad349ece7ca7aabc7b7801539d..85e23b1a8d268a742cb7f866df7e053d3c4ea627 100644
--- a/sammy/src/dex/mdex0.f
+++ b/sammy/src/dex/mdex0.f
@@ -6,7 +6,6 @@ C      use oops_common_m
       use fixedi_m, only : K2reso,
      *                     Numorr, Numrpi,
      *                     Nudwhi
-      use ifwrit_m, only : Jjjdop
       use exploc_common_m
       use array_sizes_common_m
       use oopsch_common_m, only : Nowwww, Segmen
@@ -55,7 +54,6 @@ C
       deallocate(A_Iwts)
 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 C
-      Jjjdop = 0
       CALL Write_Commons_Many
       RETURN
 C
diff --git a/sammy/src/dop/LealHwangBroadening_M.f90 b/sammy/src/dop/LealHwangBroadening_M.f90
index e9c4279309ce3dc5df2123251e89e47b0aacb6a7..4e6caf0f452d90a8bcf83ca76690b900a86ceb7b 100644
--- a/sammy/src/dop/LealHwangBroadening_M.f90
+++ b/sammy/src/dop/LealHwangBroadening_M.f90
@@ -7,8 +7,11 @@ use, intrinsic :: ISO_C_BINDING
 
 
 implicit none
-
 type, extends(FreeGasDopplerBroadening) :: LealHwangBroadening
+type(GridData)::coefGrid
+integer::Kjjjjj
+integer::K2pls1
+real(kind=8)::delT
 contains
 procedure, pass(this) :: initialize => LealHwangBroadening_initialize
 procedure, pass(this) :: broaden => LealHwangBroadening_broaden
@@ -23,16 +26,20 @@ subroutine LealHwangBroadening_initialize(this, hand, list, work)
     class(GridDataList)::list
     class(GridDataList)::work
     call FreeGasDopplerBroadening_initialize(this, hand, list, work)
+    call this%coefGrid%initialize()
+    this%Kjjjjj = 0
+    this%K2pls1 = 0
+    this%delT = 5.0d0
 end subroutine
 
 subroutine LealHwangBroadening_destroy(this)
     implicit none
     class(LealHwangBroadening) :: this
     call FreeGasDopplerBroadening_destroy(this)
+    call this%coefGrid%destroy()
 end subroutine
 
 
-
 subroutine LealHwangBroadening_broaden(this)
    class(LealHwangBroadening) :: this
 
@@ -44,6 +51,8 @@ subroutine LealHwangBroadening_broaden(this)
    call this%getData(data)
    call data%nullify()
    call data%reserve(ndatb*data%getNnnsig(), this%numPar+1)
+
+   call this%coefGrid%nullify()
 end subroutine
 
 end module LealHwangBroadening_M
diff --git a/sammy/src/dop/mdop0.f90 b/sammy/src/dop/mdop0.f90
index 7e5576c4e7cd58372e856f10acab4affe12cf1d9..536a31603eb0f8077e1ce045f04e2bf935efcaef 100644
--- a/sammy/src/dop/mdop0.f90
+++ b/sammy/src/dop/mdop0.f90
@@ -1,77 +1,50 @@
 !
 module dop_m
+  use LealHwangBroadening_M
   contains
 !
-      SUBROUTINE Samdop_0
-!
-      use over_common_m
-      use oops_common_m
-      use fixedi_m
-      use ifwrit_m
-      use exploc_common_m
-      use array_sizes_common_m
-      use oopsch_common_m
-      use cbro_common_m
-      use lbro_common_m
-      use EndfData_common_m
-      use GridData_M
+      SUBROUTINE Samdop_0(broadener)
+      use fixedi_m, only : Lllmax, numUsedPar
+      use ifwrit_m, only : Kcros, Kdebug, Ksitmp, Ksindi
+      use exploc_common_m, only : I_Iflmsc
+      use array_sizes_common_m, only : calcDataSelf
       use dop1_m
       use dop2_m
-      use AllocateFunctions_m
-      use rsl7_m
-      use mxct27_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      IMPLICIT None
 !
       type(GridData)::grid, auxGrid
-      type(GridData)::coefGrid
-      real(kind=8),allocatable,dimension(:)::A_Icoef
+      class(LealHwangBroadening)::broadener
+      integer::I2pls1, Mmmmmm
 
 !
       WRITE (6,99999)
 99999 FORMAT (' *** SAMMY-DOPPL  8 Aug 07 ***')
-      Segmen(1) = 'D'
-      Segmen(2) = 'O'
-      Segmen(3) = 'P'
-      Nowwww = 0
+
 !
 !x      CALL Initix
 !
 !
-! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
-      N = K2pls1 + 1
-      call coefGrid%initialize()
+      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()
+
 ! *** Coefgn generates Coefficients Coef for solving heat equation
 ! ***   for performing Doppler broadening
-      CALL Coefgn (coefGrid, Kjjjjj, K2pls1, Mmmmmm, I2pls1, 0)
+      CALL Coefgn (broadener%coefGrid, broadener%Kjjjjj, Mmmmmm, I2pls1)
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
-      call allocate_real_data(A_Icoef, I2pls1)
-      do i = 1, I2pls1
-        A_Icoef(i) = coefGrid%getData(i,1)
-      end do
-      call coefGrid%destroy()
-      I = Idimen (0, 0, '0, 0')
-!
-!
-      Nv = Ndatb
-      Left  = Msize - Icoef - 2*I2pls1
-      Ldatb = Ndatb
-
-      ! ensure we have an auxillary energy grid     
-      if ( expData%getLength().lt.2) then
-         write(0,*)" Missing auxillary grid"
-         stop
-      end if
 !
 !
 ! *** Dopplh performs Doppler broadening operation via Leal-Hwang method
-      CALL Dopplh(I_Iflmsc , A_Icoef  ,              &
-       Mmmmmm   , I2pls1   , Ngb      , Nss, Ldatb,  &
-       calcData, calcDataSelf)
-      deallocate(A_Icoef)
+      CALL Dopplh(broadener, I_Iflmsc ,  Mmmmmm   , I2pls1)
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 !
-      Jjjdop = 1
       CALL Write_Commons_Many
       RETURN
 !
diff --git a/sammy/src/dop/mdop1.f90 b/sammy/src/dop/mdop1.f90
index 959b6e27b9a95299fa07f263a0edb7dcb041b0e9..1f24a4717af6a474ca847786afa42ccb295d3664 100644
--- a/sammy/src/dop/mdop1.f90
+++ b/sammy/src/dop/mdop1.f90
@@ -1,64 +1,48 @@
 !
-module dop1_m
-  use convert_to_transmission_m
+module dop1_m  
+  use LealHwangBroadening_M
+  use ifwrit_m, only : Kvtemp
+  use AuxGridHelper_M, only : setAuxGridOffset
+  use DerivativeHandler_M
+  IMPLICIT None
+
+  public Dopplh
   contains
 !
 ! ____________________________________________________________________
 !
-      SUBROUTINE Dopplh (Iflmsc, Coef  ,                         &
-         Mmmmmm, M2pls1, Ngb   , Nss   , Ldatb,                  &
-         derivs, derivsSelf)
+      SUBROUTINE Dopplh (calc, Iflmsc,  Mmmmmm, M2pls1)
 !
 ! *** Purpose -- Implement Doppler broadening
 !
-      use oops_common_m
-      use fixedi_m
-      use ifwrit_m
-      use samxxx_common_m
-      use oopsch_common_m
-      use fixedr_m
-      use abro_common_m
-      use lbro_common_m
-      use hhhhhh_common_m
-      use xct2_m
-      use mxct27_m
-      use EndfData_common_m
-      use SammyGridAccess_M
-      use AuxGridHelper_M
-      use DerivativeHandler_M
-      use broad_common_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      type(SammyGridAccess)::grid
-      LOGICAL Need_isotopes
-      LOGICAL Another_Process_will_Happen
+      class(LealHwangBroadening)::calc
       integer::nauxStart
-      type(DerivativeHandler)::derivs, derivsSelf
-!
-      DIMENSION Iflmsc(*), Coef(*)
-!     DIMENSION Iflmsc(Nummsc), Coef(K2pls1)
-!
+      type(DerivativeHandler)::derivs
+!
+      integer::Iflmsc(*)
+      integer::Mmmmmm, M2pls1, Ldatb
+      real(kind=8)::Cc, Ebm, Coef
+      integer::I, Ii, Imax, Imin, Isox, K
+      integer::Kkkkkk, Kkkmin, Ldatb_M2, M, Mmax, Mmin, Mmm
+      integer::Newmin
+
+!     DIMENSION Iflmsc(Nummsc)
 !
-      Need_Isotopes = Yssmsc
-      Another_Process_will_Happen = Yresol.OR.Yssmsc.or.yaverg
 !
-      call grid%initialize()      
-      call grid%setToAuxGrid(expData)
-      call setAuxGridOffset(1) ! reset auxiallary grid offset
+      call calc%getData(derivs)
+
       call derivs%nullify()
+      Ldatb = calc%getNumEnergyUnbroadened()
 !
-      K2pls1 = M2pls1
-      Nnn = Ngb
+      calc%K2pls1 = M2pls1
       Imin = 1
-      Imax = Ldatb - K2pls1 + 1
+      Imax = Ldatb - M2pls1 + 1
       Mmin = 1
-!
-!
-      Newkkk = Kksave
+
+
 !
       Ldatb_M2 = Ldatb - Mmmmmm*2
-      Iw = 0
-      IF (Ksindi.NE.0) Iw = 1
-      call derivsSelf%nullify()
+      if (calc%hasSelf) call calc%dataSelf%nullify()
 !
 !
 !
@@ -66,30 +50,31 @@ module dop1_m
 !
       Mmax = Ldatb
 !
-      CALL Set_Negative_E (derivs, derivsSelf)
+      CALL Set_Negative_E (calc)
 !
-      DO Iso=1,derivs%getUsedIsotopes()
-         Isox = Iso
+      DO Isox=1,derivs%getUsedIsotopes()
          Mmm = Mmmmmm
          Kkkmin = 0
          Kkkkkk = 0
          nauxStart = 0
 ! ***    Loop over Energies
+
          DO I=1,Imax
             Mmm = Mmm + 1
 !
-            call derivs%reserveColumnsNs(Kkkkkk+1, derivs%getNnnsig(),numUsedPar+1)
-            IF (Ksindi.NE.0) THEN
-               call derivsSelf%reserveColumnsNs(Kkkkkk+1, derivsSelf%getNnnsig(),numUsedPar+1)
+            call derivs%reserveColumnsNs(Kkkkkk+1, derivs%getNnnsig(),calc%numPar+1)
+            IF (calc%hasSelf) THEN
+               call calc%dataSelf%reserveColumnsNs(Kkkkkk+1, calc%dataSelf%getNnnsig(),calc%numPar+1)
             END IF
 !
             Ii = I
             M  = I
 ! ***       Do the broadening integration
-            DO K=1,K2pls1
-               IF (Coef(K).NE.0.0d0.and.m.le.Mmax) THEN
-                  Cc = Coef(K) * grid%getEnergy(M, expData)
-                  CALL Multiply_by_Coef (Kkkkkk+1, derivs, derivsSelf, Cc, M, Isox)
+            DO K=1,calc%K2pls1
+               Coef = calc%coefGrid%getData(K,1)
+               IF (Coef.NE.0.0d0.and.m.le.Mmax) THEN
+                  Cc = Coef * calc%getEnergyUnbroadened(M)
+                  CALL Multiply_by_Coef (calc, Kkkkkk+1,  Cc, M, Isox)
                END IF
                M = M + 1
                IF (M.GT.Mmax) GO TO 195
@@ -97,12 +82,12 @@ module dop1_m
 !
 !
             IF (Kkkmin.EQ.0) Kkkmin = Mmm
-            Ebm = grid%getEnergy(Mmm, expData)
+            Ebm = calc%getEnergyUnbroadened(Mmm)
 
             Kkkkkk = Kkkkkk + 1
 
 ! ******    transform back to cross section (from E*Sigma)
-            CALL Transform_to_Cross ( derivs, derivsSelf, Kkkkkk, Ebm, Isox)
+            CALL Transform_to_Cross ( calc, Kkkkkk, Ebm, Isox)
 
          END DO
          GO TO 196
@@ -110,19 +95,28 @@ module dop1_m
             Imax = I
   196    CONTINUE
 !
-         IF (Ksolve.NE.2 .AND. Kvtemp.GT.0) THEN
+         IF (calc%numPar.gt.0 .AND. Kvtemp.GT.0) THEN
 ! ***       generate derivatives wrt temperature
-            CALL Temp_Deriv ( derivs, derivsSelf, Imax)
+            CALL Temp_Deriv ( calc, Imax)
          END IF
 !
       END DO
 !
-      call setAuxGridRowMax(Kkkkkk)
       nauxStart = Kkkmin
-      call dopplerOption%lealHwang%updateBroadenedOffset(nauxStart)
-      call dopplerOption%lealHwang%setLength(Kkkkkk)
-!
-      call grid%destroy()
+      ! If the broadenened data are store on the experimental
+      ! grid, the function calc%updateBroadenedOffset
+      ! does not have any effect (as we alwayas make
+      ! sure that all experimental data are calcuated).
+      ! Broadened data are on the experimental
+      ! grid if there is no more broadening
+      call calc%updateBroadenedOffset(nauxStart)
+      call calc%setLength(Kkkkkk)
+!
+      ! this function set the grid offset on the auxillary grid
+      ! regardless of whether the final broadened data
+      ! are on the experimental grid.
+      ! It is used if printing the interpolated cross section
+      ! in SAMMY output
       call setAuxGridOffset(nauxStart)
       RETURN
 !
@@ -131,35 +125,26 @@ module dop1_m
 !
 ! ____________________________________________________________________
 !
-      SUBROUTINE Set_Negative_E (derivs, derivsSelf)
+      SUBROUTINE Set_Negative_E (calc)
 ! *** Purpose -- generate the data for negative-velocity points
 ! *** NOTE    -- Integrand requires " sigma(E<0) times (sqrt(-E))**2 " for
 ! ***            which "sigma(E<0)" is negative, and "(sqrt(-E))**2" is
 ! ***            positive.  Since we Store E as negative, we'll instead
 ! ***            store "sigma(E<0)" as postive, and use "E" (which is
-! ***            negative) instead of "(sqrt(-E))**2" (which is positive).
-      use fixedi_m
-      use ifwrit_m
-      use EndfData_common_m
-      use DerivativeHandler_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      type(GridData)::grid
+! ***            negative) instead of "(sqrt(-E))**2" (which is positive).      
       real(kind=8)::ee, val
-      integer::j1,j2,ii, nLength
-      type(DerivativeHandler)::derivs, derivsSelf
-!
-      call expData%getGrid(grid, 2)  ! we ensured that it exists
-
+      class(LealHwangBroadening)::calc
+      type(DerivativeHandler)::derivs
+      integer::Ipar, Iso, Jj, Ll, N
 
-      ll = grid%getLength()
-      nLength = grid%getNumPerEnergy()
-      if (nLength.gt.1) ll = ll/nLength
+      ll = calc%getNumEnergyUnbroadened()
+      call calc%getData(derivs)
       DO Jj=1,ll
-         ee = grid%getData((Jj-1)*nLength+1, 1)
+         ee = calc%getEnergyUnbroadened(jj)
          if (ee.ge.0.0d0) exit
          DO Iso=1,derivs%getUsedIsotopes()
             ! set negative cross section to positive
-            DO Ipar=0,numUsedPar
+            DO Ipar=0,calc%numPar
                do  N=1, derivs%getNnnsig()
                  val = derivs%getDataNsOld(Jj, N, Ipar, iso)
                   if (val.ne.0.0d0) then
@@ -168,11 +153,11 @@ module dop1_m
                   end if
                 end do
 
-                IF (Ksindi.NE.0.and.N.eq.1) THEN  ! self indicated cross section only has 1 section
-                    val = derivsSelf%getDataNsOld(Jj, 1, Ipar, iso)
+                IF (calc%hasSelf) THEN  ! self indicated cross section only has 1 section
+                    val = calc%dataSelf%getDataNsOld(Jj, 1, Ipar, iso)
                     if (val.ne.0.0d0) then
                         val = -1.0d0*val
-                        call derivsSelf%addDataNsOld(Jj, 1, Ipar, iso, val)
+                        call calc%dataSelf%addDataNsOld(Jj, 1, Ipar, iso, val)
                     end if
                 end if
             END DO
@@ -185,35 +170,35 @@ module dop1_m
 !
 ! ____________________________________________________________________
 !
-      SUBROUTINE Multiply_By_Coef (irow, derivs, derivsSelf, Cc, M, Isox)
-      use fixedi_m
-      use ifwrit_m
-      use DerivativeHandler_M      
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      type(DerivativeHandler)::derivs, derivsSelf
+      SUBROUTINE Multiply_By_Coef (calc, irow, Cc, Mm, Isox)
+      class(LealHwangBroadening)::calc
+      type(DerivativeHandler)::derivs
+      real(kind=8)::Cc
+      integer::Mm, Isox
       integer::irow
+      integer::Ipar, N
+      real(kind=8)::val
 !
-      Mm = M
-      Iso = Isox
+      call calc%getData(derivs)
       DO N=1,derivs%getNnnsig()
          ! update cross section (ipar=0) as well as derivatives
-         DO Ipar=0,numUsedPar
-           val = derivs%getDataNsOld(Mm, n, Ipar, Iso)*Cc
-           val = val + derivs%getDataNs(irow, N, Ipar, Iso)
+         DO Ipar=0,calc%numPar
+           val = derivs%getDataNsOld(Mm, n, Ipar, Isox)*Cc
+           val = val + derivs%getDataNs(irow, N, Ipar, Isox)
            if (val.eq.0.0d0) cycle
-           call derivs%addDataNs(irow, n, Ipar, Iso, val)
+           call derivs%addDataNs(irow, n, Ipar, Isox, val)
          end do
       END DO
 !
-      IF (Ksindi.NE.0) THEN
+      IF (calc%hasSelf) THEN
           ! update cross section (ipar=0) as well as derivatives
-          DO N=1,derivsSelf%getNnnsig()
+          DO N=1,calc%dataSelf%getNnnsig()
              ! update cross section (ipar=0) as well as derivatives
-             DO Ipar=0,numUsedPar
-               val = derivsSelf%getDataNsOld(Mm, n, Ipar, Iso)*Cc
-               val = val + derivsSelf%getDataNs(irow, N, Ipar, Iso)
+             DO Ipar=0,calc%numPar
+               val = calc%dataSelf%getDataNsOld(Mm, n, Ipar, Isox)*Cc
+               val = val +calc%dataSelf%getDataNs(irow, N, Ipar, Isox)
                if (val.eq.0.0d0) cycle
-               call derivsSelf%addDataNs(irow, n, Ipar, Iso, val)
+               call calc%dataSelf%addDataNs(irow, n, Ipar, Isox, val)
              end do
           END DO
       END IF
@@ -224,14 +209,9 @@ module dop1_m
 !
 ! ____________________________________________________________________
 !
-      Subroutine Temp_Deriv ( derivs, derivsSelf, Imax)
-      use fixedi_m
-      use ifwrit_m
-      use abro_common_m
-      use DerivativeHandler_M
-      IMPLICIT None
-!
-      type(DerivativeHandler)::derivs, derivsSelf
+      Subroutine Temp_Deriv ( calc, Imax)
+      class(LealHwangBroadening)::calc
+      type(DerivativeHandler)::derivs
       integer::Imax
       integer::Kvt
       real(kind=8)::val1, val2, val3, resVal, Sig0,delt6
@@ -239,7 +219,8 @@ module dop1_m
 
       Sig0 = 0.0d0
       Kvt = Kvtemp
-      delt6 = delt * 6.0d0
+      delt6 = calc%delT * 6.0d0
+      call calc%getData(derivs)
       DO Iso=1,derivs%getUsedIsotopes()
          DO N=1,derivs%getNnnsig()
             val1 = derivs%getDataNs(1, 1, 0, iso)
@@ -260,20 +241,20 @@ module dop1_m
          END DO
       END DO
 !
-      IF (Ksindi.NE.0) THEN
+      IF (calc%hasSelf) THEN
          DO Iso=1,derivs%getUsedIsotopes()
-            val1 = derivsSelf%getData(1,0,Iso)
-            val2 = derivsSelf%getData(2,0,Iso)
+            val1 = calc%dataSelf%getData(1,0,Iso)
+            val2 = calc%dataSelf%getData(2,0,Iso)
             resVal = (Sig0 -  val1*2.0d0 + val2) / Delt6
-            call derivsSelf%addData(1, Kvt, Iso, resVal)
+            call calc%dataSelf%addData(1, Kvt, Iso, resVal)
          END DO
          DO I=2,Imax-1
             DO Iso=1,derivs%getUsedIsotopes()
-               val1 = derivsSelf%getData(I-1, 0, Iso)
-               val2 = derivsSelf%getData(I,   0, Iso)
-               val3 = derivsSelf%getData(I+1, 0, Iso)
+               val1 = calc%dataSelf%getData(I-1, 0, Iso)
+               val2 = calc%dataSelf%getData(I,   0, Iso)
+               val3 = calc%dataSelf%getData(I+1, 0, Iso)
                resVal = (val1 -  val2*2.0d0 + val3) / Delt6
-               call derivsSelf%addData(1, Kvt, Iso, resVal) ! Todo: Is this right? This is what it was
+               call calc%dataSelf%addData(1, Kvt, Iso, resVal) ! Todo: Is this right? This is what it was
             END DO
          END DO
       END IF
@@ -283,17 +264,17 @@ module dop1_m
 !
 ! ____________________________________________________________________
 !
-      SUBROUTINE Transform_to_Cross ( derivs, derivsSelf, irow,  Ebm, Isox)
-      use fixedi_m
-      use ifwrit_m
-      use DerivativeHandler_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      type(DerivativeHandler)::derivs, derivsSelf
+      SUBROUTINE Transform_to_Cross (calc, irow,  Ebm, Isox)
+      class(LealHwangBroadening)::calc
+      type(DerivativeHandler)::derivs
+      integer::irow, Isox
+      real(kind=8)::Ebm
       real(kind=8)::val
+      integer::Ipar, Iso, N
 
-      Iso = Isox
+      call calc%getData(derivs)
       DO N=1,derivs%getNnnsig()
-         do ipar = 0, numUsedPar
+         do ipar = 0, calc%numPar
            val = derivs%getDataNs(irow, N, Ipar, isox)
            if( val.eq.0.0d0) cycle
            val = val/Ebm
@@ -301,13 +282,13 @@ module dop1_m
          end do
       END DO
 !
-      IF (Ksindi.NE.0) THEN         
-         DO N=1,derivsSelf%getNnnsig()
-            do ipar = 0, numUsedPar
-              val = derivsSelf%getDataNs(irow, N, Ipar, isox)
+      IF (calc%hasSelf) THEN
+         DO N=1,calc%dataSelf%getNnnsig()
+            do ipar = 0, calc%numPar
+              val = calc%dataSelf%getDataNs(irow, N, Ipar, isox)
               if( val.eq.0.0d0) cycle
               val = val/Ebm
-              call derivsSelf%addDataNs(irow, N, Ipar, isox, val)
+              call calc%dataSelf%addDataNs(irow, N, Ipar, isox, val)
             end do
          END DO
       END IF
diff --git a/sammy/src/dop/mdop2.f90 b/sammy/src/dop/mdop2.f90
index 290e00f2d4ef8807015e91bdc949ea2ebd32d718..9407a07f2aa59fbc7d463c2b7da48142f5b469ec 100755
--- a/sammy/src/dop/mdop2.f90
+++ b/sammy/src/dop/mdop2.f90
@@ -1,21 +1,28 @@
 module dop2_m
   use GridData_M
+  implicit none
+
+  public Coefgn
 
   contains
 !
 !
 ! __________________________________________________________________
 !
-      SUBROUTINE Coefgn (Coef, Jjjjjj, J2pls1, Mmmmmm, M2pls1, Kkkxxx)
+      SUBROUTINE Coefgn (Coef, Jjjjjj, Mmmmmm, M2pls1)
 !
 ! *** revised extensively July 1996 to increase speed & accuracy
 !
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-       type(GridData)::Coef
-      DATA Two /2.0d0/, Four /4.0D0/
+      type(GridData)::Coef
+      integer:: Mmmmmm
+      integer::Jjjjjj,  M2pls1
+      real(kind=8),parameter::Two=2.0d0, Four=4.0D0
+      real(kind=8)::Ans, C, Cay, Cay2, Fac
+      real(kind=8)::Ratio, Ratio2, tmp, X, Y
+      integer::I, Ii, Itimes, M, Max, Mj, Mm, Mstart, Mstmin, Mstop
+      integer::Mstopx, N, Nj, Np1, Mwhich
 !
       N = Jjjjjj
-      J2pls1 = 2*N + 1
 !
       Cay = Four
       Cay2 = Cay + Two
@@ -93,12 +100,7 @@ module dop2_m
    70 CONTINUE
       Np1 = Mmmmmm + 1
       M2pls1 = Np1 + Mmmmmm
-!
-      IF (Kkkxxx.EQ.1) THEN
-         WRITE (06,10000) Jjjjjj, Mmmmmm, Mstmin, Mstopx
-10000    FORMAT (' J, M, Mstart, Mstop =', 4I10)
-      RETURN
-      END IF
+
 !
 ! *** rearrange storage of Coef ... 1 to m+1 now goes to m+1 to 2m+1
       C = Coef%getData(Np1, 1)
@@ -134,9 +136,12 @@ module dop2_m
 ! *** method:  directly if N < 128
 ! ***          via asymptotic expAnsion for ln(N) if N > 128
 !
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      DATA Dln2pi/1.8378770664125945d0/, One /1.0d0/
-!                / ln ( 2 pi)
+      real(kind=8)::Cay, Cay2, Ratio, Ratio2, Ans
+      integer::N, J, M
+      real(kind=8),parameter::Dln2pi=1.8378770664125945d0, One=1.0d0
+      !                / ln ( 2 pi)
+      integer::Jj, L, Mj, Mm, Nm, Nn
+      real(kind=8)::X,Y
 !
       X = dFLOAT(N+3)**2 + 3.0d0*dFLOAT(J)**2
       X = 2.0d0*dSQRT(X) - dFLOAT(N+3*J+9)
diff --git a/sammy/src/end/msamxx.f b/sammy/src/end/msamxx.f
index bfc06b4b1e575f7da2cc7d33e16e257a11454a63..c32d28afda0e2588ccdb4fce3fc528b4fea6550d 100755
--- a/sammy/src/end/msamxx.f
+++ b/sammy/src/end/msamxx.f
@@ -111,7 +111,7 @@ C
       ELSE IF (Ix.EQ.2) THEN
 C ***    For call from end of SAMMAS, or from DOPUSHx or others
          Tabsol = Tiniti
-C
+
       ELSE
 C     ELSE IF (Ix.GT.3 .OR. Ix.LT.0)
          STOP '[STOP in Timer in end/msamxx.f]'
@@ -123,8 +123,8 @@ C
       IF (Ix.NE.3) THEN
 C
 cq         WRITE (21,99999) (Segmen(I),I=1,3), Trel
-         WRITE ( line,99999) (Segmen(I),I=1,3), Trel
-99999    FORMAT (40X, 'CPU for SAMMY-', 3A1, ' =', 0PF9.2, ' sec')
+         WRITE ( line,99999)  Trel
+99999    FORMAT (40X, 'CPU for this section is ', 0PF9.2, ' sec')
          call printStdoutData(line)
 
 C
diff --git a/sammy/src/fin/mfin1.f90 b/sammy/src/fin/mfin1.f90
index 894b4f4a03959c51186e2f01b525656f72440fa1..ee9f28ff39fdb80ec500780a7dad553e711d8f8e 100644
--- a/sammy/src/fin/mfin1.f90
+++ b/sammy/src/fin/mfin1.f90
@@ -954,6 +954,7 @@ module fin1
                   Temp = val
                   if (associated(dopplerOption%broadener)) then
                      call dopplerOption%broadener%setTemperature(temp)
+                     dopplerOption%wantBroaden = .true.
                   end if
                ELSE IF (I.EQ.3) THEN
                   Thick = val
diff --git a/sammy/src/inp/minp06.f b/sammy/src/inp/minp06.f
index 4931451bc6884ade4e75f1934597863c1300c81b..590b96b328daae7f12a9c33b1d9f588e249deb55 100644
--- a/sammy/src/inp/minp06.f
+++ b/sammy/src/inp/minp06.f
@@ -271,7 +271,7 @@ C     DIMENSION Parbrd(Numbrd), Iflbrd(Numbrd), Bcf(Ncf), Cf2(Ncf), Cf(Ncf)
       integer::Iflbrd(*)
       real(kind=8)::Eeemin, Aw, Deltab
       real(kind=8),allocatable,dimension(:)::Cf
-      real(kind=8)::Elowbr
+      real(kind=8)::Elowbr, Delttt
 C
       real(kind=8),parameter::Zero=0.0d0
       real(kind=8)::A, B, Ao2x, Ax, Dist1, Do2x, E1, Gaussian, Xm
@@ -370,9 +370,11 @@ C
       Ao2 = Zero
       Bo2 = Zero
       Co2 = Zero
+      dopplerOption%wantBroaden = .false.
 C
       IF (Temp.NE.Zero .AND. Aw.NE.Zero) THEN
          Dopple = dSQRT(Boltzm*Aneutr*Temp/Aw)
+         dopplerOption%wantBroaden = .true.
       END IF
 C
 C
@@ -447,6 +449,9 @@ C ***    changed from 4.0d0 June 1997
          dopplerOption%freeGas%Elowbr = Elowbr
       end if
       call dopplerOption%broadener%setTemperature(temp)
+      if (associated(dopplerOption%lealHwang) ) then
+          dopplerOption%lealHwang%delT = Delttt
+      end if
 C
 C     SET STARTING VALUE FOR  EST IN SHFTGE ROUTINE
 C
@@ -489,6 +494,7 @@ C *** stuff on CARD SET 5
       Co2    = Zero
       Do2    = Zero
       Dopple = Zero
+      dopplerOption%wantBroaden = .false.
       Nobrd = 1
       RETURN
       END
diff --git a/sammy/src/int/mint0.f b/sammy/src/int/mint0.f
index b3d70340552c54354da034b44f64bd390d72d511..148df684ebb8d13f5ca2ce873ebd42947c9f7879 100644
--- a/sammy/src/int/mint0.f
+++ b/sammy/src/int/mint0.f
@@ -1,6 +1,6 @@
 C
 C
-      SUBROUTINE Samint_0(debugOut, haveMultScatter, Jwwwww)
+      SUBROUTINE Samint_0(debugOut, haveMultScatter, Jwwwww, Jjjdop)
 C
 C *** Purpose -- Print cross sections etc
 C
@@ -20,7 +20,7 @@ C
       CHARACTER*10 Sss
       character(len=80)::line
       logical::debugOut, haveMultScatter
-      integer::Jwwwww
+      integer::Jwwwww, Jjjdop
 C
 C
           WRITE (line, 99999)
@@ -35,7 +35,7 @@ C
 C
 C *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMINT
       CALL Estint (Nblmax, Iv, debugOut, haveMultScatter)
-      I = Idimen (Istarting, -1, 'Istarting, -1')
+      I = Idimen (Istarting, -1, 'Istarting, -1')    
 C
 C
 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
diff --git a/sammy/src/mso/mmso0.f b/sammy/src/mso/mmso0.f
index d1b0221484bcf1dfb613cf1b107a9e487e3df9bc..1fd1aed0c283e69f2146c275fd624fa26ae5b20d 100644
--- a/sammy/src/mso/mmso0.f
+++ b/sammy/src/mso/mmso0.f
@@ -13,7 +13,7 @@ C
      *                     Kkkrsl, Ktheta,
      *                     Ntheta, Nudwhi, Numrpi, Numorr,
      *                     Numder
-      use ifwrit_m, only : Jjjdop, Kssmsc
+      use ifwrit_m, only : Kssmsc
       use exploc_common_m
       use array_sizes_common_m
       use oopsch_common_m, only : Nowwww, Segmen
@@ -224,7 +224,6 @@ C
         Ktheta = jkn2 ! remove
         call multScat%setNumTheta(jkn3 )
         Ntheta = jkn3 ! remove
-      Jjjdop = 0
       CALL Write_Commons_Many
       RETURN
       END
diff --git a/sammy/src/old/mold1.f b/sammy/src/old/mold1.f
index cfc7639be1f94bd50e3547f87b6c817d31e525c5..af6736351fbcca172f002e1a15cc20264327e1a6 100644
--- a/sammy/src/old/mold1.f
+++ b/sammy/src/old/mold1.f
@@ -70,6 +70,9 @@ C
      &              dSQRT(Boltzm*A_Mass_Small*Temp/
      &                resParData%getMassForIsotope(Iso))
                IF (Dopple.LT.Dopwid(Iso)) Dopple = Dopwid(Iso)
+               if (associated(dopplerOption%broadener)) then
+                  dopplerOption%wantBroaden = .true.
+               end if
             END DO
          END IF
          IF (Sitemp.NE.Zero) THEN
diff --git a/sammy/src/orr/morr0.f90 b/sammy/src/orr/morr0.f90
index eb2599b08a7231c843a9447370b6b04ab2ed1805..bdd2b742cd4d2ba1b07c767216f91d8a6bd94650 100644
--- a/sammy/src/orr/morr0.f90
+++ b/sammy/src/orr/morr0.f90
@@ -8,7 +8,6 @@ module orr_m
 ! *** Purpose -- ORR resolution function calculation
 !
       use fixedi_m, only : Numorr
-      use ifwrit_m, only : Jjjdop
       use fixedr_m, only : Dist
       use EndfData_common_m, only : expData
       use DopplerAndResolutionBroadener_M
@@ -75,7 +74,6 @@ module orr_m
       call work%destroy()
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 !
-      Jjjdop = 0
       CALL Write_Commons_Many
       RETURN
 !
diff --git a/sammy/src/orr/morr1.f90 b/sammy/src/orr/morr1.f90
index 9258a87adef1665608e8d6ff4ecc8b696db062ff..5b23b941104d97646a625373f362e242f0a10a94 100644
--- a/sammy/src/orr/morr1.f90
+++ b/sammy/src/orr/morr1.f90
@@ -11,7 +11,7 @@ module orr1_m
 ! *** Purpose -- Form Resolution-broadened cross section and derivatives 
 !
       use fixedi_m, only : Numbgf, Numnbk, Numorr, numUsedPar
-      use ifwrit_m, only : Jjjdop, Kdebug
+      use ifwrit_m, only : Kdebug
       use brdd_common_m, only : Ipk, Ipnts, Iup, Kc
       use DopplerAndResolutionBroadener_M
       use mxct27_m
@@ -72,7 +72,7 @@ module orr1_m
          IF (Ipnts.LE.5) THEN
 !
 ! ********* no integration possible
-            IF (Jjjdop.NE.1) THEN
+            IF (dopplerOption%bType.ne.1) THEN
                Now = Now + 1
                do Ipar = 0, numUsedPar
                  do ns = 1, derivs%getNnnsig()
diff --git a/sammy/src/par/mpar04.f90 b/sammy/src/par/mpar04.f90
index 234c2ba4291280e226eb17444c11ad272a728e67..43372a92328451733e95ce9de67df923873196a7 100644
--- a/sammy/src/par/mpar04.f90
+++ b/sammy/src/par/mpar04.f90
@@ -359,6 +359,9 @@ module par4_m
          Sitemp = A
          Sithck = C
          Dopple = dSQRT(Boltzm*  Temp*Aneutr/Aaawww)
+         if (associated(dopplerOption%broadener)) then
+            dopplerOption%wantBroaden = .true.
+         end if
          Dosind = dSQRT(Boltzm*Sitemp*Aneutr/Aaawww)
          N = N + 2
 !
diff --git a/sammy/src/ref/mrfs3.f b/sammy/src/ref/mrfs3.f
index ba71150c7380d5526f82575c41e62b8117b832eb..f7d10bb00060d71e6851933f812831801368073c 100644
--- a/sammy/src/ref/mrfs3.f
+++ b/sammy/src/ref/mrfs3.f
@@ -79,7 +79,6 @@ C *** CARD SET 5   Broadening information
       Deltal = 0.0d0
       Deltae = 0.0d0
       Delta6 = 0.0d0
-      Delttt = 0.0d0
 C
 C *** CARD SET 6 (IF IT'S PRESENT)
       Deltab = 1.0d0
diff --git a/sammy/src/ref/mwrt1.f b/sammy/src/ref/mwrt1.f
index d09dbb4843ab060271dc60587a2b54d893cffa44..0f8cbb226668e94d80d5bf6eb8c98ee28a9e39e8 100644
--- a/sammy/src/ref/mwrt1.f
+++ b/sammy/src/ref/mwrt1.f
@@ -41,7 +41,7 @@ C
       type(RMatParticlePair)::pair
       type(RMatChannelParams)::channel
       type(SammyChannelInfo)::channelInfo
-      real(kind=8)::Elowbr
+      real(kind=8)::Elowbr, Delttt
 C
       OPEN (UNIT=11, FILE='SAMMY.INP', STATUS='unknown',
      *     FORM='formatted')
@@ -77,6 +77,7 @@ C *** Card Set 4.5   Energies at which to plot resolution fUnction
 10450 FORMAT (8F10.3)
 C
 C *** CARD SET 5   Broadening information
+      Delttt = 0.0d0
       IF (Ib.EQ.1) then
          Elowbr = 0
          if (associated(dopplerOption%freeGas)) then
@@ -84,6 +85,9 @@ C *** CARD SET 5   Broadening information
          else if (
      *    associated(dopplerOption%highEnergyFreeGas)) then
             Elowbr = dopplerOption%highEnergyFreeGas%Elowbr
+         else if (
+     *      associated(dopplerOption%lealHwang)) then
+            Delttt =  dopplerOption%lealHwang%DelT
          end if
          WRITE (11,10500) Temp, Dist, Deltal, Deltae, Deltag,
      *   Delttt, Elowbr
diff --git a/sammy/src/rpi/mrpi0.f90 b/sammy/src/rpi/mrpi0.f90
index 8bc3bda1c2359705ac6800d2e3f7ded9a6fdc082..97692cc8e3695ce7fec2b38ace57bc5eceea13a9 100644
--- a/sammy/src/rpi/mrpi0.f90
+++ b/sammy/src/rpi/mrpi0.f90
@@ -7,7 +7,6 @@ module rpi_m
       SUBROUTINE Samrpi_0
 !
       use fixedi_m, only : Medrpi, Mmmrpi
-      use ifwrit_m, only : Jjjdop
       use exploc_common_m
       use array_sizes_common_m
       use oopsch_common_m, only : Nowwww, Segmen
@@ -79,7 +78,6 @@ module rpi_m
       deallocate(A_Ixxxrp)
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 !
-      Jjjdop = 0
       CALL Write_Commons_Many
       RETURN
 !
diff --git a/sammy/src/rpi/mrpi1.f90 b/sammy/src/rpi/mrpi1.f90
index 02a2602ac15b4158365110da83a9ba1825621a58..2b0d5c4458cedea1838939fdf160c79798dd988b 100644
--- a/sammy/src/rpi/mrpi1.f90
+++ b/sammy/src/rpi/mrpi1.f90
@@ -12,7 +12,7 @@ module rpi1_m
 !
       use fixedi_m, only : numUsedPar,        &
                            Nfprpi, Nnnrpi, Numnbk, Numbgf
-      use ifwrit_m, only : Jjjdop, Kdebug, Ksolve
+      use ifwrit_m, only : Kdebug, Ksolve
       use brdd_common_m, only : Ipk, Ipnts, Iup, Kc
       use rpijnk_common_m
       use rpires_common_m
@@ -84,7 +84,7 @@ module rpi1_m
          IF (Ipnts.LE.5) THEN
 !
 ! ********* No integration possible
-            IF (Jjjdop.NE.1) THEN
+            IF (dopplerOption%bType.ne.1) THEN
                Now = Now + 1
                do Ipar = 0, numUsedPar
                  do ns = 1, derivs%getNnnsig()
diff --git a/sammy/src/rsl/mrsl0.f90 b/sammy/src/rsl/mrsl0.f90
index 8c6ba514a2a9d63a5765e64601aceb5232d7b84b..eaca6b2f2944c2475797de77db07e13bb2269c76 100644
--- a/sammy/src/rsl/mrsl0.f90
+++ b/sammy/src/rsl/mrsl0.f90
@@ -7,7 +7,6 @@ module rsl_m
 !
       use fixedi_m, only : K2reso, Kkkdex, Nudwhi,  &
                            Numorr, Numrpi
-      use ifwrit_m, only : Jjjdop
       use exploc_common_m
       use array_sizes_common_m
       use oopsch_common_m, only : Nowwww, Segmen
@@ -69,7 +68,6 @@ module rsl_m
       deallocate(A_Iwts)
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 !
-      Jjjdop = 0
       CALL Write_Commons_Many
 
       call eGrid%destroy()
diff --git a/sammy/src/rsl/mrsl1.f90 b/sammy/src/rsl/mrsl1.f90
index 3e86a34123320bdd31d43e936679e5e9a067d1f1..3155888fff1ac4511efa42e0b2cf5cc94ac22448 100644
--- a/sammy/src/rsl/mrsl1.f90
+++ b/sammy/src/rsl/mrsl1.f90
@@ -15,9 +15,9 @@ module rsl1_m
 ! ***            derivatives
 !
       use fixedi_m, only : K2reso,Numbgf, Numnbk, numUsedPar
-      use ifwrit_m, only : Jjjdop, Kdebug, Kjdele, Kjdell, Ksolve, Ndat
+      use ifwrit_m, only : Kdebug, Kjdele, Kjdell, Ksolve, Ndat
       use fixedr_m
-      use broad_common_m, only : Bo2, Co2, Dell00, Deltae,Dell00, Dele00, Dele11, Dele22, Dell11, deltal
+      use broad_common_m, only : Bo2, Co2, Dell00, Deltae,Dell00, Dele00, Dele11, Dele22, Dell11, deltal, dopplerOption
       use brdd_common_m, only : Ipk, Ipnts, Iup, Kc
       use constn_common_m, only : Sm2
       use Wdsint_m, only : Wdsint
@@ -112,7 +112,7 @@ module rsl1_m
 !
          IF (Ipnts.LE.5) THEN
 ! ********* No integration is possible
-            IF (Jjjdop.NE.1) THEN
+            IF (dopplerOption%bType.ne.1) THEN
                Now = Now + 1
                do Ipar = 0, numUsedPar
                  do ns = 1, derivs%getNnnsig()
diff --git a/sammy/src/rst/mrst0.f90 b/sammy/src/rst/mrst0.f90
index 5181d7cece11aa642a4f5edd8d3bf92d5bf4edb0..4c367a4a8d91bf43062bedba1a057e1f5346372d 100644
--- a/sammy/src/rst/mrst0.f90
+++ b/sammy/src/rst/mrst0.f90
@@ -4,7 +4,7 @@ module rst_m
 !
       SUBROUTINE Samrst_0
 !
-      use ifwrit_m, only : Jjjdop, Ndatb
+      use ifwrit_m, only : Ndatb
       use exploc_common_m
       use oopsch_common_m, only : Nowwww, Segmen
       use rst1_m
@@ -38,7 +38,6 @@ module rst_m
       deallocate(A_IWts)
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 !
-      Jjjdop = 0
       CALL Write_Commons_Few
       RETURN
 !
diff --git a/sammy/src/ssm/mssm00.f90 b/sammy/src/ssm/mssm00.f90
index d29652af0e18e9862f0418ff456982ea5f077b8c..ea3aff41dc7eabbc55e1f765043996e4be9cf253 100644
--- a/sammy/src/ssm/mssm00.f90
+++ b/sammy/src/ssm/mssm00.f90
@@ -17,7 +17,6 @@ module ssm_m
       use fixedi_m, only : Jtheta, K2reso, Kkkdex, Kkkrsl, &
                            Ktheta, Ntheta, Nudwhi,  &
                            Numder, Numorr, Numrpi
-      use ifwrit_m, only : Jjjdop
       use exploc_common_m, only  : I_Ixciso
       use oopsch_common_m, only : Nowwww, Segmen
       use logic_ssm_common_m
@@ -124,7 +123,6 @@ module ssm_m
       Ktheta = Kthhhh ! remove
       call multScat%setNumThetaNearOne( Jthhhh )
       Jtheta = Jthhhh ! remove
-      Jjjdop = 0
       CALL Write_Commons_Many
       RETURN
       end subroutine Samssm_0
diff --git a/sammy/src/the/mthe1.f90 b/sammy/src/the/mthe1.f90
index ddea116a3f9f271db8d2096ed75b0aa5b67d3c3a..8b0ee838a519295801ca0ea1db6aa649fbb5e0a0 100644
--- a/sammy/src/the/mthe1.f90
+++ b/sammy/src/the/mthe1.f90
@@ -475,12 +475,12 @@ module mthe1_m
 !
       SUBROUTINE Setthe(Nnnsig, Nnniso)
 !
-      use fixedi_m, only : Jcros, Kjjjjj,  Lllmax,   &
+      use fixedi_m, only : Jcros,  Lllmax,   &
                            Numcro, Numder, Numiso,  &
                            Jtrans
       use ifwrit_m, only : Kaverg, Kdebug, Kfinit, Kkkclq, Kcros,     &
                            Krmatx, Ksolve, Kssmsc, Ktrans, MIKE, Nfissl
-      use broad_common_m, only : Dopple      
+      use broad_common_m, only : dopplerOption
       IMPLICIT None
       integer::Nnnsig, Nnniso
       real(kind=8)::zero
@@ -527,7 +527,7 @@ module mthe1_m
       END IF
 !
       IF (Nnnsig.EQ.1 .AND. (Kssmsc.EQ.-1 .OR. Kssmsc.EQ.0) .AND.      &
-          Kaverg.NE.2 .AND. (Dopple.EQ.Zero .AND. Kjjjjj.EQ.0)) THEN
+          Kaverg.NE.2 .AND. .not.dopplerOption%wantBroaden) THEN
 ! ***    If we can add up isotopes in XCT rather than waiting for ANG or
 ! ***        FGM or SSM...
 ! ***    ie for anything other than self-shielding, multiple-scattering,
@@ -549,7 +549,7 @@ module mthe1_m
       IF (Kssmsc.EQ.-2) THEN
 ! ***    for self-shielding but not multiple-scattering (if no Doppler)
          Nnnsig = 2
-         IF (Dopple.EQ.Zero .AND. Kjjjjj.EQ.0) Nnniso = 1         
+         IF (.not.dopplerOption%wantBroaden) Nnniso = 1
       END IF
       IF (Kaverg.EQ.2) THEN
 ! ***    for Bondarenko group averages
@@ -587,12 +587,12 @@ module mthe1_m
 !
       SUBROUTINE Orgniz
 !
-      use fixedi_m, only : Jcros, Jtrans, Kjjjjj, Kkkdex, Kkkrsl,   &
+      use fixedi_m, only : Jcros, Jtrans,  Kkkdex, Kkkrsl,   &
                            Ktruet, Nudwhi, Numbgf, Numnbk, Numorr,  &
                            Numrpi
       use ifwrit_m, only : Kaverg, Kdebug, Kssmsc, Ntgrlq
       use fixedr_m, only : Dddeee
-      use broad_common_m, only : Ao2, Bo2, Co2, Dopple, Ncf
+      use broad_common_m, only : Ao2, Bo2, Co2, Ncf, dopplerOption
       use lbro_common_m, only : Debug, Yaverg, Ydoppr, Ynrmbk,  &
                                 Yresol, Yselfi, Yssmsc, Ytotrs, Ytrans
 
@@ -616,9 +616,8 @@ module mthe1_m
       Yselfi = False
       Yaverg = False
 !
-! *** IF (WANT DOPPLER BROADENING) ydoppr = True
-      IF (Dopple.NE.Zero) Ydoppr = True
-      IF (Kjjjjj.NE.0   ) Ydoppr = True
+! *** IF (WANT DOPPLER BROADENING) ydoppr = True    
+      IF (dopplerOption%wantBroaden) Ydoppr = True
 !
 ! *** if (this is transmission) ytrans = True
       IF (Jcros.EQ.1 .AND. Jtrans.EQ.1) Ytrans = True
diff --git a/sammy/src/udr/mudr0.f b/sammy/src/udr/mudr0.f
index 8bfe3776d86c374d770f65f4eb18bfa15b2252d8..aa0e71d230cb60488ead1ba6ffb5ab34cd02df9a 100644
--- a/sammy/src/udr/mudr0.f
+++ b/sammy/src/udr/mudr0.f
@@ -6,7 +6,6 @@ C
 C
       use oops_common_m, only : Msize
       use fixedi_m, only : Nudmax, Nudtim, Nudwhi, Numudr
-      use ifwrit_m, only : Jjjdop
       use exploc_common_m
       use array_sizes_common_m
       use oopsch_common_m, only : Nowwww, Segmen
@@ -80,7 +79,6 @@ C
 
 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 C
-      Jjjdop = 0
       CALL Write_Commons_Many
       RETURN
       END
diff --git a/sammy/src/udr/mudr1.f b/sammy/src/udr/mudr1.f
index e645e49a5affc886026aec4a33dc996597563d8e..ef4dcd2d3b1fdc4b7dd7d501fa2d5108a0705ee3 100644
--- a/sammy/src/udr/mudr1.f
+++ b/sammy/src/udr/mudr1.f
@@ -13,7 +13,7 @@ C
      *                     Numbgf, Numbgf,
      *                     Numnbk, Nudwhi, Ntepnt,
      *                     numUsedPar
-      use ifwrit_m, only : Jjjdop, Kdebug
+      use ifwrit_m, only : Kdebug
       use brdd_common_m, only : Ipk, Ipnts, Kc, Iup
       use mxct27_m
       use rsl3_m
@@ -80,7 +80,7 @@ C
 cx	call timer_debug('d',0)
 C
 C ********* no integration possible
-            IF (Jjjdop.NE.1) THEN
+            IF (dopplerOption%bType.ne.1) THEN
                Now = Now + 1
                do Ipar = 0, numUsedPar
                  do ns = 1, calcData%getNnnsig()