diff --git a/docs/examples/tutorial_gsw_density.ipynb b/docs/examples/tutorial_gsw_density.ipynb new file mode 100644 index 0000000000..e6eda9a03e --- /dev/null +++ b/docs/examples/tutorial_gsw_density.ipynb @@ -0,0 +1,154 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Using the `gsw` toolbox to compute density" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "This tutorial shows how to use the [`gsw` toolbox](https://teos-10.github.io/GSW-Python/) (the Gibbs SeaWater Oceanographic Toolbox of TEOS-10) within Parcels to compute density from temperature and salinity fields. The `gsw` toolbox can be installed via `conda install gsw`." + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "First, load the necessary libraries and the data:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import gsw # The package to calculate density from temperature, salinity and depth\n", + "import numpy as np\n", + "import xarray as xr\n", + "\n", + "import parcels\n", + "\n", + "# Load the CopernicusMarine data in the Agulhas region from the example_datasets\n", + "example_dataset_folder = parcels.download_example_dataset(\n", + " \"CopernicusMarine_data_for_Argo_tutorial\"\n", + ")\n", + "\n", + "ds = xr.open_mfdataset(f\"{example_dataset_folder}/*.nc\", combine=\"by_coords\")\n", + "\n", + "# TODO check how we can get good performance without loading full dataset in memory\n", + "ds.load() # load the dataset into memory\n", + "\n", + "fieldset = parcels.FieldSet.from_copernicusmarine(ds)" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "Now, define a custom Particle class that includes temperature, salinity, and density as variables, and create a ParticleSet with one particle at a known location:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "GSWParticle = parcels.Particle.add_variable(\n", + " [\n", + " parcels.Variable(\"temp\", dtype=np.float32, initial=np.nan),\n", + " parcels.Variable(\"salt\", dtype=np.float32, initial=np.nan),\n", + " parcels.Variable(\"density\", dtype=np.float32, initial=np.nan),\n", + " ]\n", + ")\n", + "\n", + "# Initiate one Argo float in the Agulhas Current\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset,\n", + " pclass=GSWParticle,\n", + " lon=[32],\n", + " lat=[-31],\n", + " z=[200],\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "Now (as the core part of this tutorial) define a custom kernel that uses the `gsw` toolbox to compute density from temperature and salinity:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "def ParcelsGSW(particles, fieldset):\n", + " particles.temp = fieldset.thetao[particles]\n", + " particles.salt = fieldset.so[particles]\n", + " pressure = gsw.p_from_z(-particles.z, particles.lat)\n", + " particles.density = gsw.density.rho(particles.salt, particles.temp, pressure)" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "Finally, run the `ParcelsGSW` Kernel for one timestep and check (for Continuous Integration purposes) that the computed density is as expected:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "pset.execute(ParcelsGSW, runtime=np.timedelta64(1, \"s\"), dt=np.timedelta64(1, \"s\"))\n", + "\n", + "np.testing.assert_allclose(pset.density, [1026.8281], rtol=1e-5)\n", + "\n", + "print(\n", + " f\"Temperature: {pset.temp[0]:.2f}, Salinity: {pset.salt[0]:.2f}, Density: {pset.density[0]:.2f}\"\n", + ")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "parcels", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pixi.toml b/pixi.toml index 672df8b8be..e78ed4880c 100644 --- a/pixi.toml +++ b/pixi.toml @@ -77,6 +77,7 @@ tests-notebooks = "pytest --nbval-lax -k 'argo' docs/examples" jupyter = "*" trajan = "*" matplotlib-base = ">=2.0.2" +gsw = "*" [feature.docs.dependencies] numpydoc = "!=1.9.0" diff --git a/src/parcels/_tutorial.py b/src/parcels/_tutorial.py index 99ceb85914..34409312cc 100644 --- a/src/parcels/_tutorial.py +++ b/src/parcels/_tutorial.py @@ -47,6 +47,7 @@ ], "CopernicusMarine_data_for_Argo_tutorial": [ "cmems_mod_glo_phy-cur_anfc_0.083deg_P1D-m_uo-vo_31.00E-33.00E_33.00S-30.00S_0.49-2225.08m_2024-01-01-2024-02-01.nc", + "cmems_mod_glo_phy-so_anfc_0.083deg_P1D-m_so_31.00E-33.00E_33.00S-30.00S_0.49-2225.08m_2024-01-01-2024-02-01.nc", "cmems_mod_glo_phy-thetao_anfc_0.083deg_P1D-m_thetao_31.00E-33.00E_33.00S-30.00S_0.49-2225.08m_2024-01-01-2024-02-01.nc", ], "DecayingMovingEddy_data": [ diff --git a/src/parcels/kernels/EOSseawaterproperties.py b/src/parcels/kernels/EOSseawaterproperties.py deleted file mode 100644 index 962b945ab4..0000000000 --- a/src/parcels/kernels/EOSseawaterproperties.py +++ /dev/null @@ -1,382 +0,0 @@ -"""Collection of pre-built eos sea water property kernels.""" - -import math - -__all__ = ["AdiabticTemperatureGradient", "PressureFromLatDepth", "PtempFromTemp", "TempFromPtemp", "UNESCODensity"] - - -def PressureFromLatDepth(particle, fieldset, time): # pragma: no cover - """ - Calculates pressure in dbars from depth in meters and latitude. - - Returns - ------- - p : array_like - pressure [db] - - References - ---------- - 1. Saunders, Peter M., 1981: Practical Conversion of Pressure to Depth. - J. Phys. Oceanogr., 11, 573-574. - doi: 10.1175/1520-0485(1981)011<0573:PCOPTD>2.0.CO;2 - """ - # Angle conversions. - deg2rad = math.pi / 180.0 - - X = math.sin(max(particle.lat * deg2rad, -1 * particle.lat * deg2rad)) - C1 = 5.92e-3 + math.pow(X, 2) * 5.25e-3 - particle.pressure = ((1 - C1) - math.pow(((math.pow((1 - C1), 2)) - (8.84e-6 * particle.depth)), 0.5)) / 4.42e-6 - - -def AdiabticTemperatureGradient(particle, fieldset, time): # pragma: no cover - """Calculates adiabatic temperature gradient as per UNESCO 1983 routines. - - - Parameters - ---------- - particle : - The particle object with the following attributes: - - S : array_like - Salinity in psu (PSS-78). - - T : array_like - Temperature in ℃ (ITS-90). - - pressure : array_like - Pressure in db. - fieldset : - The fieldset object - - Returns - ------- - array_like - Adiabatic temperature gradient in ℃ db⁻¹. - - - References - ---------- - 1. Fofonoff, P. and Millard, R.C. Jr UNESCO 1983. Algorithms for - computation of fundamental properties of seawater. UNESCO Tech. Pap. in - Mar. Sci., No. 44, 53 pp. - http://unesdoc.unesco.org/images/0005/000598/059832eb.pdf - - 2. Bryden, H. 1973. New Polynomials for thermal expansion, adiabatic - temperature gradient and potential temperature of sea water. Deep-Sea - Res. Vol20,401-408. doi:10.1016/0011-7471(73)90063-6 - - """ - s, t, pres = particle.S, particle.T, particle.pressure - - T68 = t * 1.00024 - - a = [3.5803e-5, 8.5258e-6, -6.836e-8, 6.6228e-10] - b = [1.8932e-6, -4.2393e-8] - c = [1.8741e-8, -6.7795e-10, 8.733e-12, -5.4481e-14] - d = [-1.1351e-10, 2.7759e-12] - e = [-4.6206e-13, 1.8676e-14, -2.1687e-16] - particle.adtg = ( - a[0] - + (a[1] + (a[2] + a[3] * T68) * T68) * T68 - + (b[0] + b[1] * T68) * (s - 35) - + ((c[0] + (c[1] + (c[2] + c[3] * T68) * T68) * T68) + (d[0] + d[1] * T68) * (s - 35)) * pres - + (e[0] + (e[1] + e[2] * T68) * T68) * pres * pres - ) - - -def PtempFromTemp(particle, fieldset, time): # pragma: no cover - """ - Calculates potential temperature as per UNESCO 1983 report. - - Parameters - ---------- - particle : - The particle object with the following attributes: - - S : array_like - Salinity in psu (PSS-78). - - T : array_like - Temperature in ℃ (ITS-90). - - pressure : array_like - Pressure in db. - fieldset : - The fieldset object with the following attributes: - - refpressure : array_like, optional - Reference pressure in db (default is 0). - time : float - Simulation time (not used in this function but required for consistency with other kernels). - - Returns - ------- - array_like - Potential temperature relative to reference pressure in ℃ (ITS-90). - - - References - ---------- - 1. Fofonoff, P. and Millard, R.C. Jr UNESCO 1983. Algorithms for - computation of fundamental properties of seawater. UNESCO Tech. Pap. in - Mar. Sci., No. 44, 53 pp. Eqn.(31) p.39. - http://unesdoc.unesco.org/images/0005/000598/059832eb.pdf - - 2. Bryden, H. 1973. New Polynomials for thermal expansion, adiabatic - temperature gradient and potential temperature of sea water. Deep-Sea - Res. Vol20,401-408. doi:10.1016/0011-7471(73)90063-6 - - """ - s = fieldset.psu_salinity[time, particle.depth, particle.lat, particle.lon] - t = fieldset.temperature[time, particle.depth, particle.lat, particle.lon] - pres, pr = particle.pressure, fieldset.refpressure - - # First calculate the adiabatic temperature gradient adtg - # Convert ITS-90 temperature to IPTS-68 - T68 = t * 1.00024 - - a = [3.5803e-5, 8.5258e-6, -6.836e-8, 6.6228e-10] - b = [1.8932e-6, -4.2393e-8] - c = [1.8741e-8, -6.7795e-10, 8.733e-12, -5.4481e-14] - d = [-1.1351e-10, 2.7759e-12] - e = [-4.6206e-13, 1.8676e-14, -2.1687e-16] - adtg = ( - a[0] - + (a[1] + (a[2] + a[3] * T68) * T68) * T68 - + (b[0] + b[1] * T68) * (s - 35) - + ((c[0] + (c[1] + (c[2] + c[3] * T68) * T68) * T68) + (d[0] + d[1] * T68) * (s - 35)) * pres - + (e[0] + (e[1] + e[2] * T68) * T68) * pres * pres - ) - - # Theta1. - del_P = pr - pres - del_th = del_P * adtg - th = T68 + 0.5 * del_th - q = del_th - - pprime = pres + 0.5 * del_P - adtg = ( - a[0] - + (a[1] + (a[2] + a[3] * th) * th) * th - + (b[0] + b[1] * th) * (s - 35) - + ((c[0] + (c[1] + (c[2] + c[3] * th) * th) * th) + (d[0] + d[1] * th) * (s - 35)) * pprime - + (e[0] + (e[1] + e[2] * th) * th) * pprime * pprime - ) - - # Theta2. - del_th = del_P * adtg - th = th + (1 - 1 / 2**0.5) * (del_th - q) - q = (2 - 2**0.5) * del_th + (-2 + 3 / 2**0.5) * q - - # Theta3. - adtg = ( - a[0] - + (a[1] + (a[2] + a[3] * th) * th) * th - + (b[0] + b[1] * th) * (s - 35) - + ((c[0] + (c[1] + (c[2] + c[3] * th) * th) * th) + (d[0] + d[1] * th) * (s - 35)) * pprime - + (e[0] + (e[1] + e[2] * th) * th) * pprime * pprime - ) - - del_th = del_P * adtg - th = th + (1 + 1 / 2**0.5) * (del_th - q) - q = (2 + 2**0.5) * del_th + (-2 - 3 / 2**0.5) * q - - # Theta4. - pprime = pres + del_P - adtg = ( - a[0] - + (a[1] + (a[2] + a[3] * th) * th) * th - + (b[0] + b[1] * th) * (s - 35) - + ((c[0] + (c[1] + (c[2] + c[3] * th) * th) * th) + (d[0] + d[1] * th) * (s - 35)) * pprime - + (e[0] + (e[1] + e[2] * th) * th) * pprime * pprime - ) - - del_th = del_P * adtg - particle.potemp = (th + (del_th - 2 * q) / 6) / 1.00024 - - -def TempFromPtemp(particle, fieldset, time): # pragma: no cover - """ - Calculates temperature from potential temperature at the reference - pressure PR and in situ pressure P. - - Parameters - ---------- - particle : - The particle object with the following attributes: - - S : array_like - Salinity in psu (PSS-78). - - T : array_like - Potential temperature in ℃ (ITS-90). - - pressure : array_like - Pressure in db. - fieldset : - The fieldset object with the following attributes: - - refpressure : array_like, optional - Reference pressure in db (default is 0). - time : float - Simulation time (not used in this function but required for consistency with other kernels). - - Returns - ------- - array_like - Temperature in ℃ (ITS-90). - - References - ---------- - 1. Fofonoff, P. and Millard, R.C. Jr UNESCO 1983. Algorithms for - computation of fundamental properties of seawater. UNESCO Tech. Pap. in - Mar. Sci., No. 44, 53 pp. Eqn.(31) p.39. - http://unesdoc.unesco.org/images/0005/000598/059832eb.pdf - - 2. Bryden, H. 1973. New Polynomials for thermal expansion, adiabatic - temperature gradient and potential temperature of sea water. Deep-Sea - Res. Vol20,401-408. doi:10.1016/0011-7471(73)90063-6 - - """ - s = fieldset.psu_salinity[time, particle.depth, particle.lat, particle.lon] - t = fieldset.potemperature[time, particle.depth, particle.lat, particle.lon] - pres, pr = fieldset.refpressure, particle.pressure # The order should be switched here - - # Convert ITS-90 temperature to IPTS-68 - T68 = t * 1.00024 - - a = [3.5803e-5, 8.5258e-6, -6.836e-8, 6.6228e-10] - b = [1.8932e-6, -4.2393e-8] - c = [1.8741e-8, -6.7795e-10, 8.733e-12, -5.4481e-14] - d = [-1.1351e-10, 2.7759e-12] - e = [-4.6206e-13, 1.8676e-14, -2.1687e-16] - adtg = ( - a[0] - + (a[1] + (a[2] + a[3] * T68) * T68) * T68 - + (b[0] + b[1] * T68) * (s - 35) - + ((c[0] + (c[1] + (c[2] + c[3] * T68) * T68) * T68) + (d[0] + d[1] * T68) * (s - 35)) * pres - + (e[0] + (e[1] + e[2] * T68) * T68) * pres * pres - ) - - # Theta1. - del_P = pr - pres - del_th = del_P * adtg - th = T68 + 0.5 * del_th - q = del_th - - pprime = pres + 0.5 * del_P - adtg = ( - a[0] - + (a[1] + (a[2] + a[3] * th) * th) * th - + (b[0] + b[1] * th) * (s - 35) - + ((c[0] + (c[1] + (c[2] + c[3] * th) * th) * th) + (d[0] + d[1] * th) * (s - 35)) * pprime - + (e[0] + (e[1] + e[2] * th) * th) * pprime * pprime - ) - - # Theta2. - del_th = del_P * adtg - th = th + (1 - 1 / 2**0.5) * (del_th - q) - q = (2 - 2**0.5) * del_th + (-2 + 3 / 2**0.5) * q - - # Theta3. - adtg = ( - a[0] - + (a[1] + (a[2] + a[3] * th) * th) * th - + (b[0] + b[1] * th) * (s - 35) - + ((c[0] + (c[1] + (c[2] + c[3] * th) * th) * th) + (d[0] + d[1] * th) * (s - 35)) * pprime - + (e[0] + (e[1] + e[2] * th) * th) * pprime * pprime - ) - - del_th = del_P * adtg - th = th + (1 + 1 / 2**0.5) * (del_th - q) - q = (2 + 2**0.5) * del_th + (-2 - 3 / 2**0.5) * q - - # Theta4. - pprime = pres + del_P - adtg = ( - a[0] - + (a[1] + (a[2] + a[3] * th) * th) * th - + (b[0] + b[1] * th) * (s - 35) - + ((c[0] + (c[1] + (c[2] + c[3] * th) * th) * th) + (d[0] + d[1] * th) * (s - 35)) * pprime - + (e[0] + (e[1] + e[2] * th) * th) * pprime * pprime - ) - - del_th = del_P * adtg - - particle.temp = (th + (del_th - 2 * q) / 6) / 1.00024 - - -def UNESCODensity(particle, fieldset, time): # pragma: no cover - # This is a kernel which calculates the UNESCO density - # (https://link.springer.com/content/pdf/bbm%3A978-3-319-18908-6%2F1.pdf), - # from pressure, temperature and salinity. - # density in [kg/m3] if temperature in degrees C, salinity in PSU, - # pressure in bar. - - a0 = 999.842594 - a1 = 0.06793953 - a2 = -0.009095290 - a3 = 0.0001001685 - a4 = -0.000001120083 - a5 = 0.000000006536332 - - S = fieldset.psu_salinity[time, particle.depth, particle.lat, particle.lon] # salinity - T = fieldset.cons_temperature[time, particle.depth, particle.lat, particle.lon] # temperature - P = fieldset.cons_pressure[time, particle.depth, particle.lat, particle.lon] # pressure - - rsmow = a0 + a1 * T + a2 * math.pow(T, 2) + a3 * math.pow(T, 3) + a4 * math.pow(T, 4) + a5 * math.pow(T, 5) - - b0 = 0.82449 - b1 = -0.0040899 - b2 = 0.000076438 - b3 = -0.00000082467 - b_four = 0.0000000053875 - - c0 = -0.0057246 - c1 = 0.00010227 - c2 = -0.0000016546 - - d0 = 0.00048314 - - B1 = b0 + b1 * T + b2 * math.pow(T, 2) + b3 * math.pow(T, 3) + b_four * math.pow(T, 4) - C1 = c0 + c1 * T + c2 * math.pow(T, 2) - - rho_st0 = rsmow + B1 * S + C1 * math.pow(S, 1.5) + d0 * math.pow(S, 2) - - e0 = 19652.21 - e1 = 148.4206 - e2 = -2.327105 - e3 = 0.01360477 - e4 = -0.00005155288 - - f0 = 54.6746 - f1 = -0.603459 - f2 = 0.01099870 - f3 = -0.00006167 - - g0 = 0.07944 - g1 = 0.016483 - g2 = -0.00053009 - - Kw = e0 + e1 * T + e2 * math.pow(T, 2) + e3 * math.pow(T, 3) + e4 * math.pow(T, 4) - F1 = f0 + f1 * T + f2 * math.pow(T, 2) + f3 * math.pow(T, 3) - G1 = g0 + g1 * T + g2 * math.pow(T, 2) - - K_ST0 = Kw + F1 * S + G1 * math.pow(S, 1.5) - - h0 = 3.2399 - h1 = 0.00143713 - h2 = 0.000116092 - h3 = -0.000000577905 - - i0 = 0.0022838 - i1 = -0.000010981 - i2 = -0.0000016078 - - j0 = 0.000191075 - - k0 = 0.0000850935 - k1 = -0.00000612293 - k2 = 0.000000052787 - - m0 = -0.00000099348 - m1 = 0.000000020816 - m2 = 0.00000000091697 - - Aw = h0 + h1 * T + h2 * math.pow(T, 2) + h3 * math.pow(T, 3) - A1 = Aw + (i0 + i1 * T + i2 * math.pow(T, 2)) * S + j0 * math.pow(S, 1.5) - Bw = k0 + k1 * T + k2 * math.pow(T, 2) - B2 = Bw + (m0 + m1 * T + m2 * math.pow(T, 2)) * S - - K_STp = K_ST0 + A1 * P + B2 * math.pow(T, 2) - - particle.density = rho_st0 / (1 - (P / K_STp)) diff --git a/src/parcels/kernels/TEOSseawaterdensity.py b/src/parcels/kernels/TEOSseawaterdensity.py deleted file mode 100644 index 2556326716..0000000000 --- a/src/parcels/kernels/TEOSseawaterdensity.py +++ /dev/null @@ -1,105 +0,0 @@ -"""Collection of pre-built sea water density kernels.""" - -import math - -__all__ = ["PolyTEOS10_bsq"] - - -def PolyTEOS10_bsq(particle, fieldset, time): # pragma: no cover - """Calculates density based on the polyTEOS10-bsq algorithm from Appendix A.1 and A.2 of - https://www.sciencedirect.com/science/article/pii/S1463500315000566 - requires fieldset.abs_salinity and fieldset.cons_temperature Fields in the fieldset - and a particle.density Variable in the ParticleSet - - References - ---------- - 1. Roquet, F., Madec, G., McDougall, T. J., Barker, P. M., 2014: Accurate - polynomial expressions for the density and specific volume of - seawater using the TEOS-10 standard. Ocean Modelling. - - 2. McDougall, T. J., D. R. Jackett, D. G. Wright and R. Feistel, 2003: - Accurate and computationally efficient algorithms for potential - temperature and density of seawater. Journal of Atmospheric and - Oceanic Technology, 20, 730-741. - - """ - Z = -math.fabs(particle.depth) # Z needs to be negative - SA = fieldset.abs_salinity[time, particle.depth, particle.lat, particle.lon] - CT = fieldset.cons_temperature[time, particle.depth, particle.lat, particle.lon] - - SAu = 40.0 * 35.16504 / 35.0 - CTu = 40.0 - Zu = 1e4 - deltaS = 32.0 - - zz = -Z / Zu - R00 = 4.6494977072e01 - R01 = -5.2099962525e00 - R02 = 2.2601900708e-01 - R03 = 6.4326772569e-02 - R04 = 1.5616995503e-02 - R05 = -1.7243708991e-03 - r0 = (((((R05 * zz + R04) * zz + R03) * zz + R02) * zz + R01) * zz + R00) * zz - - R000 = 8.0189615746e02 - R100 = 8.6672408165e02 - R200 = -1.7864682637e03 - R300 = 2.0375295546e03 - R400 = -1.2849161071e03 - R500 = 4.3227585684e02 - R600 = -6.0579916612e01 - R010 = 2.6010145068e01 - R110 = -6.5281885265e01 - R210 = 8.1770425108e01 - R310 = -5.6888046321e01 - R410 = 1.7681814114e01 - R510 = -1.9193502195e00 - R020 = -3.7074170417e01 - R120 = 6.1548258127e01 - R220 = -6.0362551501e01 - R320 = 2.9130021253e01 - R420 = -5.4723692739e00 - R030 = 2.1661789529e01 - R130 = -3.3449108469e01 - R230 = 1.9717078466e01 - R330 = -3.1742946532e00 - R040 = -8.3627885467e00 - R140 = 1.1311538584e01 - R240 = -5.3563304045e00 - R050 = 5.4048723791e-01 - R150 = 4.8169980163e-01 - R060 = -1.9083568888e-01 - R001 = 1.9681925209e01 - R101 = -4.2549998214e01 - R201 = 5.0774768218e01 - R301 = -3.0938076334e01 - R401 = 6.6051753097e00 - R011 = -1.3336301113e01 - R111 = -4.4870114575e00 - R211 = 5.0042598061e00 - R311 = -6.5399043664e-01 - R021 = 6.7080479603e00 - R121 = 3.5063081279e00 - R221 = -1.8795372996e00 - R031 = -2.4649669534e00 - R131 = -5.5077101279e-01 - R041 = 5.5927935970e-01 - R002 = 2.0660924175e00 - R102 = -4.9527603989e00 - R202 = 2.5019633244e00 - R012 = 2.0564311499e00 - R112 = -2.1311365518e-01 - R022 = -1.2419983026e00 - R003 = -2.3342758797e-02 - R103 = -1.8507636718e-02 - R013 = 3.7969820455e-01 - ss = math.sqrt((SA + deltaS) / SAu) - tt = CT / CTu - zz = -Z / Zu - rz3 = R013 * tt + R103 * ss + R003 - rz2 = (R022 * tt + R112 * ss + R012) * tt + (R202 * ss + R102) * ss + R002 - rz1 = (((R041 * tt + R131 * ss + R031) * tt + (R221 * ss + R121) * ss + R021) * tt + ((R311 * ss + R211) * ss + R111) * ss + R011) * tt + (((R401 * ss + R301) * ss + R201) * ss + R101) * ss + R001 # fmt: skip - rz0 = (((((R060 * tt + R150 * ss + R050) * tt + (R240 * ss + R140) * ss + R040) * tt + ((R330 * ss + R230) * ss + R130) * ss + R030) * tt + (((R420 * ss + R320) * ss + R220) * ss + R120) * ss + R020) * tt + ((((R510 * ss + R410) * ss + R310) * ss + R210) * ss + R110) * ss + R010) * tt + (((((R600 * ss + R500) * ss + R400) * ss + R300) * ss + R200) * ss + R100) * ss + R000 # fmt: skip - r = ((rz3 * zz + rz2) * zz + rz1) * zz + rz0 - - particle.density = r0 + r diff --git a/tests-v3/test_kernel_language.py b/tests-v3/test_kernel_language.py index 76b3ad1ba0..b75397e92d 100644 --- a/tests-v3/test_kernel_language.py +++ b/tests-v3/test_kernel_language.py @@ -6,19 +6,11 @@ from parcels import ( Field, - FieldSet, Kernel, Particle, ParticleSet, Variable, ) -from parcels.application_kernels.EOSseawaterproperties import ( - PressureFromLatDepth, - PtempFromTemp, - TempFromPtemp, - UNESCODensity, -) -from parcels.application_kernels.TEOSseawaterdensity import PolyTEOS10_bsq from tests.common_kernels import DoNothing from tests.utils import create_fieldset_unit_mesh @@ -307,96 +299,3 @@ def test_small_dt(dt, expectation): with expectation: pset.execute(DoNothing, dt=dt, runtime=dt * 101) assert np.allclose([p.time for p in pset], dt * 100) - - -def test_TEOSdensity_kernels(): - def generate_fieldset(xdim=2, ydim=2, zdim=2, tdim=1): - lon = np.linspace(0.0, 10.0, xdim, dtype=np.float32) - lat = np.linspace(0.0, 10.0, ydim, dtype=np.float32) - depth = np.linspace(0, 2000, zdim, dtype=np.float32) - time = np.zeros(tdim, dtype=np.float64) - U = np.ones((tdim, zdim, ydim, xdim)) - V = np.ones((tdim, zdim, ydim, xdim)) - abs_salinity = 30 * np.ones((tdim, zdim, ydim, xdim)) - cons_temperature = 10 * np.ones((tdim, zdim, ydim, xdim)) - dimensions = {"lat": lat, "lon": lon, "depth": depth, "time": time} - data = { - "U": np.array(U, dtype=np.float32), - "V": np.array(V, dtype=np.float32), - "abs_salinity": np.array(abs_salinity, dtype=np.float32), - "cons_temperature": np.array(cons_temperature, dtype=np.float32), - } - return (data, dimensions) - - data, dimensions = generate_fieldset() - fieldset = FieldSet.from_data(data, dimensions) - - DensParticle = Particle.add_variable("density", dtype=np.float32) - - pset = ParticleSet(fieldset, pclass=DensParticle, lon=5, lat=5, depth=1000) - - pset.execute(PolyTEOS10_bsq, runtime=1) - assert np.allclose(pset[0].density, 1027.45140) - - -def test_EOSseawaterproperties_kernels(): - fieldset = FieldSet.from_data( - data={"U": 0, "V": 0, "psu_salinity": 40, "temperature": 40, "potemperature": 36.89073}, - dimensions={"lat": 0, "lon": 0, "depth": 0}, - ) - fieldset.add_constant("refpressure", float(0)) - - PoTempParticle = Particle.add_variables( - [Variable("potemp", dtype=np.float32), Variable("pressure", dtype=np.float32, initial=10000)] - ) - pset = ParticleSet(fieldset, pclass=PoTempParticle, lon=5, lat=5, depth=1000) - pset.execute(PtempFromTemp, runtime=1) - assert np.allclose(pset[0].potemp, 36.89073) - - TempParticle = Particle.add_variables( - [Variable("temp", dtype=np.float32), Variable("pressure", dtype=np.float32, initial=10000)] - ) - pset = ParticleSet(fieldset, pclass=TempParticle, lon=5, lat=5, depth=1000) - pset.execute(TempFromPtemp, runtime=1) - assert np.allclose(pset[0].temp, 40) - - pset = ParticleSet(fieldset, pclass=TempParticle, lon=5, lat=30, depth=7321.45) - pset.execute(PressureFromLatDepth, runtime=1) - assert np.allclose(pset[0].pressure, 7500, atol=1e-2) - - -@pytest.mark.parametrize("pressure", [0, 10]) -def test_UNESCOdensity_kernel(pressure): - def generate_fieldset(p, xdim=2, ydim=2, zdim=2, tdim=1): - lon = np.linspace(0.0, 10.0, xdim, dtype=np.float32) - lat = np.linspace(0.0, 10.0, ydim, dtype=np.float32) - depth = np.linspace(0, 2000, zdim, dtype=np.float32) - time = np.zeros(tdim, dtype=np.float64) - U = np.ones((tdim, zdim, ydim, xdim)) - V = np.ones((tdim, zdim, ydim, xdim)) - psu_salinity = 8 * np.ones((tdim, zdim, ydim, xdim)) - cons_temperature = 10 * np.ones((tdim, zdim, ydim, xdim)) - cons_pressure = p * np.ones((tdim, zdim, ydim, xdim)) - dimensions = {"lat": lat, "lon": lon, "depth": depth, "time": time} - data = { - "U": np.array(U, dtype=np.float32), - "V": np.array(V, dtype=np.float32), - "psu_salinity": np.array(psu_salinity, dtype=np.float32), - "cons_pressure": np.array(cons_pressure, dtype=np.float32), - "cons_temperature": np.array(cons_temperature, dtype=np.float32), - } - return (data, dimensions) - - data, dimensions = generate_fieldset(pressure) - fieldset = FieldSet.from_data(data, dimensions) - - DensParticle = Particle.add_variable("density", dtype=np.float32) - - pset = ParticleSet(fieldset, pclass=DensParticle, lon=5, lat=5, depth=1000) - - pset.execute(UNESCODensity, runtime=1) - - if pressure == 0: - assert np.allclose(pset[0].density, 1005.9465) - elif pressure == 10: - assert np.allclose(pset[0].density, 1006.4179)