1. Matthew Turk
  2. grackle

Commits

Britton Smith  committed 84d69e3

Removed radiative transfer support.

  • Participants
  • Parent commits 20816bb
  • Branches default
  • Tags no-more-rad-transfer

Comments (0)

Files changed (7)

File src/clib/calculate_cooling_time.C

View file
 	gr_float *gpldl, gr_float *gphdl, gr_float *HDltea, gr_float *HDlowa,
 	gr_float *gaHIa, gr_float *gaH2a, gr_float *gaHea, gr_float *gaHpa, gr_float *gaela,
 	gr_float *gasgra, gr_int *iradfield,
-	gr_int *iradtrans, gr_float *photogamma,
 	gr_int *ih2optical, gr_int *iciecool, gr_float *ciecoa,
  	gr_int *icmbTfloor, gr_int *iClHeat,
  	gr_float *clEleFra, gr_int *clGridRank, gr_int *clGridDim,
        my_chemistry.GAHI, my_chemistry.GAH2, my_chemistry.GAHe, my_chemistry.GAHp,
        my_chemistry.GAel, my_chemistry.gas_grain,
        &my_chemistry.UVbackground_type,
-       &my_chemistry.RadiativeTransfer, my_chemistry.gammaNum, 
        &my_chemistry.h2_optical_depth_approximation, 
        &my_chemistry.cie_cooling, my_chemistry.cieco,
        &my_chemistry.cmb_temperature_floor,

File src/clib/chemistry_data.h

View file
   // Length of 1D flattened Cloudy data
   gr_int CloudyDataSize;
 
-  /******************************************************
-   *** place holder for radiative transfer parameters ***
-   *** These do nothing and should be left alone. *******
-   ******************************************************/
-
-  gr_int RadiativeTransfer;
-  gr_int RadiativeTransferCoupledRateSolver;
-  gr_int RTCoupledSolverIntermediateStep;
-  gr_int RadiativeTransferHydrogenOnly;
-  gr_float *kphHINum, *kphHeINum, *kphHeIINum, 
-    *kdissH2INum, *gammaNum;
-
 };

File src/clib/cool1d_multi.F

View file
      &                indixe, t1, t2, logtem, tdef, edot,
      &                tgas, tgasold, p2d, tdust, metallicity, rhoH, 
      &                iradtype,
-     &                iradtrans, photogamma,
      &                ih2optical, iciecool, ciecoa, cieco,
      &                icmbTfloor, iClHeat,
      &                clEleFra, clGridRank, clGridDim,
 
       integer in, jn, kn, is, ie, j, k, nratec, idim,
      &        iexpand, ih2co, ipiht, ispecies, imcool, idust, 
-     &        iradtype, iradtrans,
+     &        iradtype,
      &        imetal, igammah, ih2optical, iciecool
       real    aye, temstart, temend, z_solar,
      &        utem, uxyz, uaye, urho, utim,
       real    HM(in,jn,kn),  H2I(in,jn,kn), H2II(in,jn,kn),
      &        DI(in,jn,kn),  DII(in,jn,kn), HDI(in,jn,kn),
      &        metal(in,jn,kn)
-      real    photogamma(in,jn,kn)
       real    hyd01ka(nratec), h2k01a(nratec), vibha(nratec), 
      &        rotha(nratec), rotla(nratec), gpldla(nratec),
      &        gphdla(nratec), hdltea(nratec), hdlowa(nratec)
 
 !  Locals
 
-      integer i, j1, iter, iradfield
+      integer i, j1, iter
       real dom, qq, vibl, logtem0, logtem9, dlogtem, energy, zr
       real dt2, ttmin, comp1, comp2, scoef, acoef, HIdot,
      &     hdlte1, hdlow1, gamma2, x, fudge, fH2,
       coolunit = (uaye**5 * xbase1**2 * mh**2) / (tbase1**3 * dbase1)
       zr       = 1.d0/(aye*uaye) - 1.d0
       fudge    = 1.d0
-      iradfield = -1
 
 !     Set compton cooling coefficients (and temperature)
 
          end if
       enddo
 
-!     Photoheating from radiative transfer
-
-      if (iradtrans .eq. 1) then
-c     for radiative transfer, convert from eV/s*TimeUnits to the coolunits use
-         rtunits = ev2erg/utim/coolunit/dom
-         do i = is+1, ie+1
-            if (itmask(i)) then
-               edot(i) = edot(i) + float(ipiht) * photogamma(i,j,k) * 
-     &              rtunits * HI(i,j,k)
-c               if (photogamma(i,j,k).gt.0) then
-c                  print*, i,j,k,edot(i), photogamma(i,j,k),rtunits,dom,
-c     $                 aye,utim,coolunit,d(i,j,k),hi(i,j,k),tgas(i)
-c               endif
-               if (edot(i) .ne. edot(i)) then
-                  write(6,*) 'NaN in edot[2]: ', i, j, k, edot(i), 
-     &                 photogamma(i,j,k),HI(i,j,k), de(i,j,k), d(i,j,k), 
-     &                 e(i,j,k), p2d(i), tgas(i), dom, urho, aye, mh
-               endif
-            endif
-         enddo
-      endif
-
-
 !     --- Cloudy metal cooling and heating ---
 
       if (imcool .eq. 1) then

File src/clib/cool_multi_time.F

View file
      &                gaHIa, gaH2a, gaHea, gaHpa, gaela,
      &                gasgra, 
      &                iradtype,
-     &                iradtrans, photogamma,
      &                ih2optical, iciecool, ciecoa, 
      &                icmbTfloor, iClHeat,
      &                clEleFra, clGridRank, clGridDim,
 
       integer in, jn, kn, is, js, ks, ie, je, ke, nratec, 
      &        iexpand, ih2co, ipiht, ispecies, imetal, idim,
-     &        iradtype, iradtrans,
+     &        iradtype, 
      &        imcool, idust, igammah, ih2optical, iciecool
       real    dt, aye, temstart, temend,
      &        utem, uxyz, uaye, urho, utim,
       real    HM(in,jn,kn),  H2I(in,jn,kn), H2II(in,jn,kn),
      &        DI(in,jn,kn),  DII(in,jn,kn), HDI(in,jn,kn)
       real    metal(in,jn,kn)
-      real    photogamma(in,jn,kn)
       real    hyd01ka(nratec), h2k01a(nratec), vibha(nratec), 
      &        rotha(nratec), rotla(nratec), gpldla(nratec),
      &        gphdla(nratec), hdltea(nratec), hdlowa(nratec)
      &                indixe, t1, t2, logtem, tdef, edot,
      &                tgas, tgasold, p2d, tdust, metallicity, rhoH,
      &                iradtype,
-     &                iradtrans, photogamma, 
      &                ih2optical, iciecool, ciecoa, cieco,
      &                icmbTfloor, iClHeat,
      &                clEleFra, clGridRank, clGridDim,

File src/clib/set_default_chemistry_parameters.C

View file
   my_chemistry.CloudyCoolingGridRank          = 0;
   my_chemistry.CloudyElectronFractionFactor = 9.153959e-3; // Cloudy 07.02 abundances
 
-  /* Some inactive place holders. */
-  my_chemistry.RadiativeTransfer = 0;
-  my_chemistry.RadiativeTransferCoupledRateSolver = 0;
-  my_chemistry.RTCoupledSolverIntermediateStep = 0;
-  my_chemistry.RadiativeTransferHydrogenOnly = 0;
-  my_chemistry.kphHINum = my_chemistry.kphHeINum = 
-    my_chemistry.kphHeIINum = my_chemistry.kdissH2INum = 
-    my_chemistry.gammaNum = NULL;
-
   return my_chemistry;
 }

File src/clib/solve_chemistry.C

View file
 	gr_float *gpldl, gr_float *gphdl, gr_float *HDltea, gr_float *HDlowa,
 	gr_float *gaHIa, gr_float *gaH2a, gr_float *gaHea, gr_float *gaHpa, gr_float *gaela,
 	gr_float *gasgra, gr_int *iradtype,
-	gr_int *iradtrans, gr_int *iradcoupled, gr_int *iradstep, gr_int *ierr,
-	gr_int *irt_honly, gr_float *kphHI, gr_float *kphHeI, gr_float *kphHeII, 
-	gr_float *kdissH2I, gr_float *photogamma,
+	gr_int *ierr,
 	gr_int *ih2optical, gr_int *iciecool, gr_int *ithreebody, gr_float *ciecoa,
  	gr_int *icmbTfloor, gr_int *iClHeat,
  	gr_float *clEleFra, gr_int *clGridRank, gr_int *clGridDim,
     my_chemistry.HDlte, my_chemistry.HDlow,
     my_chemistry.GAHI, my_chemistry.GAH2, my_chemistry.GAHe, my_chemistry.GAHp,
     my_chemistry.GAel, my_chemistry.gas_grain, &my_chemistry.UVbackground_type, 
-    &my_chemistry.RadiativeTransfer, &my_chemistry.RadiativeTransferCoupledRateSolver,
-    &my_chemistry.RTCoupledSolverIntermediateStep, &ierr,
-    &my_chemistry.RadiativeTransferHydrogenOnly,
-    my_chemistry.kphHINum, my_chemistry.kphHeINum, my_chemistry.kphHeIINum, 
-    my_chemistry.kdissH2INum, my_chemistry.gammaNum,
+    &ierr,
     &my_chemistry.h2_optical_depth_approximation, &my_chemistry.cie_cooling, 
     &my_chemistry.three_body_rate, my_chemistry.cieco,
     &my_chemistry.cmb_temperature_floor,

File src/clib/solve_rate_cool.F

View file
      &                gaHIa, gaH2a, gaHea, gaHpa, gaela, 
      &                gasgra, 
      &                iradtype,
-     &                iradtrans, iradcoupled, iradstep, ierr,
-     &                irt_honly,
-     &                kphHI, kphHeI, kphHeII, kdissH2I, 
-     &                photogamma, 
+     &                ierr,
      &                ih2optical, iciecool, ithreebody, ciecoa, 
      &                icmbTfloor, iClHeat,
      &                clEleFra, clGridRank, clGridDim,
 !    HM       - H- density field
 !    H2I      - H_2 (molecular H) density field
 !    H2II     - H_2+ density field
-!    kph*     - photoionization fields
-!    gamma*   - photoheating fields
 !
 !    is,ie    - start and end indices of active region (zero based)
 !    iexpand  - comoving coordinates flag (0 = off, 1 = on)
 !    idust    - flag for H2 formation on dust grains
 !    ih2co    - flag to include H2 cooling (1 = on, 0 = off)
 !    ipiht    - flag to include photoionization heating (1 = on, 0 = off)
-!    iradtrans - flag to include radiative transfer (1 = on, 0 = off)
-!    iradcoupled - flag to indicate coupled radiative transfer
-!    iradstep - flag to indicate intermediate coupled radiative transfer timestep
 !
 !    fh       - Hydrogen mass fraction (typically 0.76)
 !    dtoh     - Deuterium to H mass ratio
 
       integer in, jn, kn, is, js, ks, ie, je, ke, nratec, 
      &        iexpand, ih2co, ipiht, ispecies, imetal, idim,
-     &        iradtype, iradtrans,
-     &        iradcoupled, iradstep, ierr, imcool, idust, 
-     &        irt_honly, igammah, ih2optical, iciecool, ithreebody,
+     &        iradtype, 
+     &        ierr, imcool, idust, 
+     &        igammah, ih2optical, iciecool, ithreebody,
      &        ndratec
       real    dt, aye, temstart, temend, gamma,
      &        utim, uxyz, uaye, urho, utem, fh, dtoh, z_solar, 
      &        u(in,jn,kn),    v(in,jn,kn),     w(in,jn,kn),
      &        metal(in,jn,kn)
 
-!  Radiation fields
-
-      real kphHI(in,jn,kn), kphHeI(in,jn,kn), kphHeII(in,jn,kn),
-     &     kdissH2I(in,jn,kn)
-      real photogamma(in,jn,kn)
-
 !  Cooling tables (coolings rates as a function of temperature)
 
       real    hyd01ka(nratec), h2k01a(nratec), vibha(nratec), 
 
 !       tolerance = 1.0e-06 * dt
 
-!       Set iteration mask to include only cells with radiation in the
-!       intermediate coupled chemistry / energy step
+!     Initialize iteration mask to true for all cells.
 
-         if (iradcoupled .eq. 1 .and. iradstep .eq. 1) then
-            do i = is+1, ie+1
-               if (kphHI(i,j,k) .gt. 0) then 
-                  itmask(i) = .true. 
-               else 
-                  itmask(i) = .false.
-               endif
-            enddo
-         endif
-
-!     Normal rate solver, but don't double-count the cells with radiation
- 
-         if (iradcoupled .eq. 1 .and. iradstep .eq. 0) then
-            do i = is+1, ie+1
-               if (kphHI(i,j,k) .gt. 0) then 
-                  itmask(i) = .false. 
-               else 
-                  itmask(i) = .true.
-               endif
-            enddo
-         endif
-
-!     No radiation timestep coupling
-
-         if (iradcoupled .eq. 0 .or. iradtrans .eq. 0) then 
-            do i = is+1, ie+1
-               itmask(i) = .true.
-            enddo
-         endif
+         do i = is+1, ie+1
+            itmask(i) = .true.
+         enddo
 
 !        Set time elapsed to zero for each cell in 1D section
 
      &                indixe, t1, t2, logtem, tdef, edot,
      &                tgas, tgasold, p2d, tdust, metallicity, rhoH, 
      &                iradtype, 
-     &                iradtrans, photogamma,
      &                ih2optical, iciecool, ciecoa, cieco,
      &                icmbTfloor, iClHeat,
      &                clEleFra, clGridRank, clGridDim,
      &                     k50, k51, k52, k53, k54, k55, k56, 
      &                     h2dust, ncrn, ncrd1, ncrd2, rhoH, 
      &                     k24shield, k25shield, k26shield,
-     &                     iradtrans, irt_honly, kphHI, kphHeI, kphHeII, 
-     &                     kdissH2I, itmask, edot, chunit, dom)
+     &                     itmask, edot, chunit, dom)
 
 !           Find timestep that keeps relative chemical changes below 10%
 
      &                     HIp, HIIp, HeIp, HeIIp, HeIIIp, dep,
      &                     HMp, H2Ip, H2IIp, DIp, DIIp, HDIp,
      &                     dedot_prev, HIdot_prev,
-     &                     iradtrans, irt_honly, kphHI, kphHeI, kphHeII, 
-     &                     kdissH2I, itmask)
+     &                     itmask)
 
 !           Add the timestep to the elapsed time for each cell and find
 !            minimum elapsed time step in this row
      &                     k50, k51, k52, k53, k54, k55, k56, 
      &                     h2dust, ncrn, ncrd1, ncrd2, rhoH, 
      &                     k24shield, k25shield, k26shield,
-     &                     iradtrans, irt_honly, kphHI, kphHeI, kphHeII, 
-     &                     kdissH2I, itmask, edot, chunit, dom)
+     &                     itmask, edot, chunit, dom)
 
 ! -------------------------------------------------------------------
 
 
 !     arguments
 
-      integer ispecies, idust, is, ie, j, k, in, jn, kn, 
-     &     iradtrans, irt_honly
+      integer ispecies, idust, is, ie, j, k, in, jn, kn
       real dedot(in), HIdot(in), dom
       double precision edot(in)
       logical itmask(in)
      &         d(in,jn,kn)
       real    HM(in,jn,kn),  H2I(in,jn,kn), H2II(in,jn,kn)
 
-!     Radiation fields
-
-      real    kphHI(in,jn,kn), kphHeI(in,jn,kn), kphHeII(in,jn,kn),
-     &        kdissH2I(in,jn,kn)
-
       double precision chunit
 
 !     Rate values
          enddo
       endif
 
-!     Add photo-ionization rates if needed
-
-      if (iradtrans .eq. 1) then
-         if (irt_honly .eq. 0) then
-            do i = is+1, ie+1
-               if (itmask(i)) then
-                  HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k)
-                  dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k)
-     &                 + kphHeI(i,j,k) * HeI(i,j,k) / 4.d0
-     &                 + kphHeII(i,j,k) *HeII(i,j,k) / 4.d0
-               endif
-            enddo
-         else
-            do i = is+1, ie+1
-               if (itmask(i)) then
-                  HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k)
-                  dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k)
-               endif
-            enddo
-         endif
-      endif
-
       return
       end
 
      &                     HIp, HIIp, HeIp, HeIIp, HeIIIp, dep,
      &                     HMp, H2Ip, H2IIp, DIp, DIIp, HDIp,
      &                     dedot_prev, HIdot_prev,
-     &                     iradtrans, irt_honly, kphHI, kphHeI, kphHeII,
-     &                     kdissH2I, itmask)
+     &                     itmask)
 c -------------------------------------------------------------------
 
       implicit NONE
 
 !     arguments
 
-      integer ispecies, idust, in, jn, kn, is, ie, j, k, 
-     &     iradtrans, irt_honly
+      integer ispecies, idust, in, jn, kn, is, ie, j, k
       real    dtit(in), dedot_prev(in), HIdot_prev(in)
       logical itmask(in)
 
       real    HM(in,jn,kn),  H2I(in,jn,kn), H2II(in,jn,kn)
       real    DI(in,jn,kn),  DII(in,jn,kn), HDI(in,jn,kn)
 
-!     Radiation fields
-
-      real    kphHI(in,jn,kn), kphHeI(in,jn,kn), kphHeII(in,jn,kn),
-     &        kdissH2I(in,jn,kn)
-
 !     Rate values
 
       real k1 (in), k2 (in), k3 (in), k4 (in), k5 (in),
             scoef  = k2(i)*HII(i,j,k)*de(i,j,k)
             acoef  = k1(i)*de(i,j,k)
      &             + k24shield(i)
-            if (iradtrans .eq. 1) acoef = acoef + kphHI(i,j,k)
             HIp(i)  = (scoef*dtit(i) + HI(i,j,k))/(1.d0 + acoef*dtit(i))
             if (HIp(i) .ne. HIp(i)) then
                write(*,*) 'HUGE HIp! :: ', i, j, k, HIp(i), HI(i,j,k),
-     $              HII(i,j,k), de(i,j,k), kphHI(i,j,k), 
+     $              HII(i,j,k), de(i,j,k), 
      $              scoef, acoef, dtit(i)
 c               ERROR_MESSAGE
             endif
 c 
             scoef  = k1(i)*HIp(i)*de(i,j,k)
      &                   + k24shield(i)*HIp(i)
-            if (iradtrans .eq. 1) scoef = scoef + kphHI(i,j,k)*HIp(i)
             acoef  = k2(i)*de (i,j,k)
             HIIp(i) = (scoef*dtit(i) + HII(i,j,k))/(1.d0 +acoef*dtit(i))
 !
                write(*,*) 'negative HIIp! :: ', i, j, k, HIIp(i), 
      $              scoef, dtit(i), HII(i,j,k), acoef,
      $              k2(i), de(i,j,k),
-     $              kphHI(i,j,k), HIp(i),
+     $              HIp(i),
      $              k24shield(i)
             endif
 
      &                 + k24shield(i)*HI(i,j,k)
      &                 + k25shield(i)*HeII(i,j,k)/4.d0
      &                 + k26shield(i)*HeI(i,j,k)/4.d0
-            if (iradtrans .eq. 1 .and. irt_honly .eq. 0) 
-     &           scoef = scoef + kphHI(i,j,k)   * HI(i,j,k)
-     &           + kphHeI(i,j,k)  * HeI(i,j,k)/4.d0
-     &           + kphHeII(i,j,k) * HeII(i,j,k)/4.d0
-            if (iradtrans .eq. 1 .and. irt_honly .eq. 1) 
-     &           scoef = scoef + kphHI(i,j,k)   * HI(i,j,k)
             acoef = -(k1(i)*HI(i,j,k)      - k2(i)*HII(i,j,k)
      &              + k3(i)*HeI(i,j,k)/4.d0  - k6(i)*HeIII(i,j,k)/4.d0
      &              + k5(i)*HeII(i,j,k)/4.d0 - k4(i)*HeII(i,j,k)/4.d0)
          scoef  = k4(i)*HeII(i,j,k)*de(i,j,k)
          acoef  = k3(i)*de(i,j,k)
      &                + k26shield(i)
-         if (iradtrans.eq.1 .and. irt_honly.eq.0) 
-     $        acoef = acoef + kphHeI(i,j,k)
          HeIp(i)   = ( scoef*dtit(i) + HeI(i,j,k) ) 
      &              / ( 1.d0 + acoef*dtit(i) )
 
          scoef  = k3(i)*HeIp(i)*de(i,j,k)
      &          + k6(i)*HeIII(i,j,k)*de(i,j,k)
      &          + k26shield(i)*HeIp(i)
-         if (iradtrans.eq.1 .and. irt_honly.eq.0) 
-     $        scoef = scoef + kphHeI(i,j,k)*HeIp(i)
          acoef  = k4(i)*de(i,j,k) + k5(i)*de(i,j,k)
      &          + k25shield(i)
-         if (iradtrans.eq.1 .and. irt_honly.eq.0) 
-     $        acoef = acoef + kphHeII(i,j,k)
          HeIIp(i)  = ( scoef*dtit(i) + HeII(i,j,k) )
      &              / ( 1.d0 + acoef*dtit(i) )
 
 
          scoef   = k5(i)*HeIIp(i)*de(i,j,k)
      &           + k25shield(i)*HeIIp(i)
-         if (iradtrans.eq.1 .and. irt_honly.eq.0) 
-     $        scoef = scoef + kphHeII(i,j,k) * HeIIp(i)
          acoef   = k6(i)*de(i,j,k)
          HeIIIp(i)  = ( scoef*dtit(i) + HeIII(i,j,k) )
      &                / ( 1.d0 + acoef*dtit(i) )
      &             + 2.d0*k18(i)* H2II(i,j,k)* de(i,j,k)/2.d0
      &             +      k19(i)* H2II(i,j,k)* HM(i,j,k)/2.d0
      &             + 2.d0*k31   * H2I(i,j,k)/2.d0
-            if (iradtrans.eq.1)
-     &           scoef = scoef + 2.d0*kdissH2I(i,j,k)* H2I(i,j,k)/2.d0
             acoef  =      k1(i) * de(i,j,k)
      &             +      k7(i) * de(i,j,k)  
      &             +      k8(i) * HM(i,j,k)
      &             +      k10(i)* H2II(i,j,k)/2.d0
      &             + 2.d0*k22(i)* HI(i,j,k)**2
      &             + k24shield(i)
-            if (iradtrans .eq. 1) acoef = acoef + kphHI(i,j,k)
 
             if (idust .gt. 0) then
                acoef = acoef + 2.d0 * h2dust(i) * rhoH(i)
      &                      ( 1. + acoef*dtit(i) )
             if (HIp(i) .ne. HIp(i)) then
                write(*,*) 'HUGE HIp! :: ', i, j, k, HIp(i), HI(i,j,k),
-     $              HII(i,j,k), de(i,j,k), H2I(i,j,k), kdissH2I(i,j,k),
-     $              kphHI(i,j,k)
+     $              HII(i,j,k), de(i,j,k), H2I(i,j,k)
             endif
 
 !          2) HII
             scoef  =    k1(i)  * HI(i,j,k) * de(i,j,k)
      &             +    k10(i) * H2II(i,j,k)*HI(i,j,k)/2.d0
      &             + k24shield(i)*HI(i,j,k)
-            if (iradtrans .eq. 1) scoef = scoef + kphHI(i,j,k)*HI(i,j,k)
             acoef  =    k2(i)  * de(i,j,k)
      &             +    k9(i)  * HI(i,j,k)
      &             +    k11(i) * H2I(i,j,k)/2.d0
      &             + k24shield(i)*HIp(i)
      &             + k25shield(i)*HeIIp(i)/4.d0
      &             + k26shield(i)*HeIp(i)/4.d0
-            if (iradtrans .eq. 1 .and. irt_honly.eq.0) 
-     &           scoef = scoef + kphHI(i,j,k)   * HIp(i)
-     &           + kphHeI(i,j,k)  * HeIp(i)/4.d0
-     &           + kphHeII(i,j,k) * HeIIp(i)/4.d0
-            if (iradtrans .eq. 1 .and. irt_honly.eq.1) 
-     &           scoef = scoef + kphHI(i,j,k)   * HIp(i)
             acoef = - (k1(i) *HI(i,j,k)    - k2(i)*HII(i,j,k)
      &              +  k3(i) *HeI(i,j,k)/4.d0 - k6(i)*HeIII(i,j,k)/4.d0
      &              +  k5(i) *HeII(i,j,k)/4.d0- k4(i)*HeII(i,j,k)/4.d0
             acoef = ( k13(i)*HI(i,j,k) + k11(i)*HII(i,j,k)
      &              + k12(i)*de(i,j,k) )
      &              + k29 + k31
-            if (iradtrans.eq.1) acoef = acoef + kdissH2I(i,j,k)
 
             if (idust .gt. 0) then
                scoef = scoef + 2.d0 * h2dust(i) * HI(i,j,k) * rhoH(i)
      &             +    k54(i) * H2I(i,j,k)/2.d0
      &             +    k56(i) * HM(i,j,k)
      &             + k24shield(i)
-            if (iradtrans .eq. 1) acoef = acoef + kphHI(i,j,k)
             DIp(i)    = ( scoef*dtit(i) + DI(i,j,k) ) / 
      &                  ( 1.d0 + acoef*dtit(i) )
 
      &            +  2.d0*k53(i) * HII(i,j,k)* HDI(i,j,k)/3.d0
      &            )
      &            + k24shield(i)*DI(i,j,k)
-            if (iradtrans .eq. 1) scoef = scoef + kphHI(i,j,k)*DI(i,j,k)
             acoef =    k2(i)  * de(i,j,k)
      &            +    k51(i) * HI(i,j,k)
      &            +    k52(i) * H2I(i,j,k)/2.d0