Releasing particles at different times output question

[pecos27]

Hello!

I am using ocean parcels to release particles from reanalysis from NEMO (c-grid) output across a section but at staggered times i.e. 1,000 particles every 5 days over the time period 1993 - 2020 in the Nordic Seas. For one year, it's the equivalent of 73,000 particles released.
This is more of an advice/recommendation question, but would it be best to store the output zarr files in files split up by the year of release?

Or is it sufficient to store all the outputs (staggered particle release and running the experiment for 27 years) in a single zarr file?

Please let me know if my question is unclear at all or needs further clarification. Thank you!

This is my code for reference (just trying to release particles every 5 days for year 2019, running for 1 year) -- it also doesn't include the kernels, functions etc:

n_particles_per_release = 1_000

lon_bds = (-6.5, -2.5)
lat_bds = (61.3, 60.3)

lon=np.random.uniform(*lon_bds, size=(n_particles_per_release, ))
lat=np.random.uniform(*lat_bds, size=(n_particles_per_release, ))
depth = np.random.uniform(650, 1100, size=(n_particles_per_release, )) 


release_times = np.datetime64('2019-01-01') + np.arange(1 * 73) * np.timedelta64(5 * 24 * 3600, 's')

from parcels import (
    FieldSet, 
    ParticleSet, 
    JITParticle, 
    AdvectionRK4_3D, 
    ParticleFile,
    Variable
)
from parcels import ParticleFile

inpath='../'
ufiles = sorted(glob(f"{inpath}/U/*.nc"))
ufiles = [f for f in ufiles if '2019' in f]

vfiles = sorted(glob(f"{inpath}/V/*.nc")) 
vfiles = [f for f in vfiles if '2019' in f]

wfiles = sorted(glob(f"{inpath}/W/*.nc")) 
wfiles = [f for f in wfiles if '2019' in f]

sfiles = sorted(glob(f"{inpath}/S/*.nc"))
sfiles = [f for f in sfiles if '2019' in f]

tfiles = sorted(glob(f"{inpath}/T/*.nc")) 
tfiles = [f for f in tfiles if '2019' in f]



variables = {'U': 'vozocrtx', 'V': 'vomecrty', 'W': 'vovecrtz', 'S': 'vosaline', 'T': 'votemper'}

filenames = {
    "U": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": ufiles}, # must use same everywhere w,files. but w depth is 0, northeast corner of T grid is lid
    "V": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": vfiles},
    "W": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": wfiles},
    "S": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": sfiles},
    "T": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": tfiles}
}


c_grid_dimensions = {
    "lon": "nav_lon",
    "lat": "nav_lat",
    "depth": "depthw",
    "time": "time_counter",
}

dimensions = {
    "U": c_grid_dimensions,
    "V": c_grid_dimensions,
    "W": c_grid_dimensions,
    "S": c_grid_dimensions,
    "T": c_grid_dimensions
}

with warnings.catch_warnings():
    warnings.simplefilter("ignore", parcels.FileWarning)
    fieldsetC = parcels.FieldSet.from_nemo(filenames, variables, dimensions,
        allow_time_extrapolation=True)
        
        
        class SampleParticle(JITParticle):
    temp = Variable('temp', dtype=np.float32, initial=-100) #fieldsetC.T
    salt = Variable('salt', dtype=np.float32, initial=-100)
    age  = Variable('age',  dtype=np.float32, initial=0)
    uvel = Variable('uvel', dtype=np.float32, initial=0)
    vvel = Variable('vvel', dtype=np.float32, initial=0)




lon_release = lon # longitude of release
lat_release = lat # latitude of release
depth_release = depth  # depth of release, meters
time_release = release_times

pset = ParticleSet.from_list(fieldset=fieldsetC,
                                pclass=SampleParticle,
                                lon=np.tile(lon_release, len(release_times)),
                                lat=np.tile(lat_release, len(release_times)),
                                depth=np.tile(depth_release, len(release_times)),
                                time=np.repeat(release_times, len(lon_release)))
                                
outputfile = ParticleFile('parcels_output.zarr',pset,timedelta(hours=6), chunks=(n_particles_per_release,1_000)) 

adv = pset.Kernel(parcels.AdvectionRK4_3D)
age = pset.Kernel(age)
sample = pset.Kernel(sampling)
sample_UV = pset.Kernel(velocity_sampling)
deleteparticle = pset.Kernel(DeleteParticle_outside_domain_beached)

kernels = adv+sample+sample_UV+age+deleteparticle 


pset.execute(kernels, runtime=timedelta(days=365), dt=-timedelta(minutes=10), output_file=outputfile) #time step must be smaller than fieldset resolution

[erikvansebille]

Hi @pecos27, thanks for your question. See https://docs.oceanparcels.org/en/latest/examples/documentation_LargeRunsOutput.html and https://docs.oceanparcels.org/en/latest/examples/documentation_advanced_zarr.html (both written by @JamiePringle; who may be able to provide more information) for some starting points on the use of Parcles zarr output.

One note I already want to make is that setting your chunks=(n_particles_per_release,1_000) may not be ideal. I don't know how 1_000 is defined, but I would recommend the first argument to be len(pset) instead of n_particles_per_release; otherwise you may end up with a lot of files in your zarr folder and a big overhead burden to create and manage these files.

[pecos27]

Hi @erikvansebille, thank you so much for your quick reply; this is really helpful!