loop over years and use lastest location

[LAWTIDEBOY]

Hello everyone !

Sorry in advance if that is already documented, I went through the documentation but did not find any clear guidance.
In order to minimize memory issues (150 years simulations...), I am trying to run shorter (yearly) experiments starting over the previous year latest positions.

I saw that there is a "write_latest_locations" function, but I remain highly unsure how to use it, as a kernel ? ... any help welcome... thank you ! :/

Here an example of the script I am trying to design, I would like to write the last position in the ZARR file, and start the next simulation from the previous last positions of the particles:

from glob import glob
import xarray as xr

from parcels import (
    AdvectionRK4_3D,
    AdvectionRK4,
    FieldSet,
    JITParticle,
    ParticleSet,
    XarrayDecodedFilter,
    logger,
    Variable,
    ErrorCode,
)

from datetime import timedelta, datetime
from operator import attrgetter
import numpy as np
import math
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from netCDF4 import Dataset
from scipy.interpolate import griddata

[xi,zi] = np.meshgrid(np.arange(-10,130,40), z, indexing='ij')
y = np.random.randn(xi.shape[0],zi.shape[1])
xi = xi.reshape(-1, 1)
zi = zi.reshape(-1, 1)
yi = np.ones(len(xi))*80
print('number of particles = ',len(xi))

variables = {'U': 'u',
             'V': 'v',
             'W': 'w'}

full_dims = {'time': 'time',
             'depth': 'depth',
             'lat': 'lat',
             'lon': 'lon'}

dimensions = {'U': full_dims, 'V': full_dims, 'W': full_dims}

def periodicBC(particle, fieldset, time):
    if particle.lon < fieldset.halo_west:
        particle_dlon += fieldset.halo_east - fieldset.halo_west
    elif particle.lon > fieldset.halo_east:
        particle_dlon -= fieldset.halo_east - fieldset.halo_west

def pushback(particle, fieldset, time):
    particle.depth = 2.5

def CheckOutOfBounds(particle, fieldset, time):
    particle.delete()

timestep = 30.5
length_experience = 365 # 150 years
output_dir = '/scratch/projects/hbk00083/model_outputs/fesom2.1_recom/EXP1_REGULAR/TRAJECTORIES/'

for year in np.arange(1950,1951):
    print(year)
    ds_parcels = xr.open_mfdataset('/scratch/projects/hbk00083/model_outputs/fesom2.1_recom/EXP1_REGULAR/*'+str(year)+'.nc')
    ds_parcels = ds_parcels.rename({'nz1': 'depth'})
    ds_parcels = ds_parcels.transpose("time", "depth", "lat", "lon")
    ds_parcels
    print('loading data structure done...')
    fieldset = FieldSet.from_xarray_dataset(ds_parcels, variables, dimensions, allow_time_extrapolation=True, chunksize='auto')
    fieldset.add_constant('halo_west', fieldset.U.grid.lon[0])
    fieldset.add_constant('halo_east', fieldset.U.grid.lon[-1])
    fieldset.add_periodic_halo(zonal=True)

    pset = ParticleSet(fieldset=fieldset, 
                       pclass=JITParticle, 
                       lon=xi, # a vector of release longitudes      
                       lat=yi, # a vector of release latitudes
                       time=datetime(year, 1, 1),
                       depth = 2.5*np.ones(len(xi))) # zi, a vector of release depths

    output_file = pset.ParticleFile(name=output_dir+"Particles_"+str(year)+".zarr", outputdt=timedelta(days=timestep))
    kernels = pset.Kernel([AdvectionRK4_3D, periodicBC])
    print('executing kernel...')
    pset.execute(kernels,
                 runtime=timedelta(days=length_experience),
                 dt=timedelta(days=timestep),# positive = forward in time, -timedelta = backward in time
                 output_file=output_file,
                 recovery={ErrorCode.ErrorThroughSurface: pushback, ErrorCode.ErrorOutOfBounds: CheckOutOfBounds})

   # TO BE CONTINUED HERE PICKING UP LATEST LOCATION Xi, Yi, Zi FOR THE NEXT LOOP ?

PS: I am using latest parcels+python version in a conda environment, linux environment.

[erikvansebille]

You can use the write_latest_locations() function after the pset.execute, so output file.write_latest_locations(set, pset.time_nextloop)

Note also that you can use ParticleSet.from_particlefile() to create a new ParticleSet from a Zarr-file; that might be useful too?

[LAWTIDEBOY]

Thank you @erikvansebille, both solutions work very nicely :)

However, I figured out that it is quite impossible to install last version >3.0 of parcels on a HPC where the base conda environment is not controlled by user. I know, out of topic, I am creating a new topic for this purpose.

[michaeldenes]

While not an elegant solution, as a workaround, could you install parcels directly from github on your hpc? Something like pip install --upgrade git+https://github.com/OceanParcels/parcels.git might work?