diff --git a/Common/derived_types.f90 b/Common/derived_types.f90 index 153e8d15eb58a01dd48b077ac9a2edfe22e35d07..b7d72fc5982e9b0b07a16694c20d278a128e473f 100644 --- a/Common/derived_types.f90 +++ b/Common/derived_types.f90 @@ -51,14 +51,14 @@ module derived_types end type eddy type snapshot - type(eddy), allocatable:: list(:) ! (number_extr) Visible - ! "eddies" at a given date. In general, there may be both - ! cyclonic and anticyclonic eddies in the list. Visible eddies - ! include "eddies" without an outermost contour, that is, only - ! the extremum is visible. The subscript value in list is the - ! "identification number" of the extremum. - - integer number_extr ! number of visible extrema + type(eddy), allocatable:: list(:) ! (number_extr) "Eddies" at a + ! given date. In general, there may be both cyclonic and + ! anticyclonic eddies in the list. Eddies include "eddies" + ! without an outermost contour, that is, only the extremum is + ! defined. The subscript value in list is the "identification + ! number" of the extremum. + + integer number_extr ! number of extrema integer, allocatable:: extr_map(:, :) ! (1 - copy:nlon + copy, nlat) if the grid is periodic in diff --git a/Common/util_eddies.py b/Common/util_eddies.py index 4beb0d64a4b180ccf6c319d9f42e2bfbd02ef40d..13a97424f725f6b7311ffd737abbf1a7d47a2c4c 100644 --- a/Common/util_eddies.py +++ b/Common/util_eddies.py @@ -230,10 +230,11 @@ class SHPC_class: return self._i_slice[date_index - self.d_init[0]] def comp_ishape(self, date, eddy_index, orientation): - """Compute the location in the shapefiles: return `(i_slice, ishape)`. + """Compute the location in the shapefiles, given the date index and + the intra-date eddy index: return `(i_slice, ishape)`. - Returns None if ishape_last is None or eddy_index is greater - than the maximum value. + Returns None for ishape if ishape_last is None or eddy_index + is greater than the maximum value. """ diff --git a/Inst_eddies/Documentation_texfol/documentation.tex b/Inst_eddies/Documentation_texfol/documentation.tex index 75824fd481812c49b6f0657a03616cf06127d173..5a87b9db98efe07fb9056717a59b1dc10b18e4e6 100644 --- a/Inst_eddies/Documentation_texfol/documentation.tex +++ b/Inst_eddies/Documentation_texfol/documentation.tex @@ -737,59 +737,52 @@ contour de vitesse maximale, 720 latitudes par 1440 longitudes. \item[longitude, latitude] Vecteurs de réels \item[ssh, u, v] Tableaux de réels, deux dimensions (nlon, nlat). Environ 4 MiB chacun. -\item[eddy] La comparaison entre \verb+out_cont+\%ssh et ssh\_extr donne la - distinction anticyclone -- cyclone mais il est quand même plus - confortable d'avoir cette information directement dans les - sorties. Pour cela, j'ajoute un champ cyclone. Mémoire utilisée par - un scalaire de ce type : +\item[eddy] La comparaison entre \verb+out_cont+\%ssh et extr\%ssh + donne la distinction anticyclone -- cyclone mais il est quand même + plus confortable d'avoir cette information directement dans les + sorties. Pour cela, j'ajoute une composante cyclone. Mémoire + utilisée par un scalaire de ce type : \begin{equation*} \mathrm{taille}(\mathtt{eddy}) - = 9\ \mathrm{mots} + \mathrm{taille}(\mathtt{out\_cont}) + = 10\ \mathrm{mots} + \mathrm{taille}(\mathtt{out\_cont}) + \mathrm{taille}(\mathtt{speed\_cont}) \end{equation*} Notons $n_o$ le nombre de points dans le contour extérieur (out\_cont\%n\_points) et $n_s$ le nombre de points dans le contour de vitesse maximale (speed\_cont\%n\_points). On obtient : - \begin{align*} - \mathrm{taille}(\mathtt{eddy}) - & = [2 (n_o + n_s) + 17]\ \mathrm{mots} \\ - & = 4 [2 (n_o + n_s) + 17]\ \mathrm{B} - \end{align*} -soit environ \np{0.2} KiB. + \begin{equation*} + \mathrm{taille}(\mathtt{eddy}) = 2 (n_o + n_s + 9)\ \mathrm{mots} + \end{equation*} + soit $8 (n_o + n_s + 9)\ \mathrm{B}$, soit environ \np{0.2} KiB en + moyenne. \item[snapshot\%list\_vis] On a besoin d'encapsuler le vecteur de type eddy dans un scalaire de type snapshot pour pouvoir considérer un - vecteur de snapshots. Notons $n_{ve}$ le nombre de tourbillons à une - date donnée (number\_vis\_extr). Mémoire utilisée par ce champ : + vecteur de snapshots. Notons $n_e$ le nombre d'extremums à une + date donnée (snapshot\%number\_extr). Mémoire utilisée par ce champ : \begin{equation*} - n_{ve} \langle \mathrm{taille}(\mathtt{eddy}) \rangle - = n_{ve} - [2 (\langle n_o \rangle + \langle n_s \rangle) + 17]\ \mathrm{mots} + n_e \langle \mathrm{taille}(\mathtt{eddy}) \rangle + = 2 n_e (\langle n_o \rangle + \langle n_s \rangle + 9)\ \mathrm{mots} \end{equation*} - soit environ 3 MiB. + soit environ 3 MiB en moyenne. \item[snapshot\%extr\_map] Notons $n_\lambda$ le nombre de longitudes, $n_\phi$ le nombre de latitudes. Taille : $n_\lambda n_\phi$ mots, soit environ 4 MiB. -\item[snapshot\%ind\_extr] Taille : $2 n_{ve}$ mots, soit environ - \np{0.1} MiB. \end{description} On a : \begin{align*} \mathrm{taille}(\mathtt{snapshot}) - & = n_{ve} \langle \mathrm{taille}(\mathtt{eddy}) \rangle - + (n_\lambda n_\phi + 2 n_{ve} + 2)\ \mathrm{mots} \\ - & = [n_{ve} (2 (\langle n_o \rangle + \langle n_s \rangle) + 19) - + n_\lambda n_\phi + 2]\ \mathrm{mots} + & = n_e \langle \mathrm{taille}(\mathtt{eddy}) \rangle + + (n_\lambda n_\phi + 1)\ \mathrm{mots} \\ + & = [2 n_e (\langle n_o \rangle + \langle n_s \rangle + 9) + + n_\lambda n_\phi + 1]\ \mathrm{mots} \end{align*} D'où la taille d'un scalaire de type snapshot : environ 7 MiB. -J'avais d'abord pensé mettre dans le type eddy les champs ind\_lon et -ind\_lat. Ce n'est pas très cohérent parce que les indices se réfèrent -au maillage global, alors que les autres champs sont des propriétés -intrinsèques. Ainsi, j'étais embêté pour définir un tourbillon à -partir d'une procédure travaillant sur une section du maillage -global. J'ai donc choisi donc de remplacer ces champs par un champ -ind\_extr dans snapshot. Et j'ajoute à la place des champs longitude, -latitude dans le type eddy. +J'ai choisi de mettre dans le type eddy les indices de l'extremum, +bien que les indices se réfèrent au maillage global, alors que les +autres champs sont des propriétés intrinsèques. Il faut donc faire +attention pour définir un tourbillon à partir d'une procédure +travaillant sur une section du maillage global. Cf. figure (\ref{fig:copy}). \begin{figure} @@ -966,7 +959,7 @@ dans le programme pour Matlab : \begin{algorithmic} \STATE appel de random\_number(level\_try) - \STATE level\_try = e\%ssh\_extr $\pm$ [4 + (2 * level\_try - + \STATE level\_try = e\%extr\%ssh $\pm$ [4 + (2 * level\_try - 1) * precision] \COMMENT{un peu de bruit sur l'amplitude, autour de 4} \end{algorithmic} @@ -1776,16 +1769,17 @@ Dans les données Aviso, les champs SSH, u et v sont sur une même grille rectangulaire longitude, latitude. Dans les sorties du modèle utilisé par Jonathan Gula, le champ $\zeta$ et les composantes de la vitesse dans le plan de projection sont chacun sur une grille -propre. Les trois grilles sont rectangulaires mais non -longitude-latitude. Discrétisation Arakawa C (Arakawa 1977 k0737, § -III.A). La grille de $\zeta$ et de l'angle $\alpha$ de la vitesse -projetée est appelée grille $\rho$. La grille $\psi$ est la grille de -la vorticité. Le couple de coordonnées dans le plan de projection est -$(\xi, \eta)$. Il faut de toutes façons se ramener à la même grille -pour passer des composantes de la vitesse dans l'espace projeté aux -composantes zonale et méridienne dans l'espace physique. Le plus -simple est, au moment de la lecture, de ramener les vitesses projetées -sur la grille de $\zeta$ et de faire la rotation. +propre. Les trois grilles sont structurées mais ne sont pas +cartésiennes en longitude et latitude. Discrétisation Arakawa C +(Arakawa 1977 k0737, § III.A). La grille de $\zeta$ et de l'angle +$\alpha$ de la vitesse projetée est appelée grille $\rho$. La grille +$\psi$ est la grille de la vorticité. Le couple de coordonnées dans le +plan de projection est $(\xi, \eta)$. Il faut de toutes façons se +ramener à la même grille pour passer des composantes de la vitesse +dans l'espace projeté aux composantes zonale et méridienne dans +l'espace physique. Le plus simple est, au moment de la lecture, de +ramener les vitesses projetées sur la grille de $\zeta$ et de faire la +rotation. La position en longitude et latitude intervient dans le programme \verb+inst_eddies+ essentiellement à travers les variables corner, @@ -1796,59 +1790,20 @@ dans le calcul de la surface des contours (appels à \verb+set_max_speed+) et dans le calcul de la vitesse moyenne (procédure \verb+mean_speed+). -Donc plutôt créer une procédure \verb+find_contours_curv+ dans -\verb+Contour_531+ ? Dans \verb+find_contours_curv+, il faut faire une -interpolation bilinéaire de $\lambda$ et $\phi$ aux points d'un -contour. Comme les points d'un contour sont par construction sur les -arêtes du maillage, l'interpolation est en fait linéaire à une seule -dimension. - Le passage $(i, j) \in \mathbb{N}^2 \mapsto (\lambda, \phi)$ est ok : c'est juste un accès direct aux tableaux venant du fichier NetCDF. Le -passage $(i, j) \in \mathbb{R}^2 \mapsto (\lambda, \phi)$ est ok aussi : -c'est une interpolation bilinéaire. Le problème est le passage -$(\lambda, \phi) \mapsto (i, j) \in \mathbb{R}^2$. - -\verb+eddy_graph+ utilise \verb+ind_extr+ et \verb+extr_map+ et ne -peut pas les reconstruire facilement à partir de la longitude et de la -latitude des extremums pour une grille non longitude -- latitude. Idée -: faire passer \verb+snapshot%ind_extr+ dans eddy et l'écrire dans -\verb+extremum.dbf+. Mais inconvénient : cf. § -\ref{sec:structures}. Donc garder \verb+snapshot%ind_extr+ et l'écrire -dans un fichier indépendant (binaire Fortran à accès direct ? HDF5 ?) -ou écrire un fichier NetCDF contenant \verb+extr_map+ ? Ajouter dans -eddy une composante polyline \verb+out_cont_no_coord+. Dans -\verb+set_max_speed+, calculer \verb+mean_speed+ sur (radius4 - 1) -$\times$ 4 contours. Mettre \verb+outside_points+ vide dans -\verb+good_contour+ pour trouver ces contours supplémentaires. Faire -la recherche avec \verb+good_contour+ dans l'espace des indices. - -type ssh contour : ssh, area, speed, contour xy, contour $\lambda -\phi$. - -type eddy -\begin{itemize} -\item extremum : coord proj x y, coord $\lambda \phi$, ssh -\item innermost level, cyclone -\item type(ssh contour) contour list\verb+(:)+ (dernier : outer) -\item indice dans contour list du contour de vitesse maximale -\item valid, delta in delta out -\end{itemize} - -Algorithme principal : -\begin{algorithmic} - \STATE lire ($\lambda$, $\phi$, ssh, u, v) - \STATE appel de set\_all\_extr - \STATE appel de set\_contours -\end{algorithmic} +passage $(i, j) \in \mathbb{R}^2 \mapsto (\lambda, \phi)$ est ok aussi +: c'est une interpolation bilinéaire. Pour un contour, comme les +points d'un contour sont par construction sur les arêtes du maillage, +l'interpolation est en fait linéaire à une seule dimension. -\verb+get_1_outerm+ : dans l'espace des indices et dans l'espace -$(\lambda, \phi)$. +Le problème est le passage +$(\lambda, \phi) \mapsto (i, j) \in \mathbb{R}^2$. -\verb+mean_speed+ : faire l'interpolation dans l'espace des -indices. (La rotation de la vitesse dans le plan de projection a été -effectuée au moment de la lecture des entrées du programme principal.) -Il faut le contour extérieur dans l'espace des indices pour pouvoir -interpoler la vitesse aux points du contour extérieur. +\verb+eddy_graph+ utilise \verb+extr%coord_proj+ et \verb+extr_map+ et +ne peut pas les reconstruire facilement à partir de la longitude et de +la latitude des extremums pour une grille non longitude -- +latitude. Idée : écrire un shapefile contenant \verb+extr%coord_proj+ +(beaucoup moins de place que \verb+extr_map+). \end{document} diff --git a/Inst_eddies/mean_speed.f90 b/Inst_eddies/mean_speed.f90 index f681d037868900069fa5e4c5fae0f2e35eaf565f..5daf3fcad233efc1ba570782bd5e3c50fcfce283 100644 --- a/Inst_eddies/mean_speed.f90 +++ b/Inst_eddies/mean_speed.f90 @@ -7,11 +7,13 @@ contains real function mean_speed(u, v, p, points_xy, center, corner) ! Interpolates the wind at each point of the input polyline and - ! computes the mean azimuthal speed at interpolation points. We - ! assume (and do not check) that the longitude of the center has - ! been brought modulo 2 pi in the longitude range of the polyline - ! `p` and that the polyline is in the longitude range of the grid - ! corresponding to corner and step. + ! computes the mean azimuthal speed at interpolation + ! points. Interpolation is in projection space. + + ! We assume (and do not check) that the longitude of the center + ! has been brought modulo 2 pi in the longitude range of the + ! polyline `p` and that the polyline is in the longitude range of + ! the grid corresponding to corner and step. use, intrinsic:: iso_fortran_env, only: error_unit @@ -27,7 +29,7 @@ contains ! radians. Should be closed. real, intent(in):: points_xy(:, :) ! (2, p%n_points) - ! Points given by projection coordinates. + ! Points of the contour `p` given by projection coordinates. real, intent(in):: center(:) ! (2) ! Longitude and latitude, in rad. Azimuthal speed is computed with diff --git a/Inst_eddies/nearby_extr.f90 b/Inst_eddies/nearby_extr.f90 index c028d1f7af338f78d4f10740801a5da1754147b8..30c0aa6bb93c8373bba7b1c5960e8501b98cfc42 100644 --- a/Inst_eddies/nearby_extr.f90 +++ b/Inst_eddies/nearby_extr.f90 @@ -19,8 +19,8 @@ contains ! extremum. 0 at other points. type(eddy), intent(in):: list(:) - ! Visible eddies at a given date. We need components extr%coord, - ! cyclone and out_cont to be defined. + ! Eddies at a given date. We need components extr%coord, cyclone + ! and out_cont to be defined. integer, intent(in):: i ! identifying number of the target extremum diff --git a/Overlap/candidate_overlap.f90 b/Overlap/candidate_overlap.f90 index e94405a0eaec79d82ac91d29ebe6531df580a304..7f5bb58cb0853cdcc3bfee216399204d66b4041b 100644 --- a/Overlap/candidate_overlap.f90 +++ b/Overlap/candidate_overlap.f90 @@ -19,7 +19,7 @@ contains ! extremum. 0 at other points. type(eddy), intent(in):: list(:) - ! Visible eddies at a given date. We need component delta_in to be + ! Eddies at a given date. We need component delta_in to be ! defined. Arriving in this subroutine, list%delta_in could be <= ! delta or huge(0).