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 5952322271..a9860bf35d 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,10 @@
 			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."
@@ -62,7 +66,39 @@
 		<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>
@@ -328,11 +364,35 @@
 			 description="flag to monitor if the particle was transfered"
 		/>
 <!-- TEMP FIELDS THAT NEED MOVED TO NonHalo below after verifying interpolation is correct -->
-		<var name="temperature" type="real" dimensions="nParticles Time" units="kg m^{-3}"
+		<var name="particleTemperature" type="real" dimensions="nParticles Time" units="degrees Celsius"
 			 description="sampled temperature for particle"
 		/>
-		<var name="salinity" type="real" dimensions="nParticles Time" units="kg m^{-3}"
+		<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>
 <!--
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 7ddfe41a95..ea43daf262 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
@@ -284,6 +284,7 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
       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, &
@@ -309,14 +310,61 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
       integer :: err_tmp
 
       ! sampling variables
-      integer, pointer :: index_temperature, index_salinity
+      ! active tracers
+      integer, pointer :: index_temperature, &
+                          index_salinity
       real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers
 
       logical, pointer :: config_AM_lagrPartTrack_sample_temperature
-      real (kind=RKIND), pointer :: 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 :: 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
 
@@ -338,9 +386,20 @@ 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.
       else
@@ -401,6 +460,17 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
         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)
@@ -411,6 +481,7 @@ 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)
@@ -506,8 +577,20 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
           call mpas_pool_get_array(particle % haloDataPool, 'currentCellGlobalID', currentCellGlobalID)
 
           ! sampling variables
-          call mpas_pool_get_array(particle % haloDataPool, 'temperature', temperature)
-          call mpas_pool_get_array(particle % haloDataPool, 'salinity', salinity)
+          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)
@@ -734,67 +817,9 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
           call mpas_timer_stop("get_validated_cell_idLPT")
 #endif
 
-          ! ensure that buoyancySurface particles have the correct zLevelParticle (interpolated based on their location)
-          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
-            iLevel = 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(iLevel, 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
-            if (verticalTreatment == 3) then !('passiveFloat')
-              ! constrain vertical level to be within dz/2 of ssh or bottom depth
-              zLevelParticle = min(max(zMid(maxLevelCell(iCell),iCell), zLevelParticle), zMid(1,iCell))
-            end if
-            iLevel = mpas_get_vertical_id(maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell))
-
-            if (verticalTreatment == 1) then !('indexLevel')
-              ! assign indexLevel particle vertical position to dz/2 location in vertical cell midpoint
-              zLevelParticle = zMid(indexLevel, iCell)
-              iLevel = indexLevel
-            end if
-#ifdef MPAS_DEBUG
-            call mpas_timer_stop("mpas_get_vertical_idLPT")
-#endif
-          end if
-
-#ifdef MPAS_DEBUG
-          call mpas_timer_start("verifyParticleDepthInDomainLPT")
-#endif
-          !verify particles are still in domain !{{{
-#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
+          call get_vertical_level_info(iLevel, zLevelParticle, &
+            verticalTreatment, maxLevelCell, iCell, buoyancyInterp, buoyancyTimeInterp, &
+            zMid, indexLevel, ssh, bottomDepth, indexToParticleID)
 
            !! sampling setup (if needed)
            !if (&
@@ -808,17 +833,218 @@ subroutine ocn_compute_lagrangian_particle_tracking(domain, timeLevel, err)!{{{
 
            if ( associated(activeTracers) ) then
              if (config_AM_lagrPartTrack_sample_temperature) then
-              !call mpas_log_write(' testing output for temperature')
-              call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
-                activeTracers(index_temperature,:,iCell), temperature)
-              !call mpas_log_write('done testing output for temperature')
+              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, &
+                  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')
-              call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
-                activeTracers(index_salinity,:,iCell), salinity)
-              !call mpas_log_write('done testing output for salinity')
+              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, &
+                  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, &
+                     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, &
+                     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, &
+                     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, &
+                     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, &
+                     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, &
+                     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, &
+                     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, &
+                     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
 
@@ -1774,6 +2000,7 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
       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
@@ -1796,15 +2023,54 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
       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 :: temperature
+      real (kind=RKIND), pointer :: particleTemperature
       logical, pointer :: config_AM_lagrPartTrack_sample_salinity
-      real (kind=RKIND), pointer :: 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
 
@@ -1843,9 +2109,21 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
         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)
@@ -1889,8 +2167,19 @@ 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, 'temperature', temperature)
-          call mpas_pool_get_array(particle % haloDataPool, 'salinity', salinity)
+          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)
@@ -1920,22 +2209,9 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
           call mpas_timer_stop("get_validated_cell_idLPT")
 #endif
 
-          if(verticalTreatment == 4) then  !('buoyancySurface') !{{{
-            iLevel = mpas_get_vertical_id(maxLevelCell(iCell), buoyancyParticle, buoyancyTimeInterp(:,iCell))
-
-            ! interpolate 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
-            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
+          call get_vertical_level_info(iLevel, zLevelParticle, &
+            verticalTreatment, maxLevelCell, iCell, buoyancyParticle, buoyancyTimeInterp, &
+            zMid, indexLevel)
 
           ! store global cell ID
           currentCellGlobalID = indexToCellID(iCell)
@@ -1965,17 +2241,218 @@ subroutine initialize_particle_properties(domain, timeLevel, err)!{{{
           !}}}
           if ( associated(activeTracers) ) then
             if (config_AM_lagrPartTrack_sample_temperature) then
-             !call mpas_log_write(' testing output for temperature')
-             call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
-               activeTracers(index_temperature,:,iCell), temperature)
-             !call mpas_log_write('done testing output for temperature')
+             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, &
+                   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')
-             call interp_cell_scalars(iLevel, maxLevelCell(iCell), -zLevelParticle, -zMid(:,iCell), &
-               activeTracers(index_salinity,:,iCell), salinity)
-             !call mpas_log_write('done testing output for salinity')
+             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, &
+                   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, &
+                    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, &
+                    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, &
+                     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, &
+                     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, &
+                     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, &
+                     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, &
+                     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, &
+                     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
 
@@ -2498,6 +2975,414 @@ 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
+
+#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, &
+       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
+
+     !-----------------------------------------------------------------
+     ! output variables
+     !-----------------------------------------------------------------
+     real (kind=RKIND), pointer :: particleScalar                     !< interpolated scalar value for the particle
+
+     !-----------------------------------------------------------------
+     ! local variables
+     !-----------------------------------------------------------------
+     integer :: iLevel                                       !< vertical level index
+     real (kind=RKIND) :: zLevelParticle                     !< vertical level z location
+     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
+
+     ! 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
+
+     ! 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
+       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
+     !}}}
+
+
+#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
+     !-----------------------------------------------------------------
+
+     ! {{{
+     ! cellsOnVertex assumes an "empty" neighbor (edge case) is 0
+     do aCell = 1, vertexDegree
+       ! only handles the case where the vertex is an interior vertex, otherwise treat as cell-based value
+       ! (no horizontal interpolation)
+       if (vertexID .eq. -1 .or. &
+         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
+
+         !call mpas_log_write('Using column interpolation for $i because it is on an edge', intArgs = (/iCell/))
+
+         ! 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
+       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)
+
+!#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
@@ -2798,7 +3683,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( 'Warning: 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)
@@ -2850,6 +3735,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
 
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/make_particle_file.py b/testing_and_setup/compass/utility_scripts/LIGHTparticles/make_particle_file.py
index 9ba5586150..70536ad97c 100644
--- a/testing_and_setup/compass/utility_scripts/LIGHTparticles/make_particle_file.py
+++ b/testing_and_setup/compass/utility_scripts/LIGHTparticles/make_particle_file.py
@@ -340,19 +340,41 @@ def build_isopycnal_particles(xCell, yCell, zCell, buoysurf, afilter): #{{{
     return Particles(x, y, z, cellindices, 'buoyancySurface', buoypart=buoypart, buoysurf=buoysurf, spatialfilter=afilter) #}}}
 
 
-def build_passive_floats(xCell, yCell, zCell, f_init, nvertlevels, afilter): #{{{
+def build_passive_floats(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')
-    zlevel = -np.kron(np.linspace(0,1,nvertlevels+2)[1:-1], f_init.variables['bottomDepth'][:])
+    if vertseedtype == 'linear':
+        wgts = np.linspace(0, 1, nvertlevels+2)[1:-1]
+    elif vertseedtype == 'log':
+        wgts = np.geomspace(1/(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'][:])
     cellindices = np.tile(np.arange(len(xCell)), (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(xCell, yCell, zCell, afilter): #{{{
 
@@ -364,7 +386,7 @@ def build_surface_floats(xCell, yCell, zCell, afilter): #{{{
     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): #{{{
+def build_particle_file(f_init, f_name, f_decomp, types, spatialfilter, buoySurf, nVertLevels, downsample, vertseedtype): #{{{
 
     xCell, yCell, zCell = cell_centers(f_init, downsample)
 
@@ -373,7 +395,7 @@ def build_particle_file(f_init, f_name, f_decomp, types, spatialfilter, buoySurf
     if 'buoyancy' in types or 'all' in types:
         particlelist.append(build_isopycnal_particles(xCell, yCell, zCell, buoySurf, spatialfilter))
     if 'passive' in types or 'all' in types:
-        particlelist.append(build_passive_floats(xCell, yCell, zCell, f_init, nVertLevels, spatialfilter ))
+        particlelist.append(build_passive_floats(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:
@@ -415,6 +437,10 @@ def build_particle_file(f_init, f_name, f_decomp, types, spatialfilter, buoySurf
             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",
@@ -454,7 +480,7 @@ def build_particle_file(f_init, f_name, f_decomp, types, spatialfilter, buoySurf
         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))
+                int(args.downsample), args.vertseedtype)
         print('Done building particle file')
     else:
         print('Remapping particles...')
diff --git a/testing_and_setup/compass/utility_scripts/LIGHTparticles/streams.ocean b/testing_and_setup/compass/utility_scripts/LIGHTparticles/streams.ocean
index 5674660df7..272cc2d81e 100644
--- a/testing_and_setup/compass/utility_scripts/LIGHTparticles/streams.ocean
+++ b/testing_and_setup/compass/utility_scripts/LIGHTparticles/streams.ocean
@@ -68,7 +68,7 @@
 
 <stream name="lagrPartTrackOutput"
         runtime_format="single_file"
-        filename_interval="01-00-00_00:00:00"
+        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"
@@ -82,6 +82,9 @@
     <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"/>