diff --git a/src/core_ocean/analysis_members/Registry_lagrangian_particle_tracking.xml b/src/core_ocean/analysis_members/Registry_lagrangian_particle_tracking.xml
index 63b797d774..911456c0cf 100644
--- a/src/core_ocean/analysis_members/Registry_lagrangian_particle_tracking.xml
+++ b/src/core_ocean/analysis_members/Registry_lagrangian_particle_tracking.xml
@@ -55,6 +55,50 @@
 			description="Specify whether particles should be reset when they enter the resetInsideRegionMaskValue1 mask."
 			possible_values=".false. or .true."
 		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_horizontal_interp" type="logical" default_value=".true." units="unitless"
+			description="If true, particles horizontally interpolate sample quantities."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_temperature" type="logical" default_value=".true." units="unitless"
+			description="If true, particles sample temperature."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_salinity" type="logical" default_value=".true." units="unitless"
+			description="If true, particles sample salinity."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_DIC" type="logical" default_value=".false." units="unitless"
+			description="If true, particles sample DIC."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_ALK" type="logical" default_value=".false." units="unitless"
+			description="If true, particles sample ALK."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_PO4" type="logical" default_value=".false." units="unitless"
+			description="If true, particles sample PO4."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_NO3" type="logical" default_value=".false." units="unitless"
+			description="If true, particles sample NO3."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_SiO3" type="logical" default_value=".false." units="unitless"
+			description="If true, particles sample SiO3."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_NH4" type="logical" default_value=".false." units="unitless"
+			description="If true, particles sample NH4."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_Fe" type="logical" default_value=".false." units="unitless"
+			description="If true, particles sample Fe."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_AM_lagrPartTrack_sample_O2" type="logical" default_value=".false." units="unitless"
+			description="If true, particles sample O2."
+			possible_values=".true. or .false."
+		/>
 	</nml_record>
 
 	<dims>
@@ -84,9 +128,12 @@
 			<var name="indexToParticleID"/>
 			<var name="currentBlock"/>
 			<var name="currentCell"/>
+			<var name="currentCellGlobalID"/>
 			<var name="xParticle"/>
 			<var name="yParticle"/>
 			<var name="zParticle"/>
+			<var name="lonParticle"/>
+			<var name="latParticle"/>
 			<var name="zLevelParticle"/>
 			<var name="xParticleReset"/>
 			<var name="yParticleReset"/>
@@ -248,7 +295,10 @@
 			 description="current block a particle is on"
 		/>
 		<var name="currentCell" type="integer" dimensions="nParticles Time" units="unitless"
-			 description="current cell a particle is on"
+			 description="current local cell a particle is on"
+		/>
+		<var name="currentCellGlobalID" type="integer" dimensions="nParticles Time" units="unitless"
+			 description="current global cell a particle is on"
 		/>
 		<var name="indexToParticleID" type="integer" dimensions="nParticles" units="unitless"
 			 description="designates global ID for a particle"
@@ -262,6 +312,12 @@
 		<var name="zParticle" type="real" dimensions="nParticles Time" units="m"
 			 description="z location of horizontal particle position"
 		/>
+		<var name="lonParticle" type="real" dimensions="nParticles Time" units="m"
+			 description="longitude location of horizontal particle position"
+		/>
+		<var name="latParticle" type="real" dimensions="nParticles Time" units="m"
+			 description="latitude location of horizontal particle position"
+		/>
 		<var name="zLevelParticle" type="real" dimensions="nParticles Time" units="m"
 			 description="z-level for vertical elevation of particle position"
 		/>
@@ -307,6 +363,37 @@
 		<var name="transfered" type ="integer" dimensions="nParticles Time" units="unitless"
 			 description="flag to monitor if the particle was transfered"
 		/>
+<!-- TEMP FIELDS THAT NEED MOVED TO NonHalo below after verifying interpolation is correct -->
+		<var name="particleTemperature" type="real" dimensions="nParticles Time" units="degrees Celsius"
+			 description="sampled temperature for particle"
+		/>
+		<var name="particleSalinity" type="real" dimensions="nParticles Time" units="grams salt per kilogram seawater"
+			 description="sampled salinity for particle"
+		/>
+		<var name="particleDIC" type="real" dimensions="nParticles Time" units="mmol m^{-3}"
+			 description="sampled dissolved inorganic carbon for particle"
+		/>
+		<var name="particleALK" type="real" dimensions="nParticles Time" units="meq m^{-3}"
+			 description="sampled alkalinity for particle"
+		/>
+		<var name="particlePO4" type="real" dimensions="nParticles Time" units="mmol m^{-3}"
+			 description="sampled dissolved inorganic phosphate for particle"
+		/>
+		<var name="particleNO3" type="real" dimensions="nParticles Time" units="meq m^{-3}"
+			 description="sampled dissolved inorganic nitrate for particle"
+		/>
+		<var name="particleSiO3" type="real" dimensions="nParticles Time" units="mmol m^{-3}"
+			 description="sampled dissolved inorganic silicate for particle"
+		/>
+		<var name="particleNH4" type="real" dimensions="nParticles Time" units="meq m^{-3}"
+			 description="sampled dissolved inorganic ammonia for particle"
+		/>
+		<var name="particleFe" type="real" dimensions="nParticles Time" units="mmol m^{-3}"
+			 description="sampled dissolved inorganic iron for particle"
+		/>
+		<var name="particleO2" type="real" dimensions="nParticles Time" units="meq m^{-3}"
+			 description="sampled dissolved oxygen for particle"
+		/>
 	</var_struct>
 <!--
 	<var_struct name="lagrPartTrackNonHalo" time_levs="1" packages="lagrPartTrackAMPKG">
diff --git a/src/core_ocean/analysis_members/mpas_ocn_lagrangian_particle_tracking.F b/src/core_ocean/analysis_members/mpas_ocn_lagrangian_particle_tracking.F
index 608f099c65..33f364a456 100644
--- a/src/core_ocean/analysis_members/mpas_ocn_lagrangian_particle_tracking.F
+++ b/src/core_ocean/analysis_members/mpas_ocn_lagrangian_particle_tracking.F
@@ -19,8 +19,8 @@
 !-----------------------------------------------------------------------
 #define COMMA ,
 #ifdef MPAS_DEBUG
-#define LIGHT_DEBUG_WRITE(Z)   print *, Z
-#define LIGHT_DEBUG_ALL_WRITE(Z)  print *, Z
+#define LIGHT_DEBUG_WRITE(Z) !  print *, Z
+#define LIGHT_DEBUG_ALL_WRITE(Z) ! print *, Z
 #else
 #define LIGHT_DEBUG_WRITE(Z)  ! print *, Z
 #define LIGHT_DEBUG_ALL_WRITE(Z) ! print *, Z
@@ -253,19 +253,21 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
       type (mpas_pool_type), pointer :: meshPool
       type (mpas_pool_type), pointer :: scratchPool
       type (mpas_pool_type), pointer :: diagnosticsPool
+      type (mpas_pool_type), pointer :: tracersPool
       type (mpas_pool_type), pointer :: lagrPartTrackFieldsPool, lagrPartTrackCellsPool, lagrPartTrackScratchPool
       integer, dimension(:), pointer :: cellOwnerBlock
-      integer, pointer :: currentBlock, ioBlock, indexToParticleID, transfered
+      integer, pointer :: currentBlock, ioBlock, indexToParticleID, transfered, currentCellGlobalID
       type (mpas_particle_list_type), pointer :: particlelist
       type (mpas_particle_type), pointer :: particle
 
       real (kind=RKIND), dimension(3) :: particlePosition, particleVelocity
       real (kind=RKIND), pointer :: xParticle, yParticle, zParticle, buoyancyParticle !, &
 !        lonVel, latVel, sumU, sumV, sumUU, sumUV, sumVV
+      real (kind=RKIND), pointer :: lonParticle, latParticle
       real (kind=RKIND), dimension(3) :: xSubStep, diffSubStep, diffParticlePosition
       real (kind=RKIND), pointer :: zLevelParticle
       real (kind=RKIND), dimension(:,:), pointer :: zTop, vertVelocityTop, zMid, areaBArray
-      real (kind=RKIND), dimension(:), pointer :: bottomDepth
+      real (kind=RKIND), dimension(:), pointer :: bottomDepth, ssh
       type (field2DReal), pointer :: normalVelocity, uVertexVelocity, vVertexVelocity, wVertexVelocity, layerThickness
 
       real (kind=RKIND), dimension(:,:), pointer :: uVertexVelocityArray, vVertexVelocityArray, wVertexVelocityArray, &
@@ -275,13 +277,14 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
       real(kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
 
       integer, dimension(:,:), pointer :: verticesOnCell, boundaryVertex
-      integer, dimension(:), pointer :: maxLevelCell
+      integer, dimension(:), pointer :: maxLevelCell, indexToCellID
 
       integer theVertex, iLevel, iLevelBuoyancy, aVertex, &
         nSteps, timeStep, subStep, subStepOrder, timeInterpOrder, aTimeLevel, nCellVertices, blockProc, arrayIndex
       integer, pointer :: nCells, nVertLevels, iCell
       integer, dimension(:), pointer :: nCellVerticesArray
       integer, dimension(:,:), pointer :: cellsOnCell
+      integer, dimension(:,:), pointer :: cellsOnVertex
       logical, pointer :: onSphere
       logical :: config_AM_lagrPartTrack_reset_particles
       character (len=StrKIND), pointer :: config_AM_lagrPartTrack_reset_criteria, &
@@ -306,6 +309,66 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
       logical :: resetParticle, resetParticleAny
       integer :: err_tmp
 
+      ! sampling variables
+      ! active tracers
+      integer, pointer :: index_temperature, &
+                          index_salinity
+      real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers
+
+      logical, pointer :: config_AM_lagrPartTrack_sample_temperature
+      logical, pointer :: config_AM_lagrPartTrack_sample_horizontal_interp
+      real (kind=RKIND), pointer :: particleTemperature
+
+      logical, pointer :: config_AM_lagrPartTrack_sample_salinity
+      real (kind=RKIND), pointer :: particleSalinity
+
+      ! BGC tracers
+      logical, pointer :: config_use_ecosysTracers
+      integer, pointer :: index_DIC, &
+                          index_ALK, &
+                          index_PO4, &
+                          index_NO3, &
+                          index_SiO3, &
+                          index_NH4, &
+                          index_Fe, &
+                          index_O2
+      real (kind=RKIND), dimension(:,:,:), pointer :: ecosysTracers
+
+      ! DIC
+      logical, pointer :: config_AM_lagrPartTrack_sample_DIC
+      real (kind=RKIND), pointer :: particleDIC
+
+      ! ALK
+      logical, pointer :: config_AM_lagrPartTrack_sample_ALK
+      real (kind=RKIND), pointer :: particleALK
+
+      ! PO4
+      logical, pointer :: config_AM_lagrPartTrack_sample_PO4
+      real (kind=RKIND), pointer :: particlePO4
+
+      ! NO3
+      logical, pointer :: config_AM_lagrPartTrack_sample_NO3
+      real (kind=RKIND), pointer :: particleNO3
+
+      ! SiO3
+      logical, pointer :: config_AM_lagrPartTrack_sample_SiO3
+      real (kind=RKIND), pointer :: particleSiO3
+
+      ! NH4
+      logical, pointer :: config_AM_lagrPartTrack_sample_NH4
+      real (kind=RKIND), pointer :: particleNH4
+
+      ! Fe
+      logical, pointer :: config_AM_lagrPartTrack_sample_Fe
+      real (kind=RKIND), pointer :: particleFe
+
+      ! O2
+      logical, pointer :: config_AM_lagrPartTrack_sample_O2
+      real (kind=RKIND), pointer :: particleO2
+
+      REAL (kind=RKIND), parameter :: eps=1e-12
+
+
       err = 0
 
       !! don't do compute for debugging purposes
@@ -323,6 +386,19 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
       call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_compute_interval', config_AM_lagrPartTrack_compute_interval)
       call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_output_stream', config_AM_lagrPartTrack_output_stream)
       call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_timeIntegration', timeIntegration)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_horizontal_interp', &
+                                             config_AM_lagrPartTrack_sample_horizontal_interp)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_temperature', config_AM_lagrPartTrack_sample_temperature)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_salinity', config_AM_lagrPartTrack_sample_salinity)
+      call mpas_pool_get_config(ocnConfigs, 'config_use_ecosysTracers', config_use_ecosysTracers)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_DIC', config_AM_lagrPartTrack_sample_DIC)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_ALK', config_AM_lagrPartTrack_sample_ALK)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_PO4', config_AM_lagrPartTrack_sample_PO4)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_NO3', config_AM_lagrPartTrack_sample_NO3)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_SiO3', config_AM_lagrPartTrack_sample_SiO3)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_NH4', config_AM_lagrPartTrack_sample_NH4)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_Fe', config_AM_lagrPartTrack_sample_Fe)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_O2', config_AM_lagrPartTrack_sample_O2)
 
       if (trim(config_AM_lagrPartTrack_reset_criteria) == 'none') then
         config_AM_lagrPartTrack_reset_particles = .False.
@@ -351,6 +427,7 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
         call mpas_pool_get_subpool(block % structs, 'scratch', scratchPool)
         call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool)
+        call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
         call mpas_pool_get_subpool(block % structs, 'lagrPartTrackFields', lagrPartTrackFieldsPool)
         call mpas_pool_get_subpool(block % structs, 'lagrPartTrackScratch', lagrPartTrackScratchPool)
         call mpas_pool_get_subpool(block % structs, 'lagrPartTrackCells', lagrPartTrackCellsPool)
@@ -374,10 +451,27 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
         ! now there is a time level so selection of the correct time level appears to be tricky
         ! the issue is large, nearly NAN normalVelocities with timeLevel=2
         call mpas_pool_get_field(statePool, 'normalVelocity', normalVelocity, timeLevel=timeLevel)
+        call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel=timeLevel)
         ! potentially a costly function call, but ensures halos are ok
         !call mpas_dmpar_exch_halo_field(normalVelocity)
         call mpas_pool_get_field(statePool, 'layerThickness', layerThickness, timeLevel=timeLevel)
 
+        ! sampling variables
+        call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1)
+        call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature)
+        call mpas_pool_get_dimension(tracersPool, 'index_salinity', index_salinity)
+        if (config_use_ecosysTracers) then
+          call mpas_pool_get_array(tracersPool, 'ecosysTracers', ecosysTracers, timeLevel)
+          call mpas_pool_get_dimension(tracersPool, 'index_DIC', index_DIC)
+          call mpas_pool_get_dimension(tracersPool, 'index_ALK', index_ALK)
+          call mpas_pool_get_dimension(tracersPool, 'index_PO4', index_PO4)
+          call mpas_pool_get_dimension(tracersPool, 'index_NO3', index_NO3)
+          call mpas_pool_get_dimension(tracersPool, 'index_SiO3', index_SiO3)
+          call mpas_pool_get_dimension(tracersPool, 'index_NH4', index_NH4)
+          call mpas_pool_get_dimension(tracersPool, 'index_Fe', index_Fe)
+          call mpas_pool_get_dimension(tracersPool, 'index_O2', index_O2)
+        end if
+
         call mpas_pool_get_array(meshPool, 'verticesOnCell', verticesOnCell)
         call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
         call mpas_pool_get_array(meshPool, 'boundaryVertex', boundaryVertex)
@@ -387,11 +481,14 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
         call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
         call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nCellVerticesArray)
         call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
+        call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex)
+        call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID)
 
         call mpas_pool_get_config(meshPool, 'on_a_sphere', onSphere)
         call mpas_pool_get_config(meshPool, 'sphere_radius', sphereRadius)
         !}}}
 
+        ! super-cycling:
         ! get the appropriate time scale for particle integration, e.g., dtSim
         if (trim(config_AM_lagrPartTrack_compute_interval) == 'dt') then
           ! from time step value
@@ -459,6 +556,12 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
           call mpas_pool_get_array(particle % haloDataPool, 'xParticle', xParticle)
           call mpas_pool_get_array(particle % haloDataPool, 'yParticle', yParticle)
           call mpas_pool_get_array(particle % haloDataPool, 'zParticle', zParticle)
+          call mpas_pool_get_array(particle % haloDataPool, 'lonParticle', lonParticle)
+          call mpas_pool_get_array(particle % haloDataPool, 'latParticle', latParticle)
+
+          ! compute particle lat/lon position
+          call convert_latlon_from_xyz(latParticle, lonParticle, xParticle, yParticle, zParticle)
+
           particlePosition(1) = xParticle
           particlePosition(2) = yParticle
           particlePosition(3) = zParticle
@@ -471,6 +574,22 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
           call mpas_pool_get_array(particle % haloDataPool, 'buoyancyParticle', buoyancyParticle)
           call mpas_pool_get_array(particle % haloDataPool, 'indexToParticleID', indexToParticleID)
           call mpas_pool_get_array(particle % haloDataPool, 'currentCell', iCell)
+          call mpas_pool_get_array(particle % haloDataPool, 'currentCellGlobalID', currentCellGlobalID)
+
+          ! sampling variables
+          call mpas_pool_get_array(particle % haloDataPool, 'particleTemperature', particleTemperature)
+          call mpas_pool_get_array(particle % haloDataPool, 'particleSalinity', particleSalinity)
+
+          if (config_use_ecosysTracers) then
+            call mpas_pool_get_array(particle % haloDataPool, 'particleDIC', particleDIC)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleALK', particleALK)
+            call mpas_pool_get_array(particle % haloDataPool, 'particlePO4', particlePO4)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleNO3', particleNO3)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleSiO3', particleSiO3)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleNH4', particleNH4)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleFe', particleFe)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleO2', particleO2)
+          end if
 
 !          call mpas_pool_get_array(particle % haloDataPool, 'lonVel', lonVel)
 !          call mpas_pool_get_array(particle % haloDataPool, 'latVel', latVel)
@@ -610,7 +729,8 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
 #ifdef MPAS_DEBUG
                 call mpas_timer_start("mpas_get_vertical_idLPT")
 #endif
-                iLevel = mpas_get_vertical_id(maxLevelCell(iCell), zSubStep, zMid(:,iCell))
+                ! must have inputs into mpas_get_vertical_id be positive to get right top/bottom
+                iLevel = mpas_get_vertical_id(maxLevelCell(iCell), -zSubStep, -zMid(:,iCell))
 #ifdef MPAS_DEBUG
                 call mpas_timer_stop("mpas_get_vertical_idLPT")
 #endif
@@ -679,7 +799,7 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
 #endif
           !}}}
 
-          !!!!!!!!!! PERFORM SAMPLING (VELOCITY, TEMP, SALINITY, ETC) !!!!!!!!!!
+          !!!!!!!!!! PERFORM SAMPLING (VELOCITY, TEMP, SALINITY, ETC) AND ERROR CHECKING!!!!!!!!!!
           ! this could be done just before the output anyway
 
           ! need iCell computed for final position
@@ -697,52 +817,259 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
           call mpas_timer_stop("get_validated_cell_idLPT")
 #endif
 
-          if(verticalTreatment == 4) then  !('buoyancySurface') !{{{
-            !! determine index level (don't need validated version because that will "fix" values which we may not want
-            !! however, if the particle's target buoyancy surface doesn't exist then we will need to
-            !! ensure that vertical location is valid, placing particles outside of range of zMid back inside
-            !! we don't validate this because we want the code to fail, at least initially, in the case that
-            !! the particle is in a cell that does not have the proper buoyancy target because this implies
-            !! that the buoyancy tracking mode has completely failed.
-            !! need to make sure it is validated for buoyancy particles
-            iLevelBuoyancy = mpas_get_vertical_id(maxLevelCell(iCell), buoyancyInterp, buoyancyTimeInterp(:,iCell))
-            ! interpolate the scalars now (assumes that scalar value is constant within a particular cell)
-            call interp_cell_scalars(iLevelBuoyancy, maxLevelCell(iCell), buoyancyInterp, buoyancyTimeInterp(:,iCell), &
-              zMid(:,iCell), zLevelParticle)
-            !}}}
-          else
-            ! make sure final zLevelParticle is ok so that it can't extent past zMid range
-#ifdef MPAS_DEBUG
-            call mpas_timer_start("mpas_get_vertical_idLPT")
-#endif
-            iLevel = mpas_get_vertical_id(maxLevelCell(iCell), zLevelParticle, zMid(:,iCell))
-#ifdef MPAS_DEBUG
-            call mpas_timer_stop("mpas_get_vertical_idLPT")
-#endif
-          end if
-
-          ! compute necessary information for autocorrelation !{{{
-          ! get the updated velocity
-          ! ensure that buoyancy is fixed for each run
-          buoyancyInterp = buoyancyParticle
-          ! we just need the last part of the velocity field interpolation because we are at the end of the timestep
-          timeInterpOrder = 2
-          timeCoeff(1) = 0.0_RKIND
-          timeCoeff(2) = 1.0_RKIND
-          ! return interpolated horizontal velocity "particleVelocity" and vertical velocity "particleVelocityVert"
-          ! noting we use the final positions
-#ifdef MPAS_DEBUG
-          call mpas_timer_start("velocity_time_interpolationLPT")
-#endif
-          call velocity_time_interpolation(particleVelocity, particleVelocityVert, &
-            diagnosticsPool, lagrPartTrackFieldsPool, &
-            timeInterpOrder, timeCoeff, iCell, iLevel, buoyancyInterp, maxLevelCell, &
-            verticalTreatment, indexLevel, nCellVertices, verticesOnCell, boundaryVertex, &
-            particlePosition, zLevelParticle, zMid, zTop, vertVelocityTop, xVertex, yVertex, zVertex, &
-            meshPool, areaBArray)
-#ifdef MPAS_DEBUG
-          call mpas_timer_stop("velocity_time_interpolationLPT")
-#endif
+          call get_vertical_level_info(iLevel, zLevelParticle, &
+            verticalTreatment, maxLevelCell, iCell, buoyancyInterp, buoyancyTimeInterp, &
+            zMid, indexLevel, ssh, bottomDepth, indexToParticleID)
+
+           !! sampling setup (if needed)
+           !if (&
+           !  config_AM_lagrPartTrack_sample_temperature .or. &
+           !  config_AM_lagrPartTrack_sample_salinity) then
+
+           !end if
+
+           ! store global cell ID
+           currentCellGlobalID = indexToCellID(iCell)
+
+           if ( associated(activeTracers) ) then
+             if (config_AM_lagrPartTrack_sample_temperature) then
+              !call mpas_log_write(' testing output for temperature')
+              if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                call interp_scalars(particleTemperature, activeTracers(index_temperature,:,:), &
+                  verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                  particlePosition(1),particlePosition(2),particlePosition(3), &
+                  xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                  yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                  zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                  xCell, yCell, zCell, &
+                  meshPool, &
+                  cellsOnVertex, &
+                  verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                  zLevelParticle, zMid, indexLevel)
+              else
+                call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                  activeTracers(index_temperature,:,iCell), particleTemperature)
+              end if
+              !call mpas_log_write('done testing output for temperature')
+             end if
+
+             if (config_AM_lagrPartTrack_sample_salinity) then
+              !call mpas_log_write(' testing output for salinity')
+              if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                call interp_scalars(particleSalinity, activeTracers(index_salinity,:,:), &
+                  verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                  particlePosition(1),particlePosition(2),particlePosition(3), &
+                  xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                  yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                  zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                  xCell, yCell, zCell, &
+                  meshPool, &
+                  cellsOnVertex, &
+                  verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                  zLevelParticle, zMid, indexLevel)
+              else
+                call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                  activeTracers(index_salinity,:,iCell), particleSalinity)
+              end if
+              !call mpas_log_write('done testing output for salinity')
+             end if
+           end if
+
+           if (config_use_ecosysTracers) then
+             call mpas_pool_get_array(tracersPool, 'ecosysTracers', ecosysTracers, timeLevel)
+             if ( associated(ecosysTracers) ) then
+               if (config_AM_lagrPartTrack_sample_DIC) then
+                 !call mpas_log_write('testing output for DIC')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleDIC, ecosysTracers(index_DIC,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_DIC,:,iCell), particleDIC)
+                 end if
+                 !call mpas_log_write('done testing output for DIC')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_ALK) then
+                 !call mpas_log_write('testing output for ALK')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleALK, ecosysTracers(index_ALK,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_ALK,:,iCell), particleALK)
+                 end if
+                 !call mpas_log_write('done testing output for ALK')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_PO4) then
+                 !call mpas_log_write('testing output for PO4')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particlePO4, ecosysTracers(index_PO4,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_PO4,:,iCell), particlePO4)
+                 end if
+                 !call mpas_log_write('done testing output for PO4')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_NO3) then
+                 !call mpas_log_write('testing output for NO3')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleNO3, ecosysTracers(index_NO3,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_NO3,:,iCell), particleNO3)
+                 end if
+                 !call mpas_log_write('done testing output for NO3')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_SiO3) then
+                 !call mpas_log_write('testing output for SiO3')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleSiO3, ecosysTracers(index_SiO3,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_SiO3,:,iCell), particleSiO3)
+                 end if
+                 !call mpas_log_write('done testing output for SiO3')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_NH4) then
+                 !call mpas_log_write('testing output for NH4')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleNH4, ecosysTracers(index_NH4,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_NH4,:,iCell), particleNH4)
+                 end if
+                 !call mpas_log_write('done testing output for NH4')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_Fe) then
+                 !call mpas_log_write('testing output for Fe')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleFe, ecosysTracers(index_Fe,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_Fe,:,iCell), particleFe)
+                 end if
+                 !call mpas_log_write('done testing output for Fe')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_O2) then
+                 !call mpas_log_write('testing output for O2')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleO2, ecosysTracers(index_O2,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_O2,:,iCell), particleO2)
+                 end if
+                 !call mpas_log_write('done testing output for O2')
+               end if
+             end if
+           end if
+
+!          ! compute necessary information for autocorrelation !{{{
+!          ! get the updated velocity
+!          ! ensure that buoyancy is fixed for each run
+!          buoyancyInterp = buoyancyParticle
+!          ! we just need the last part of the velocity field interpolation because we are at the end of the timestep
+!          timeInterpOrder = 2
+!          timeCoeff(1) = 0.0_RKIND
+!          timeCoeff(2) = 1.0_RKIND
+!          ! return interpolated horizontal velocity "particleVelocity" and vertical velocity "particleVelocityVert"
+!          ! noting we use the final positions
+!#ifdef MPAS_DEBUG
+!          call mpas_timer_start("velocity_time_interpolationLPT")
+!#endif
+!          call velocity_time_interpolation(particleVelocity, particleVelocityVert, &
+!            diagnosticsPool, lagrPartTrackFieldsPool, &
+!            timeInterpOrder, timeCoeff, iCell, iLevel, buoyancyInterp, maxLevelCell, &
+!            verticalTreatment, indexLevel, nCellVertices, verticesOnCell, boundaryVertex, &
+!            particlePosition, zLevelParticle, zMid, zTop, vertVelocityTop, xVertex, yVertex, zVertex, &
+!            meshPool, areaBArray)
+!#ifdef MPAS_DEBUG
+!          call mpas_timer_stop("velocity_time_interpolationLPT")
+!#endif
           ! convert horizontal velocity to lat/lon velocity
 
 !          ! store velocity for use in computing normalized autocorrelation offline
@@ -774,7 +1101,11 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
           yParticle = particlePosition(2)
           zParticle = particlePosition(3)
 
+          ! compute particle lat/lon position
+          call convert_latlon_from_xyz(latParticle, lonParticle, xParticle, yParticle, zParticle)
+
           !!!!!!!!!! PASS PARTICLES FROM PROCESSOR TO PROCESSOR !!!!!!!!!!
+          ! note, assumes that integration above didn't move particle past halo
           !{{{
           ! 1. determine if iCell is on halo (just set each particle's currentBlock to the correct currentBlock
           ! 2. determine owning block in halo, update particle's currentBlock
@@ -1648,6 +1979,7 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
       real (kind=RKIND), dimension(3) :: particlePosition, particleVelocity
       real (kind=RKIND), pointer :: xParticle, yParticle, zParticle, lonVel, latVel, buoyancyParticle !, &
 !        lonVel, latVel, sumU, sumV, sumUU, sumUV, sumVV
+      real (kind=RKIND), pointer :: lonParticle, latParticle
       real (kind=RKIND), dimension(3) :: xSubStep, diffSubStep, diffParticlePosition
       real (kind=RKIND), pointer :: zLevelParticle
       real (kind=RKIND), dimension(:,:), pointer :: zTop, vertVelocityTop, zMid, areaBArray
@@ -1661,13 +1993,14 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
       real(kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
 
       integer, dimension(:,:), pointer :: verticesOnCell, boundaryVertex
-      integer, dimension(:), pointer :: maxLevelCell
+      integer, dimension(:), pointer :: maxLevelCell, indexToCellID
 
       integer theVertex, iLevel, iLevelBuoyancy, aVertex, &
         nSteps, timeStep, subStep, subStepOrder, timeInterpOrder, aTimeLevel, nCellVertices, blockProc, arrayIndex
       integer, pointer :: nCells, nVertLevels
       integer, dimension(:), pointer :: nCellVerticesArray
       integer, dimension(:,:), pointer :: cellsOnCell
+      integer, dimension(:,:), pointer :: cellsOnVertex
       logical, dimension(:,:), pointer :: ioProcRecvList
       logical, dimension(:), pointer :: ioProcSendList
       logical, pointer :: onSphere
@@ -1685,10 +2018,60 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
       real (kind=RKIND) :: diffSubStepVert, diffParticlePositionVert, particleVelocityVert, verticalVelocityInterp
       real (kind=RKIND) :: buoyancyInterp
       real (kind=RKIND), pointer :: sphereRadius
-      integer, pointer :: verticalTreatment, indexLevel, filterNum, iCell
+      integer, pointer :: verticalTreatment, indexLevel, filterNum, iCell, currentCellGlobalID
+
+      type (mpas_pool_type), pointer :: tracersPool
+      real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers
+      integer, pointer :: index_temperature, index_salinity
+      logical, pointer :: config_AM_lagrPartTrack_sample_horizontal_interp
+      logical, pointer :: config_AM_lagrPartTrack_sample_temperature
+      real (kind=RKIND), pointer :: particleTemperature
+      logical, pointer :: config_AM_lagrPartTrack_sample_salinity
+      real (kind=RKIND), pointer :: particleSalinity
+
+      real (kind=RKIND), dimension(:,:,:), pointer :: ecosysTracers
+      logical, pointer :: config_use_ecosysTracers
+      integer, pointer :: index_DIC, &
+                          index_ALK, &
+                          index_PO4, &
+                          index_NO3, &
+                          index_SiO3, &
+                          index_NH4, &
+                          index_Fe, &
+                          index_O2
+      logical, pointer :: config_AM_lagrPartTrack_sample_DIC
+      real (kind=RKIND), pointer :: particleDIC
+      logical, pointer :: config_AM_lagrPartTrack_sample_ALK
+      real (kind=RKIND), pointer :: particleALK
+      logical, pointer :: config_AM_lagrPartTrack_sample_PO4
+      real (kind=RKIND), pointer :: particlePO4
+      logical, pointer :: config_AM_lagrPartTrack_sample_NO3
+      real (kind=RKIND), pointer :: particleNO3
+      logical, pointer :: config_AM_lagrPartTrack_sample_SiO3
+      real (kind=RKIND), pointer :: particleSiO3
+      logical, pointer :: config_AM_lagrPartTrack_sample_NH4
+      real (kind=RKIND), pointer :: particleNH4
+      logical, pointer :: config_AM_lagrPartTrack_sample_Fe
+      real (kind=RKIND), pointer :: particleFe
+      logical, pointer :: config_AM_lagrPartTrack_sample_O2
+      real (kind=RKIND), pointer :: particleO2
 
       err = 0
 
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_horizontal_interp', &
+                                             config_AM_lagrPartTrack_sample_horizontal_interp)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_temperature', config_AM_lagrPartTrack_sample_temperature)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_salinity', config_AM_lagrPartTrack_sample_salinity)
+      call mpas_pool_get_config(ocnConfigs, 'config_use_ecosysTracers', config_use_ecosysTracers)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_DIC', config_AM_lagrPartTrack_sample_DIC)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_ALK', config_AM_lagrPartTrack_sample_ALK)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_PO4', config_AM_lagrPartTrack_sample_PO4)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_NO3', config_AM_lagrPartTrack_sample_NO3)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_SiO3', config_AM_lagrPartTrack_sample_SiO3)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_NH4', config_AM_lagrPartTrack_sample_NH4)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_Fe', config_AM_lagrPartTrack_sample_Fe)
+      call mpas_pool_get_config(ocnConfigs, 'config_AM_lagrPartTrack_sample_O2', config_AM_lagrPartTrack_sample_O2)
+
       dminfo = domain % dminfo
 
       block => domain % blocklist
@@ -1701,6 +2084,7 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
         call mpas_pool_get_subpool(block % structs, 'lagrPartTrackFields', lagrPartTrackFieldsPool)
         call mpas_pool_get_subpool(block % structs, 'lagrPartTrackScratch', lagrPartTrackScratchPool)
         call mpas_pool_get_subpool(block % structs, 'lagrPartTrackCells', lagrPartTrackCellsPool)
+        call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
 
         ! particlelist should be stored in the structs pool probably (need to seriously rework the code!)
         particlelist => block % particlelist
@@ -1711,6 +2095,7 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
         call mpas_pool_get_array(meshPool, 'xVertex', xVertex)
         call mpas_pool_get_array(meshPool, 'yVertex', yVertex)
         call mpas_pool_get_array(meshPool, 'zVertex', zVertex)
+        call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID)
         call mpas_pool_get_array(lagrPartTrackCellsPool, 'wachspressAreaB', areaBArray)
         call mpas_pool_get_array(diagnosticsPool, 'zTop', zTop)
         call mpas_pool_get_array(diagnosticsPool, 'zMid', zMid)
@@ -1721,8 +2106,24 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
         call mpas_dmpar_exch_halo_field(normalVelocity)
         call mpas_pool_get_field(statePool, 'layerThickness', layerThickness, timeLevel=timeLevel)
 
+        call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, timeLevel=timeLevel)
+        call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature)
+        call mpas_pool_get_dimension(tracersPool, 'index_salinity', index_salinity)
+        if (config_use_ecosysTracers) then
+          call mpas_pool_get_array(tracersPool, 'ecosysTracers', ecosysTracers, timeLevel=timeLevel)
+          call mpas_pool_get_dimension(tracersPool, 'index_DIC', index_DIC)
+          call mpas_pool_get_dimension(tracersPool, 'index_ALK', index_ALK)
+          call mpas_pool_get_dimension(tracersPool, 'index_PO4', index_PO4)
+          call mpas_pool_get_dimension(tracersPool, 'index_NO3', index_NO3)
+          call mpas_pool_get_dimension(tracersPool, 'index_SiO3', index_SiO3)
+          call mpas_pool_get_dimension(tracersPool, 'index_NH4', index_NH4)
+          call mpas_pool_get_dimension(tracersPool, 'index_Fe', index_Fe)
+          call mpas_pool_get_dimension(tracersPool, 'index_O2', index_O2)
+        end if
+
         call mpas_pool_get_array(meshPool, 'verticesOnCell', verticesOnCell)
         call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
+        call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex)
         call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
         call mpas_pool_get_array(meshPool, 'boundaryVertex', boundaryVertex)
         call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
@@ -1745,13 +2146,19 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
           call mpas_pool_get_array(particle % haloDataPool, 'xParticle', xParticle)
           call mpas_pool_get_array(particle % haloDataPool, 'yParticle', yParticle)
           call mpas_pool_get_array(particle % haloDataPool, 'zParticle', zParticle)
+          call mpas_pool_get_array(particle % haloDataPool, 'lonParticle', lonParticle)
+          call mpas_pool_get_array(particle % haloDataPool, 'latParticle', latParticle)
           call mpas_pool_get_array(particle % haloDataPool, 'currentCell', iCell)
+          call mpas_pool_get_array(particle % haloDataPool, 'currentCellGlobalID', currentCellGlobalID)
           iCell = -1
 
           particlePosition(1) = xParticle
           particlePosition(2) = yParticle
           particlePosition(3) = zParticle
 
+          ! compute particle lat/lon position
+          call convert_latlon_from_xyz(latParticle, lonParticle, xParticle, yParticle, zParticle)
+
           call mpas_pool_get_array(particle % haloDataPool, 'zLevelParticle', zLevelParticle)
 
           call mpas_pool_get_array(particle % haloDataPool, 'verticalTreatment', verticalTreatment)
@@ -1760,6 +2167,20 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
           call mpas_pool_get_array(particle % haloDataPool, 'buoyancyParticle', buoyancyParticle)
           call mpas_pool_get_array(particle % haloDataPool, 'indexToParticleID', indexToParticleID)
 
+          call mpas_pool_get_array(particle % haloDataPool, 'particleTemperature', particleTemperature)
+          call mpas_pool_get_array(particle % haloDataPool, 'particleSalinity', particleSalinity)
+          if (config_use_ecosysTracers) then
+            call mpas_pool_get_array(particle % haloDataPool, 'particleDIC', particleDIC)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleALK', particleALK)
+            call mpas_pool_get_array(particle % haloDataPool, 'particlePO4', particlePO4)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleNO3', particleNO3)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleSiO3', particleSiO3)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleNH4', particleNH4)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleFe', particleFe)
+            call mpas_pool_get_array(particle % haloDataPool, 'particleO2', particleO2)
+          end if
+
+
 !          call mpas_pool_get_array(particle % haloDataPool, 'lonVel', lonVel)
 !          call mpas_pool_get_array(particle % haloDataPool, 'latVel', latVel)
 !          call mpas_pool_get_array(particle % haloDataPool, 'sumU', sumU)
@@ -1788,17 +2209,12 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
           call mpas_timer_stop("get_validated_cell_idLPT")
 #endif
 
-          if(verticalTreatment == 4) then  !('buoyancySurface') !{{{
-            iLevelBuoyancy = mpas_get_vertical_id(maxLevelCell(iCell), buoyancyParticle, buoyancyTimeInterp(:,iCell))
+          call get_vertical_level_info(iLevel, zLevelParticle, &
+            verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+            zMid, indexLevel)
 
-            ! interpolate the scalars now (assumes that scalar value is constant within a particular cell)
-            call interp_cell_scalars(iLevelBuoyancy, maxLevelCell(iCell), buoyancyParticle, buoyancyTimeInterp(:,iCell), &
-              zMid(:,iCell), zLevelParticle)
-          else
-            ! make sure final zLevelParticle is ok so that it can't extent past zMid range
-
-            iLevel = mpas_get_vertical_id(maxLevelCell(iCell), zLevelParticle, zMid(:,iCell))
-          end if
+          ! store global cell ID
+          currentCellGlobalID = indexToCellID(iCell)
 
           ! compute necessary information for autocorrelation !{{{
           ! get the updated velocity
@@ -1823,6 +2239,222 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
 !          ! store velocity for use in computing normalized autocorrelation offline
 !          call mpas_convert_xyz_velocity_to_latlon(lonVel, latVel, particlePosition, particleVelocity)
           !}}}
+          if ( associated(activeTracers) ) then
+            if (config_AM_lagrPartTrack_sample_temperature) then
+             !call mpas_log_write(' testing output for temperature')
+             if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+               call interp_scalars(particleTemperature, activeTracers(index_temperature,:,:), &
+                   verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                   particlePosition(1),particlePosition(2),particlePosition(3), &
+                   xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                   yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                   zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                   xCell, yCell, zCell, &
+                   meshPool, &
+                   cellsOnVertex, &
+                   verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                   zLevelParticle, zMid, indexLevel)
+             else
+               call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                 activeTracers(index_temperature,:,iCell), particleTemperature)
+             end if
+             !call mpas_log_write('done testing output for temperature')
+            end if
+
+            if (config_AM_lagrPartTrack_sample_salinity) then
+             !call mpas_log_write(' testing output for salinity')
+             if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+               call interp_scalars(particleSalinity, activeTracers(index_salinity,:,:), &
+                   verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                   particlePosition(1),particlePosition(2),particlePosition(3), &
+                   xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                   yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                   zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                   xCell, yCell, zCell, &
+                   meshPool, &
+                   cellsOnVertex, &
+                   verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                   zLevelParticle, zMid, indexLevel)
+             else
+               call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                 activeTracers(index_salinity,:,iCell), particleSalinity)
+             end if
+             !call mpas_log_write('done testing output for salinity')
+            end if
+          end if
+
+          if (config_use_ecosysTracers) then
+            call mpas_pool_get_array(tracersPool, 'ecosysTracers', ecosysTracers, timeLevel)
+            if ( associated(ecosysTracers) ) then
+              if (config_AM_lagrPartTrack_sample_DIC) then
+                !call mpas_log_write('testing output for DIC')
+                if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                  call interp_scalars(particleDIC, ecosysTracers(index_DIC,:,:), &
+                    verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                    particlePosition(1),particlePosition(2),particlePosition(3), &
+                    xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                    yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                    zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                    xCell, yCell, zCell, &
+                    meshPool, &
+                    cellsOnVertex, &
+                    verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                    zLevelParticle, zMid, indexLevel)
+                else
+                  call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                    ecosysTracers(index_DIC,:,iCell), particleDIC)
+                end if
+                !call mpas_log_write('done testing output for DIC')
+              end if
+
+              if (config_AM_lagrPartTrack_sample_ALK) then
+                !call mpas_log_write('testing output for ALK')
+                if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                  call interp_scalars(particleALK, ecosysTracers(index_ALK,:,:), &
+                    verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                    particlePosition(1),particlePosition(2),particlePosition(3), &
+                    xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                    yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                    zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                    xCell, yCell, zCell, &
+                    meshPool, &
+                    cellsOnVertex, &
+                    verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                    zLevelParticle, zMid, indexLevel)
+                else
+                  call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                    ecosysTracers(index_ALK,:,iCell), particleALK)
+                end if
+                !call mpas_log_write('done testing output for ALK')
+              end if
+
+               if (config_AM_lagrPartTrack_sample_PO4) then
+                 !call mpas_log_write('testing output for PO4')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particlePO4, ecosysTracers(index_PO4,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_PO4,:,iCell), particlePO4)
+                 end if
+                 !call mpas_log_write('done testing output for PO4')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_NO3) then
+                 !call mpas_log_write('testing output for NO3')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleNO3, ecosysTracers(index_NO3,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_NO3,:,iCell), particleNO3)
+                 end if
+                 !call mpas_log_write('done testing output for NO3')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_SiO3) then
+                 !call mpas_log_write('testing output for SiO3')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleSiO3, ecosysTracers(index_SiO3,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_SiO3,:,iCell), particleSiO3)
+                 end if
+                 !call mpas_log_write('done testing output for SiO3')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_NH4) then
+                 !call mpas_log_write('testing output for NH4')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleNH4, ecosysTracers(index_NH4,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_NH4,:,iCell), particleNH4)
+                 end if
+                 !call mpas_log_write('done testing output for NH4')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_Fe) then
+                 !call mpas_log_write('testing output for Fe')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleFe, ecosysTracers(index_Fe,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_Fe,:,iCell), particleFe)
+                 end if
+                 !call mpas_log_write('done testing output for Fe')
+               end if
+
+               if (config_AM_lagrPartTrack_sample_O2) then
+                 !call mpas_log_write('testing output for O2')
+                 if (config_AM_lagrPartTrack_sample_horizontal_interp) then
+                   call interp_scalars(particleO2, ecosysTracers(index_O2,:,:), &
+                     verticesOnCell(:,iCell), nCellVerticesArray(iCell), &
+                     particlePosition(1),particlePosition(2),particlePosition(3), &
+                     xVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     yVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     zVertex(verticesOnCell(1:nCellVerticesArray(iCell),iCell)), &
+                     xCell, yCell, zCell, &
+                     meshPool, &
+                     cellsOnVertex, &
+                     verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+                     zLevelParticle, zMid, indexLevel)
+                 else
+                   call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+                     ecosysTracers(index_O2,:,iCell), particleO2)
+                 end if
+                 !call mpas_log_write('done testing output for O2')
+               end if
+            end if
+          end if
 
           ! get next particle to process on the list
           particlelist => particlelist % next
@@ -1885,7 +2517,7 @@ subroutine particle_vertical_treatment(verticalTreatment, indexLevel, nCellVerti
      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
-     integer :: aVertex, theVertex
+     integer :: aVertex, theVertex, iLevelTmp
 
      vertVelocityInterp = 0.0_RKIND
 
@@ -1905,7 +2537,7 @@ subroutine particle_vertical_treatment(verticalTreatment, indexLevel, nCellVerti
 
          !! ensure that the boundary condition is enforced
          !call zero_boundary_nodal_values(nCellVertices, verticesOnCell, &
-         !  boundaryVertex, uvUcell, uvVcell, uvWcell)
+         !  boundaryVertex, uvCell)
 
          ! no vertical motion, just horizontal motion
          return
@@ -1920,7 +2552,7 @@ subroutine particle_vertical_treatment(verticalTreatment, indexLevel, nCellVerti
 
          !! ensure that there is zero nodal velocity on the boundary
          !call zero_boundary_nodal_values(nCellVertices, verticesOnCell, &
-         !  boundaryVertex, uvUcell, uvVcell, uvWcell)
+         !  boundaryVertex, uvCell)
 
          ! no vertical motion, just horizontal motion
          return
@@ -1928,18 +2560,20 @@ subroutine particle_vertical_treatment(verticalTreatment, indexLevel, nCellVerti
 
        case (3) verticalTreatmentCase !('passiveFloat') !{{{
 
-         ! interpolate the vertical velocity
+         ! interpolate the vertical velocity (can be different than zMid...,
+         ! which is assumed input)
+         iLevelTmp = mpas_get_vertical_id(nVertLevels, -zLoc, -zTop)
          vertVelocityInterp = interp_vert_velocity_to_zlevel( &
-           iLevel, zLoc, zTop, vertVelocityTop)
+           iLevelTmp, -zLoc, -zTop, vertVelocityTop)
 
          ! interpolate the horizontal velocity based on z-levels
          call interp_nodal_vectors(ncellvertices, verticesoncell, &
            ilevel, nVertLevels, zLoc, zmid, &
-           uvertexvelocity, vvertexvelocity, wvertexvelocity, uvCell)
+           uVertexVelocity, vVertexVelocity, wVertexVelocity, uvCell)
 
-         !! ensure that there is zero nodal velocity on the boundary
+         ! ensure that there is zero nodal velocity on the boundary
          !call zero_boundary_nodal_values(nCellVertices, verticesOnCell, &
-         !  boundaryVertex, uvUcell, uvVcell, uvWcell)
+         !  boundaryVertex, uvCell)
 
          return
          !}}}
@@ -1950,11 +2584,11 @@ subroutine particle_vertical_treatment(verticalTreatment, indexLevel, nCellVerti
          ! interpolate the horizontal velocity
          call interp_nodal_vectors(ncellvertices, verticesoncell, &
            ilevel, nVertLevels, phiInterp, phiMid, &
-           uvertexvelocity, vvertexvelocity, wvertexvelocity, uvCell)
+           uVertexVelocity, vVertexVelocity, wVertexVelocity, uvCell)
 
          ! ensure that there is zero nodal velocity on the boundary
          !call zero_boundary_nodal_values(nCellVertices, verticesOnCell, &
-         !  boundaryVertex, uvUcell, uvVcell, uvWcell)
+         !  boundaryVertex, uvCell)
 
          ! no vertical motion, just horizontal motion
          return
@@ -2047,9 +2681,8 @@ subroutine velocity_time_interpolation(particleVelocity, particleVelocityVert, &
     call mpas_pool_get_config(meshPool, 'x_period', x_period)
     call mpas_pool_get_config(meshPool, 'y_period', y_period)
 
-    ! allocations for particular cell !{{{
+    ! allocations for particular cell
     allocate(vertCoords(3,nCellVertices), uvCell(3,nCellVertices), areaB(nCellVertices))
-    !}}}
 
     ! get horizontal vertex locations (noting that there may be a
     ! bit of error because the particle could be at the top
@@ -2342,6 +2975,446 @@ subroutine particle_horizontal_movement(meshPool, pParticle, dpParticle) !{{{
 
    end subroutine particle_horizontal_movement!}}}
 
+!***********************************************************************
+!
+!  routine get_vertical_level_info
+!
+!> \brief   Returns iLevel and zLevelParticle information for vertical location
+!> \author  Phillip Wolfram
+!> \date    11/13/2018
+!> \details
+!>  This routine computes the iLevel and zLevelParticle associated with a particular
+!>  verticalTreatment.
+!
+!-----------------------------------------------------------------------
+   subroutine get_vertical_level_info(iLevel, zLevelParticle, &
+       verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+       zMid, indexLevel, ssh, bottomDepth, indexToParticleID)  !{{{
+
+     implicit none
+
+     !-----------------------------------------------------------------
+     ! input variables
+     !-----------------------------------------------------------------
+     integer, intent (in) :: verticalTreatment                        !< verticalTreatment used for particles
+     integer, dimension(:), pointer :: maxLevelCell                   !< maximum level for each horizontal cell
+     integer, intent(in) :: iCell                                     !< index for horizontal cell location
+     real (kind=RKIND), intent(in) :: buoyancyParticle                !< buoyancy value for particle
+     real (kind=RKIND), dimension(:,:), pointer :: buoyancyTimeInterp !< intepolated buoyancy value
+     real (kind=RKIND), dimension(:,:), pointer :: zMid               !< mid depth z value for each cell
+     integer, intent(in) :: indexLevel                                !< index level for particle vertical location
+     real (kind=RKIND), dimension(:), pointer, optional :: bottomDepth, ssh  !< error checking variables
+     integer, intent(in), pointer, optional :: indexToParticleID      !< index to particle ID
+
+
+     !-----------------------------------------------------------------
+     ! input/output variables
+     !-----------------------------------------------------------------
+
+     integer, intent(inout) :: iLevel                                 !< vertical level / cell of phiInterp
+     real (kind=RKIND), intent(inout) :: zLevelParticle                     !< vertical level for particle
+
+     !-----------------------------------------------------------------
+     ! local variables
+     !-----------------------------------------------------------------
+     REAL (kind=RKIND), parameter :: eps=1e-12
+
+
+     ! compute zLevelParticle and iLevel ! {{{
+     if(verticalTreatment == 4) then  !('buoyancySurface')
+       iLevel = mpas_get_vertical_id(maxLevelCell(iCell), buoyancyParticle, buoyancyTimeInterp(:,iCell))
+
+       ! get zLevelParticle by interpolating the scalars now (assumes that
+       ! scalar value is constant within a particular cell)
+       call interp_cell_scalars(iLevel, maxLevelCell(iCell), buoyancyParticle, buoyancyTimeInterp(:,iCell), &
+         zMid(:,iCell), zLevelParticle)
+     else if (verticalTreatment == 1) then !('indexLevel')
+       zLevelParticle = zMid(indexLevel, iCell)
+       iLevel = indexLevel
+     else if (verticalTreatment == 3 .or. verticalTreatment == 2) then !('passiveFloat' or 'fixedZLevel')
+       ! make sure final zLevelParticle is ok so that it can't extent past zMid range
+       ! note, zlevel can be modified near top / bottom of ocean to keep particle in bounds
+       zLevelParticle = min(max(zMid(maxLevelCell(iCell),iCell), zLevelParticle), zMid(1,iCell))
+       iLevel = mpas_get_vertical_id(maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell))
+       !call mpas_log_write('iLevel = $i', intArgs=(/iLevel/))
+     else
+       call mpas_log_write('Vertical treatment is not defined!', MPAS_LOG_CRIT)
+     end if
+     ! }}}
+
+     !verify particles are still in domain !{{{
+     if ( present(ssh) .and. present(bottomDepth) .and. present(indexToParticleID) ) then
+#ifdef MPAS_DEBUG
+       call mpas_timer_start("verifyParticleDepthInDomainLPT")
+#endif
+
+      ! debug !{{{
+!#ifdef MPAS_DEBUG
+!       ! check and verify that free floating particle hasn't gone out of boundary
+!       if (iLevel == -1) then
+!         call mpas_log_write(' Warning: Particle $i location of $r is higher than zMid of $r', &
+!           realArgs=(/ zLevelParticle, zMid(1, iCell) /), intArgs=(/iCell/))
+!       end if
+!       if (iLevel == 0) then
+!         call mpas_log_write(' Warning: Particle $i location of $r is deeper than zMid at bottom $r', &
+!           realArgs=(/ zLevelParticle, zMid(maxLevelCell(iCell),iCell) /), intArgs=(/iCell/))
+!       end if
+!#endif
+      !}}}
+
+       ! verify that it is less than SSH and greater than bottomDepth
+       if (zLevelParticle + bottomDepth(iCell) < eps) then
+         call mpas_log_write(' Error: $i ($i) Particle of type $i height $r is lower than bottom depth $r', &
+           realArgs=(/ zLevelParticle, -bottomDepth(iCell) /), intArgs=(/iCell, indexToParticleID, verticalTreatment/),&
+           messageType=MPAS_LOG_ERR)
+       end if
+       if (zLevelParticle - ssh(iCell) > eps) then
+         call mpas_log_write(' Error: Particle $i height $r is higher than ssh of $r', &
+           realArgs=(/ zLevelParticle, ssh(iCell) /), intArgs=(/iCell/), &
+           messageType=MPAS_LOG_ERR)
+       end if
+#ifdef MPAS_DEBUG
+       call mpas_timer_stop("verifyParticleDepthInDomainLPT")
+#endif
+     end if
+     !}}}
+
+   end subroutine get_vertical_level_info !}}}
+
+
+!***********************************************************************
+!
+!  routine interp_scalars
+!
+!> \brief   Interpolate scalars horizontal with linear interpolation and
+!>          based on a z level criteria (z-level, buoyancy, etc)
+!> \author  Phillip Wolfram
+!> \date    11/13/2018
+!> \details
+!>  This routine interpolates cell scalars linearly horizontally and veritcally
+!>  depending upon a criteria such as z-level, buoyancy, etc.
+!
+!-----------------------------------------------------------------------
+   subroutine interp_scalars(particleScalar, tracerarray, & !{{{
+       verticesOnTheCell, nVertices, &
+       xp,yp,zp, &
+       xv,yv,zv, &
+       xc,yc,zc, &
+       meshPool, &
+       cellsOnVertex, &
+       verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+       zLevelParticleIn, zMid, indexLevel)
+
+     implicit none
+
+     !-----------------------------------------------------------------
+     ! input variables
+     !-----------------------------------------------------------------
+     integer, intent (in) :: verticalTreatment                        !< verticalTreatment used for particles
+     integer, dimension(:), pointer :: maxLevelCell                   !< maximum level for each horizontal cell
+     real (kind=RKIND), intent(in) :: buoyancyParticle                !< buoyancy value for particle
+     real (kind=RKIND), dimension(:,:), pointer :: buoyancyTimeInterp !< intepolated buoyancy value
+     real (kind=RKIND), dimension(:,:), pointer :: zMid               !< mid depth z value for each cell
+     integer, intent(in) :: indexLevel                                !< index level for particle vertical location
+     integer, dimension(:), intent(in) :: verticesOnTheCell           !< vertex indices on cell
+     integer, dimension(:,:), pointer, intent(in) :: cellsOnVertex    !< cells ajoined to vertex
+     integer, intent(in) :: nVertices                                 !< number of vertices on the cell
+     real (kind=RKIND), intent(in) :: xp,yp,zp                        !< particle location
+     real (kind=RKIND), dimension(:), intent(in) :: xv,yv,zv          !< spatial location of cell vertices
+     real (kind=RKIND), dimension(:), intent(in) :: xc,yc,zc          !< spatial location of cell centers
+     type (mpas_pool_type), intent(in), pointer :: meshPool           !< meshPool pointer
+     real (kind=RKIND), dimension(:,:), intent(in) :: tracerarray     !< tracer array to be interpolated
+     real (kind=RKIND), pointer :: zLevelParticleIn                   !< vertical level z location
+
+     !-----------------------------------------------------------------
+     ! output variables
+     !-----------------------------------------------------------------
+     real (kind=RKIND), pointer :: particleScalar                     !< interpolated scalar value for the particle
+
+     !-----------------------------------------------------------------
+     ! local variables
+     !-----------------------------------------------------------------
+     real (kind=RKIND) :: zLevelParticle                     !< vertical level z location
+     integer :: iLevel                                       !< vertical level index
+     integer :: iCell                                        !< local cell index for horizontal cell location
+     integer :: triangleIndex                                !< index of triangle for interpolation
+     real (kind=RKIND), dimension(3) :: scalarTriangle       !< interpolated values on triangle
+     integer :: iVertex                                      !< vertex index variable
+     real (kind=RKIND), dimension(3, nVertices) :: xVertex   !< coordinates of vertex
+     logical, dimension(nVertices) :: interiorVertex         !< determines if vertex is on the interior
+     real (kind=RKIND), dimension(3) :: xPoint               !< coordinates of transformed point for nearest neighbor calc
+     real (kind=RKIND), dimension(3) :: lambda               !< barycentric weights for triangle
+
+     real (kind=RKIND) :: pointRadius   ! magnitude of a point radius
+     logical, pointer :: on_a_sphere    ! flag designating if we are on a sphere
+     logical, pointer :: is_periodic    ! flag designating if periodicity is important
+     real(kind=RKIND), pointer :: x_period, y_period
+     integer, parameter :: vertexDegree = 3                  !< max number of cells on Vertex
+     real(kind=RKIND) :: dx, dy, dz, r2, r2Min
+     real (kind=RKIND), dimension(3,3) :: vertexTriangle
+     integer :: aVertex, vertexID, aCell, aPoint, localCell
+     integer, pointer :: nCells
+     real(kind=RKIND), parameter :: eps = 1.0e-12_RKIND
+     integer :: ierr = 0
+     logical :: boundaryInterp
+
+     ! assume that we cannot do spatial interpolation
+     vertexID = -1
+
+     !-----------------------------------------------------------------
+     ! find the vertexID of the polygon containing the point aPoint
+     !-----------------------------------------------------------------
+
+     ! {{{
+     ! for each vertex assess whether it is an interior vertex
+     ! interior vertex until proven otherwise
+     interiorVertex(:) = .true.
+     do aPoint = 1, nVertices
+       do aVertex = 1, vertexDegree
+         !call mpas_log_write('CellsOnVertex($i, $i) = $i', intArgs = (/aPoint, aVertex, cellsOnVertex(aPoint, aVertex)/))
+         if (cellsOnVertex(aVertex, verticesOnTheCell(aPoint)) == 0) then
+           interiorVertex(aPoint) = .false.
+         end if
+       end do
+     end do
+     !! hard code as exterior vertex for error checking / debugging
+     !interiorVertex(:) = .false.
+
+     ! prepare vertex array xVertex needed to find closest vertex to point xPoint ! {{{
+     call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere)
+     call mpas_pool_get_config(meshPool, 'is_periodic', is_periodic)
+     call mpas_pool_get_config(meshPool, 'x_period', x_period)
+     call mpas_pool_get_config(meshPool, 'y_period', y_period)
+     call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
+
+     ! check to see if point is in vertex polygon
+     xPoint(1) = xp
+     xPoint(2) = yp
+     xPoint(3) = zp
+
+     ! for each vertex assess whether it is valid to use
+     do aVertex = 1, nVertices
+       if (interiorVertex(aVertex)) then
+         do aCell = 1, vertexDegree
+
+           localCell = cellsOnVertex(aCell, verticesOnTheCell(aVertex))
+
+           ! get coordinates of vertex for triangle for barycentric interpolation prepared
+           if ( (.not. on_a_sphere) .and. is_periodic) then
+             vertexTriangle(aCell, 1) = mpas_fix_periodicity(xc(localCell), xp, x_period)
+             vertexTriangle(aCell, 2) = mpas_fix_periodicity(yc(localCell), yp, y_period)
+             vertexTriangle(aCell, 3) = zc(localCell)
+           else
+             vertexTriangle(aCell, 1) = xc(localCell)
+             vertexTriangle(aCell, 2) = yc(localCell)
+             vertexTriangle(aCell, 3) = zc(localCell)
+           end if
+         end do
+
+         ! check to make sure the point xPoint is in the vertexTriangle
+         if ( mpas_point_in_polygon(xPoint, vertexTriangle, on_a_sphere) ) then
+           vertexID = aVertex
+         end if
+         ! otherwise default to "boundary" or edge case
+       end if
+     end do
+
+     !! normalize 3D points on sphere
+     !if ( on_a_sphere .and. .not. is_periodic) then
+     !  pointRadius = sqrt(xp*xp + yp*yp + zp*zp)
+     !  xPoint = xPoint/pointRadius
+     !  vertexTriangle = vertexTriangle/pointRadius
+     !end if
+     !}}}
+
+     ! debugging output !{{{
+!#ifdef MPAS_DEBUG
+!     ! check to make sure nearest vertex found makes sense
+!     call mpas_log_write('vertexID = $i, xPoint = $r $r $r', intArgs = (/ vertexID /), &
+!       realArgs = (/ xPoint(1), xPoint(2), xPoint(3) /))
+!     do aVertex = 1, nVertices
+!       call mpas_log_write('$i xVertex = $r $r $r$', intArgs = (/ aVertex /), &
+!         realArgs = (/ &
+!         xVertex(1, aVertex), &
+!         xVertex(2, aVertex), &
+!         xVertex(3, aVertex)  &
+!         /))
+!     end do
+!#endif
+     !}}}
+     !}}}
+     ! }}}
+
+     !-----------------------------------------------------------------
+     ! for each vertex, determine how to interpolate from neighbooring
+     ! cell values if an edge-case; if a boundary just do vertical
+     ! interpolation
+     !-----------------------------------------------------------------
+
+     ! {{{
+     ! cellsOnVertex assumes an "empty" neighbor (edge case) is 0
+     do aCell = 1, vertexDegree
+
+       boundaryInterp = .false.
+
+       ! only handles the case where the vertex is an exterior vertex, otherwise treat as cell-based value
+       ! (no horizontal interpolation)
+       if (vertexID .eq. -1) then
+         boundaryInterp = .true.
+       else if (cellsOnVertex(aCell, verticesOnTheCell(vertexID)) == 0 .or. &
+                cellsOnVertex(aCell, verticesOnTheCell(vertexID)) == nCells+1 .or. &
+                (cellsOnVertex(1, verticesOnTheCell(vertexID)) == cellsOnVertex(2, verticesOnTheCell(vertexID))) .or. &
+                (cellsOnVertex(2, verticesOnTheCell(vertexID)) == cellsOnVertex(3, verticesOnTheCell(vertexID))) .or. &
+                (cellsOnVertex(3, verticesOnTheCell(vertexID)) == cellsOnVertex(1, verticesOnTheCell(vertexID)))) then
+         boundaryInterp = .true.
+       end if
+
+       if (boundaryInterp) then
+
+         !call mpas_log_write('Using column interpolation for $i because it is on an edge', intArgs = (/iCell/))
+
+         ! note that zLevelParticle corresponding to adjacent cells will be a
+         ! little different on boundary, hence the need to use this temporary variable
+         ! so interpolation near  boundaries interpolates on the surface
+         ! of the boundary (e.g., to make sure interpolated values are within interior
+         ! away from dz/2 of column boundary)
+         zLevelParticle = zLevelParticleIn
+
+         ! in the case of any edge, just use cell value and return it
+         call get_vertical_level_info(iLevel, zLevelParticle, &
+           verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+           zMid, indexLevel)
+
+         call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
+           tracerarray(:,iCell), particleScalar)
+
+         return
+       end if
+     end do
+     ! }}}
+
+     !-----------------------------------------------------------------
+     ! perform barycentric interpolation horizontal from
+     ! vertical interpolants
+     !-----------------------------------------------------------------
+
+     ! {{{
+     ! all cellsOnVertex(verticesOnTheCell(vertexID), 1:3) are real cells, for each of triangle vertices
+     ! get cell values for the vertical interpolant for each point of triangle for interpolation
+     vertexTriangle = 0.0_RKIND
+     do aCell = 1, vertexDegree
+
+       localCell = cellsOnVertex(aCell, verticesOnTheCell(vertexID))
+
+       ! perform vertical interpolation for each cell
+       ! return iLevel, zLevelParticle for vertical interpolation of the scalar
+       ! note that zLevelParticle corresponding to adjacent cells will be a
+       ! little different on boundary, hence the need to use this temporary variable
+       ! so interpolation near  boundaries interpolates on the surface
+       ! of the boundary (e.g., to make sure interpolated values are within interior
+       ! away from dz/2 of column boundary)
+       zLevelParticle = zLevelParticleIn
+       call get_vertical_level_info(iLevel, zLevelParticle, &
+         verticalTreatment, maxLevelCell, localCell, &
+         buoyancyParticle, buoyancyTimeInterp, &
+         zMid, indexLevel)
+
+       ! perform vertical interpolation for the scalar
+       call interp_cell_scalars(iLevel, maxLevelCell(localCell), &
+         -zLevelParticle, -zMid(:,localCell), &
+         tracerarray(:,localCell), scalarTriangle(aCell))
+
+       ! get coordinates of vertex for triangle for barycentric interpolation prepared
+       if ( (.not. on_a_sphere) .and. is_periodic) then
+         vertexTriangle(1, aCell) = mpas_fix_periodicity(xc(localCell), xp, x_period)
+         vertexTriangle(2, aCell) = mpas_fix_periodicity(yc(localCell), yp, y_period)
+         vertexTriangle(3, aCell) = zc(localCell)
+       else
+         vertexTriangle(1, aCell) = xc(localCell)
+         vertexTriangle(2, aCell) = yc(localCell)
+         vertexTriangle(3, aCell) = zc(localCell)
+       end if
+     end do
+
+     ! perform barycentric interpolation in the horizontal
+     xPoint(1) = xp
+     xPoint(2) = yp
+     xPoint(3) = zp
+
+     !! normalize 3D points on sphere
+     !if ( on_a_sphere .and. .not. is_periodic) then
+     !  pointRadius = sqrt(xp*xp + yp*yp + zp*zp)
+     !  xPoint = xPoint/pointRadius
+     !  vertexTriangle = vertexTriangle/pointRadius
+     !end if
+
+     !!! this wachspress routine is designed for
+     !!! polygons with more than 3 sides (making it have higher cost too) ...
+     lambda = mpas_wachspress_coordinates(vertexDegree, vertexTriangle, xPoint, meshPool)
+
+     ! more robust (and less expensive) barycentric interpolation approach,
+     ! assuming a PURE interpolation (which cannot be guaranteed here)
+     !call mpas_calculate_barycentric_weights(xPoint, &
+     !  vertexTriangle(:,1), vertexTriangle(:,2), vertexTriangle(:,3), &
+     !  meshPool, &
+     !  lambda(1), lambda(2), lambda(3), ierr)
+
+! debugging output !{{{
+#ifdef MPAS_DEBUG
+     ! ensure interpolation is obtained
+     if (minval(lambda) < -eps .or. maxval(lambda) > (1.0_RKIND + eps) .or. &
+         abs(1.0_RKIND - sum(lambda)) > eps ) then
+       call mpas_log_write('Not a valid barycentric interpolation: ' // &
+         'lambda = $r $r $r vertexID = $i cellsOnVertex(:, verticesOnTheCell(vertexID)) = $i $i $i', MPAS_LOG_ERR, &
+         realArgs=(/lambda(1), lambda(2), lambda(3)/), &
+         intArgs=(/vertexID, &
+         cellsOnVertex(1, verticesOnTheCell(vertexID)), &
+         cellsOnVertex(2, verticesOnTheCell(vertexID)), &
+         cellsOnVertex(3, verticesOnTheCell(vertexID))/))
+     end if
+#endif
+!}}}
+
+     particleScalar = scalarTriangle(1) * lambda(1) + &
+                      scalarTriangle(2) * lambda(2) + &
+                      scalarTriangle(3) * lambda(3)
+     !}}}
+
+     ! debugging output !{{{
+!#ifdef MPAS_DEBUG
+!     ! check to make sure nearest vertex found makes sense
+!     call mpas_log_write('vertexID = $i, xPoint = $r $r $r', intArgs = (/ vertexID /), &
+!       realArgs = (/ xPoint(1), xPoint(2), xPoint(3) /))
+!     do aVertex = 1, nVertices
+!       call mpas_log_write('$i xVertex = $r $r $r$', intArgs = (/ aVertex /), &
+!         realArgs = (/ &
+!         xVertex(1, aVertex), &
+!         xVertex(2, aVertex), &
+!         xVertex(3, aVertex) /))
+!     end do
+!     do aCell = 1, vertexDegree
+!       localCell = cellsOnVertex(aCell, verticesOnTheCell(vertexID))
+!       call mpas_log_write(' cell id $i $r $r $r', intArgs=(/ localCell/), &
+!         realArgs = (/&
+!         vertexTriangle(1, aCell), &
+!         vertexTriangle(2, aCell), &
+!         vertexTriangle(3, aCell)/))
+!     end do
+!     call mpas_log_write('scalar Triangle = $r $r $r', realArgs = &
+!       (/scalarTriangle(1), &
+!         scalarTriangle(2), &
+!         scalarTriangle(3)/))
+!     call mpas_log_write('lambda = $r $r $r', realArgs = &
+!       (/lambda(1), &
+!         lambda(2), &
+!         lambda(3)/))
+!     call mpas_log_write('ierr = $i', intArgs=(/ierr/))
+!#endif
+     !}}}
+
+   end subroutine interp_scalars !}}}
+
+
 !***********************************************************************
 !
 !  routine interp_cell_scalars
@@ -2365,7 +3438,7 @@ subroutine interp_cell_scalars(iLevel, nVertLevels, zInterp, zVals, & !{{{
      real (kind=RKIND), dimension(:), intent(in) :: zVals             !< scalar values (x) on cell for interpolant
      integer, intent(in) :: iLevel                                    !< vertical level / cell of phiInterp
      integer, intent(in) :: nVertLevels                               !< number of vertical levels
-     real (kind=RKIND), intent(in) :: zInterp                       !< location to interpolate
+     real (kind=RKIND), intent(in) :: zInterp                         !< location to interpolate
      real (kind=RKIND), dimension(:), intent(in) :: phiVals           !< values at elevation of cell middle (where vertex
                                                                       !< velocities are defined)
 
@@ -2581,7 +3654,7 @@ subroutine zero_boundary_nodal_values(nCellVertices, verticesOnCell, & !{{{
      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
-     real (kind=RKIND), dimension(:,:), intent(out) :: uvCell  !< components of vertex velocity (vertically selected)
+     real (kind=RKIND), dimension(:,:), intent(inout) :: uvCell  !< components of vertex velocity (vertically selected)
 
      !-----------------------------------------------------------------
      ! local variables
@@ -2642,7 +3715,7 @@ subroutine get_bounding_indices(iLow, iHigh, phiInterp, phiVals, iLevel, nVertLe
      ! check to make sure points are in range
      ! optimization point: smarter algorithm won't have to call this ever
      if(.not.((phiInterp <= phiVals(iHigh)) .and. (phiInterp >= phiVals(iLow)))) then
-       !call mpas_log_write( 'fast interpolation failed, trying general, brute force search for interpolation bounds')
+       !call mpas_log_write( 'Warning: fast interpolation failed, trying general, brute force search for interpolation bounds')
        !print *,  'iLow = ', iLow, ' iHigh = ', iHigh
        !print *,  'phiInterp = ', phiInterp , ' phiLow =',  phiVals(iLow), ' phiHigh = ', phiVals(iHigh)
        call get_bounding_indices_brute_force(nVertLevels, phiInterp, phiVals, iLow, iHigh)
@@ -2694,6 +3767,8 @@ subroutine get_bounding_indices_brute_force(nVertLevels, phiInterp, phiVals, iLo
          iLow  = aLevel+1
          iHigh = aLevel
          exit
+       else
+         call mpas_log_write('Value for interpolation was not found!', MPAS_LOG_CRIT)
        end if
      end do
 
@@ -2702,15 +3777,14 @@ subroutine get_bounding_indices_brute_force(nVertLevels, phiInterp, phiVals, iLo
        ! we are ok
        LIGHT_DEBUG_ALL_WRITE('brute force interpolation successful')
        LIGHT_DEBUG_ALL_WRITE('iLow = ' COMMA iLow COMMA ' iHigh = ' COMMA iHigh)
-       write(message, *) 'phiInterp = ', phiInterp, ' phiLow =', phiVals(iLow), ' phiHigh = ', phiVals(iHigh)
+       LIGHT_DEBUG_ALL_WRITE('phiInterp = ' COMMA phiInterp COMMA ' phiLow =' COMMA phiVals(iLow) COMMA ' phiHigh = ' COMMA phiVals(iHigh))
        LIGHT_DEBUG_ALL_WRITE(message)
      else
-       print *,  'brute force interpolation failed with phiInterp = ', phiInterp, ' phiLow = ', phiVals(iLow), &
-                        ' phiHigh = ', phiVals(iHigh)
+       LIGHT_DEBUG_ALL_WRITE('brute force interpolation failed with phiInterp = ' COMMA phiInterp COMMA ' phiLow = ' COMMA phiVals(iLow) COMMA ' phiHigh = ' COMMA phiVals(iHigh))
        LIGHT_DEBUG_ALL_WRITE(' phiVals = ' COMMA phiVals(1:nVertLevels))
      end if
 
-     call mpas_log_write( 'Warning!: brute force interpolation used, boundary condition may be wrong!')
+     LIGHT_DEBUG_WRITE('Warning!: brute force interpolation used, boundary condition may be wrong!')
 #endif
 
    end subroutine get_bounding_indices_brute_force!}}}
@@ -2763,7 +3837,9 @@ real (kind=RKIND) function interp_vert_velocity_to_zlevel( & !{{{
          alpha = (zSubStep - zTop(iLevel))/(zTop(iLevel+1)- zTop(iLevel))
 #ifdef MPAS_DEBUG
        else
-         call mpas_log_write( 'Error with vertical velocity interpolation!')
+         call mpas_log_write('Error with vertical velocity interpolation! ' // &
+           'iLevel = $i, zTop(iLevel) = $r, zTop(iLevel+1)=$r zTop(iLevel+2) = $r zSubstep=$r', &
+           realArgs =(/zTop(iLevel), zTop(iLevel+1), zTop(iLevel+2), zSubStep/), intArgs=(/iLevel/), messageType=MPAS_LOG_CRIT)
 #endif
        end if
        interp_vert_velocity_to_zlevel = alpha * vertVelocityTop(iLevel)+ (1.0_RKIND - alpha) * vertVelocityTop(iLevel+1)
@@ -2801,6 +3877,48 @@ subroutine time_interp_field(basePool, timeInterpOrder, timeCoeff, field, fieldn
      end do
 
    end subroutine time_interp_field !}}}
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  convert_latlon_from_xyz
+!
+!> \brief
+!> \author Adrian K. Turner and Phillip J. Wolfram, LANL
+!> \date 2013-2014, 10/08/2018
+!> \details
+!>  Adapted from seaice_latlon_from_xyz, with
+!>          1. assume xyz points define r
+!>          2. shift to be in range of 0 to 2*pi as in mesh tools
+!>
+!> Note, should be moved into framework.
+!
+!-----------------------------------------------------------------------
+  subroutine convert_latlon_from_xyz(lat, lon, x, y, z) !{{{
+    implicit none
+
+    ! given xyz coordinates determine the latitude and longitude
+    real(kind=RKIND), intent(out) :: &
+         lat, & !< Output: latitude of point
+         lon    !< Output: longitude of point
+
+    real(kind=RKIND), intent(in) :: &
+         x,   & !< Input: x position of point
+         y,   & !< Input: y position of point
+         z      !< Input: z position of point
+
+    real(kind=RKIND) :: r      !< Computed: radius of earth
+
+    r = sqrt(x*x + y*y + z*z)
+    lon = atan2(y, x)
+    lat = asin(z/r)
+
+    ! ensure range in 0, 2*pi
+    if (lon < 0.0_RKIND) then
+      lon = 2*pii + lon
+    end if
+
+  end subroutine convert_latlon_from_xyz !}}}
+
 !}}}
 
 end module ocn_lagrangian_particle_tracking
diff --git a/src/core_ocean/analysis_members/mpas_ocn_particle_list.F b/src/core_ocean/analysis_members/mpas_ocn_particle_list.F
index fa63c4e8ff..ac694adeba 100644
--- a/src/core_ocean/analysis_members/mpas_ocn_particle_list.F
+++ b/src/core_ocean/analysis_members/mpas_ocn_particle_list.F
@@ -18,8 +18,8 @@
 !-----------------------------------------------------------------------
 #define COMMA ,
 #ifdef MPAS_DEBUG
-#define LIGHT_DEBUG_WRITE(Z)   print *, Z
-#define LIGHT_DEBUG_ALL_WRITE(Z)  print *, Z
+#define LIGHT_DEBUG_WRITE(Z) ! print *, Z
+#define LIGHT_DEBUG_ALL_WRITE(Z) ! print *, Z
 #else
 #define LIGHT_DEBUG_WRITE(Z)  ! print *, Z
 #define LIGHT_DEBUG_ALL_WRITE(Z) ! print *, Z
@@ -3478,7 +3478,7 @@ subroutine communicate_particle_halo_data(domain, procNeighs, nPartSend, nPartRe
 #ifdef _MPI
      !call MPI_Barrier(domain % dminfo % comm)
 #endif
-     
+
      numReals = 0
      numInts = 0
      call mpas_pool_get_subpool(domain % blocklist % structs, 'lagrPartTrackHalo', lagrPartTrackPool)
@@ -3486,9 +3486,9 @@ subroutine communicate_particle_halo_data(domain, procNeighs, nPartSend, nPartRe
      do while(mpas_pool_get_next_member(lagrPartTrackPool, dimItr))
        if (dimItr % memberType == MPAS_POOL_FIELD) then
          !print *,  'transfering ', trim(dimItr % memberName)
-         if (dimItr % dataType == MPAS_POOL_REAL) then 
+         if (dimItr % dataType == MPAS_POOL_REAL) then
            numReals = numReals + 1
-         elseif (dimItr % dataType == MPAS_POOL_INTEGER) then 
+         elseif (dimItr % dataType == MPAS_POOL_INTEGER) then
            numInts = numInts + 1
          end if
        else
@@ -3496,13 +3496,13 @@ subroutine communicate_particle_halo_data(domain, procNeighs, nPartSend, nPartRe
          write(message, *) "Different type expected in registry for key ", trim(dimItr % memberName), &
                            " of type", dimItr % memberType, "in halo data for communication-- don't know what to do!"
          LIGHT_DEBUG_ALL_WRITE(message)
-       end if 
+       end if
      end do
 
 ! NOTE, code assumes that memberName is consistent between processors in terms of pool interation
 !#ifdef mpas_debug
 !     ! verify that membername order is the same between processors
-!     call mpas_log_write('printing particle pool contents (verify that membername consistent between processors):') 
+!     call mpas_log_write('printing particle pool contents (verify that membername consistent between processors):')
 !     call mpas_pool_print_summary(lagrPartTrackPool, mpas_pool_field, .true.)
 !#endif
 
@@ -3524,7 +3524,7 @@ subroutine communicate_particle_halo_data(domain, procNeighs, nPartSend, nPartRe
      allocate(array1DRealSend(numProcs), array1DRealRecv(numProcs))
      ! recv
      j = 1
-     do i=1,numProcs ! {{{ 
+     do i=1,numProcs ! {{{
        if(nPartRecv(i) > 0) then
          allocate(array1DRealRecv(i)%val(nPartRecv(i)*numReals))
 #ifdef _MPI
@@ -3607,12 +3607,12 @@ subroutine communicate_particle_halo_data(domain, procNeighs, nPartSend, nPartRe
 
      deallocate(array1DRealSend, array1DRealRecv)
      !!!! end REALS !!!! ! }}}
-     
+
      !!!!! begin INTS !!!! ! {{{
      allocate(array1DIntSend(numProcs), array1DIntRecv(numProcs))
      ! recv
      j = 1
-     do i=1,numProcs ! {{{ 
+     do i=1,numProcs ! {{{
        if(nPartRecv(i) > 0) then
          allocate(array1DIntRecv(i)%val(nPartRecv(i)*numInts))
 #ifdef _MPI
@@ -3695,7 +3695,7 @@ subroutine communicate_particle_halo_data(domain, procNeighs, nPartSend, nPartRe
 
      deallocate(array1DIntSend, array1DIntRecv)
      !!!! end INTS !!!! ! }}}
-     
+
      deallocate(recvRequestID, sendRequestID)
      LIGHT_DEBUG_WRITE('Finished primary MPI communication for halo')
 
diff --git a/testing_and_setup/compass/ocean/soma/32km/default/config_analysis.xml b/testing_and_setup/compass/ocean/soma/32km/default/config_analysis.xml
new file mode 100644
index 0000000000..b0cffc078e
--- /dev/null
+++ b/testing_and_setup/compass/ocean/soma/32km/default/config_analysis.xml
@@ -0,0 +1,11 @@
+<?xml version="1.0"?>
+<config case="analysis">
+	<add_link source_path="mpas_model" source="testing_and_setup/compass/ocean/soma/32km/default/analysis/check_particle_sampling.py" dest="check_particle_sampling.py"/>
+	<add_link source="../forward/output/output.0001-01-01_00.00.00.nc" dest="output.nc"/>
+	<add_link source="../forward/analysis_members/lagrPartTrack.0001-01-01_00.00.00.nc" dest="particleoutput.nc"/>
+
+	<run_script name="run.py">
+		<step executable="./check_particle_sampling.py">
+		</step>
+	</run_script>
+</config>
diff --git a/testing_and_setup/compass/ocean/soma/32km/default/config_driver.xml b/testing_and_setup/compass/ocean/soma/32km/default/config_driver.xml
index 046265dd2b..6a0588b7c1 100644
--- a/testing_and_setup/compass/ocean/soma/32km/default/config_driver.xml
+++ b/testing_and_setup/compass/ocean/soma/32km/default/config_driver.xml
@@ -8,12 +8,17 @@
 	<case name="forward">
 		<step executable="./run.py" quiet="true" pre_message=" * Running forward" post_message="     Complete"/>
 	</case>
+	<case name="analysis">
+		<step executable="./run.py" quiet="true" pre_message=" * Running analysis" post_message="     Complete"/>
+	</case>
 	<validation>
 		<compare_fields file1="forward/analysis_members/lagrPartTrack.0001-01-01_00.00.00.nc">
 			<template file="LIGHT_comparison.xml" path_base="script_core_dir" path="templates/validations"/>
 			<field name="currentCell" l1_norm="0.0" l2_norm="0.0" linf_norm="0.0"/>
 			<field name="transfered" l1_norm="0.0" l2_norm="0.0" linf_norm="0.0"/>
 			<field name="numTimesReset" l1_norm="0.0" l2_norm="0.0" linf_norm="0.0"/>
+			<field name="particleTemperature" l1_norm="0.0" l2_norm="0.0" linf_norm="0.0"/>
+			<field name="particleSalinity" l1_norm="0.0" l2_norm="0.0" linf_norm="0.0"/>
 		</compare_fields>
 		<compare_fields file1="forward/output/output.0001-01-01_00.00.00.nc">
 			<field name="temperature" l1_norm="0.0" l2_norm="0.0" linf_norm="0.0"/>
diff --git a/testing_and_setup/compass/ocean/soma/32km/default/config_forward.xml b/testing_and_setup/compass/ocean/soma/32km/default/config_forward.xml
index 4c3db84f45..2e5b3a23af 100644
--- a/testing_and_setup/compass/ocean/soma/32km/default/config_forward.xml
+++ b/testing_and_setup/compass/ocean/soma/32km/default/config_forward.xml
@@ -35,8 +35,14 @@
 		<template file="lagrangian_particle_tracking.xml" path_base="script_core_dir" path="templates/analysis_members"/>
 		<template file="soma_template.xml" path_base="script_configuration_dir"/>
 		<stream name="output">
-                        <attribute name="output_interval">0000_00:01:30</attribute>
-                </stream>
+			<attribute name="output_interval">0000_00:01:30</attribute>
+		</stream>
+		<stream name="lagrPartTrackOutput">
+			<add_contents>
+				<member name="particleTemperature" type="var"/>
+				<member name="particleSalinity" type="var"/>
+			</add_contents>
+		</stream>
 	</streams>
 
 	<run_script name="run.py">
diff --git a/testing_and_setup/compass/ocean/soma/analysis/check_particle_sampling.py b/testing_and_setup/compass/ocean/soma/analysis/check_particle_sampling.py
new file mode 100755
index 0000000000..18f927e6fd
--- /dev/null
+++ b/testing_and_setup/compass/ocean/soma/analysis/check_particle_sampling.py
@@ -0,0 +1,58 @@
+#!/usr/bin/env python
+"""
+    Verifies that temperature and salinity profiles for SOMA case from particle
+    interpolations are consistent with initialized profiles.
+
+    Phillip J. Wolfram
+    04/22/2019
+"""
+import xarray as xr
+import numpy as np
+import matplotlib.pyplot as plt
+
+ds = xr.open_mfdataset('particleoutput.nc')
+
+# obtain particle data
+zLevelParticle = ds.zLevelParticle.values
+particleTemperature = ds.particleTemperature.values
+particleSalinity = ds.particleSalinity.values
+
+# compute analytical formulations
+
+# hard coded factors from config files
+gamma = 1.0/1250.0 # psu /m
+s0 = 34 # psu
+alpha=0.25 # kg / m^3
+beta=0.05
+drho = 5 # kg/m^3
+zt = 300 # m
+t0 = 15 # deg C
+hshelf = 100 # m
+h0 = 2400 # m
+z = np.linspace(0, zLevelParticle.min(), 30)
+
+# compute profiles
+salinity = s0 - gamma*z
+temperature = ((1-beta)*drho*np.tanh(z/zt) + beta*drho*(z/(h0+hshelf)))/alpha + t0
+
+# temperature comparison
+plt.figure()
+plt.plot(particleTemperature.ravel(), zLevelParticle.ravel(), '.',
+         label='Particle value')
+plt.plot(temperature, z, '-', label='Initialized scalar value')
+plt.ylabel('Depth (m)')
+plt.xlabel('Temperature $^\circ$C')
+plt.title('Particle temperature comparison')
+plt.savefig('particle_temperature_comparison.png')
+
+
+# salinity comparison
+plt.figure()
+plt.plot(particleSalinity.ravel(), zLevelParticle.ravel(), '.',
+         label='Particle value')
+plt.plot(salinity, z, '-', label='Initialized scalar value')
+plt.title('Particle salinity comparison')
+plt.ylabel('Depth (m)')
+plt.xlabel('Salinity psu')
+plt.savefig('particle_salinity_comparison.png')
+
diff --git a/testing_and_setup/compass/ocean/soma/soma_analysis_members.xml b/testing_and_setup/compass/ocean/soma/soma_analysis_members.xml
index 052672dfde..93d9c32155 100644
--- a/testing_and_setup/compass/ocean/soma/soma_analysis_members.xml
+++ b/testing_and_setup/compass/ocean/soma/soma_analysis_members.xml
@@ -43,6 +43,6 @@
 		<stream name="lagrPartTrackInput">
 			<attribute name="filename_template">particles.nc</attribute>
 		</stream>
-  </streams>
+	</streams>
 
 </template>
diff --git a/testing_and_setup/compass/ocean/ziso/20km/default/config_forward.xml b/testing_and_setup/compass/ocean/ziso/20km/default/config_forward.xml
index 299e3d44da..af04163b38 100644
--- a/testing_and_setup/compass/ocean/ziso/20km/default/config_forward.xml
+++ b/testing_and_setup/compass/ocean/ziso/20km/default/config_forward.xml
@@ -5,11 +5,7 @@
 	<add_link source="../init_step2/mesh.nc" dest="mesh.nc"/>
 	<add_link source="../init_step2/graph.info" dest="graph.info"/>
 
-	<get_file dest_path="initial_condition_database" file_name="make_particle_resets.151222.py">
-		<mirror protocol="wget" url="http://oceans11.lanl.gov/mpas_data/initial_condition_database"/>
-	</get_file>
-
-	<add_link source_path="initial_condition_database" source="make_particle_resets.151222.py" dest="make_particle_resets.py"/>
+  <add_link source_path="utility_scripts" source="/LIGHTparticles/make_particle_file.py" dest="make_particles.py"/>
 
 	<add_executable source="model" dest="ocean_model"/>
 	<add_executable source="metis" dest="metis"/>
@@ -23,6 +19,8 @@
 		<option name="config_pio_stride">4</option>
 		<option name="config_write_output_on_startup">.false.</option>
 		<option name="config_mom_del4">5.0e10</option>
+		<option name="config_AM_lagrPartTrack_sample_temperature">.true.</option>
+		<option name="config_AM_lagrPartTrack_sample_salinity">.true.</option>
 	</namelist>
 
 	<streams name="streams.ocean" keep="immutable" mode="forward">
@@ -39,20 +37,31 @@
 		<template file="lagrangian_particle_tracking.xml" path_base="script_core_dir" path="templates/analysis_members"/>
 		<template file="ziso_template.xml" path_base="script_configuration_dir"/>
 		<stream name="output">
-                        <attribute name="output_interval">0000_00:01:30</attribute>
-                </stream>
+			<attribute name="output_interval">0000_00:01:30</attribute>
+	</stream>
+		<stream name="lagrPartTrackOutput">
+			<add_contents>
+				<member name="particleTemperature" type="var"/>
+				<member name="particleSalinity" type="var"/>
+			</add_contents>
+	</stream>
 	</streams>
-
 	<run_script name="run.py">
 		<step executable="./metis">
 			<argument flag="graph.info">4</argument>
 		</step>
 		<step executable="python">
-			<argument flag="">make_particle_resets.py</argument>
-			<argument flag="">mesh.nc</argument>
+			<argument flag="">make_particles.py</argument>
+			<argument flag="">-i</argument>
+			<argument flag="">init.nc</argument>
+			<argument flag="">-o</argument>
 			<argument flag="">particles.nc</argument>
-			<argument flag="">graph.info.part.4</argument>
-			<argument flag="">11</argument>
+			<argument flag="">-g</argument>
+			<argument flag="">graph.info.part.</argument>
+			<argument flag="">-p</argument>
+			<argument flag="">4</argument>
+			<argument flag="">--types</argument>
+			<argument flag="">["buoyancy"]</argument>
 		</step>
 		<model_run procs="4" threads="1" namelist="namelist.ocean" streams="streams.ocean"/>
 	</run_script>
diff --git a/testing_and_setup/compass/utility_scripts/LIGHTparticles/animated_particles.py b/testing_and_setup/compass/utility_scripts/LIGHTparticles/animated_particles.py
new file mode 100644
index 0000000000..9b2bef1578
--- /dev/null
+++ b/testing_and_setup/compass/utility_scripts/LIGHTparticles/animated_particles.py
@@ -0,0 +1,133 @@
+#! /usr/bin/env python
+"""
+    Take an input netcdf database and extract particle information, converting to
+    vtu format for input into ParaView.
+
+    Usage exampe:
+    ./animated_particles.py -f free_floating.nc
+
+    produces
+
+    particles_free_floating*.vtu
+    particles_free_floating.pvd
+
+    Phillip Wolfram
+    LANL
+    07/15/2014
+
+"""
+from pyevtk.vtk import VtkGroup
+from pyevtk.hl import pointsToVTK
+import numpy as np
+
+import netCDF4
+import os
+import pdb
+
+def cleanup(data):
+    if type(data) == np.ma.core.MaskedArray:
+        data.data[data.mask] = data.fill_value
+        data = data.data
+    return data
+
+def build_particle_file(fname_in, fname_outbase):
+    """
+    Take existing netCDF4 input f_in and produce
+    vtu files for f_outbase with all the particle data
+
+    Phillip Wolfram
+    LANL
+    07/15/2014
+    """
+    ## initialize the starting points
+    fullpath = fname_in.split(os.sep)
+
+    # don't want to place this in starting directory
+    #g = VtkGroup(os.sep.join(fullpath[:-1]) + '/' + fname_outbase + '_' + os.path.splitext(fullpath[-1])[0] )
+    g = VtkGroup(fname_outbase + '_' + os.path.splitext(fullpath[-1])[0] )
+
+    # load the netCDF database
+
+    ## get time
+    #time = 0.0
+
+
+    # open the database
+    f_in = netCDF4.Dataset(fname_in, 'r')
+
+    # num particles and data frames
+    Npoints = len(f_in.dimensions['nParticles'])
+    Ntime = len(f_in.dimensions['Time'])
+
+    # initialize empty dictionary
+    particleData = dict()
+
+    for t in np.arange(Ntime):
+
+        # get the points locations
+        #x = f_in.variables['xParticle'][t] #+ f_in.variables['zLevelParticle'][t]
+        #y = f_in.variables['yParticle'][t] #+ f_in.variables['zLevelParticle'][t]
+        #z = f_in.variables['zParticle'][t] #+ f_in.variables['zLevelParticle'][t]
+        x = f_in.variables['xParticle'][t]
+        y = f_in.variables['yParticle'][t]
+        z = f_in.variables['zParticle'][t]
+        print t
+        print np.vstack((x,y,z)).T
+
+        particleType = f_in.variables['xParticle'].dimensions
+        particleTypeStatic = f_in.variables['indexToParticleID'].dimensions
+
+        # get particle data
+        particleData.clear()
+        for v in f_in.variables:
+            dim = f_in.variables[v].dimensions
+            if dim == particleType:
+                particleData[str(v)] = cleanup(f_in.variables[v][t])
+            if dim == particleTypeStatic:
+                particleData[str(v)] = cleanup(f_in.variables[v][:])
+
+        #print particleData
+
+        # set the new time (just use index for simplicity for now)
+        time = 2*t
+        #time = ''.join(f_in.variables['xtime'][t,0,:]).replace(' ','')
+
+        # file name
+        #f_out = fname_outbase + '_' + os.path.splitext(fname_in)[0] + '_' + str(t)
+        f_out = fname_outbase + '_' + os.path.splitext(fullpath[-1])[0] + '_' + str(t)
+        # output file
+        pointsToVTK(f_out, x, y, z, data = particleData)
+        # make sure the file is added to the time vector
+        g.addFile(filepath = f_out + '.vtu', sim_time = time)
+
+
+    # save the animation file
+    g.save()
+
+
+if __name__ == "__main__":
+    from optparse import OptionParser
+
+    # Get command line parameters
+    parser = OptionParser()
+    parser.add_option("-f", "--file", dest="inputfilename",
+            help="file to open for appending \
+                    particle data 'particle_' extension",
+                    metavar="FILE")
+
+    parser.add_option("-o", "--output", dest="outputfilename",
+            help="output file name base for export to *.vtu \
+                    file which stores particle data",
+                    metavar="FILE")
+
+    options, args = parser.parse_args()
+
+    if not options.inputfilename:
+        parser.error("Input filename is a required input.")
+
+    if not options.outputfilename:
+        parser.error("Output filename base for vtu files is required input.")
+
+    build_particle_file(options.inputfilename, options.outputfilename)
+
+# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python
diff --git a/testing_and_setup/compass/utility_scripts/LIGHTparticles/make_particle_file.py b/testing_and_setup/compass/utility_scripts/LIGHTparticles/make_particle_file.py
new file mode 100644
index 0000000000..735d4d8300
--- /dev/null
+++ b/testing_and_setup/compass/utility_scripts/LIGHTparticles/make_particle_file.py
@@ -0,0 +1,538 @@
+#!/usr/bin/env python
+"""
+File writes a particle input dataset for use in MPAS-O / E3SM.
+
+Example usage is
+
+    python make_particle_file.py -i init.nc -g graph.info.part.6 \
+            -o particles.nc -p 6 --spatialfilter SouthernOceanXYZ \
+            --downsample=1
+
+Phillip J. Wolfram
+Last Modified: 08/03/2018
+"""
+
+import netCDF4
+import numpy as np
+from scipy import spatial
+from scipy import sparse
+from pyamg.classical import split, interpolate as amginterp
+
+verticaltreatments = {'indexLevel':1, 'fixedZLevel': 2, 'passiveFloat': 3, 'buoyancySurface': 4, 'argoFloat': 5}
+defaults = {'dt': 300, 'resettime': 1.0*24.0*60.0*60.0}
+
+
+def use_defaults(name, val): #{{{
+    if (val is not None) or (val is not np.nan):
+        return val
+    else:
+        return defaults[name] #}}}
+
+
+def ensure_shape(start, new): #{{{
+    if isinstance(new, (int, float)):
+        new = new*np.ones_like(start)
+    return new #}}}
+
+
+def southern_ocean_only_xyz(x, y, z, maxNorth=-45.0): #{{{
+    sq = np.sqrt(x**2 + y**2 + z**2)
+    lat = np.arcsin(z / sq)
+    #lon = np.arctan2(y, x)
+    ok = np.pi/180.0*maxNorth
+    ids = lat < ok
+    return ids #}}}
+
+
+def southern_ocean_only_planar(x, y, z, maxy=1000.0*1e3): #{{{
+    ids = y < maxy
+    return ids #}}}
+
+
+def remap_particles(fin, fpart, fdecomp): #{{{
+    """
+    Remap particles onto a new grid decomposition.
+
+    Load in particle positions, locations of grid cell centers, and decomposition
+    corresponding to fin.
+
+    The goal is to update particle field currentBlock to comply with the new grid
+    as defined by fin and decomp.  NOTE: FIN AND FDECOMP MUST BE COMPATIBLE!
+
+    We assume that all particles will be within the domain such that a nearest
+    neighbor search is sufficient to make the remap.
+
+    Phillip Wolfram
+    LANL
+    Origin: 08/19/2014, Updated: 07/13/2018
+    """
+    # load the files
+    f_in = netCDF4.Dataset(fin, 'r')
+    f_part = netCDF4.Dataset(fpart,'r+')
+
+    # get the particle data
+    xpart = f_part.variables['xParticle']
+    ypart = f_part.variables['yParticle']
+    zpart = f_part.variables['zParticle']
+    currentBlock = f_part.variables['currentBlock']
+    try:
+        currentCell = f_part.variables['currentCell']
+        currentCellGlobalID = f_part.variables['currentCellGlobalID']
+    except:
+        currentCell = f_part.createVariable('currentCell', 'i', ('nParticles'))
+        currentCellGlobalID = f_part.createVariable('currentCellGlobalID', 'i', ('nParticles'))
+
+    # get the cell positions
+    xcell = f_in.variables['xCell']
+    ycell = f_in.variables['yCell']
+    zcell = f_in.variables['zCell']
+
+    # build the spatial tree
+    tree = spatial.cKDTree(np.vstack((xcell,ycell,zcell)).T)
+
+    # get nearest cell for each particle
+    dvEdge = f_in.variables['dvEdge']
+    maxdist = 2.0*max(dvEdge[:])
+    _, cellIndices = tree.query(np.vstack((xpart,ypart,zpart)).T,distance_upper_bound=maxdist,k=1)
+
+    # load the decomposition (apply to latest time step)
+    decomp = np.genfromtxt(fdecomp)
+    currentBlock[-1,:] = decomp[cellIndices]
+    currentCell[-1,:] = -1
+    currentCellGlobalID[-1,:] = cellIndices + 1
+
+    # close the files
+    f_in.close()
+    f_part.close()
+
+    #}}}
+
+
+def downsample_points(x, y, z, tri, nsplit): #{{{
+    """
+    Downsample points using algebraic multigrid splitting.
+
+    Note, currently assumes that all points on grid are equidistant, which does
+    a numeric (not area-weighted) downsampling.
+
+    Phillip Wolfram
+    LANL
+    Origin: 03/09/2015, Updated: 01/14/2019
+    """
+    # reference on cleanest way to do this calculation:
+    # https://www.mathworks.com/matlabcentral/answers/
+    # 369143-how-to-do-delaunay-triangulation-and-return-an-adjacency-matrix
+
+    # allocate the memory
+    Np = x.shape[0]
+    A = sparse.lil_matrix((Np, Np))
+
+    # cleanup impartial cells (don't include the triangles on boundary)
+    tri = tri[np.logical_not(np.any(tri==-1, axis=1)), :]
+
+    # handle one direction for triangles
+    A[tri[:,0], tri[:,1]] = 1
+    A[tri[:,1], tri[:,2]] = 1
+    A[tri[:,2], tri[:,0]] = 1
+
+    # handle other direction (bi-directional graph)
+    A[tri[:,1], tri[:,0]] = 1
+    A[tri[:,2], tri[:,1]] = 1
+    A[tri[:,0], tri[:,2]] = 1
+
+    A = A.tocsr()
+
+    Cpts = np.arange(Np)
+    # Grab root-nodes (i.e., Coarse / Fine splitting)
+    for ii in np.arange(nsplit):
+        splitting = split.PMIS(A)
+        # convert to index for subsetting particles
+        Cpts = Cpts[np.asarray(splitting, dtype=bool)]
+
+        if ii < nsplit - 1:
+            P = amginterp.direct_interpolation(A, A, splitting)
+            R = P.T.tocsr()
+            A = R * A * P
+
+    return Cpts, x[Cpts], y[Cpts], z[Cpts] #}}}
+
+class Particles(): #{{{
+
+    def __init__(self, x, y, z, cellindices, verticaltreatment, dt=np.nan, zlevel=np.nan,
+            indexlevel=np.nan, buoypart=np.nan, buoysurf=None, spatialfilter=None,
+            resettime=np.nan, xreset=np.nan, yreset=np.nan, zreset=np.nan, zlevelreset=np.nan): #{{{
+
+        # start with all the indicies and restrict
+        ids = np.ones((len(x)), dtype=bool)
+        if type(spatialfilter) is str:
+            spatialfilter = [spatialfilter]
+        if spatialfilter:
+            if np.max(['SouthernOceanXYZ' == afilter for afilter in spatialfilter]):
+                ids = np.logical_and(ids, southern_ocean_only_xyz(x,y,z))
+            if np.max(['SouthernOceanPlanar' == afilter for afilter in spatialfilter]):
+                ids = np.logical_and(ids, southern_ocean_only_planar(x,y,z))
+
+        self.x = x[ids]
+        self.y = y[ids]
+        self.z = z[ids]
+        self.verticaltreatment = ensure_shape(self.x, verticaltreatments[verticaltreatment])
+        self.nparticles = len(self.x)
+
+        self.dt = dt
+
+        # 3D passive floats
+        self.zlevel = ensure_shape(x, zlevel)[ids]
+
+        # isopycnal floats
+        if buoysurf is not None:
+            self.buoysurf = buoysurf
+        self.buoypart = ensure_shape(x, buoypart)[ids]
+        self.cellindices = cellindices[ids]
+        self.cellGlobalID = cellindices[ids]
+
+        # index level following floats
+        self.indexlevel = ensure_shape(x, indexlevel)[ids]
+
+        # reset features
+        self.resettime = ensure_shape(x, resettime)[ids]
+        self.xreset = ensure_shape(x, xreset)[ids]
+        self.yreset = ensure_shape(x, yreset)[ids]
+        self.zreset = ensure_shape(x, zreset)[ids]
+        self.zlevelreset = ensure_shape(x, zlevelreset)[ids]
+
+        return #}}}
+
+    def compute_lat_lon(self): #{{{
+        """
+        Ripped out whole-sale from latlon_coordinate_transforms.py
+        PJW 01/15/2019
+        """
+
+        x = self.x
+        y = self.y
+        z = self.z
+
+        self.latParticle = np.arcsin(z / np.sqrt(x ** 2 + y ** 2 + z ** 2))
+        self.lonParticle = np.arctan2(y, x)
+
+        return #}}}
+#}}}
+
+
+class ParticleList(): #{{{
+
+    def __init__(self, particlelist): #{{{
+        self.particlelist = particlelist #}}}
+
+
+    def aggregate(self): #{{{
+        self.len()
+
+        # buoyancysurf
+        buoysurf = np.array([])
+        for alist in self.particlelist:
+            if 'buoysurf' in dir(alist):
+                buoysurf = np.unique(np.setdiff1d(np.append(buoysurf, alist.buoysurf), None))
+        if len(buoysurf) > 0:
+            self.buoysurf = np.asarray(buoysurf, dtype='f8')
+        else:
+            self.buoysurf = None
+
+        return #}}}
+
+
+    def __getattr__(self, name): #{{{
+        # __getattr__ ensures self.x is concatenated properly
+        return self.concatenate(name) #}}}
+
+
+    def concatenate(self, varname): #{{{
+        var = getattr(self.particlelist[0], varname)
+        for alist in self.particlelist[1:]:
+            var = np.append(var, getattr(alist, varname))
+        return var #}}}
+
+
+    def append(particlelist): #{{{
+        self.particlelist.append(particlelist[:]) #}}}
+
+
+    def len(self): #{{{
+        self.nparticles = 0
+        for alist in self.particlelist:
+            self.nparticles += alist.nparticles
+
+        return self.nparticles #}}}
+
+    # probably a cleaner way to have this "fall through" to the particle instances themselves,
+    # but didn't have time to sort this all out so this isn't general for now
+    def compute_lat_lon(self): #{{{
+        for alist in self.particlelist:
+            alist.compute_lat_lon()
+        return #}}}
+
+
+    def write(self, f_name, f_decomp):  #{{{
+
+        decomp = np.genfromtxt(f_decomp)
+
+        self.aggregate()
+        assert max(decomp) < self.nparticles, \
+                'Number of particles must be larger than decomposition!'
+
+        f_out = netCDF4.Dataset(f_name, 'w',format='NETCDF3_64BIT_OFFSET')
+
+        f_out.createDimension('Time')
+        f_out.createDimension('nParticles', self.nparticles)
+
+        f_out.createVariable('xParticle', 'f8', ('Time','nParticles'))
+        f_out.createVariable('yParticle', 'f8', ('Time','nParticles'))
+        f_out.createVariable('zParticle', 'f8', ('Time','nParticles'))
+        f_out.createVariable('lonParticle', 'f8', ('Time','nParticles'))
+        f_out.createVariable('latParticle', 'f8', ('Time','nParticles'))
+        f_out.createVariable('zLevelParticle', 'f8', ('Time','nParticles'))
+        f_out.createVariable('dtParticle', 'f8', ('Time','nParticles'))
+        f_out.createVariable('buoyancyParticle', 'f8', ('Time','nParticles'))
+        f_out.createVariable('currentBlock', 'i', ('Time', 'nParticles'))
+        f_out.createVariable('currentCell', 'i', ('Time', 'nParticles'))
+        f_out.createVariable('currentCellGlobalID', 'i', ('Time', 'nParticles'))
+        f_out.createVariable('indexToParticleID', 'i', ('nParticles'))
+        f_out.createVariable('verticalTreatment', 'i', ('Time','nParticles'))
+        f_out.createVariable('indexLevel', 'i', ('Time','nParticles'))
+        f_out.createVariable('resetTime', 'i', ('nParticles'))
+        f_out.createVariable('currentBlockReset', 'i', ('nParticles'))
+        f_out.createVariable('currentCellReset', 'i', ('nParticles'))
+        f_out.createVariable('xParticleReset', 'f8', ('nParticles'))
+        f_out.createVariable('yParticleReset', 'f8', ('nParticles'))
+        f_out.createVariable('zParticleReset', 'f8', ('nParticles'))
+        f_out.createVariable('zLevelParticleReset', 'f8', ('nParticles'))
+
+        f_out.variables['xParticle'][0,:] = self.x
+        f_out.variables['yParticle'][0,:] = self.y
+        f_out.variables['zParticle'][0,:] = self.z
+
+        self.compute_lat_lon()
+        f_out.variables['lonParticle'][0,:] = self.lonParticle
+        f_out.variables['latParticle'][0,:] = self.latParticle
+
+        f_out.variables['verticalTreatment'][0,:] = self.verticaltreatment
+
+        f_out.variables['zLevelParticle'][0,:] = self.zlevel
+
+        if self.buoysurf is not None and len(self.buoysurf) > 0:
+            f_out.createDimension('nBuoyancySurfaces', len(self.buoysurf))
+            f_out.createVariable('buoyancySurfaceValues', 'f8', ('nBuoyancySurfaces'))
+            f_out.variables['buoyancyParticle'][0,:] = self.buoypart
+            f_out.variables['buoyancySurfaceValues'][:] = self.buoysurf
+
+        f_out.variables['dtParticle'][0,:] = defaults['dt']
+        # assume single-processor mode for now
+        f_out.variables['currentBlock'][:] = 0
+        f_out.variables['resetTime'][:] = defaults['resettime'] # reset each day
+        f_out.variables['indexLevel'][:] = 1
+        f_out.variables['indexToParticleID'][:] = np.arange(self.nparticles)
+
+        # resets
+        f_out.variables['currentBlock'][0,:] = decomp[self.cellindices]
+        f_out.variables['currentBlockReset'][:] = decomp[self.cellindices]
+        f_out.variables['currentCell'][0,:] = -1
+        f_out.variables['currentCellGlobalID'][0,:] = self.cellGlobalID + 1
+        f_out.variables['currentCellReset'][:] = -1
+        f_out.variables['xParticleReset'][:] = f_out.variables['xParticle'][0,:]
+        f_out.variables['yParticleReset'][:] = f_out.variables['yParticle'][0,:]
+        f_out.variables['zParticleReset'][:] = f_out.variables['zParticle'][0,:]
+        f_out.variables['zLevelParticleReset'][:] = f_out.variables['zLevelParticle'][0,:]
+
+        f_out.close()
+        return #}}}
+#}}}
+
+
+def get_cell_coords(f_init): #{{{
+    return f_init.variables['xCell'][:], \
+           f_init.variables['yCell'][:], \
+           f_init.variables['zCell'][:] #}}}
+
+
+def expand_nlevels(x, n): #{{{
+    return np.tile(x, (n)) #}}}
+
+
+def cell_centers(f_init, downsample): #{{{
+
+    f_init = netCDF4.Dataset(f_init,'r')
+    nparticles = len(f_init.dimensions['nCells'])
+    xCell, yCell, zCell = get_cell_coords(f_init)
+    if downsample:
+        tri = f_init.variables['cellsOnVertex'][:,:] - 1
+        cpts, xCell, yCell, zCell = downsample_points(xCell, yCell, zCell, tri, downsample)
+    else:
+        cpts = np.arange(len(xCell))
+    f_init.close()
+
+    return cpts, xCell, yCell, zCell  #}}}
+
+
+def build_isopycnal_particles(cpts, xCell, yCell, zCell, buoysurf, afilter): #{{{
+
+    nparticles = len(xCell)
+    nbuoysurf = buoysurf.shape[0]
+
+    x = expand_nlevels(xCell, nbuoysurf)
+    y = expand_nlevels(yCell, nbuoysurf)
+    z = expand_nlevels(zCell, nbuoysurf)
+
+    buoypart = (np.tile(buoysurf,(nparticles,1))).reshape(nparticles*nbuoysurf,order='F').copy()
+    cellindices = np.tile(cpts, (nbuoysurf))
+
+    return Particles(x, y, z, cellindices, 'buoyancySurface', buoypart=buoypart, buoysurf=buoysurf, spatialfilter=afilter) #}}}
+
+
+def build_passive_floats(cpts, xCell, yCell, zCell, f_init, nvertlevels, afilter, vertseedtype): #{{{
+
+    x = expand_nlevels(xCell, nvertlevels)
+    y = expand_nlevels(yCell, nvertlevels)
+    z = expand_nlevels(zCell, nvertlevels)
+    f_init = netCDF4.Dataset(f_init,'r')
+    if vertseedtype == 'linear':
+        wgts = np.linspace(0, 1, nvertlevels+2)[1:-1]
+    elif vertseedtype == 'log':
+        wgts = np.geomspace(1.0/(nvertlevels-1), 1, nvertlevels+1)[0:-1]
+    elif vertseedtype == 'denseCenter':
+        wgts = dense_center_seeding(nvertlevels)
+    else:
+        raise ValueError("Must designate `vertseedtype` as one of the following: ['linear', 'log', 'denseCenter']")
+    zlevel = -np.kron(wgts, f_init.variables['bottomDepth'][cpts])
+    cellindices = np.tile(cpts, (nvertlevels))
+    f_init.close()
+
+    return Particles(x, y, z, cellindices, 'passiveFloat', zlevel=zlevel, spatialfilter=afilter) #}}}
+
+def dense_center_seeding(nVert): #{{{
+    """
+    Distributes passive floats with 50% of them occurring between 40% and 60%
+    of the bottom depth.
+    """
+    nMid = np.ceil((1/2) * nVert)
+    nRem = nVert - nMid
+    if nRem % 2 != 0:
+        nMid += 1
+        nRem -= 1
+    upper = np.linspace(0, 0.4, (int(nRem) // 2) + 1)
+    center = np.linspace(0.4, 0.6, int(nMid) + 2)
+    lower = np.linspace(0.6, 1, (int(nRem) // 2) + 1)
+    c_wgts = np.concatenate([upper[1:], center[1:-1], lower[0:-1]])
+    return c_wgts
+
+def build_surface_floats(cpts, xCell, yCell, zCell, afilter): #{{{
+
+    x = expand_nlevels(xCell, 1)
+    y = expand_nlevels(yCell, 1)
+    z = expand_nlevels(zCell, 1)
+    cellindices = cpts
+
+    return Particles(x, y, z, cellindices, 'indexLevel', indexlevel=1, zlevel=0, spatialfilter=afilter) #}}}
+
+
+def build_particle_file(f_init, f_name, f_decomp, types, spatialfilter, buoySurf, nVertLevels, downsample, vertseedtype): #{{{
+
+    cpts, xCell, yCell, zCell = cell_centers(f_init, downsample)
+
+    # build particles
+    particlelist = []
+    if 'buoyancy' in types or 'all' in types:
+        particlelist.append(build_isopycnal_particles(cpts, xCell, yCell, zCell, buoySurf, spatialfilter))
+    if 'passive' in types or 'all' in types:
+        particlelist.append(build_passive_floats(cpts, xCell, yCell, zCell, f_init, nVertLevels, spatialfilter, vertseedtype))
+    # apply surface particles everywhere to ensure that LIGHT works
+    # (allow for some load-imbalance for filters)
+    if 'surface' in types or 'all' in types:
+        particlelist.append(build_surface_floats(cpts, xCell, yCell, zCell, spatialfilter))
+
+    # write particles to disk
+    ParticleList(particlelist).write(f_name, f_decomp)
+
+    return #}}}
+
+
+if __name__ == "__main__":
+    import argparse
+    import os
+
+    parser = argparse.ArgumentParser(
+            description=__doc__, formatter_class=argparse.RawTextHelpFormatter)
+
+    parser.add_argument("-i", "--init", dest="init",
+            help="Name of netCDF init file",
+            default="init.nc",
+            metavar="PATH/INIT_NAME.nc")
+    parser.add_argument("-g", "--graph", dest="graph",
+            default='graph.info.part.',
+            help="Path / name of graph file of form */*.info.part.",
+            metavar="PATH/graph.info.part.")
+    parser.add_argument("-o", "--particlefile", dest="particles",
+            default="particles.nc",
+            help="Path / name of netCDF particle file",
+            metavar="PATH/particles.nc")
+    parser.add_argument("-p", "--procs", dest="procs",
+            help="Number of processors",
+            metavar="INT")
+    parser.add_argument("-t", "--types", dest="types",
+            help="Types of particles",
+            default="all",
+            metavar="One or more of ['buoyancy', 'passive', 'surface', 'all']")
+    parser.add_argument("--nvertlevels", dest="nvertlevels",
+            default=10,
+            help="Number of vertical levels for passive, 3D floats",
+            metavar="INT")
+    parser.add_argument("--vertseedtype", dest="vertseedtype",
+            default="linear",
+            help="Method for seeding in the vertical",
+            metavar="One of ['linear', 'log', 'denseCenter']")
+    parser.add_argument("--nbuoysurf", dest="nbuoysurf",
+            default=11,
+            help="Number of buoyancy surfaces for isopycnally-constrained particles",
+            metavar="INT")
+    parser.add_argument("--potdensmin", dest="potdensmin",
+            default=1028.5,
+            help="Minimum value of potential density surface for isopycnally-constrained particles",
+            metavar="INT")
+    parser.add_argument("--potdensmax", dest="potdensmax",
+            default=1030.0,
+            help="Maximum value of potential density surface for isopycnally-constrained particles",
+            metavar="INT")
+    parser.add_argument("--spatialfilter", dest="spatialfilter",
+            default=None,
+            help="Apply a certain type of spatial filter, e.g., ['SouthernOceanPlanar', 'SouthernOceanXYZ']",
+            metavar="STRING")
+    parser.add_argument("--remap", dest="remap",
+            action="store_true",
+            help="Remap particle file based on input mesh and decomposition.")
+    parser.add_argument("-d", "--downsample", dest="downsample",
+            metavar="INT", default=0,
+            help="Downsample particle positions using AMG a number of times.")
+
+    args = parser.parse_args()
+
+    if not '.info.part.' in args.graph:
+        OSError('Graph file processor count is inconsistent with processors specified!')
+    if not ('.' + str(args.procs)) in args.graph:
+        args.graph = args.graph + str(args.procs)
+
+    if not os.path.exists(args.init):
+        raise OSError('Init file {} not found.'.format(args.init))
+    if not os.path.exists(args.graph):
+        raise OSError('Graph file {} not found.'.format(args.graph))
+
+    if not args.remap:
+        print('Building particle file...')
+        build_particle_file(args.init, args.particles, args.graph, args.types, args.spatialfilter,
+                np.linspace(args.potdensmin, args.potdensmax, int(args.nbuoysurf)), int(args.nvertlevels),
+                int(args.downsample), args.vertseedtype)
+        print('Done building particle file')
+    else:
+        print('Remapping particles...')
+        remap_particles(args.init, args.particles, args.graph)
+        print('Done remapping particles')
+
+# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python
diff --git a/testing_and_setup/compass/utility_scripts/LIGHTparticles/streams.ocean b/testing_and_setup/compass/utility_scripts/LIGHTparticles/streams.ocean
new file mode 100644
index 0000000000..272cc2d81e
--- /dev/null
+++ b/testing_and_setup/compass/utility_scripts/LIGHTparticles/streams.ocean
@@ -0,0 +1,108 @@
+<stream name="lagrPartTrackInput"
+        filename_template="/lcrc/group/acme/pwolfram/acme_scratch/20180509.GMPAS-IAF.T62_oQU240.anvil.particles/run/particles.nc"
+        input_interval="initial_only"
+        packages="lagrPartTrackAMPKG"
+        type="input"
+        io_type="pnetcdf"
+        immutable="true">
+
+	<var name="indexToParticleID"/>
+	<var name="currentBlock"/>
+	<var name="currentCell"/>
+	<var name="xParticle"/>
+	<var name="yParticle"/>
+	<var name="zParticle"/>
+	<var name="zLevelParticle"/>
+	<var name="xParticleReset"/>
+	<var name="yParticleReset"/>
+	<var name="zParticleReset"/>
+	<var name="zLevelParticleReset"/>
+	<var name="currentCellReset"/>
+	<var name="currentBlockReset"/>
+	<var name="timeSinceReset"/>
+	<var name="resetTime"/>
+	<var name="numTimesReset"/>
+	<var name="buoyancyParticle"/>
+	<var name="verticalTreatment"/>
+	<var name="dtParticle"/>
+	<var name="indexLevel"/>
+	<var name="transfered"/>
+	<var name="buoyancySurfaceValues"/>
+</stream>
+
+<stream name="lagrPartTrackRestart"
+        filename_interval="output_interval"
+        clobber_mode="truncate"
+        output_interval="stream:restart:output_interval"
+        filename_template="mpaso.rst.am.lagrangianParticleTrackingRestart.$Y-$M-$D_$S.nc"
+        input_interval="initial_only"
+        packages="lagrPartTrackAMPKG"
+        type="input;output"
+        io_type="pnetcdf"
+        immutable="true">
+
+    <var name="xtime"/>
+    <var name="indexToParticleID"/>
+    <var name="currentBlock"/>
+    <var name="currentCell"/>
+    <var name="xParticle"/>
+    <var name="yParticle"/>
+    <var name="zParticle"/>
+    <var name="zLevelParticle"/>
+    <var name="xParticleReset"/>
+    <var name="yParticleReset"/>
+    <var name="zParticleReset"/>
+    <var name="zLevelParticleReset"/>
+    <var name="currentCellReset"/>
+    <var name="currentBlockReset"/>
+    <var name="timeSinceReset"/>
+    <var name="resetTime"/>
+    <var name="numTimesReset"/>
+    <var name="buoyancyParticle"/>
+    <var name="verticalTreatment"/>
+    <var name="dtParticle"/>
+    <var name="indexLevel"/>
+    <var name="transfered"/>
+    <var name="buoyancySurfaceValues"/>
+</stream>
+
+<stream name="lagrPartTrackOutput"
+        runtime_format="single_file"
+        filename_interval="00-01-00_00:00:00"
+        clobber_mode="truncate"
+        reference_time="0001-01-01_00:00:00"
+        output_interval="00-00-02_00:00:00"
+        filename_template="mpaso.hist.am.lagrPartTrack.$Y-$M-$D_$h.$m.$s.nc"
+        packages="lagrPartTrackAMPKG"
+        io_type="pnetcdf"
+        type="output">
+
+    <var name="boundaryVertex"/>
+    <var name="xtime"/>
+    <var name="indexToParticleID"/>
+    <var name="currentBlock"/>
+    <var name="currentCell"/>
+    <var name="currentCellGlobalID"/>
+    <var name="lonParticle"/>
+    <var name="latParticle"/>
+    <var name="xParticle"/>
+    <var name="yParticle"/>
+    <var name="zParticle"/>
+    <var name="zLevelParticle"/>
+    <var name="xParticleReset"/>
+    <var name="yParticleReset"/>
+    <var name="zParticleReset"/>
+    <var name="zLevelParticleReset"/>
+    <var name="currentCellReset"/>
+    <var name="currentBlockReset"/>
+    <var name="timeSinceReset"/>
+    <var name="resetTime"/>
+    <var name="numTimesReset"/>
+    <var name="buoyancyParticle"/>
+    <var name="verticalTreatment"/>
+    <var name="dtParticle"/>
+    <var name="indexLevel"/>
+    <var name="cellOwnerBlock"/>
+    <var name="transfered"/>
+    <var name="buoyancySurfaceValues"/>
+</stream>
diff --git a/testing_and_setup/compass/utility_scripts/LIGHTparticles/user_nl_mpaso b/testing_and_setup/compass/utility_scripts/LIGHTparticles/user_nl_mpaso
new file mode 100644
index 0000000000..3a5992dc35
--- /dev/null
+++ b/testing_and_setup/compass/utility_scripts/LIGHTparticles/user_nl_mpaso
@@ -0,0 +1,32 @@
+!----------------------------------------------------------------------------------
+! Users should add all user specific namelist changes after these comments
+! in the form of
+!   namelist_var = new_namelist_value
+! *** EXCEPT FOR ***
+! 1. DO NOT CHANGE config_start_time, config_run_duration, config_stop_time,
+!    config_do_restart, config_Restart_timestamp_filename, config_calendar_type,
+!    config_set_restingThickness_to_IC, config_alter_ICs_for_pbcs
+! 2. To preview the namelists, invoke $CASEROOT preview-namelists and look at
+!    $CASEROOT/CaseDocs/mpaso_in
+!
+!
+! Users cannot use this file to configure streams. In order to configure streams,
+! write a stream configuration file, and place it in
+! SourceMods/src.mpaso/streams.ocean_forward
+!----------------------------------------------------------------------------------
+config_write_output_on_startup = .true.
+&am_lagrparttrack
+ config_am_lagrparttrack_compute_interval = 'dt'
+ config_am_lagrparttrack_compute_on_startup = .false.
+ config_am_lagrparttrack_enable = .true.
+ config_am_lagrparttrack_filter_number = 0
+ config_am_lagrparttrack_input_stream = 'lagrPartTrackInput'
+ config_am_lagrparttrack_output_stream = 'lagrPartTrackOutput'
+ config_am_lagrparttrack_restart_stream = 'lagrPartTrackRestart'
+ config_am_lagrparttrack_region_stream = 'none'
+ config_am_lagrparttrack_reset_criteria = 'none'
+ config_am_lagrparttrack_reset_global_timestamp = '0029_23:59:59'
+ config_am_lagrparttrack_reset_if_inside_region = .false.
+ config_am_lagrparttrack_reset_if_outside_region = .false.
+ config_am_lagrparttrack_write_on_startup = .false.
+/