Ma tâche est de calculer le tourbillon potentiel isentropique sur une grille.Toutes les données (ens ECMWF) sont données sur une grille de pression hybride sigma.Puisque je ne trouve vraiment pas de référence, je suppose (!) que le tourbillon relatif (vo, id = 138) est évalué au niveau du modèle et est donc inutile pour moi.

Mon approche (pour chaque point de la grille):

1. Calculer le gradient thêta, en utilisant le trois dimensions du modèle, convertissant les distances verticales et horizontales en mètres.

2. Définissez deux vecteurs en utilisant chacun trois contraintes:

   • Le produit scalaire avec un gradient thêta disparaît .
   • Le premier vecteur n'a pas de méridien, le second pas de composante zonale.
   • La distance au point de grille horizontal adjacent définit la magnitude de la composante méridionale / zonale.
3. Ces deux vecteurs couvrent maintenant un plan tangentiel à la surface isentropique au point de grille courant. Ils pointent vers des positions décalées verticalement uniquement des voisins horizontaux du point de grille actuel.
4. Les valeurs de u et v sont interpolées - encore une fois, en utilisant des mètres (hauteur géopotentielle) - aux emplacements donnés par les vecteurs tangents.
5. La longueur géométrique de ces vecteurs tangents est déterminée et utilisée avec les valeurs interpolées de u et v pour calculer la vorticité relative isentropique.
6. La stabilité statique est calculée en utilisant les valeurs de pression et thêta au point courant et au voisin vertical.
7. En utilisant l'accélération gravitationnelle et le facteur de conversion de 10 ^ -6, j'utilise l'équation donnée dans le célèbre livre de Wallace & Hobbs pour Potential Vorticity: $ g (f + \ eta _ {\ theta}) \ frac {\ partial \ theta} {\ partial p} $

Mon problème:

Le résultat ne ressemble pas à ce que j'attends Les valeurs varient d'environ -30 à +50 pvu dans toute l'atmosphère, alors que la surface de 1,5 pvu ne représente clairement pas la tropopause. Cela ressemble plus à un accident. J'utilise "omp" pour augmenter la vitesse de calcul, mais les choses se ressemblent si je ne l'utilise pas aussi, donc cela ne peut pas être une confusion d'index.

Ma question:

Qu'est-ce qui ne va pas? Pouvez-vous repérer des erreurs systématiques? Peut-être que je suis devenu aveugle avec le temps. Savez-vous comment le tourbillon relatif est évalué dans le modèle ECMWF? (J'ai essayé, vraiment ...)

Pour être plus précis: Ma méthode est d'être implémentée en C ++ dans un environnement déjà existant pour visualiser les données NWP en 3D en temps réel. Les champs de données de base - donnés en hybride coordonnées sigma - sont lues, traitées, écrites dans un champ de données dérivées, interpolées et rendues sur les niveaux de pression déclenchés par l'interaction de l'utilisateur.Le programme doit posséder les méthodes de traitement à lui seul, c'est pourquoi je ne peux pas utiliser de packages Python, par exemple .

Voici le code source de mon approche. Si quelque chose n'est pas clair, veuillez demander! Peut-être que quelqu'un repèrera une erreur.

  // VORTICITÉ POTENTIELLE ISENTROPIQUE // Grilles d'entrée: 0 t, 1 q, 2 u, 3 v, 4 sfcPhi, 5 sfc t_d, 6 sfc t, 7 sfc p, 8 sfc u, 9 sfc velse if ((derivedVarName == "Isentropic PV (an)") || (derivedVarName == "Isentropic PV (fc)") || ( dérivéVarName == "Isentropic PV (ens)")) {unsigned int k, j, i; int n, l; #pragma omp parallel {#pragma omp pour private (k, j, i, n, l) for (k = 1; k < derivedGrid->getNumLevels () - 1; k ++) for (j = 0; j < derivedGrid->getNumLats () - 1; j ++) for (i = 0; i < derivedGrid->getNumLons () - 1; i ++) {// Point de grille actuel. double T000 = inputGrids.at (0) ->getValue (k, j, i); double p000 = inputGrids.at (0) ->getPressure (k, j, i); double q000 = inputGrids.at (1) ->getValue (k, j, i); double thêta000 = potTemp (T000, p000); // Voisin supérieur. double T100 = inputGrids.at (0) ->getValue (k-1, j, i); double p100 = inputGrids.at (0) ->getPressure (k-1, j, i); double q100 = inputGrids.at (1) ->getValue (k-1, j, i);
double thêta100 = potTemp (T100, p100); // voisin x. double theta001 = potTemp (inputGrids.at (0) ->getValue (k, j, i + 1), inputGrids.at (0) ->getPressure (k, j, i + 1)); // voisin y. double theta010 = potTemp (inputGrids.at (0) ->getValue (k, j + 1, i), inputGrids.at (0) ->getPressure (k, j + 1, i)); // (1) Calcule le gradient thêta au point courant, en utilisant l'épaisseur géopotentielle dZ [m] pour la verticale, // la distance du grand cercle [m] pour les distances horizontales. // Calcul de l'épaisseur géopotentielle entre le point courant et le voisin supérieur. double T_v = (virtualTempFromSpecHum (T000, q000) + virtualTempFromSpecHum (T100, q100)) / 2 .; double dPhi = geopotThickness (T_v, p100, p000); double dZ = geopotToZ (dPhi); // Calcul du gradient. double * grad = gradient (theta000, inputGrids.at (0) ->getEastInterfaceLon (i), inputGrids.at (0) ->getNorthInterfaceLat (j), // Point de grille actuel. theta001, inputGrids.at (0) -AgiliLpronp +1), // voisin xlon. Theta010, inputGrids.at (0) ->getNorthInterfaceLat (j + 1), // voisin-ylat. Theta100, dZ // voisin Z.); // Composants de dégradé. double xgrad = grad [2]; double ygrad = grad [1]; double Zgrad = grad [0]; // (2) Définissez les vecteurs tangents et calculez leurs composantes Z. // Définition: // 1) Le produit scalaire de vecteurs tangents avec gradient s'annule. // 2) Un vecteur tangent dans le plan y-z, l'autre dans le plan x-z.
// 3) La longueur des composants horizontaux est égale à l'espacement entre les points de la grille. double radius_m = MetConstants :: EARTH_RADIUS_km * 1000 .; // Vecteur tangent dans le plan y-z. double dx = gcDistance_deg (inputGrids.at (0) ->getEastInterfaceLon (i), inputGrids.at (0) ->getNorthInterfaceLat (j), inputGrids.at (0) ->getEast (0) ->getEast (0) ->getEast (entrée) >getNorthInterfaceLat (j), radius_m); // Composant Z. double dZ_x = -dx * xgrad / Zgrad; // Vecteur tangent dans le plan x-z. double dy = gcDistance_deg (inputGrids.at (0) ->getEastInterfaceLon (i), inputGrids.at (0) ->getNorthInterfaceLat (j), inputGrids.atorth (0) ->getNorthInterfaceLat (j), inputGrids. j + 1), rayon_m); // Composant Z. double dZ_y = -dy * ygrad / Zgrad; // Les vecteurs tangents pointent désormais du point de grille actuel vers des emplacements uniquement // déplacés verticalement depuis le voisin x ou y, respectivement. // (3) Recherche des niveaux dont les valeurs Z renferment des valeurs de hauteur géopotentielle et // interpole u et v à ces valeurs. // La hauteur géopotentielle du point de grille actuel est calculée dans la section // coordonnée x et réutilisée dans la section coordonnée y. // double coordonnée x vdx; // Valeur interpolée de v. Double Z_x; // Valeur d'interpolation de Z. // (3a) Calcule la hauteur géopotentielle du point de grille actuel et ajoute la composante Z // du vecteur tangent du plan x-Z.
// PV: 0 t, 1 q, 2 u, 3 v, 4 sfcPhi, 5 sfc t_d, 6 sfc t, 7 sfc p // Géopotentiel de surface. double geopotSurf_J = inputGrids.at (4) ->getValue (0, j, i); // En supposant une température comprise entre -45 et +60 degrés Celsius. // Raphido: Seul le sfc T_v est utilisé. Utilisez plutôt la valeur moyenne, également dans la méthode de hauteur géopotentielle. // Température virtuelle de surface. double pSurf_hPa = inputGrids.at (7) ->getValue (0, j, i) / 100 .; double virtualTempSurf_K = virtualTempFromDewPoint (inputGrids.at (5) ->getValue (0, j, i), pSurf_hPa, inputGrids.at (6) ->getValue (0, j, i)); // Épaisseur géopotentielle de la couche de surface. double dPhiSurfLayer = geopotThickness (virtualTempSurf_K, inputGrids.at (0) ->getPressure (inputGrids.at (0) ->getNumLevels () - 1, j, i), pSurf_hPa); double geopotential000_J = 0 .; double virtualTemp_av_K = 0 .; double dPhi_J = 0 .; // Pour chaque point, boucle du niveau de pression de surface au niveau de pression actuel (k). for (l = derivedGrid->getNumLevels () - 2; l > k-1; l--) {// Calcule la température virtuelle moyenne entre le niveau de pression actuel // et le niveau inférieur. virtualTemp_av_K = (virtualTempFromSpecHum (inputGrids.at (0) ->getValue (l, j, i), inputGrids.at (1) ->getValue (l, j, i))
+ virtualTempFromSpecHum (inputGrids.at (0) ->getValue (l + 1, j, i), inputGrids.at (1) ->getValue (l + 1, j, i))) / 2 .; // Calcule l'épaisseur actuelle. dPhi_J = geopotThickness (virtualTemp_av_K, inputGrids.at (0) ->getPressure (l, j, i), inputGrids.at (0) ->getPressure (l + 1, j, i)); // Faites la somme des différentiels. geopotential000_J + = dPhi_J; } // Ajout de l'épaisseur géopotentielle de la couche de surface geopotential000_J + = dPhiSurfLayer; // Ajout du géopotentiel de surface geopotential000_J + = geopotSurf_J; // Conversion en hauteur géopotentielle en divisant par g_0. Z_x = geopotToZ (geopotential000_J); // Ajouter une contribution de dégradé. Z_x + = dZ_x; // (3b) Parcourez les niveaux du modèle pour trouver les points de grille adjacents entourant Z_x. // À partir du niveau le plus bas, calculez la hauteur géopotentielle des points de grille adjacents. for (n = derivedGrid->getNumLevels () - 1; n > 1; n--) {// Géopotentiel de surface. double geopotSurf_J = inputGrids.at (4) ->getValue (0, j, i + 1); // En supposant une température comprise entre -45 et +60 degrés Celsius // Température virtuelle de surface. double pSurf_hPa = inputGrids.at (7) ->getValue (0, j, i + 1) / 100 .; double virtualTempSurf_K = virtualTempFromDewPoint (inputGrids.at (5) ->getValue (0, j, i + 1),
pSurf_hPa, inputGrids.at (6) ->getValue (0, j, i + 1)); // Épaisseur géopotentielle de la couche de surface double dPhiSurfLayer = geopotThickness (virtualTempSurf_K, inputGrids.at (0) ->getPressure (inputGrids.at (0) ->getNumLevels () - 1, j, i + 1), pSurf_hPa); // Voisin inférieur. double geopotential_J_lower = geopotSurf_J; double Zx_bot = geopotToZ (geopotential_J_lower); // Voisin supérieur. double geopotential_J_upper = geopotential_J_lower + dPhiSurfLayer; double Zx_top = geopotToZ (geopotential_J_upper); // Vérifie si les points de grille adjacents englobent Z_x. // Si c'est le cas, vérifiez si la surface est impliquée ou non, puis interpolez. if ((Zx_bot < Z_x) && (Z_x < Zx_top)) {// Niveau de surface impliqué. if (n == derivedGrid->getNumLevels () - 1) {double v_top = inputGrids.at (3) ->getValue (inputGrids.at (0) ->getNumLevels () - 1, j, i + 1); double v_bot = inputGrids.at (9) ->getValue (0, j, i + 1); // Valeur interpolée de v. Vdx = v_bot + (v_top - v_bot) * (Z_x - Zx_bot) / (Zx_top - Zx_bot); Pause; } // Niveau de surface non impliqué. autre
{double v_top = inputGrids.at (3) ->getValue (n-1, j, i + 1); double v_bot = inputGrids.at (3) ->getValue (n, j, i + 1); // Valeur interpolée de v. Vdx = v_bot + (v_top - v_bot) * (Z_x - Zx_bot) / (Zx_top - Zx_bot); Pause; }} // Z_x non inclus, monte d'un niveau. else {geopotential_J_lower = geopotential_J_upper; Zx_bot = geopotToZ (geopotential_J_lower); double virtualTemp_av_K = 0 .; double dPhi_J = 0 .; // Calcul de la température virtuelle moyenne entre le niveau de pression actuel // et le niveau inférieur. virtualTemp_av_K = (virtualTempFromSpecHum (inputGrids.at (0) ->getValue (n, j, i + 1), inputGrids.at (1) ->getValue (n, j, i + 1)) + virtualTempFromSpecH (entrée virtuelle) ->getValue (n-1, j, i + 1), inputGrids.at (1) ->getValue (n-1, j, i + 1))) / 2 .; // Calcule l'épaisseur actuelle. dPhi_J = geopotThickness (virtualTemp_av_K, inputGrids.at (0) ->getPressure (n-1, j, i + 1), inputGrids.at (0) ->getPressure (n, j, i + 1)); geopotential_J_upper + = dPhi_J; Zx_top = geopotToZ (geopotential_J_upper); }}
// coordonnée y double udy; // Valeur interpolée de u. double Z_y; // Valeur d'interpolation de Z. // (3c) Calcule la hauteur géopotentielle du point de grille actuel et ajoute la composante Z // du vecteur tangent du plan y-Z // Géopotentiel au point de grille actuel déjà connu. // Conversion en hauteur géopotentielle. Z_y = geopotToZ (geopotential000_J); // Ajouter une contribution de dégradé. Z_y + = dZ_y; // (3d) Parcourez les niveaux du modèle pour trouver des points de grille adjacents entourant Z_y. // À partir du niveau le plus bas, calculez la hauteur géopotentielle des points de grille adjacents. for (n = derivedGrid->getNumLevels () - 1; n > 1; n--) {// Géopotentiel de surface. double geopotSurf_J = inputGrids.at (4) ->getValue (0, j + 1, i); // En supposant une température comprise entre -45 et +60 degrés Celsius. // Température virtuelle de surface. double pSurf_hPa = inputGrids.at (7) ->getValue (0, j + 1, i) / 100 .; double virtualTempSurf_K = virtualTempFromDewPoint (inputGrids.at (5) ->getValue (0, j + 1, i), pSurf_hPa, inputGrids.at (6) ->getValue (0, j + 1, i)); // Double épaisseur géopotentielle de la couche de surface dPhiSurfLayer = geopotThickness (virtualTempSurf_K, inputGrids.at (0) ->getPressure (inputGrids.at (0) ->getNumLevels () - 1, j + 1, i), pSurf_hP // Voisin inférieur.
double geopotential_J_lower = geopotSurf_J; double Zy_bot = geopotToZ (geopotential_J_lower); // Voisin supérieur. double geopotential_J_upper = geopotential_J_lower + dPhiSurfLayer; double Zy_top = geopotToZ (geopotential_J_upper); // Vérifie si les points de grille adjacents englobent Z_x. // Si c'est le cas, vérifiez si la surface est impliquée ou non, puis interpolez. if ((Zy_bot < Z_y) && (Z_y < Zy_top)) {// Niveau de surface impliqué. if (n == derivedGrid->getNumLevels () - 1) {double u_top = inputGrids.at (2) ->getValue (inputGrids.at (0) ->getNumLevels () - 1, j + 1, i); double u_bot = inputGrids.at (8) ->getValue (0, j + 1, i); // Valeur interpolée de v. Udy = u_bot + (u_top - u_bot) * (Z_y - Zy_bot) / (Zy_top - Zy_bot); Pause; } // Niveau de surface non impliqué. else {double u_top = inputGrids.at (2) ->getValue (n-1, j + 1, i); double u_bot = inputGrids.at (2) ->getValue (n, j + 1, i); // Valeur interpolée de u. udy = u_bot + (u_top - u_bot) * (Z_y - Zy_bot) / (Zy_top - Zy_bot); Pause; }} // Z_y non inclus, monte d'un niveau. autre                        {
geopotential_J_lower = geopotential_J_upper; Zy_bot = geopotToZ (geopotential_J_lower); double virtualTemp_av_K = 0 .; double dPhi_J = 0 .; // Calcul de la température virtuelle moyenne entre le niveau de pression actuel // et le niveau inférieur. virtualTemp_av_K = (virtualTempFromSpecHum (inputGrids.at (0) ->getValue (n, j + 1, i), inputGrids.at (1) ->getValue (n, j + 1, i)) + virtualTempFromSpecH. (entrée virtuelle 0) ->getValue (n-1, j + 1, i), inputGrids.at (1) ->getValue (n-1, j + 1, i))) / 2 .; // Calcule l'épaisseur actuelle. dPhi_J = geopotThickness (virtualTemp_av_K, inputGrids.at (0) ->getPressure (n-1, j + 1, i), inputGrids.at (0) ->getPressure (n, j + 1, i)); geopotential_J_upper + = dPhi_J; Zy_top = geopotToZ (geopotential_J_upper); }} // OLD: calcule la longueur géométrique des vecteurs tangents. // NOUVEAU: Ecrire une nouvelle méthode en utilisant les dZ déjà calculés. double dx_Z = vecLengthZ (dx, 0, dZ_x); double dy_Z = vecLengthZ (0, dy, dZ_y); // Calculer les différences de vent double du = udy - inputGrids.at (2) ->getValue (k, j, i); double dv = vdx - inputGrids.at (3) ->getValue (k, j, i); // Calcul de la vorticité relative isentropique. double vo_rel = (du / dy_Z) - (dv / dx_Z); // Calcul de la stabilité statique. double dTheta = theta100 - theta000;
double dp = 100. * (p100 - p000); // Calcul du tourbillon planétaire. // vo_plan = 2. * w_0 * sin (2. * M_PI * lat / 360.); double w_0 = MetConstants :: EARTH_ROTATION_RATE; double lat000 = inputGrids.at (0) ->getNorthInterfaceLat (j); double vo_plan = 2. * w_0 * sin (2. * M_PI * lat000 / 360.); // Calculer PV [pvu] double facteur = MetConstants :: PVU_FACTOR; double g = MetConstants :: GRAVITY_ACCELERATION; double PV_pvu = facteur * g * (vo_rel + vo_plan) * (dTheta / dp); dérivésGrid->setValue (k, j, i, PV_pvu); }}}