diff --git a/unstructured/M3Dmodules.f90 b/unstructured/M3Dmodules.f90
index 79cfe966b..e484e08a2 100644
--- a/unstructured/M3Dmodules.f90
+++ b/unstructured/M3Dmodules.f90
@@ -459,6 +459,33 @@ module basic
 
   type(spline1d) :: q_spline
 
+contains
+
+  !   Added 1/1/2016 to get consistency between 2D,3D,Cyl,Tor
+  !   BCL moved 12/6/2019
+  subroutine tpi_factors(tpifac,tpirzero)
+    use math
+    implicit none
+    real, intent(out) :: tpifac, tpirzero
+    if(nplanes.eq.1) then
+       if(itor.eq.1) then
+          tpifac = 1.
+          tpirzero = 1.
+       else
+          tpifac = 1./rzero
+          tpirzero = 1.
+       endif
+    else
+       if(itor.eq.1) then
+          tpifac = twopi
+          tpirzero = twopi
+       else
+          tpifac = twopi
+          tpirzero = twopi*rzero
+       endif
+    endif
+  end subroutine tpi_factors
+
 end module basic
 
 module arrays
diff --git a/unstructured/adapt.f90 b/unstructured/adapt.f90
index 348e7cd73..41827ad78 100644
--- a/unstructured/adapt.f90
+++ b/unstructured/adapt.f90
@@ -233,8 +233,7 @@ subroutine adapt_by_psi
     i_control%err_p_old = 0.
     n_control%err_i = 0.
     n_control%err_p_old = 0.
-    if(myrank.eq.0 .and. iprint.ge.2) print *, "  transport coefficients"
-    call define_transport_coefficients
+    call define_transport_coefficients(1)
     call derived_quantities(1)
     !ke_previous = ekin
   end subroutine adapt_by_psi
@@ -406,11 +405,11 @@ subroutine adapt_by_error
        n_control%err_i = 0.
        n_control%err_p_old = 0.
        call reset_scalars
-       if(myrank.eq.0 .and. iprint.ge.2) print *, "  transport coefficients"
-       call define_transport_coefficients
        if(eqsubtract.eq.1) then
-         call derived_quantities(0)
+          call define_transport_coefficients(0)
+          call derived_quantities(0)
        end if
+       call define_transport_coefficients(1)
        call derived_quantities(1)
        meshAdapted =1
     end if
diff --git a/unstructured/diagnostics.f90 b/unstructured/diagnostics.f90
index 12936696d..5464578e6 100644
--- a/unstructured/diagnostics.f90
+++ b/unstructured/diagnostics.f90
@@ -80,9 +80,6 @@ module diagnostics
   real, dimension(iflux_loops_max) :: flux_loop_val
   integer, dimension(iflux_loops_max) :: flux_loop_itri
 
-  vectype, dimension(MAX_PTS) :: rhop79
-  vectype, dimension(MAX_PTS) :: Lorentz_pel
-
 contains
 
   ! ======================================================================
@@ -164,7 +161,6 @@ end subroutine distribute_timings
   ! resets diagnostic energy and scalar quantities to zero
   ! ======================================================================
   subroutine reset_scalars()
-    use pellet
 
     implicit none
 
@@ -235,10 +231,6 @@ subroutine reset_scalars()
     nfluxd = 0.
     nfluxv = 0.
     nsource = 0.
-    if(ipellet_abl.gt.0) then
-       nsource_pel = 0. ! this is an array
-       temp_pel = 0.    ! this is an array
-    end if
 
     bwb2 = 0.
 
@@ -272,7 +264,6 @@ end subroutine reset_scalars
   ! ======================================================================
   subroutine distribute_scalars()    
     use basic
-    use pellet
 
     implicit none
 
@@ -281,7 +272,6 @@ subroutine distribute_scalars()
     integer, parameter :: num_scalars = 73
     integer :: ier
     double precision, dimension(num_scalars) :: temp, temp2
-    double precision, allocatable  :: ptemp(:)
 
     ! Allreduce energy terms
     if(maxrank .gt. 1) then
@@ -437,19 +427,6 @@ subroutine distribute_scalars()
        wall_force_n0_z_halo = temp2(72)
        helicity        = temp2(73)
 
-       if(ipellet_abl.gt.0) then
-          allocate(ptemp(npellets))
-
-          ptemp = nsource_pel
-          call mpi_allreduce(ptemp, nsource_pel, npellets, MPI_DOUBLE_PRECISION,  &
-                             MPI_SUM, MPI_COMM_WORLD, ier)
-          ptemp = temp_pel
-          call mpi_allreduce(ptemp, temp_pel, npellets, MPI_DOUBLE_PRECISION,  &
-                            MPI_SUM, MPI_COMM_WORLD, ier)
-
-          deallocate(ptemp)
-       end if
-
     endif
 
   end subroutine distribute_scalars
@@ -626,31 +603,6 @@ subroutine second(tcpu)
   end subroutine second
 
 
-!   Added 1/1/2016 to get consistency between 2D,3D,Cyl,Tor
-subroutine tpi_factors(tpifac,tpirzero)
-  use basic
-  use math
-  implicit none
-  real, intent(out) :: tpifac, tpirzero
-  if(nplanes.eq.1) then
-     if(itor.eq.1) then
-        tpifac = 1.
-        tpirzero = 1.
-     else
-        tpifac = 1./rzero
-        tpirzero = 1.
-     endif
-  else
-     if(itor.eq.1) then
-        tpifac = twopi
-        tpirzero = twopi
-     else
-        tpifac = twopi
-        tpirzero = twopi*rzero
-     endif
-  endif
-end subroutine tpi_factors
-
 ! ======================================================================
 ! calculate scalars
 ! -----------------
@@ -669,7 +621,6 @@ subroutine calculate_scalars()
   use boundary_conditions
   use math
   use gyroviscosity
-  use pellet
 
   implicit none
  
@@ -755,13 +706,9 @@ subroutine calculate_scalars()
   call finalize(field_vec)
 
   numelms = local_elements()
-  
-  if(ipellet.ne.0) call calculate_Lor_vol
 
-  ! BCL Warning: nsource_pel and temp_pel are now vectors
-  !              this compiles, but may break at runtime for OpenMP (OMP=1)
 !$OMP PARALLEL DO PRIVATE(mr,dum1,ier,is_edge,n,iedge,idim,izone,izonedim,i) &
-!$OMP& REDUCTION(+:ekinp,ekinpd,ekinph,ekint,ekintd,ekinth,ekin3,ekin3d,ekin3h,wallcur,emagp,emagpd,emagph,emagt,emagtd,emagth,emag3,area,parea,totcur,pcur,m_iz,tflux,pflux,tvor,volume,pvol,totden,pden,totrad,linerad,bremrad,ionrad,reckrad,recprad,totre,nsource,epotg,tmom,pmom,bwb2,efluxp,efluxt,efluxs,efluxk,tau_em,tau_sol,tau_com,tau_visc,tau_gyro,tau_parvisc,nfluxd,nfluxv,xray_signal,Lor_vol,nsource_pel,temp_pel,wall_force_n0_x,wall_force_n0_y,wall_force_n0_z,wall_force_n1_x,wall_force_n1_y,wall_force_n1_z,totne,w_pe,pcur_co,pcur_sn,m_iz_co,m_iz_sn,w_m,w_p,wall_force_n0_x_halo,wall_force_n0_z_halo,helicity)
+!$OMP& REDUCTION(+:ekinp,ekinpd,ekinph,ekint,ekintd,ekinth,ekin3,ekin3d,ekin3h,wallcur,emagp,emagpd,emagph,emagt,emagtd,emagth,emag3,area,parea,totcur,pcur,m_iz,tflux,pflux,tvor,volume,pvol,totden,pden,totrad,linerad,bremrad,ionrad,reckrad,recprad,totre,nsource,epotg,tmom,pmom,bwb2,efluxp,efluxt,efluxs,efluxk,tau_em,tau_sol,tau_com,tau_visc,tau_gyro,tau_parvisc,nfluxd,nfluxv,xray_signal,wall_force_n0_x,wall_force_n0_y,wall_force_n0_z,wall_force_n1_x,wall_force_n1_y,wall_force_n1_z,totne,w_pe,pcur_co,pcur_sn,m_iz_co,m_iz_sn,w_m,w_p,wall_force_n0_x_halo,wall_force_n0_z_halo,helicity)
   do itri=1,numelms
 
      !call zonfac(itri, izone, izonedim)
@@ -912,25 +859,8 @@ subroutine calculate_scalars()
 #endif
 
      ! particle source
-     if(idens.eq.1) then        
-        nsource = nsource - twopi*int1(sig79)/tpifac
+     if(idens.eq.1) nsource = nsource - twopi*int1(sig79)/tpifac
      
-        ! Pellet radius and density/temperature at the pellet surface
-        if(ipellet_abl.gt.0) then
-           do ip=1,npellets
-              if(r_p(ip).ge.1e-8) then
-                 ! weight density/temp by pellet distribution (normalized)
-                 temp79a = pellet_distribution(ip, x_79, phi_79, z_79, real(pt79(:,OP_1)), 1)
-                 nsource_pel(ip) = nsource_pel(ip) + twopi*int2(net79(:,OP_1),temp79a)/tpifac
-                 temp_pel(ip) = temp_pel(ip) + twopi*int2(pet79(:,OP_1)/net79(:,OP_1),temp79a)*p0_norm/(1.6022e-12*n0_norm*tpifac)
-              else
-                 nsource_pel(ip) = 0.
-                 temp_pel(ip) = 0.
-              end if
-           end do
-       endif
-     endif
-
      ! gravitational potential energy
      epotg = epotg + grav_pot()
 
@@ -1020,8 +950,6 @@ subroutine calculate_scalars()
 
   call distribute_scalars
 
-  if(ipellet_abl.gt.0) call calculate_ablation
-
   ekin = ekinp + ekint + ekin3
   emag = emagp + emagt + emag3
   ekind = ekinpd + ekintd + ekin3d
@@ -1067,73 +995,10 @@ subroutine calculate_scalars()
      print *, "  Ionization loss = ", ionrad
      print *, "  Recombination radiation (kinetic) = ", reckrad
      print *, "  Recombination radiation (potential) = ", recprad
-     if(ipellet_abl.gt.0 .and. iprint.ge.2) then
-        do ip=1,npellets
-           print *, "  Pellet #", ip
-           print *, "    particles injected = ",pellet_rate(ip)*dt*(n0_norm*l0_norm**3)
-           print *, "    radius (in cm) = ", r_p(ip)*l0_norm
-           print *, "    local electron temperature (in eV) = ", temp_pel(ip)
-           print *, "    local electron density (in ne14) = ", nsource_pel(ip)
-           print *, "    rpdot (in cm/s) = ", rpdot(ip)*l0_norm/t0_norm
-           print *, "    Lor_vol = ", Lor_vol(ip)
-           print *, "    R position: ", pellet_r(ip)*l0_norm
-           print *, "    phi position: ", pellet_phi(ip)
-           print *, "    Z position: ", pellet_z(ip)*l0_norm
-        end do
-     endif
   endif
 
 end subroutine calculate_scalars
 
-
-subroutine calculate_Lor_vol()
-
-  use basic
-  use mesh_mod
-  use m3dc1_nint
-  use math
-  use pellet
-
-  implicit none
- 
-  include 'mpif.h'
-
-  integer :: itri, numelms, ier
-  integer :: is_edge(3)  ! is inode on boundary
-  real :: tpifac,tpirzero
-  integer :: izone, izonedim
-  real, allocatable :: temp(:)
-  integer :: ip
-
-  call tpi_factors(tpifac,tpirzero)
-
-  numelms = local_elements()
-
-  allocate(temp(npellets))
-  temp = 0.
-
-  do itri=1,numelms
-
-     call m3dc1_ent_getgeomclass(2, itri-1,izonedim,izone)
-     if(izone.ne.1) cycle
-     call define_element_quadrature(itri, int_pts_diag, int_pts_tor)
-     call define_fields(itri, FIELD_P, 0, 0)
-
-     ! perform volume integral of pellet cloud (without normalization)
-     do ip=1,npellets
-        temp79a  = pellet_distribution(ip, x_79, phi_79, z_79, real(pt79(:,OP_1)), 0)
-        temp(ip) = temp(ip) + twopi*int1(temp79a)/tpifac
-     end do
-
-  end do
-
-  call mpi_allreduce(temp, Lor_vol, npellets, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ier )
-
-  deallocate(temp)
-
-end subroutine calculate_Lor_vol
-
-
 !======================================================================
 ! magnetic_region
 ! ~~~~~~~~~~~~~~~
diff --git a/unstructured/hdf5_output.f90 b/unstructured/hdf5_output.f90
index 0d52d9bc3..eefce77eb 100644
--- a/unstructured/hdf5_output.f90
+++ b/unstructured/hdf5_output.f90
@@ -552,11 +552,13 @@ subroutine output_scalar(parent_id, name, value, t, error)
     integer(HSIZE_T) :: chunk_size(1) = (/ 100 /)
     integer(HSIZE_T) :: dims(1)
     integer(HSIZE_T) :: maxdims(1)
-    integer(HSIZE_T), parameter :: local_dims(1) = (/ 1 /)
-    integer(HSIZE_T), dimension(1,1) :: coord
-    integer(SIZE_T), parameter :: num_elements = 1
+    integer(HSIZE_T) :: local_dims(1)
+    integer(HSIZE_T), allocatable :: coord(:,:)
+    integer(SIZE_T)  :: num_elements
     integer(HID_T) :: memspace, filespace, dset_id, p_id, plist_id
-    real :: values(1)
+    integer(HSIZE_T) :: curdim
+    integer :: i
+    real, allocatable :: values(:)
     logical :: exists
 
 #ifdef USETAU
@@ -566,12 +568,12 @@ subroutine output_scalar(parent_id, name, value, t, error)
 
     dims(1) = t+1
     maxdims(1) = H5S_UNLIMITED_F
-    values(1) = value
-    coord(1,1) = t + 1
-    
+
+    !@if(myrank.eq.0 .and. iprint.ge.2) print *, "   ", name
     call h5lexists_f(parent_id, name, exists, error)
 
     if(.not.exists) then
+       num_elements = 1
        call h5screate_simple_f(1, dims, filespace, error, maxdims)
        call h5pcreate_f(H5P_DATASET_CREATE_F, p_id, error)
        call h5pset_chunk_f(p_id, 1, chunk_size, error)
@@ -584,20 +586,35 @@ subroutine output_scalar(parent_id, name, value, t, error)
        end if
        call h5pclose_f(p_id, error)
        call h5sclose_f(filespace, error)
+       curdim = dims(1)
     else
        call h5dopen_f(parent_id, name, dset_id, error)
-       call h5dset_extent_f(dset_id, dims, error)
+       call h5dget_space_f(dset_id, filespace, error)
+       call h5sget_simple_extent_npoints_f(filespace, curdim, error)
+       call h5sclose_f(filespace, error)
+       if(dims(1).gt.curdim) then
+          num_elements = 1
+          call h5dset_extent_f(dset_id, dims, error)
+          curdim = dims(1)
+       else
+          num_elements = curdim - dims(1) + 1
+       end if
     endif
 
+    local_dims(1) = num_elements
+    allocate(values(num_elements))
+    allocate(coord(1,num_elements))
+    values = sqrt(-1.)
+    values(1) = value
+    coord(1,:) = (/ (i, i = dims(1), curdim) /)
+
     if(myrank.eq.0) then
        call h5screate_simple_f(1, local_dims, memspace, error)
        call h5dget_space_f(dset_id, filespace, error)
-       call h5sselect_elements_f(filespace, H5S_SELECT_SET_F, 1, &
-            num_elements, coord, error)
-       
        call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error)
        call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_INDEPENDENT_F, error)
-
+       call h5sselect_elements_f(filespace, H5S_SELECT_SET_F, 1, &
+            num_elements, coord, error)
        call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, values, local_dims, error, &
             file_space_id=filespace, mem_space_id=memspace, xfer_prp=plist_id)
 
@@ -608,6 +625,9 @@ subroutine output_scalar(parent_id, name, value, t, error)
     endif
     call h5dclose_f(dset_id, error)
 
+    deallocate(values)
+    deallocate(coord)
+
   end subroutine output_scalar
 
   ! read_scalar
@@ -873,7 +893,6 @@ subroutine read_1dextendarr(parent_id, name, value, NMAX, t, error)
     integer(HSIZE_T) :: dims(2), maxdims(2), local_dims(2), off(2)
     integer(SIZE_T) :: num_elements
     integer(HID_T) :: memspace, filespace, dset_id, p_id, plist_id
-    logical :: exists
 
 #ifdef USETAU
     integer :: dummy     ! this is necessary to prevent TAU from
diff --git a/unstructured/init_conds.f90 b/unstructured/init_conds.f90
index b098dd725..2ce17b526 100644
--- a/unstructured/init_conds.f90
+++ b/unstructured/init_conds.f90
@@ -514,7 +514,6 @@ subroutine initial_conditions()
   use init_common
   use kprad_m3dc1
   use pellet
-  use diagnostics
   use cylinder
 
   implicit none
diff --git a/unstructured/makefile b/unstructured/makefile
index 23feecc6a..08802e8b6 100644
--- a/unstructured/makefile
+++ b/unstructured/makefile
@@ -143,7 +143,7 @@ OBJS := fftw_fortran.o read_namelist.o gsl_wrapper.o \
 	iterdb.o read_gyro.o read_neo.o radiation.o \
 	element.o $(V_OBJ) field.o nintegrate_mod.o \
 	M3Dmodules.o resistive_wall.o kprad.o \
-	 m3dc1_nint.o boundary.o gyroviscosity.o bootstrap.o \
+	m3dc1_nint.o boundary.o gyroviscosity.o bootstrap.o \
 	metricterms_new.o two_fluid.o harned_mikic.o biharmonic.o \
 	electric_field.o pellet.o parallel_heat_flux.o \
 	temperature_plots.o electrostatic_potential.o \
@@ -160,8 +160,7 @@ OBJS := fftw_fortran.o read_namelist.o gsl_wrapper.o \
         init_eqdsk.o init_dskbal.o init_jsolver.o \
 	init_3dwave.o init_3ddiffusion.o init_frs.o init_ftz.o init_eigen.o \
 	init_intkink.o init_lz.o init_kstar.o init_basicq.o init_cyl.o \
-	transport.o \
-	auxiliary_fields.o  \
+	auxiliary_fields.o transport.o \
 	time_step_split.o time_step_unsplit.o \
         hdf5_output.o particle.o   \
 	time_step.o output.o restart_hdf5.o \
diff --git a/unstructured/newpar.f90 b/unstructured/newpar.f90
index 04570b33c..6e159af33 100644
--- a/unstructured/newpar.f90
+++ b/unstructured/newpar.f90
@@ -276,8 +276,7 @@ Program Reducedquintic
   if(ntime.eq.0 .or. (ntime.eq.ntime0 .and. eqsubtract.eq.1)) then
 
      if(eqsubtract.eq.1) then
-        if(myrank.eq.0 .and. iprint.ge.2) print *, "  transport coefficients"
-        call define_transport_coefficients
+        call define_transport_coefficients(0)
         call derived_quantities(0)
         if(iwrite_aux_vars.eq.1) call calculate_auxiliary_fields(0)
      end if
@@ -296,8 +295,7 @@ Program Reducedquintic
 
   ! Calculate all quantities derived from basic fields
   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  if(myrank.eq.0 .and. iprint.ge.2) print *, "  transport coefficients"
-  call define_transport_coefficients
+  call define_transport_coefficients(1)
   call derived_quantities(1)
 
 
@@ -678,8 +676,7 @@ end subroutine smooth_fields
 ! ======================================================================
 ! derived_quantities
 ! ~~~~~~~~~~~~~~~~~~
-! calculates all derived quantities, including auxiliary fields
-! and scalars
+! calculates all derived quantities, including scalars
 ! ======================================================================
 subroutine derived_quantities(ilin)
   use basic
@@ -699,6 +696,14 @@ subroutine derived_quantities(ilin)
 
   vectype :: temp
 
+  ! Update density and temperature
+  if(myrank.eq.0 .and. iprint.ge.2) print *, "  Calculate electron density and temperatures"
+  call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract)
+  if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) &
+       call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), &
+                                   pe_field(ilin), p_field(ilin), ne_field(ilin), &
+                                   den_field(ilin),eqsubtract)
+  
   ! Find lcfs
   ! ~~~~~~~~~
   if(myrank.eq.0 .and. iprint.ge.2) print *, "  finding lcfs"
@@ -739,25 +744,14 @@ subroutine derived_quantities(ilin)
     endif
   else
     if(myrank.eq.0 .and. iprint.ge.1) then
-      write(*,'(A,2e12.4)') ' no temperatue maximum found near ',xmag,zmag
+      write(*,'(A,2e12.4)') ' no temperature maximum found near ',xmag,zmag
     endif
   endif
 
-  ! Electron density
-  call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract)
-
-
   ! Define auxiliary fields
   ! ~~~~~~~~~~~~~~~~~~~~~~~
   if(myrank.eq.0 .and. itimer.eq.1) call second(tstart)
 
-  if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) then
-     if(myrank.eq.0 .and. iprint.ge.2) print *, "  temperatures"
-     call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), &
-          pe_field(ilin), p_field(ilin), ne_field(ilin), den_field(ilin), &
-          eqsubtract)
-  end if
-
   !   toroidal current
   if(myrank.eq.0 .and. iprint.ge.2) print *, "  toroidal current"
   if(inocurrent_tor.eq.1) then
diff --git a/unstructured/pellet.f90 b/unstructured/pellet.f90
index 80f615103..4189a6438 100644
--- a/unstructured/pellet.f90
+++ b/unstructured/pellet.f90
@@ -42,7 +42,7 @@ module pellet
   subroutine pellet_init()
     use basic
     use read_ascii
-!    use diagnostics
+
     implicit none
     character(LEN=10), parameter :: pellet_filename = 'pellet.dat'
 
@@ -114,10 +114,56 @@ subroutine pellet_init()
 
   end subroutine pellet_init
 
+  subroutine calculate_Lor_vol()
+
+    use basic
+    use mesh_mod
+    use m3dc1_nint
+    use math
+
+    implicit none
+
+    include 'mpif.h'
+
+    integer :: itri, numelms, ier
+    real :: tpifac,tpirzero
+    integer :: izone, izonedim
+    real, allocatable :: temp(:)
+    integer :: ip
+
+    call tpi_factors(tpifac,tpirzero)
+
+    numelms = local_elements()
+
+    allocate(temp(npellets))
+    temp = 0.
+
+    do itri=1,numelms
+
+       call m3dc1_ent_getgeomclass(2, itri-1,izonedim,izone)
+       if(izone.ne.1) cycle
+       call define_element_quadrature(itri, int_pts_diag, int_pts_tor)
+       call define_fields(itri, FIELD_P, 0, 0)
+
+       ! perform volume integral of pellet cloud (without normalization)
+       do ip=1,npellets
+          temp79a  = pellet_distribution(ip, x_79, phi_79, z_79, real(pt79(:,OP_1)), 0)
+          temp(ip) = temp(ip) + twopi*int1(temp79a)/tpifac
+       end do
+
+    end do
+
+    call mpi_allreduce(temp, Lor_vol, npellets, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ier )
+
+    deallocate(temp)
+
+  end subroutine calculate_Lor_vol
+
+
   vectype elemental function pellet_distribution(ip, r, phi, z, pres, inorm)
     use math
     use basic
-!    use diagnostics
+
     implicit none
     integer, intent(in) :: ip
     real, intent(in) :: r, phi, z, pres
@@ -215,7 +261,7 @@ end function pellet_distribution
 
   subroutine pellet_advance
     use basic
-!    use diagnostics
+
     implicit none
     real, allocatable :: x(:), y(:)
 
diff --git a/unstructured/time_step.f90 b/unstructured/time_step.f90
index 17d69414e..f851cda32 100644
--- a/unstructured/time_step.f90
+++ b/unstructured/time_step.f90
@@ -107,6 +107,11 @@ subroutine onestep
      calc_matrices = 0
   endif
 
+  if(ipellet.ne.0) then
+     if(myrank.eq.0 .and. iprint.ge.2) print *, "  define pellet source"
+     call define_pellet_source(1)
+  end if
+
   ! Advance impurity charge states
   if(myrank.eq.0 .and. itimer.eq.1) call second(tstart)
   call kprad_ionize(dt)
@@ -116,7 +121,10 @@ subroutine onestep
   endif
 
   if(myrank.eq.0 .and. iprint.ge.2) print *, "  transport coefficients"
-  call define_transport_coefficients
+  call define_transport_coefficients(1)
+
+  if(myrank.eq.0 .and. iprint.ge.2) print *, "  sources"
+  call define_sources(1)
 
   ! start of loop to repeat timestep if max iterations exceeded in 3D
   do irepeat = 1, max_repeat
@@ -177,8 +185,6 @@ subroutine onestep
   dtold = dt
   if(ntime.gt.1 .and. linear.eq.0) call variable_timestep
 
-  call pellet_advance
-
   call runaway_advance
 
   ! copy time advance vectors to field data
@@ -187,9 +193,6 @@ subroutine onestep
 ! if(eqsubtract.eq.0) call subtract_axi    !DEBUG
   call export_time_advance_vectors
 
-  ! Calculate all quantities derived from basic fields
-  call derived_quantities(1)
-
   ! Advect impurity charge states
   if(myrank.eq.0 .and. itimer.eq.1) call second(tstart)
   call kprad_advect(dt)
@@ -198,6 +201,10 @@ subroutine onestep
      t_kprad = t_kprad + tend - tstart
   endif
 
+  ! Calculate all quantities derived from basic fields
+  call derived_quantities(1)
+
+
   ! Conserve toroidal flux
   if(iconstflux.eq.1 .and. numvar.ge.2) then
      call conserve_flux
diff --git a/unstructured/time_step_split.f90 b/unstructured/time_step_split.f90
index 764e6d75e..4e93d08f7 100644
--- a/unstructured/time_step_split.f90
+++ b/unstructured/time_step_split.f90
@@ -864,16 +864,17 @@ subroutine step_split(calc_matrices)
      call destroy_vector(temp)
      call destroy_vector(temp2)
 
-     if(irecalc_eta.eq.1) then
-        call export_time_advance_vectors_split
-        call define_transport_coefficients
-     end if
-
-     ! Electron temperature
+     ! Electron density
      neold_vec = ne_vec
      call calculate_ne(1, den_v, ne_v, eqsubtract)
      ne_field(1) = ne_v    ! This is needed so that boundary_te works properly
      den_field(1) = den_v
+
+     if(irecalc_eta.eq.1) then
+        call export_time_advance_vectors_split
+        call define_transport_coefficients(1)
+     end if
+
   endif    ! on idens=1
      
   !
@@ -1103,7 +1104,7 @@ subroutine step_split(calc_matrices)
      call matvecmult(d2_mat,phi_vec,b2_phi)
      call add(b1_phi, b2_phi)
 
-     ! Inculde density terms
+     ! Include density terms
      if(idens.eq.1) then
         call matvecmult(r42_mat,ne_vec,b2_phi)
         call add(b1_phi, b2_phi)
@@ -1204,7 +1205,7 @@ subroutine step_split(calc_matrices)
         call export_time_advance_vectors_split
         ! redefine transport coefficients with new den/pe values
         call lcfs(psi_field(1))
-        call define_transport_coefficients
+        call define_transport_coefficients(1)
         ! revert fields to old values
         phi_vec = b2_phi
         
diff --git a/unstructured/transport.f90 b/unstructured/transport.f90
index 57ff68eb0..b1c728882 100644
--- a/unstructured/transport.f90
+++ b/unstructured/transport.f90
@@ -1104,7 +1104,7 @@ end function bs_func
 
 ! define_transport_coefficients
 ! =============================
-subroutine define_transport_coefficients()
+subroutine define_transport_coefficients(ilin)
 
   use basic
   use arrays
@@ -1115,21 +1115,20 @@ subroutine define_transport_coefficients()
   use pellet
   use diagnostics
   use kprad_m3dc1
+  use auxiliary_fields
 
   implicit none
 
   include 'mpif.h'
 
+  integer, intent(in) :: ilin
   integer :: itri, izone
   integer :: numelms, def_fields,ier
 
   logical, save :: first_time = .true.
-  logical :: solve_sigma, solve_kappa, solve_visc, solve_resistivity, &
-       solve_visc_e, solve_q, solve_totrad, solve_linerad, solve_bremrad, &
-       solve_ionrad, solve_reckrad, solve_recprad, solve_cd, solve_f, &
-       solve_fp
+  logical :: solve_kappa, solve_visc, solve_resistivity, solve_visc_e
 
-  integer, parameter :: num_scalars = 15
+  integer, parameter :: num_scalars = 4
   integer, dimension(num_scalars) :: temp, temp2
   vectype, dimension(dofs_per_element) :: dofs
 
@@ -1137,15 +1136,186 @@ subroutine define_transport_coefficients()
   if((linear.eq.1).and.(.not.first_time)) return
   first_time = .false.
 
-  if(myrank.eq.0 .and. iprint.ge.1) &
-       print *, "Calculating transport coefficients"
+  if(myrank.eq.0 .and. iprint.ge.1) print *, "Calculating transport coefficients"
+
+  if(myrank.eq.0 .and. iprint.ge.2) print *, "  Calculate electron density and temperatures"
+  call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract)
+  if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) &
+       call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), &
+                                   pe_field(ilin), p_field(ilin), ne_field(ilin), &
+                                   den_field(ilin),eqsubtract)
 
   ! which transport coefficients need matrix solve
   solve_resistivity = .false.
   solve_visc = .false.
   solve_kappa = .false.
-  solve_sigma = .false.
   solve_visc_e = .false.
+
+  ! clear variables
+  resistivity_field = 0.
+  kappa_field = 0.
+  visc_field = 0.
+  if(ibootstrap.ne.0) visc_e_field = 0.
+
+  call finalize(field0_vec)
+  call finalize(field_vec)
+
+  ! specify which primitive fields are to be evalulated
+  def_fields = FIELD_N + FIELD_PE + FIELD_P + FIELD_PSI + FIELD_I + FIELD_B2I
+  if(itemp.ge.1) def_fields = def_fields + FIELD_TE
+  if(iresfunc.eq.2 .or. iresfunc.eq.3 .or. iresfunc.eq.4) &
+       def_fields = def_fields + FIELD_ETA
+  if(ikappafunc.eq.5) def_fields = def_fields + FIELD_KAP
+  if(ivisfunc.eq.3) def_fields = def_fields + FIELD_MU
+
+  if(myrank.eq.0 .and. iprint.ge.2) print *, '  defining...'
+
+  ! Calculate RHS
+  numelms = local_elements()
+!$OMP PARALLEL DO &
+!$OMP& PRIVATE(dofs)
+  do itri=1,numelms
+
+     call define_element_quadrature(itri, int_pts_aux, 5)
+     call define_fields(itri, def_fields, 1, linear)
+
+     call get_zone(itri, izone)
+
+
+     dofs = resistivity_func()
+     if(.not.solve_resistivity) solve_resistivity = any(dofs.ne.0.)
+
+!$OMP CRITICAL
+     if(solve_resistivity) &
+          call vector_insert_block(resistivity_field%vec,itri,1,dofs,VEC_ADD)
+!$OMP END CRITICAL
+
+     dofs = kappa_func()
+     if(.not.solve_kappa) solve_kappa = any(dofs.ne.0.)
+!$OMP CRITICAL
+     if(solve_kappa) &
+          call vector_insert_block(kappa_field%vec,itri,1,dofs,VEC_ADD)
+!$OMP END CRITICAL
+
+     dofs = viscosity_func()
+     if(.not.solve_visc) solve_visc = any(dofs.ne.0.)
+!$OMP CRITICAL
+     if(solve_visc) &
+          call vector_insert_block(visc_field%vec,itri,1,dofs,VEC_ADD)
+!$OMP END CRITICAL
+
+     if(ibootstrap.ne.0) then
+        dofs = electron_viscosity_func()
+        if(.not.solve_visc_e) solve_visc_e = any(dofs.ne.0.)
+!$OMP CRITICAL
+        if(solve_visc_e) &
+             call vector_insert_block(visc_e_field%vec,itri,1,dofs,VEC_ADD)
+!$OMP END CRITICAL
+     end if
+  end do
+!$OMP END PARALLEL DO
+
+
+  ! Solve all the variables that have been defined
+  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+  ! make sure all processes agree on what needs to be solved
+  if(maxrank.gt.1) then
+     temp = 0
+     temp2 = 0
+     if(solve_resistivity) temp(1) = 1
+     if(solve_kappa)       temp(2) = 1
+     if(solve_visc)        temp(3) = 1
+     if(solve_visc_e)      temp(4) = 1
+
+     call mpi_allreduce(temp, temp2, num_scalars, MPI_INTEGER, &
+          MPI_MAX, MPI_COMM_WORLD, ier)
+
+     solve_resistivity = temp2(1).eq.1
+     solve_kappa       = temp2(2).eq.1
+     solve_visc        = temp2(3).eq.1
+     solve_visc_e      = temp2(4).eq.1
+  end if
+
+  if(myrank.eq.0 .and. iprint.ge.1) print *, ' solving...'
+
+  if(solve_resistivity) then
+     if(myrank.eq.0 .and. iprint.ge.1) print *, '  resistivity'
+     call newvar_solve(resistivity_field%vec, mass_mat_lhs)
+  end if
+
+  if(solve_kappa) then
+     if(myrank.eq.0 .and. iprint.ge.1) print *, '  kappa'
+     call newvar_solve(kappa_field%vec, mass_mat_lhs)
+  endif
+
+  if(solve_visc) then
+     if(myrank.eq.0 .and. iprint.ge.1) print *, '  viscosity'
+     call newvar_solve(visc_field%vec, mass_mat_lhs)
+  endif
+
+  if(solve_visc_e) then
+     if(myrank.eq.0 .and. iprint.ge.1) print *, '  electron viscosity'
+     call newvar_solve(visc_e_field%vec, mass_mat_lhs)
+  endif
+
+  ! the "compressible" viscosity is the same as the "incompressible"
+  ! viscosity up to a constant
+  visc_c_field = visc_field
+
+  ! add in constant components
+  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
+  call add(resistivity_field, etar*eta_fac)
+  call add(visc_field, amu)
+  call add(visc_c_field, amuc)
+
+  if(myrank.eq.0 .and. iprint.ge.2) print *, 'done define_transport_coefficients'
+
+end subroutine define_transport_coefficients
+
+
+! define_sources
+! =============================
+subroutine define_sources(ilin)
+
+  use basic
+  use arrays
+  use m3dc1_nint
+  use newvar_mod
+  use sparse
+  use neutral_beam
+  use pellet
+  use diagnostics
+  use kprad_m3dc1
+  use auxiliary_fields
+
+  implicit none
+
+  include 'mpif.h'
+
+  integer, intent(in) :: ilin
+  integer :: itri, izone
+  integer :: numelms, def_fields,ier
+
+  logical :: solve_sigma, solve_q, solve_totrad, solve_linerad, solve_bremrad, &
+             solve_ionrad, solve_reckrad, solve_recprad, solve_cd, solve_f, solve_fp
+
+  integer, parameter :: num_scalars = 11
+  integer, dimension(num_scalars) :: temp, temp2
+  vectype, dimension(dofs_per_element) :: dofs
+
+  if(myrank.eq.0 .and. iprint.ge.1) print *, "Calculating sources"
+
+  if(myrank.eq.0 .and. iprint.ge.2) print *, "  Calculate electron density and temperatures"
+  call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract)
+  if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) &
+       call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), &
+                                   pe_field(ilin), p_field(ilin), ne_field(ilin), &
+                                   den_field(ilin),eqsubtract)
+
+
+  ! which transport coefficients need matrix solve
+  solve_sigma = .false.
   solve_f = .false.
   solve_q = .false.
   solve_totrad = .false.
@@ -1158,10 +1328,6 @@ subroutine define_transport_coefficients()
   solve_fp = .false.
 
   ! clear variables
-  resistivity_field = 0.
-  kappa_field = 0.
-
-  visc_field = 0.
   if(density_source) sigma_field = 0.  
   if(momentum_source) Fphi_field = 0.
   if(heat_source) Q_field = 0.
@@ -1174,7 +1340,6 @@ subroutine define_transport_coefficients()
      Recprad_field = 0.
   end if
   if(icd_source .gt. 0) cd_field = 0.
-  if(ibootstrap.ne.0) visc_e_field = 0.
   if(ipforce.gt.0) pforce_field = 0.
   if(ipforce.gt.0) pmach_field = 0.
 
@@ -1194,11 +1359,6 @@ subroutine define_transport_coefficients()
   if(iarc_source.ne.0) def_fields = def_fields + FIELD_WALL
 
   if(myrank.eq.0 .and. iprint.ge.2) print *, '  defining...'
-
-  if(ipellet.ne.0) then
-     ! make sure normalization for pellet_distribution defined
-     call calculate_Lor_vol
-  end if
   
   ! Calculate RHS
   numelms = local_elements()
@@ -1211,23 +1371,6 @@ subroutine define_transport_coefficients()
 
      call get_zone(itri, izone)
 
-
-     dofs = resistivity_func()
-     if(.not.solve_resistivity) solve_resistivity = any(dofs.ne.0.)
-
-!$OMP CRITICAL
-     if(solve_resistivity) &
-          call vector_insert_block(resistivity_field%vec,itri,1,dofs,VEC_ADD)
-!$OMP END CRITICAL
-
-     dofs = kappa_func()
-     if(.not.solve_kappa) solve_kappa = any(dofs.ne.0.)
-!$OMP CRITICAL
-     if(solve_kappa) &
-          call vector_insert_block(kappa_field%vec,itri,1,dofs,VEC_ADD)
-!$OMP END CRITICAL
-
-
      if(density_source) then
         dofs = sigma_func(izone)
         if(.not.solve_sigma) solve_sigma = any(dofs.ne.0.)
@@ -1237,13 +1380,6 @@ subroutine define_transport_coefficients()
 !$OMP END CRITICAL
      end if
 
-     dofs = viscosity_func()
-     if(.not.solve_visc) solve_visc = any(dofs.ne.0.)
-!$OMP CRITICAL
-     if(solve_visc) &
-          call vector_insert_block(visc_field%vec,itri,1,dofs,VEC_ADD)
-!$OMP END CRITICAL
-
      if(momentum_source) then 
         dofs = force_func(izone)
         if(.not.solve_f) solve_f = any(dofs.ne.0.)
@@ -1328,14 +1464,6 @@ subroutine define_transport_coefficients()
              call vector_insert_block(cd_field%vec,itri,1,dofs,VEC_ADD)
      end if
 
-     if(ibootstrap.ne.0) then
-        dofs = electron_viscosity_func()
-        if(.not.solve_visc_e) solve_visc_e = any(dofs.ne.0.)
-!$OMP CRITICAL
-        if(solve_visc_e) &
-             call vector_insert_block(visc_e_field%vec,itri,1,dofs,VEC_ADD)
-!$OMP END CRITICAL
-     end if
   end do
 !$OMP END PARALLEL DO
 
@@ -1347,69 +1475,41 @@ subroutine define_transport_coefficients()
   if(maxrank.gt.1) then 
      temp = 0
      temp2 = 0
-     if(solve_resistivity) temp(1) = 1
-     if(solve_kappa)       temp(2) = 1
-     if(solve_sigma)       temp(3) = 1
-     if(solve_visc)        temp(4) = 1
-     if(solve_visc_e)      temp(5) = 1
-     if(solve_f)           temp(6) = 1
-     if(solve_q)           temp(7) = 1
-     if(solve_fp)          temp(8) = 1
-     if(solve_cd)          temp(9) = 1
-     if(solve_totrad)      temp(10) = 1
-     if(solve_linerad)      temp(11) = 1
-     if(solve_bremrad)     temp(12) = 1
-     if(solve_ionrad)      temp(13) = 1
-     if(solve_reckrad)     temp(14) = 1
-     if(solve_recprad)     temp(15) = 1
+     if(solve_sigma)       temp(1) = 1
+     if(solve_f)           temp(2) = 1
+     if(solve_q)           temp(3) = 1
+     if(solve_fp)          temp(4) = 1
+     if(solve_cd)          temp(5) = 1
+     if(solve_totrad)      temp(6) = 1
+     if(solve_linerad)     temp(7) = 1
+     if(solve_bremrad)     temp(8) = 1
+     if(solve_ionrad)      temp(9) = 1
+     if(solve_reckrad)     temp(10) = 1
+     if(solve_recprad)     temp(11) = 1
 
      call mpi_allreduce(temp, temp2, num_scalars, MPI_INTEGER, &
           MPI_MAX, MPI_COMM_WORLD, ier)
 
-     solve_resistivity = temp2(1).eq.1
-     solve_kappa       = temp2(2).eq.1
-     solve_sigma       = temp2(3).eq.1
-     solve_visc        = temp2(4).eq.1
-     solve_visc_e      = temp2(5).eq.1
-     solve_f           = temp2(6).eq.1
-     solve_q           = temp2(7).eq.1
-     solve_fp          = temp2(8).eq.1
-     solve_cd          = temp2(9).eq.1
-     solve_totrad      = temp2(10).eq.1
-     solve_linerad      = temp2(11).eq.1
-     solve_bremrad     = temp2(12).eq.1
-     solve_ionrad      = temp2(13).eq.1
-     solve_reckrad     = temp2(14).eq.1
-     solve_recprad     = temp2(15).eq.1
+     solve_sigma       = temp2(1).eq.1
+     solve_f           = temp2(2).eq.1
+     solve_q           = temp2(3).eq.1
+     solve_fp          = temp2(4).eq.1
+     solve_cd          = temp2(5).eq.1
+     solve_totrad      = temp2(6).eq.1
+     solve_linerad     = temp2(7).eq.1
+     solve_bremrad     = temp2(8).eq.1
+     solve_ionrad      = temp2(9).eq.1
+     solve_reckrad     = temp2(10).eq.1
+     solve_recprad     = temp2(11).eq.1
   end if
 
   if(myrank.eq.0 .and. iprint.ge.1) print *, ' solving...'
 
-  if(solve_resistivity) then
-     if(myrank.eq.0 .and. iprint.ge.1) print *, '  resistivity'
-     call newvar_solve(resistivity_field%vec, mass_mat_lhs)
-  end if
-
-  if(solve_kappa) then
-     if(myrank.eq.0 .and. iprint.ge.1) print *, '  kappa'
-     call newvar_solve(kappa_field%vec, mass_mat_lhs)
-  endif
-
-  if(solve_sigma) then
+   if(solve_sigma) then
      if(myrank.eq.0 .and. iprint.ge.1) print *, '  sigma'
      call newvar_solve(sigma_field%vec, mass_mat_lhs_dc)
   endif
 
-  if(solve_visc) then
-     if(myrank.eq.0 .and. iprint.ge.1) print *, '  viscosity'
-     call newvar_solve(visc_field%vec, mass_mat_lhs)
-  endif
-
-  if(solve_visc_e) then
-     if(myrank.eq.0 .and. iprint.ge.1) print *, '  electron viscosity'
-     call newvar_solve(visc_e_field%vec, mass_mat_lhs)
-  endif
-
   if(solve_f) then
      if(myrank.eq.0 .and. iprint.ge.1) print *, '  fphi'
      call newvar_solve(Fphi_field%vec, mass_mat_lhs_dc)
@@ -1458,27 +1558,53 @@ subroutine define_transport_coefficients()
   if(solve_fp) then
      if(myrank.eq.0 .and. iprint.ge.1) print *, '  pforce'
      call newvar_solve(pforce_field%vec, mass_mat_lhs)
+  endif
 
-     if(myrank.eq.0 .and. iprint.ge.1) print *, '  pmach'
-     call newvar_solve(pmach_field%vec, mass_mat_lhs)
+  if(myrank.eq.0 .and. iprint.ge.2) &
+       print *, 'done define_sources'
 
-  endif
+end subroutine define_sources
 
 
-  ! the "compressible" viscosity is the same as the "incompressible"
-  ! viscosity up to a constant
-  visc_c_field = visc_field
+! define pellet source
+! =============================
+subroutine define_pellet_source(ilin)
+
+  use basic
+  use mesh_mod
+  use arrays
+  use m3dc1_nint
+  use metricterms_new
+  use kprad_m3dc1
+  use pellet
+  use math
+  use auxiliary_fields
+
+  implicit none
+
+  include 'mpif.h'
+
+  integer, intent(in) :: ilin    ! 0 for equilibrium fields, 1 for perturbed
+  integer :: itri, numelms, def_fields, ier
+  integer :: ip
+  integer :: izone, izonedim
+  real :: tpifac, tpirzero
+  double precision, allocatable  :: ptemp(:)
 
+  if(myrank.eq.0 .and. iprint.ge.1) &
+       print *, "begin define_pellet_source"
   
-  ! add in constant components
-  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
-  call add(resistivity_field, etar*eta_fac)
-  call add(visc_field, amu)
-  call add(visc_c_field, amuc)
+  if(ntime.gt.1) then
+     if(myrank.eq.0 .and. iprint.ge.2) print *, "  Advance pellet position"
+     call pellet_advance
+  end if
 
+  ! Define normalization for pellet_distribution
+  call calculate_Lor_vol
 
   ! Read LP data
   if(iread_lp_source.eq.1) then
+     if(myrank.eq.0 .and. iprint.ge.2) print *, "  Read LP source"
      call read_lp_source('cloud.txt', ier)
      if(ier.ne.0) then
         if(myrank.eq.0) print *, 'Error reading LP source ', 'cloud.txt'
@@ -1486,9 +1612,83 @@ subroutine define_transport_coefficients()
      end if
   end if
 
-  if(myrank.eq.0 .and. iprint.ge.2) &
-       print *, 'done define_transport_coefficients'
+  if(myrank.eq.0 .and. iprint.ge.2) print *, "  Calculate electron density and temperatures"
+  call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract)
+  if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) &
+       call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), &
+                                   pe_field(ilin), p_field(ilin), ne_field(ilin), &
+                                   den_field(ilin),eqsubtract)
+
+  if(ipellet_abl.gt.0) then
+     if(myrank.eq.0 .and. iprint.ge.2) print *, "  Calculate ablation inputs"
+
+     call tpi_factors(tpifac,tpirzero)
+     nsource_pel = 0.
+     temp_pel = 0.
+
+     def_fields = FIELD_P + FIELD_N
+     call finalize(field0_vec)
+     call finalize(field_vec)
+     numelms = local_elements()
+
+!$OMP PARALLEL DO PRIVATE(itri,ier,izone,izonedim,ip) &
+!$OMP& REDUCTION(+:nsource_pel,temp_pel)
+     do itri=1,numelms
+
+        call m3dc1_ent_getgeomclass(2, itri-1,izonedim,izone)
+        if(izone.ne.1) cycle
+        call define_element_quadrature(itri, int_pts_diag, int_pts_tor)
+        call define_fields(itri, def_fields, 0, 0)
+
+        ! Pellet radius and density/temperature at the pellet surface
+        do ip=1,npellets
+           if(r_p(ip).ge.1e-8) then
+              ! weight density/temp by pellet distribution (normalized)
+              temp79a = pellet_distribution(ip, x_79, phi_79, z_79, real(pt79(:,OP_1)), 1)
+              nsource_pel(ip) = nsource_pel(ip) + twopi*int2(net79(:,OP_1),temp79a)/tpifac
+              temp_pel(ip) = temp_pel(ip) + twopi*int2(pet79(:,OP_1)/net79(:,OP_1),temp79a)*p0_norm/(1.6022e-12*n0_norm*tpifac)
+           else
+              nsource_pel(ip) = 0.
+              temp_pel(ip) = 0.
+           end if
+        end do
+     end do
+!$OMP END PARALLEL DO
+
+     allocate(ptemp(npellets))
+     ptemp = nsource_pel
+     call mpi_allreduce(ptemp, nsource_pel, npellets, MPI_DOUBLE_PRECISION,  &
+                        MPI_SUM, MPI_COMM_WORLD, ier)
+     ptemp = temp_pel
+     call mpi_allreduce(ptemp, temp_pel, npellets, MPI_DOUBLE_PRECISION,  &
+                        MPI_SUM, MPI_COMM_WORLD, ier)
+     deallocate(ptemp)
+
+     if(myrank.eq.0 .and. iprint.ge.2) print *, "  Calculate ablation rate"
+     call calculate_ablation
+
+  end if
+
+  if(myrank.eq.0 .and. iprint.ge.2) then
+     do ip=1,npellets
+        print *, "  Pellet #", ip
+        print *, "    particles injected = ",pellet_rate(ip)*dt*(n0_norm*l0_norm**3)
+        print *, "    Lor_vol = ", Lor_vol(ip)
+        print *, "    R position: ", pellet_r(ip)*l0_norm
+        print *, "    phi position: ", pellet_phi(ip)
+        print *, "    Z position: ", pellet_z(ip)*l0_norm
+        if(ipellet_abl.gt.0) then
+           print *, "    radius (in cm) = ", r_p(ip)*l0_norm
+           print *, "    local electron temperature (in eV) = ", temp_pel(ip)
+           print *, "    local electron density (in ne14) = ", nsource_pel(ip)
+           print *, "    rpdot (in cm/s) = ", rpdot(ip)*l0_norm/t0_norm
+        end if
+     end do
+  end if
+
+  if(myrank.eq.0 .and. iprint.ge.1) print *, 'done define_pellet_source'
+
+end subroutine define_pellet_source
 
-end subroutine define_transport_coefficients
 
 end module transport_coefficients