Pb with sampling temperature with parcels.FieldSet.from_croco

[soizicg]

Question

Hello,

I have been testing the new implementation of croco and it's really good to read ROMS file. Advection of particles seems to behave really well, thanks for adding that option into Parcels!

However it looks like some features are not completely working yet. For my application I need to sample the temperature at the particle position, but it sends me the following error:

FieldOutOfBoundSurfaceError: Field sampled out-of-bound at the surface, at (depth=-44.6170337677002, lat=53.955386881397935, lon=-4.461946337284217)
 67%|██████▋   | 173400.0/259200.0 [00:57<00:01, 70961.93it/s]

In this test (code below), the temperature is not sampled before 2 days of simulation, and that's when the error is produced (67% of the 3 days simulation). The particle is kept at a depth of 50 m and is relocated at at the seabed (considered as the depth z = -0.9 * h) if the bathymetry h is shallower at the particle's position. In the AdvectionRK4_cst3D_CROCO kernel, I used fieldset.UVW because fieldset.UV doesn't work yet with a croco fieldset.
Note that the following code works if I use a ScipyParticle class instead of JITParticle, but the running time is very long for a 70 days simulation. I was wondering if there is a way to ensure it's also working with JITParticle?

import xarray as xr
from datetime import timedelta as delta
import parcels


class AgeParticle(parcels.JITParticle): 
    age = parcels.Variable("age", initial=-600/3600)  # Variable 'age' is added with initial value 0.
    T = parcels.Variable("T")  

# Create a custom kernel which keeps track of the particle age (hours)
def AgingKernel(particle, fieldset, time):
    particle.age += particle.dt / 3600

def SampleT(particle, fieldset, time):
    if particle.age > 48:
        h = fieldset.H[0, 0, particle.lat,  particle.lon]
        if h <= -fieldset.cst_depth:
            depth = -0.9 * h
        else:
            depth = fieldset.cst_depth     
        particle.depth = depth
        
        particle.T = fieldset.T[time, particle.depth, particle.lat, particle.lon, particle]
    
def AdvectionRK4_cst3D_CROCO(particle, fieldset, time):
    """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity.
    This kernel assumes the vertical velocity is the 'w' field from CROCO output and works on sigma-layers.
    """
    h1 = fieldset.H[0, 0, particle.lat, particle.lon]
    if h1 <= -fieldset.cst_depth:
        depth1 = -0.9 * h1
    else:
        depth1 = fieldset.cst_depth    
    (u1, v1, w1) = fieldset.UVW[time, depth1, particle.lat, particle.lon, particle]
    lon1 = particle.lon + u1 * 0.5 * particle.dt
    lat1 = particle.lat + v1 * 0.5 * particle.dt
    
    h2 = fieldset.H[0, 0, lat1, lon1]
    if h2 <= -fieldset.cst_depth:
        depth2 = -0.9 * h2
    else:
        depth2 = fieldset.cst_depth
    (u2, v2, w2) = fieldset.UVW[time + 0.5 * particle.dt, depth2, lat1, lon1, particle]
    lon2 = particle.lon + u2 * 0.5 * particle.dt
    lat2 = particle.lat + v2 * 0.5 * particle.dt
    
    h3 = fieldset.H[0, 0, lat2, lon2]
    if h3 <= -fieldset.cst_depth:
        depth3 = -0.9 * h3
    else:
        depth3 = fieldset.cst_depth        
    (u3, v3, w3) = fieldset.UVW[time + 0.5 * particle.dt, depth3, lat2, lon2, particle]
    lon3 = particle.lon + u3 * particle.dt
    lat3 = particle.lat + v3 * particle.dt

    h4 = fieldset.H[0, 0, lat3, lon3]
    if h4 <= -fieldset.cst_depth:
        depth4 = -0.9 * h4
    else:
        depth4 = fieldset.cst_depth      
    (u4, v4, w4) = fieldset.UVW[time + particle.dt, depth4, lat3, lon3, particle]

    particle_dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * particle.dt  # noqa
    particle_dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * particle.dt  # noqa
    h5 = fieldset.H[0, 0, particle.lat + particle_dlat, particle.lon + particle_dlon]
    if h5 <= -fieldset.cst_depth:
        depth5 = -0.9 * h5
    else:
        depth5 = fieldset.cst_depth     
    particle.depth = depth5
    

                
# forcing_filename = "../forcing/may2019_all_oceanparcels.txt" 
forcing_filename = "../forcing/may2019_test.txt" 
netcdf_files = open(forcing_filename,'r').read().split('\n')


variables = {"U": "u", "V": "v", "W": "w0", "H": "h", "Zeta": "zeta0", "Cs_w": "Cs_w", "T": "temp"}

lon_rho = "lon_psi"  # Note, this would be "lon_rho" for a dataset on a spherical grid
lat_rho = "lat_psi"  # Note ,this would be "lat_rho" for a dataset on a spherical grid

dimensions = {
    "U": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
    "V": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
    "W": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
    "H": {"lon": lon_rho, "lat": lat_rho},
    "Zeta": {"lon": lon_rho, "lat": lat_rho, "time": "time"},
    "Cs_w": {"depth": "s_w"},
    "T": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
}

fieldset = parcels.FieldSet.from_croco(
    netcdf_files,
    variables,
    dimensions,
    hc=10,  # Note this stretching parameter is needed for the vertical grid
    allow_time_extrapolation=True,  # Note, this is only needed for this specific example dataset, that has only one snapshot
    mesh="spherical",  # Note, this is only needed for this specific example dataset, that has been created on a 'flat' mesh (i.e. in km instead of in degrees)
)

fieldset.add_constant("cst_depth", -50)

    
pset = parcels.ParticleSet(
    fieldset, # the fields that the particleset uses
    pclass=AgeParticle, # define the type of particle
    lon=-4.7740,  # release longitude
    lat=54.0097, # release latitude 
    depth=-50, # release latitude
    time=0 # release time
)      

output_filename = "../output/croco3D_10m_forum"
outputfile = pset.ParticleFile(
    name=(output_filename+".zarr"), 
    outputdt=3600
)

pset.execute(
    [AdvectionRK4_cst3D_CROCO, SampleT, AgingKernel],
    runtime=delta(days=3),
    dt=600,
    output_file=outputfile,
)

[erikvansebille]

Hi @soizicg, thanks for using Parcels and particularly your enthusiastic testing of the CROCO version!

I had a quick look through your code, and one thing I noticed is that you use particle.depth = depth5 at the end of your AdvectionRK4_cst3D_CROCO Kernel.

In general, it's better to change particle_ddepth (see also this part of the Kernel loop tutorial), so for example

particle_ddepth += depth5 - particle.depth

Does that help?

[soizicg]

Thanks @erikvansebille for your suggestion, it is a much cleaner way to code indeed, but it doesn't solve the issue.

I think the source of the issue might lie with the sigma coordinates. When I had a look into the code of Parcels it looks like an extra function croco_from_z_to_sigma helps to convert the vertical coordinates.
The error I get:
FieldOutOfBoundSurfaceError: Field sampled out-of-bound at the surface, at (depth=-44.6170337677002, lat=53.955386881397935, lon=-4.461946337284217)
suggests that the particle went through the surface even if the depth is z=-44, so clearly in the water column. So I wonder if the depth coordinates are not properly converted for JITParticles with croco fields?

[erikvansebille]

AH, it could indeed be a bug in the JIT implementation. Can you run the pest.execute() with delete_cfiles=False and then upload the generated C-file? Perhaps we see an obvious issue by inspecting the raw C-file?

[soizicg]

Can you tell me where I find the C-file? Is it the file in parcels/pycache directory of my parcels directory?

[erikvansebille]

If you added delete_cfiles=False to the pset.execute() call, you should see an output message that the C-file is stored in a certain location on your hard-drive. Can you find it in that location?

[soizicg]

Thanks for the directions, I have found them and attach them:

cfiles.zip

[soizicg]

@erikvansebille, thanks for pointing me to the cfiles, I have found the issue.
It was not a problem with Parcels. I was experimenting with another use of croco, and my modifications caused the issue...

Apologies for that, but thank you I have learn a lot!

[erikvansebille]

Ah good; so what was the issue then? Would it be useful to briefly explain that here, so that others don't make the same mistake?