import warnings
import jax.numpy as jnp
import numpy as np
import pandas as pd
from tracmass.particle import ParticleState
from tracmass.gcm import FluxDataset, GridFields
[docs]
def ijk(
i1: int,
i2: int,
j1: int,
j2: int,
k1: int,
k2: int,
grid: GridFields | None = None,
use_subset_coordinates: bool = False,
**sub_kw,
) -> ParticleState:
"""Seed particles in every cell of an i/j/k index box.
Parameters
----------
i1, i2, j1, j2, k1, k2 : int
1-based TRACMASS cell index ranges (inclusive). By default *i* and
*j* are absolute, full-grid indices -- the frame of saved
trajectories -- whatever ``subset`` the grid was loaded with.
grid : GridFields or None, optional
The run grid. Its ``i_offset`` / ``j_offset`` convert absolute
indices to the subset-relative ones particles use while
integrating, cells outside the grid are dropped, and land cells
(``kmt == 0``) are skipped. Without a grid the indices are used
as given (no offset, no masking).
use_subset_coordinates : bool, optional
If True, *i*/*j* are already subset-relative (cell 1 is the
subset's first cell) and no offset is applied. Default False.
From settings: ``use_subset_coordinates`` in ``[env.seed]``.
**sub_kw
Sub-cell options (``sub_dist``, ``i_sub_num``, ...).
Returns
-------
ParticleState
Particles in subset-relative indices.
Raises
------
ValueError
If *grid* is given and no cell of the box lies inside it.
Examples
--------
>>> # glorys subset i1 = 2600, j1 = 500: absolute cells 2611..2621
>>> state = seed.ijk(2611, 2621, 511, 521, 1, 1, grid=grid)
>>> # the same cells in subset coordinates
>>> state = seed.ijk(11, 21, 11, 21, 1, 1, grid=grid,
... use_subset_coordinates=True)
"""
iarr, jarr, karr = np.mgrid[i1 : i2 + 1, j1 : j2 + 1, k1 : k2 + 1]
ivec, jvec, kvec = iarr.flatten(), jarr.flatten(), karr.flatten()
if grid is not None:
if not use_subset_coordinates:
ivec = ivec - int(grid.i_offset)
jvec = jvec - int(grid.j_offset)
inside = (ivec >= 1) & (ivec <= grid.imt) & (jvec >= 1) & (jvec <= grid.jmt)
if ivec.size and not inside.any():
frame = "subset" if use_subset_coordinates else "absolute"
raise ValueError(
f"seed.ijk: box i {i1}..{i2}, j {j1}..{j2} ({frame} "
f"coordinates) lies outside the grid, which spans absolute "
f"i {grid.i_offset + 1}..{grid.i_offset + grid.imt}, "
f"j {grid.j_offset + 1}..{grid.j_offset + grid.jmt}"
)
ivec, jvec, kvec = ivec[inside], jvec[inside], kvec[inside]
# i/j are 1-based TRACMASS cell indices; kmt is indexed 0-based.
mask = np.asarray(grid.kmt).astype(bool)[ivec - 1, jvec - 1]
ivec, jvec, kvec = ivec[mask], jvec[mask], kvec[mask]
return _expand_sub(ivec, jvec, kvec, **sub_kw)
[docs]
def from_settings(grid: GridFields | None = None, cfg=None) -> ParticleState:
"""Seed particles as described by the ``[env.seed]`` settings table.
Parameters
----------
grid : GridFields or None, optional
Passed on to the seeding function (``ijk`` uses it to skip land).
cfg : mapping or None, optional
The seed table. Default: ``seed`` from the active
``tracmass.config.settings`` environment.
Returns
-------
ParticleState
Raises
------
ValueError
If there is no seed table, or its ``approach`` is not a key of
:data:`SETTINGS_APPROACHES`.
Notes
-----
The table names an ``approach`` and gives its arguments in a ``box``
sub-table; every other key (``sub_dist``, ``i_sub_num``, ...) is passed
on as a sub-cell seeding option::
[omi.seed]
approach = "ijk"
sub_dist = "rect_xyz"
i_sub_num = 40
j_sub_num = 40
k_sub_num = 1
box = { i1 = 10, i2 = 20, j1 = 10, j2 = 20, k1 = 1, k2 = 1 }
``i``/``j``/``k`` in ``box`` are 1-based TRACMASS cell indices, with
``i``/``j`` absolute (full-grid) unless the table sets
``use_subset_coordinates = true`` (see :func:`ijk`). With
``approach = "latlon_box"`` the box is geographic instead (see
:func:`latlon_box`)::
box = { lat1 = -64, lat2 = -62, lon1 = -70, lon2 = -68, k1 = 1, k2 = 1 }
and ``approach = "latlon_point"`` takes ``lat``, ``lon``, ``dlat``
(and optionally ``dlon``, ``k1``, ``k2``). ``approach =
"latlon_list"`` seeds one particle per point, from a table file
(``seed_file = "points.csv"`` with ``lat``/``lon`` columns) or inline
``lat = [...]`` / ``lon = [...]`` lists (see :func:`latlon_list`).
Examples
--------
>>> tracmass.settings.setenv("omi")
>>> state = seed.from_settings(grid)
"""
if cfg is None:
from tracmass.config import settings
cfg = settings.get("seed", None)
if not cfg:
raise ValueError(
"seed.from_settings: no [env.seed] table in the active settings"
)
kw = {str(k).lower(): v for k, v in dict(cfg).items()}
approach = kw.pop("approach", "ijk")
kw.pop("seed_freq", None) # repeat seeding, read by the integrators
try:
fn = SETTINGS_APPROACHES[approach]
except KeyError:
raise ValueError(
f"seed.from_settings: unknown approach {approach!r}; "
f"expected one of {sorted(SETTINGS_APPROACHES)}"
) from None
box = {str(k).lower(): v for k, v in dict(kw.pop("box", None) or {}).items()}
use_subset = bool(kw.pop("use_subset_coordinates", False))
if fn is ijk:
kw["use_subset_coordinates"] = use_subset
return fn(**box, **kw, grid=grid)
def _grid_lonlat(grid: GridFields, gcm: str | None):
"""Cell-centre ``(lon, lat)`` of *grid* from the GCM adapter, ``(imt, jmt)``."""
import importlib
if gcm is None:
from tracmass.config import gcm_name
gcm = gcm_name()
if gcm is None:
raise ValueError(
"seed: lat/lon seeding needs the GCM; pass gcm=... or set "
"[env.gcm] name"
)
module = importlib.import_module(f"tracmass.gcm.{gcm}")
if not hasattr(module, "lonlat_grid"):
raise ValueError(
f"seed: GCM {gcm!r} has no lonlat_grid; lat/lon seeding needs an "
f"adapter with a geographic grid (glorys, ecco2, scrippsgog)"
)
return module.lonlat_grid(grid)
[docs]
def latlon_box(
lat1: float,
lat2: float,
lon1: float,
lon2: float,
kpos: int = 1,
grid: GridFields | None = None,
*,
k1: int | None = None,
k2: int | None = None,
gcm: str | None = None,
**sub_kw,
) -> ParticleState:
"""Seed every ocean cell whose centre lies in a lat/lon box.
Parameters
----------
lat1, lat2 : float
Latitude range in degrees (inclusive).
lon1, lon2 : float
Longitude range in degrees (inclusive), eastward from *lon1* to
*lon2*: any longitude convention works, and ``lon1 > lon2`` gives
a box across the dateline (e.g. ``170`` → ``-170``).
kpos : int, optional
1-based level to seed when *k1* is not given. Default 1.
grid : GridFields
The run grid (required): selects the window of the GCM's lon/lat
(for a ``subset``) and its ``kmt`` masks out land.
k1, k2 : int or None, optional
1-based level range (inclusive), as in :func:`ijk`; *k2* defaults
to *k1*. Levels below the local bottom (``k > kmt``) are skipped.
gcm : str or None, optional
Adapter providing ``lonlat_grid``. Default: ``[env.gcm] name``.
**sub_kw
Sub-cell options (``sub_dist``, ``i_sub_num``, ...), see :func:`ijk`.
Returns
-------
ParticleState
Raises
------
ValueError
If *grid* is missing, the GCM has no geographic grid, or the box
contains no ocean cell of the grid.
Examples
--------
>>> state = seed.latlon_box(-64, -62, -70, -68, grid=grid, k1=1,
... i_sub_num=10, j_sub_num=10)
From settings::
[omi.seed]
approach = "latlon_box"
i_sub_num = 10
box = { lat1 = -64, lat2 = -62, lon1 = -70, lon2 = -68, k1 = 1, k2 = 1 }
"""
if grid is None:
raise ValueError("seed.latlon_box: grid is required")
lon, lat = _grid_lonlat(grid, gcm)
width = (lon2 - lon1) % 360.0
in_box = (
(lat >= lat1)
& (lat <= lat2)
& ((lon - lon1) % 360.0 <= (360.0 if lon2 - lon1 >= 360.0 else width))
)
if k1 is None:
k1 = kpos
if k2 is None:
k2 = k1
kmt = np.asarray(grid.kmt)
ivec, jvec, kvec = [], [], []
for k in range(k1, k2 + 1):
ii, jj = np.nonzero(in_box & (kmt >= k))
ivec.append(ii + 1)
jvec.append(jj + 1)
kvec.append(np.full(ii.shape, k))
ivec, jvec, kvec = (np.concatenate(v) for v in (ivec, jvec, kvec))
if ivec.size == 0:
raise ValueError(
f"seed.latlon_box: no ocean cells in lat {lat1}..{lat2}, "
f"lon {lon1}..{lon2} (levels {k1}..{k2}); the grid spans lat "
f"{np.nanmin(lat):.1f}..{np.nanmax(lat):.1f}, lon "
f"{np.nanmin(lon):.1f}..{np.nanmax(lon):.1f}"
)
return _expand_sub(ivec, jvec, kvec, **sub_kw)
[docs]
def latlon_point(
lat: float,
lon: float,
dlat: float,
dlon: float | None = None,
kpos: int = 1,
grid: GridFields | None = None,
**kw,
) -> ParticleState:
"""Seed the ocean cells within *dlat* / *dlon* degrees of a point.
A :func:`latlon_box` of ``lat ± dlat`` by ``lon ± dlon`` (*dlon*
defaults to *dlat*); other keywords (``k1``, ``k2``, ``gcm``, sub-cell
options) are passed on.
"""
dlon = dlat if dlon is None else dlon
return latlon_box(
lat - dlat, lat + dlat, lon - dlon, lon + dlon, kpos, grid, **kw
)
[docs]
def latlon_list(
df: pd.DataFrame | str | None = None,
grid: GridFields | None = None,
*,
seed_file: str | None = None,
lat: list[float] | None = None,
lon: list[float] | None = None,
lat_col: str = "lat",
lon_col: str = "lon",
k_col: str = "k",
kpos: int = 1,
gcm: str | None = None,
**_ignored,
) -> ParticleState:
"""Seed one particle at each lat/lon position of a table.
Each row becomes a particle at the exact fractional grid position of
its point (not snapped to the cell centre), found with
:func:`tracmass.trajs.coords.lonlat_to_xy` on the adapter's
``lonlat_grid(grid)``. A point on land (or outside the grid) is moved
to the centre of the nearest ocean cell.
Parameters
----------
df : pd.DataFrame or str or None
Table with *lat_col* / *lon_col* columns (and optionally *k_col*,
the 1-based level of each point). A str is read as a file, like
*seed_file*.
grid : GridFields
The run grid (required): selects the GCM's lon/lat window (for a
``subset``) and its ``kmt`` masks out land.
seed_file : str or None, optional
Path to a ``.csv`` or ``.parquet`` file to read *df* from.
lat, lon : list of float or None, optional
Positions given directly instead of *df* / *seed_file*.
lat_col, lon_col, k_col : str, optional
Column names. Defaults ``"lat"``, ``"lon"``, ``"k"``.
kpos : int, optional
1-based level for every particle when there is no *k_col* column.
Default 1.
gcm : str or None, optional
Adapter providing ``lonlat_grid``. Default: ``[env.gcm] name``.
**_ignored
Sub-cell options from a settings table are accepted and ignored:
each point is seeded exactly once.
Returns
-------
ParticleState
Particles in subset-relative coordinates, in row order (rows whose
level lies below the local bottom, ``k > kmt``, are dropped).
Raises
------
ValueError
If *grid* or the positions are missing, a column is absent, or no
point lies in an ocean cell deep enough for its level.
Examples
--------
>>> df = pd.DataFrame({"lat": [-63.0, -62.5], "lon": [-69.0, -68.2]})
>>> state = seed.latlon_list(df, grid=grid)
From settings::
[omi.seed]
approach = "latlon_list"
seed_file = "stations.csv" # columns lat, lon (and optionally k)
or with the positions inline::
[omi.seed]
approach = "latlon_list"
lat = [-63.0, -62.5]
lon = [-69.0, -68.2]
"""
from tracmass.trajs.coords import lonlat_to_xy
if grid is None:
raise ValueError("seed.latlon_list: grid is required")
if isinstance(df, str):
df, seed_file = None, df
if df is None and seed_file is not None:
df = (
pd.read_parquet(seed_file)
if str(seed_file).endswith(".parquet")
else pd.read_csv(seed_file)
)
if df is None:
if lat is None or lon is None:
raise ValueError(
"seed.latlon_list: give a DataFrame, a seed_file, or lat and lon"
)
df = pd.DataFrame({lat_col: np.atleast_1d(lat), lon_col: np.atleast_1d(lon)})
missing = [c for c in (lat_col, lon_col) if c not in df.columns]
if missing:
raise ValueError(
f"seed.latlon_list: missing column(s) {missing}; "
f"the table has {list(df.columns)}"
)
lat_v = np.asarray(df[lat_col], dtype=float)
lon_v = np.asarray(df[lon_col], dtype=float)
if k_col in df.columns:
k_v = np.asarray(df[k_col], dtype=int)
else:
k_v = np.full(lat_v.shape, int(kpos))
lon_arr, lat_arr = _grid_lonlat(grid, gcm)
kmt = np.asarray(grid.kmt)
x, y = lonlat_to_xy(lon_v, lat_v, (lon_arr, lat_arr, kmt))
x, y = np.asarray(x, dtype=float), np.asarray(y, dtype=float)
i0, j0 = np.floor(x).astype(int), np.floor(y).astype(int)
keep = (k_v >= 1) & (k_v <= kmt[i0, j0])
if not keep.any():
raise ValueError(
f"seed.latlon_list: none of the {lat_v.size} points lies in an "
f"ocean cell deep enough for its level"
)
if not keep.all():
warnings.warn(
f"seed.latlon_list: dropped {int((~keep).sum())} of {keep.size} "
f"points whose level is below the local bottom (k > kmt)",
stacklevel=2,
)
# Cell k spans [k-1, k]; place the particle at the level's centre.
return _expand_sub(
x[keep], y[keep], k_v[keep] - 0.5, sub_dist="already_xyz"
)
#: Seeding approaches usable from a ``[env.seed]`` settings table, by name.
#: Each takes the table's ``box`` entries and sub-cell keys as keyword
#: arguments plus ``grid=``.
SETTINGS_APPROACHES = {
"ijk": ijk,
"latlon_box": latlon_box,
"latlon_point": latlon_point,
"latlon_list": latlon_list,
}
[docs]
def mask(mask: np.ndarray, kpos: int = 1, **sub_kw) -> ParticleState:
"""Seed particles at every True cell of a raw ``(j, i)``-ordered mask.
``mask`` must be in the grid's *native* row-major orientation — the
same orientation grid fields come in straight from xarray/NetCDF (e.g.
``seagrid.ecco2.grid_definition().kmt.values`` or ``seagrid.glorys``'s raw
``kmt``): latitude (``j``) first, longitude (``i``) second. This is
the *opposite* of ``GridFields.kmt``'s ``(imt, jmt)`` convention — do
not pass ``grid.kmt`` here directly; use :func:`ijk` for that.
Parameters
----------
mask : np.ndarray
Boolean array, shape ``(jmt, imt)`` for a 2-D mask or
``(jmt, imt, kmt)`` for a 3-D one. True marks a cell to seed.
kpos : int, optional
1-based vertical level for seeded particles when ``mask`` is 2-D.
Ignored (and superseded by the mask's own k-axis) when ``mask``
is 3-D. Default 1.
**sub_kw
Forwarded to :func:`_expand_sub` (e.g. ``sub_dist``, ``i_sub_num``).
Returns
-------
ParticleState
One particle per True cell, placed at the cell centre, at level
``kpos`` (2-D mask) or the cell's own k (3-D mask).
Warns
-----
UserWarning
If ``mask.shape[0] > mask.shape[1]`` (more j than i). Ocean grids
are normally near-global in longitude but span at most the poles
in latitude, so ``jmt > imt`` usually means ``mask`` was passed
already transposed to ``(imt, jmt)`` instead of the expected raw
``(jmt, imt)``.
Examples
--------
>>> kmt = seagrid.ecco2.grid_definition().kmt.values # (jmt, imt)
>>> state = seed.mask(kmt > 0, kpos=1)
"""
if mask.shape[0] > mask.shape[1]:
warnings.warn(
f"seed.mask() expects a raw (jmt, imt) grid (j first), but got "
f"shape {mask.shape} with more j than i "
f"({mask.shape[0]} > {mask.shape[1]}). Ocean grids are normally "
"wider in longitude (i) than latitude (j) -- this may be an "
"already-transposed (imt, jmt) array.",
stacklevel=2,
)
if mask.ndim == 2:
jmt, imt = mask.shape
jgrid, igrid = np.mgrid[1 : jmt + 1, 1 : imt + 1]
ivec, jvec = igrid[mask], jgrid[mask]
kvec = np.full(ivec.shape, kpos)
else:
jmt, imt, kmt = mask.shape
jgrid, igrid, kgrid = np.mgrid[1 : jmt + 1, 1 : imt + 1, 1 : kmt + 1]
ivec, jvec, kvec = igrid[mask], jgrid[mask], kgrid[mask]
return _expand_sub(ivec, jvec, kvec, **sub_kw)
def from_df(
df: pd.DataFrame,
i_col_name: str = "i",
j_col_name: str = "j",
k_col_name: str = "k",
flux: FluxDataset | None = None,
**sub_kw,
) -> ParticleState:
ivec, jvec, kvec = df[i_col_name].data, df[j_col_name].data, df[k_col_name].data
if flux is not None:
mask = flux.kmt.astype(bool).values[ivec, jvec]
ivec, jvec, kvec = ivec[mask], jvec[mask], kvec[mask]
return _expand_sub(ivec, jvec, kvec, **sub_kw)
def from_df_xyz(
df: pd.DataFrame,
i_col_name: str = "i",
j_col_name: str = "j",
k_col_name: str = "k",
flux: FluxDataset | None = None,
**sub_kw,
) -> ParticleState:
sub_kw["sub_dist"] = "already_xyz"
return from_df(
df=df,
i_col_name=i_col_name,
j_col_name=j_col_name,
k_col_name=k_col_name,
flux=flux,
**sub_kw,
)
def _expand_sub(
ivec: np.array, jvec: np.array, kvec: np.array, **sub_kw
) -> ParticleState:
sub_dist = sub_kw.get("sub_dist", "rect_xyz")
if sub_dist == "rect_xyz":
def expand(vec1, vec2, vec3, sub_num):
subvec = np.linspace(-1, 0, sub_num + 2)[1:-1]
vec1 = np.hstack([vec1 + pos for pos in subvec])
vec2 = np.hstack([vec2 for pos in subvec])
vec3 = np.hstack([vec3 for pos in subvec])
return vec1, vec2, vec3
x0, y0, z0 = expand(ivec, jvec, kvec, sub_kw.get("i_sub_num", 1))
y0, z0, x0 = expand(y0, z0, x0, sub_kw.get("j_sub_num", 1))
z0, x0, y0 = expand(z0, x0, y0, sub_kw.get("k_sub_num", 1))
elif sub_dist == "gaus_xyz":
mean = [0, 0]
cov = [[1, 0], [0, 100]]
rng = np.random.default_rng()
def norm(vec):
vec = vec / vec.max() / 2
vec = vec / np.abs(vec.min()) / 2
return vec[(vec < 0.5) & (vec > -0.5)][:5000]
def expand(vec, sub_num):
n_part = len(vec)
sub_x, sub_y = rng.multivariate_normal(mean, cov, n_part + 20).T
sub_x, sub_y = norm(sub_x), norm(sub_y)
x0 = expand(ivec, sub_kw.get("tot_sub_num", 1))
y0 = expand(jvec, sub_kw.get("tot_sub_num", 1))
z0 = expand(kvec, sub_kw.get("tot_sub_num", 1))
elif sub_dist == "already_xyz":
x0, y0, z0 = ivec, jvec, kvec
else:
# Cell ia spans [ia-1, ia]; place particle at cell centre ia-0.5.
x0, y0, z0 = ivec - 0.5, jvec - 0.5, kvec - 0.5
ia = jnp.floor(x0).astype(jnp.int32) + 1
ja = jnp.floor(y0).astype(jnp.int32) + 1
ka = jnp.floor(z0).astype(jnp.int32) + 1
n_particles = len(ia)
state = ParticleState(
x=x0,
y=y0,
z=z0,
ia=ia,
ja=ja,
ka=ka,
tt=jnp.zeros(n_particles),
ts=jnp.zeros(n_particles),
t0=jnp.zeros(n_particles),
subvol=jnp.ones(n_particles) * 86400.0, # tseas as subvol for ts tracking
niter=jnp.zeros(n_particles, dtype=jnp.int32),
active=jnp.ones(n_particles, dtype=jnp.bool_),
reason=jnp.zeros(n_particles, dtype=jnp.int16),
)
return state
"""
! loneparticles
IF (loneparticle>0 .AND. log_level>0) THEN
PRINT '(A,I7)','WARNING! Loneparticle is set to : ', loneparticle
ELSEIF (log_level>0) THEN
PRINT '(A,I9)',' Cells used : ', nsdMax
END IF
SUBROUTINE seed()
! --------------------------------------------------
!
! Purpose:
! Seeding particles in TRACMASS.
!
!
! Method:
! Populate the position seeding array in trj and nrj.
!
! --------------------------------------------------
REAL(DP) :: vol
IF (log_level >= 5) THEN
PRINT*,' entering seed '
END IF
findTime: DO jsd=1,nsdTim
IF (seed_tim(jsd) == ints) THEN
itim = seed_tim(jsd)
EXIT findTime
ELSE IF (seed_tim(jsd) /= ints .AND. jsd == nsdTim) THEN
itim = -1
RETURN
END IF
END DO findTime
!Loop over the seed size (nsdMax)
startLoop: DO jsd=1,nsdMax
IF ( (seedTime == 1 .OR. seedTime == 2) .AND. &
& (ints /= itim) ) THEN
CYCLE startLoop
ELSE IF (seedTime /= 1 .AND. seedTime /= 2) THEN
PRINT*,'timeStart =',seedTime,' is not a valid configuration!'
STOP
END IF
iist = seed_ijk(jsd,1)
ijst = seed_ijk(jsd,2)
ikst = seed_ijk(jsd,3)
CALL reverse()
isec = seed_set(jsd,1)
idir = seed_set(jsd,2)
! Skip the loop
IF (iist < 1) CYCLE startLoop
IF (ijst < 1) CYCLE startLoop
IF (ikst < 1) CYCLE startLoop
IF (iist > imt) CYCLE startLoop
IF (ijst > jmt) CYCLE startLoop
IF (ikst > km) CYCLE startLoop
IF (kmt(iist,ijst) == 0) CYCLE startLoop
vol = 0.
kb = ikst
jb = ijst
ib = iist
ibm = ib - 1
IF (ibm == 0) ibm = imt
SELECT CASE(isec)
CASE(1) ! Through eastern meridional-vertical surface
vol = uflux(iist,ijst,ikst,nsm)
CASE (2) ! Through northern zonal-vertical surface
vol = vflux(iist,ijst,ikst,nsm)
CASE (3) ! Through upper zonal-meridional surface
vol = 1.
! Vertical
CALL vertvel (ib,ibm,jb,kb)
#if defined w_explicit
vol = wflux(ib,jb,kb,nsm)
#elif w_2dim
vol = 1.
#else
vol = wflux(kb,nsm)
#endif
END SELECT
IF ((nff*idir*vol <= 0.d0 .AND. idir /= 0 ) .OR. (vol == 0.)) CYCLE startLoop
! Volume/mass transport needs to be positive
vol = ABS(vol)
! Number of trajectories for box (iist,ijst,ikst)
SELECT CASE (nqua)
CASE (1) ! partQuant particles per seed gridcell
num = INT(partQuant)
CASE (2) ! particles reflect mass transport at seeding.
! set by partQuant
num = INT(vol/partQuant)
IF (MOD(vol,partQuant) /= 0) num = num + 1
! For num>4, num will be assigned to a square number
IF (FLOAT(INT(SQRT(FLOAT(num)))) /= SQRT(FLOAT(num)) .AND. num>4) THEN
num = (INT(SQRT(FLOAT(num))) + 1)**2
END IF
CASE (3) ! particle reflects air/water mass/volume at seeding
vol = dzt(ib,jb,kb,1)
num = INT(vol/partQuant)
END SELECT
IF (num == 0) num = 1
! Subvol definition
IF (num == 1) THEN
ijt = 1
ikt = 1
ELSE
CALL split_grid()
END IF
#if w_2dim
! Makes no difference in depth
ijt = num; ikt = 1
#endif
subvol = vol / DBLE (ijt*ikt)
IF (subvol == 0.d0) THEN
PRINT *, ' Transport of particle is zero!!!'
PRINT *, ' vol :', vol
PRINT *, ' uflux :', uflux (ib, jb, kb,nsm)
PRINT *, ' vflux :', vflux (ib, jb, kb,nsm)
PRINT *, ' subvol :', subvol
STOP
ENDIF
! Determine start position for each particle
! --------------------------------------------------
ijjLoop: DO jjt=1,ijt
kkkLoop: DO jkt=1,ikt
IF(x1>IMT .or. x1<0 .or. y1>JMT .or. y1<0 .or. z1>KM .or. z1<0 ) THEN
PRINT *, ' Particles outside the domain !!!'
PRINT *, ' ib:', ib, 'ibm:',ib-1
PRINT *, ' jb:', jb, ' kb:',kb
PRINT *, ' [x1,y1,z1] = ',x1,y1,z1
STOP
END IF
! If tracers are activated, initialise only those defined between the limits imposed by tracer0min and tracermax0
IF (l_tracers) THEN
DO itrac = 1, numtracers
! Depending on the dimensions
ktracer = kb
IF (tracers(itrac)%dimension == '2D') ktracer = 1
! Compute the value of the tracer
IF (isec == 1) THEN ! Zonal wall
iistp = iist + 1
IF (iist == IMT) iistp = 1
tracervalue(itrac) = 0.5*tracers(itrac)%data(iist,jb,ktracer,nsm) &
+ 0.5*tracers(itrac)%data(iistp,jb,ktracer,nsm)
ELSE IF (isec == 2) THEN ! Meridional wall
tracervalue(itrac) = 0.5*tracers(itrac)%data(ib,ijst,ktracer,nsm) &
+ 0.5*tracers(itrac)%data(ib,ijst+1,ktracer,nsm)
ELSE IF (isec == 3) THEN ! Horizontal wall
IF (ikst == KM) THEN
tracervalue(itrac) = tracers(itrac)%data(ib,jb,ikst,nsm)
ELSE
tracervalue(itrac) = 0.5*tracers(itrac)%data(ib,jb,ikst,nsm) &
+ 0.5*tracers(itrac)%data(ib,jb,ikst+1,nsm)
END IF
END IF
tracerbinvalue(itrac,2) = tracerbin(tracervalue(itrac),itrac)
IF ((tracervalue(itrac)<tracer0min(itrac)) .OR. (tracervalue(itrac)>tracer0max(itrac))) THEN
CYCLE kkkLoop
END IF
END DO
END IF
! Only one particle for diagnistics purposes
IF ((loneparticle>0) .and. (ntrac.ne.loneparticle)) THEN
trajectories(ntrac)%active = .FALSE.
CYCLE kkkLoop
END IF
! Rerun/streamfunction options
IF ((l_rerun .EQV..TRUE.) .AND. (trajectories(ntrac)%lbas==-1)) THEN
trajectories(ntrac)%active = .FALSE.
CYCLE kkkLoop
END IF
! tt - time, fractions of ints
! ts - time [s] rel to start
ts = DBLE (ints-1)
tt = ts * tseas
tss = 0.d0
! Boxface
IF (isec ==1 .AND. nff*idir>0) boxface = 1
IF (isec ==1 .AND. nff*idir<0) boxface = 2
IF (isec ==2 .AND. nff*idir>0) boxface = 3
IF (isec ==2 .AND. nff*idir<0) boxface = 4
IF (isec ==3 .AND. nff*idir>0) boxface = 5
IF (isec ==3 .AND. nff*idir<0) boxface = 6
!Save initial particle position
IF(log_level >= 3) THEN
PRINT*,' write initial trajectory position '
END IF
CALL write_data('ini')
CALL write_data('run')
END DO kkkLoop
END DO ijjLoop
END DO startLoop
IF (log_level >= 5) THEN
PRINT*,' leaving seed'
END IF
END SUBROUTINE seed
"""
