How to do vertical nearest extrapolation of tracer fields all the way up to the surface (on the NEMO grid?)

[willirath]

I try to sample tracer field on the NEMO grid all the way to the surface. For adding the tracer variable, I largely follow the NEMO 3D tutorial and add the tracer fields (temperature in my case) as an additional fieldset.

Fully reproducible "minimal" example: https://nbviewer.org/gist/willirath/88689100f3ce573d50c551cf9c4e116f#Parcels-Experiment

The essential bits are

# all above is identical to the tutorial
fieldset = FieldSet.from_nemo(filenames, variables, dimensions)

temperature_field = Field.from_netcdf(
    filenames=tfiles,
    variable=("temperature", "votemper"),
    dimensions={"lon": "nav_lon", "lat": "nav_lat", "depth": "deptht", "time": "time_counter"},
    interp_method="nearest",
)

fieldset.add_field(temperature_field, "temperature")

and

class TemperatureParticle(JITParticle):
    temperature = Variable("temperature", initial=0)

def SampleTemperature(particle, fieldset, time):
    particle.temperature = fieldset.temperature[particle]

pset = ParticleSet.from_line(
    fieldset=fieldset,
    pclass=TemperatureParticle,
    size=1_000,
    start=(-80.0, 30.0),
    finish=(-55.0, 35.0),
    # depth=float(fieldset.temperature.grid.depth[0] + 0.1),  # this is working
    # depth=float(fieldset.temperature.grid.depth[0] + 0.0),  # this is working
    depth=float(fieldset.temperature.grid.depth[0] - 0.1),  # this is NOT working
)

This works as expected and picks the correct temperature values using nearest-neighbor interpolation if the vertical position of the particles is below the center of the uppermost tracer grid cell but fails with a ThroughSurfaceError if the particles go into the upper half of the uppermost grid box.

Is there a way to extend the tracer field to the surface?

[willirath]

(A possible hack could be to add a 2d temperature field directly at the surface? But I didn't try that yet.)

[erikvansebille]

Thanks for reporting, @willirath. I think that the new #1402 provides a nice, clean solution to this problem. In #1402 (which will become the basis of Parcels v3, to be launched this autumn), we do away with the RecoveryKernels system, allowing a much more fine-grained control of what to do when interpolation throws an error.

In this new code, you have access to the particle.state in the Kernels, so you can write a small Kernel that interpolates the right velocity if the particle throws a ThroughSurfaceError. For example (not tested):

def InterpolateNearSurfaceTemperature(particle, fieldset, time):
    if particle.state == StatusCode.ErrorThroughSurface:
        if particle.depth > 0:  # location still in the ocean 
            particle.temperature = fieldset.temperature[particle.time, mindepth, particle.lat, particle.lon]
            particle.state = StatusCode.Success
        else:
            particle.delete()  # assuming you want to delete particles in the air

where mindepth is the minimum depth in the model (e.g. 0.5?).

This of course has a slight overhead because it requires one more interpolation, but it is done in the main kernel loop itself so should be fast.

For more info, see also the new tutorial at https://github.com/OceanParcels/parcels/blob/vectorial_summing_zarr_approach/docs/examples/tutorial_kernelloop.ipynb that I wrote today

I'd be very keen to hear if this indeed works for you; as we're developing Parcels v3.0, we need these kinds of real world tests!