| blob | b986f9a5b9399ce0efc99c6e67dcb526b7111d13 |
1 % TEST GRAVITY WAVES %
2 % propagation visqueuse totale ou partielle
4 % pour le moment, il n'y a pas distinction entre les deux axes horizontaux,
5 % c'est du 3D planaire disons
7 % {{{ pré-processing
9 clear all;
10 close all;
12 % flags
14 has_viscosity = 1;
16 % atmospheric data from USSA76
17 load modele1
19 % modele1:
20 % altitude (0 - 500km)
21 % density
22 % dérivée de la density
23 % dérivée du logarithme de la density
25 % tests à résolution variable
26 res = 1 % 10, 100, etc.
27 binf = 1
28 bsup = 9999
29 Z = modele1(binf:res:bsup,1);
30 rho = modele1(binf:res:bsup,2);
31 drho = modele1(binf:res:bsup,3);
32 dlgrho = modele1(binf:res:bsup,4);
34 xx = size(Z); % 9999
35 n_mod = xx(1) % 9999 le nombre de modes (autant que de tranches
36 % d'altitudes ?)
37 I(1:n_mod) = 1;
39 % calcul de la dérivée spatiale verticale
40 % TODO à vectoriser ^^
41 for j = 1:n_mod-1
42 dz(j) = (Z(j+1)-Z(j))*1.e3; % pas spatial en mètres, utilisé dans les
43 end % dérivées en dz j'imagine
44 dz(n_mod) = dz(n_mod-1); % FIXME ajouté pour la version vectorisée, bof
45 dz = dz'; % multiplication matricielle plus bas
47 % conditions limites sur la vitesse horizontale : pour le calcul du gradient de vent
48 u0b = -50; % vitesse en bas (bottom)
49 u0t = -52; % vitesse en haut (top)
50 % en l'absence de profil des vents externes, on s'en crée un linéaire, basique
51 % coefficient d'un gradient linéaire des vents, égal à d(u0)/dz
52 alpha = (u0t-u0b)/((Z(n_mod)-Z(1))*1.e3);
53 % FIXME gruik pour le gradient linéaire
54 alpha_local = alpha.*ones(n_mod,1);
55 % profil des vitesses moyennes qui en découle
56 u0 = Z.*(1.e3*alpha) + u0b;
58 M = 5.9768e+24; % masse de la Terre
59 G = 6.67e-11; % constante gravitationnelle
60 R = 6378.e3; % rayon de la Terre
62 %----------------------------------
63 % omega, kx et ky à la main
64 %----------------------------------
66 T = [50 100 150].*60 % trois périodes de temps : 3000, 6000, 9000 secondes
67 wn = 2*pi./T; % les trois fréquences associées...
69 h = 2500; % hauteur de ?
70 g = 9.8;
71 c = sqrt(g*h); % vitesse de phase des ondes tsunamigéniques
72 % elle est imposée par continuité à l'interface, c'est
73 % une condition initiale/limite si on peut dire
74 kx = wn./c; % nombres d'ondes horizontaux
75 ky = kx; % id.
77 %----------------------------------
78 % itérateurs spatio-temporels
80 nlat = max(size(ky)); % en fait, on a fait en sorte que nlon = nlat...
81 nlon = max(size(kx));
82 nt = max(size(wn));
84 %----------------------------------
85 if has_viscosity
86 % viscosité cinématique
87 % FIXME à récupérer comme data du modèle d'atmosphère
88 visc = 1.3e-5; % valeur moyenne pour commencer
89 end
91 %----------------------------------
92 % Calculs principaux
93 % L'objectif est de calculer les modes pour les IGW, ie. la valeur des paramètres
94 % des équations du système à chaque altitude (« mode ») pour laquelle est fixé
95 % un kz approprié, sous la contrainte des conditions limites. La condition limite
96 % bottom fait office de condition initiale, couplée à une condition sur w, une
97 % condition sur P (continuité à l'interface).
98 %
99 % On considère trois périodes de temps, donc trois fréquences wn, donc trois
100 % nombre d'ondes k à travers la vitesse de phase imposée comme condition au limite
101 % implicite par continuité à l'interface. Il y a donc trois cas étudiés (ik).
102 % Pour chaque cas, on analyse les trois fréquences (iw). La boucle sur les iw
103 % résout le propagateur, sous les conditions initiales/limites.
104 %----------------------------------
106 load BestView2
108 %----------------------------------
110 % pré-processing }}}
112 % {{{ boucle sur les nombres d'onde (sur les longeurs d'onde), ici trois.
113 % une figure par mode ik
114 for ik = 1:nlon % (nlon-1)/2 + 2:nlon
115 ik
116 drawnow; % DRAWNOW causes figure windows and their children to update and
117 % flushes the system event queue. Any callbacks generated by incoming
118 % events - e.g. mouse or key events - will be dispatched before
119 % DRAWNOW returns.
121 k = kx(ik) % nombre d'onde horizontal courant
123 figure(2*ik); % FIGURE(H) makes H the current figure, forces it to become visible,
124 % and raises it above all other figures on the screen. If Figure H
125 % does not exist, and H is an integer, a new figure is created with
126 % handle H.
127 set(2*ik, 'position', [44 250 500 750]);
128 set(gcf, 'Color', 'w');
129 set(gca, 'LineWidth', [2])
130 set(gca, 'fontsize', 18);
131 view(VV);
132 title(['Hwater = ', num2str(h), ' m and kx = ', num2str(k), ' rad/m'], 'fontsize', 18);
133 % xlabel('omega w (rad/s)', 'fontsize', 18);
134 xlabel('period T (min)', 'fontsize', 18);
135 ylabel('Vertical V_r_e_a_l (m/s)', 'fontsize', 18);
136 zlabel('altitude (km)', 'fontsize', 18);
137 grid on;
138 hold on;
140 figure(2*ik+1);
141 set(2*ik+1, 'position', [44 250 500 750]);
142 set(gcf, 'Color', 'w');
143 set(gca, 'LineWidth', [2])
144 set(gca, 'fontsize', 18);
145 view(VV);
146 title(['Hwater = ', num2str(h), ' m and kx = ', num2str(k), ' rad/m'], 'fontsize', 18);
147 % xlabel('omega (rad/s)', 'fontsize', 18);
148 xlabel('period T (min)', 'fontsize', 18);
149 ylabel('Vertical V_i_m_a_g (m/s)', 'fontsize', 18);
150 zlabel('altitude (km)', 'fontsize', 18);
151 grid on;
152 hold on;
154 % FIXME relou ces figures :)
155 close(2*ik);
156 close(2*ik+1);
158 % }}}
160 % {{{ boucle sur les trois fréquences par mode, sachant qu'un mode est une altitude ici
161 % avec nt = max(size(wn)) == 3 ici
162 for iw = 1:nt %(nt-1)/2 + 2:nt
164 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
165 % vectorisation en K*X = B
166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
168 % {{{ constantes et utilitaires
170 N2 = -g.*dlgrho;
171 N2over2g = (N2./(2*g)).^2;
172 omega = wn(iw)
173 BigOmega = omega - k.*u0; % FIXME k à vectoriser quand je passerai en full 3D
174 kxy2 = k^2; % FIXME no full 3D
175 phase = k - omega; % (pas unitaires) FIXME à modifier pour la full 3D
176 expniphase = exp(-i*phase); % FIXME vérifier cette formule
177 kz2 = k.*k.*(N2./omega./omega - 1) - N2over2g;
179 % pour le cas visqueux
180 BigGamma = BigOmega + i*visc*kxy2;
181 BigGammadz = BigGamma.*dz;
182 BigGammadz2 = BigGammadz.*dz;
183 phi = 2.*visc.*sqrt(rho);
184 phiOverBigGammadz2 = phi./BigGammadz2;
186 % }}}
188 % {{{ résolution vectorielle en différences finies
190 % {{{ coefficients de propagation non visqueuse, hors des conditions initiales
192 c1 = -(-k.*alpha.*2.*dz./BigOmega + dlgrho.*dz);
193 c2 = -i.*k.*k.*2.*dz./BigOmega;
194 c3 = -(-dlgrho.*dz);
195 c4 = -(i.*2.*dz.*(BigOmega + dlgrho.*g./BigOmega)); % }}}
197 % {{{ coefficients de propagation visqueuse, hors des conditions initiales
199 % coefficients pour la perturbation de la vitesse horizontale
200 cu1 = 1 + i.*phiOverBigGammadz2;
201 cu2 = i./2.*phiOverBigGammadz2;
202 cu3 = cu2;
203 cu4 = k./BigGamma;
204 for j = 3:n_mod-1
205 cu5(j) = (i*(u0(j+1) - u0(j-1)))/(2*BigGammadz(j));
206 end
207 % TODO FIXME à terme, il faudra peut-être adopter cette technique du step-forward d'orde 1 pour les conditions au sommet pour d'autres coefficients ?
208 cu5(n_mod) = (i*(u0(n_mod) - u0(n_mod-1)))/(BigGammadz(n_mod));
209 cu5 = cu5';
211 % coefficients pour la perturbation de la vitesse verticale
212 cw1 = -(dz.*dlgrho - 2.*k.*alpha./BigGamma);
213 cw2 = 2.*i.*kxy2.*dz./BigGamma;
214 cw3 = -(k.*phiOverBigGammadz2);
215 cw4 = 2.*k.*phiOverBigGammadz2;
216 cw5 = -(phiOverBigGammadz2);
218 % coefficients pour la perturbation de pression
219 cp1 = dz.*dlgrho;
220 cp2 = 2.*i.*(BigOmega + g./BigOmega.*dlgrho + i.*visc.*kxy2) + 2.*phi./dz;
221 cp3 = -phi./dz;
222 cp4 = -cp3; % }}}
224 % {{{ conditions initiales et coefficients associés
226 % {{{ C.I. pures (1)
228 % C.I.1 pour les perturbations de vitesse verticale et de pression
229 % on se donne une perturbation unitaire, dans l'optique de l'étude des modes
231 w1 = complex(1,0);
232 w(1) = w1;
234 kz2 = k*k*(-g/omega/omega*dlgrho(1) - 1) - 0.25/6.4e7;
235 N2 = -g*dlgrho(1);
236 if N2(1) >= omega*omega
237 % cas d'une onde propagative
238 p(1) = i/(k*k)*(alpha_local(1)*k - i*sqrt(kz2(1))*(omega - k*u0(1)) - 0.5*(omega - k*u0(1))*dlgrho(1))*w1;
239 else
240 % cas d'une onde non propagative
241 p(1) = i/(k*k)*(alpha_local(1)*k - sqrt(kz2(1))*(omega - k*u0(1)) - 0.5*(omega - k*u0(1))*dlgrho(1))*w1;
242 end
244 % C.I.1 pour la perturbation de la vitesse horizontale
245 % FIXME je trouve un facteur 1/rho0 pour u(j)
246 % u(j) = 1/BigOmega(j) * (k*p(j) - i*alpha*w(j));
247 % d'où les coefficients :
248 ui1 = -i*alpha_local(1)*w(1);
249 ui2 = k/BigOmega(1);
250 % par ailleurs, ces valeurs explicites sont nécessaires :
251 % C.I.1 perturbation horizontale
252 u(1) = 1/BigOmega(1) * (k*p(1) - i*alpha*w(1));
254 % C.I. pures (1) }}}
256 % {{{ C.I. au premier pas de temps (2)
258 b1 = -(-i*k*dz(1) * (1 / BigOmega(1) * (-i*alpha)) + 0.5*dlgrho(1)*dz(1) + 1);
259 b2 = -(-i*k*dz(1) * (1 / BigOmega(1) * (k)));
260 b3 = -(-0.5*dlgrho(1)*dz(1) + 1);
261 b4 = -(i*BigOmega(1)*dz(1) + i/BigOmega(1)*dlgrho(1)*dz(1)*g);
262 ui3 = -i*alpha_local(2)/BigOmega(2);
263 ui4 = k/BigOmega(2);
265 % C.I. au premier pas de temps (2) }}}
267 % conditions initiales et coefficients associés }}}
269 % {{{ construction de la matrice creuse
271 % construction des diagonales i.e. des vecteurs colonnes alternés, pour spdiags()
273 % attention, il faut faire une permutation circulaire* pour les diagonales d'indices impairs, de façon à avoir le bon coefficient sur la première ligne. Ici, permutations d'ordre 3, selon le numéro de colonne
275 % construction : deux lignes, puis lecture globale (colonne par colonne, donc en alternant les valeurs -- ex. si une ligne de 1 et une ligne de 0, donnera une suite de 1 et 0 alternés)
276 % {{{ FIXME diagonales pour le cas non visqueux
277 %C13 = [c1'; c3'];
278 %C2 = [zeros(1, size(c2,1)); c2']; % *
279 %C4 = [c4'; zeros(1, size(c4,1))]; % *
281 %diagC13 = C13(:);
282 %diagC2 = C2(:);
283 %% FIXME pb de signe ici ! en attendant, je mets un moins pour retrouver
284 %% des parties imaginaires positives...
285 %diagC4 = -C4(:);
287 %% matrice des conditions limites
288 %cl0 = [1 0; 0 1];
289 %cl1 = [b1 b2 1 0; b4 b3 0 1];
290 %%cl = [cl0 zeros(size(cl0,1), n_mod-size(cl0,1)) ; cl1 zeros(size(cl1,1), n_mod-size(cl1,1))]
291 %n_mod2 = 2*n_mod;
292 %cl = [cl0 zeros(2, n_mod2-2) ; cl1 zeros(2, n_mod2-4)];
294 %% construction de la matrice creuse des coefficients du système linéaire pour
295 %% l'itération en différences finies
296 %% partie utile de la matrice diagonale à laquelle sont ajoutées les 2*2
297 %% conditions limites/initiales
299 %plop = ones(n_mod2, 1);
300 %preK = spdiags([-plop diagC4 diagC13 diagC2 plop], 0:4, n_mod2, n_mod2);
301 %preK = preK(1:n_mod2 - 4, :);
302 %K = [cl; preK];
303 %%Kf = full(K);
304 %%Kf(1:8,1:8)
305 %%diagC2(1)
306 %%diagC2(2)
307 %%diagC2(3)
308 %%diagC2(4)
309 %%Kf(n_mod2-7:end,n_mod2-7:end)
310 %%size(Kf)
311 % }}}
313 % {{{ diagonales pour le cas visqueux
315 % helpers
316 dimvisc = size(cu1, 1);
317 myones = ones(1, dimvisc);
318 myzeros = zeros(1,dimvisc);
320 % diagonales alternées (les * dénotent celles affectées d'une permutation circulaire)
321 diag1 = [cw5'; cw4'; myzeros]; % *
322 diag1 = diag1(:);
324 diag2 = [cu2'; -myones; -myones];
325 diag2 = diag2(:);
327 diag3 = [myzeros myzeros myzeros]; % *
328 diag3 = diag3';
330 diag4 = [cp2'; myzeros; cw4']; %*
331 diag4 = diag4(:);
333 diag5 = [cu1'; cw1'; cp2'];
334 diag5 = diag5(:);
336 diag6 = [cw2'; myzeros; cu5']; % *
337 diag6 = diag6(:);
339 diag7 = [cp3'; cu4'; cw3']; % *
340 diag7 = diag7(:);
342 diag8 = [cu3'; myones; myones];
343 diag8 = diag8(:);
345 % cas visqueux }}}
347 % {{{ matrices des conditions initiales
349 ci0 = [1 ui1 ui2; ...
350 0 1 0; ...
351 0 0 1];
353 ci1 = [0 0 0 1 ui3 ui4; ...
354 0 b1 b2 0 1 0; ...
355 0 b4 b3 0 0 1];
357 n_mod3 = 3*n_mod;
358 ci = [ci0 zeros(3, n_mod3-3) ; ci1 zeros(3, n_mod3-6)];
360 % }}}
362 % {{{ construction effective de la matrice creuse
364 % construction de la matrice creuse des coefficients du système linéaire pour
365 % l'itération en différences finies
366 % partie utile de la matrice diagonale à laquelle sont ajoutées les 2*2
367 % conditions limites/initiales
369 preK = spdiags([diag1 diag2 diag3 diag4 diag5 diag6 diag7 diag8], ...
370 -1:6, n_mod3, n_mod3);
371 preK = preK(1:n_mod3 - 6, :);
372 K = [ci; preK];
373 disp('K :')
374 size(K)
375 %Kf = full(K);
376 %Kf(1:8,1:8)
377 %diagC2(1)
378 %diagC2(2)
379 %diagC2(3)
380 %diagC2(4)
381 %Kf(n_mod2-7:end,n_mod2-7:end)
382 %size(Kf)
384 % construction de la matrice creuse }}}
386 % {{{ construction du second membre
388 % on construit le vecteur B des seconds membres (SM), avec :
389 % - une partie conditions initiales
390 % - TODO éventuellement, une partie de propagation non visqueuse (SM = 0)
391 % ce qui nécessite de modifier K !
392 % - une partie de propagation visqueuse avec les SM appropriés
394 % {{{ calculs des seconds membres visqueux
396 % TODO il faut éventuellement calculer un linspace des u0 à partir du gradient linéaire, si on utilise pas de profil de vents. Pour les profils de vents, s'il n'y a pas le même nombre de pas en espace, il faut également faire des linspaces entre les points de donnée (interpolations linéaires)
397 for j = 6:n_mod-1
398 sm_tmp = (phi(j)/(BigGamma(j)*dz(j)))*expniphase*(u0(j+1) - 2*u0(j) + u0(j-1));
399 end
400 sm_tmp(n_mod) = (phi(n_mod)/(BigGamma(n_mod)*dz(n_mod-1)))*expniphase*(u0(n_mod) - 2*u0(n_mod-1) + u0(n_mod-2));
402 sm_u = i.*sm_tmp./2;
403 sm_w = k.*sm_tmp;
404 sm_p = zeros(1,n_mod-6);
406 % vecteur des valeurs alternées
407 sm = [sm_u'; sm_w'; sm_p'];
408 sm = sm(:);
410 disp('sm :')
411 size(sm)
413 % }}}
415 B = [u(1); w(1); p(1); ... % conditions initiales
416 0; 0; 0; ... % le premier pas est non visqueux
417 sm ... % la propagation visqueuse implique des termes non nuls,
418 ]; % tous calculés a priori avec les données du problème
420 disp('B :')
421 size(B)
423 % construction du second membre }}}
425 % {{{ résolution du système linéaire K*X = B
427 %%%%%%%%%%%%%%%%%%%%
428 U = K\B; % beautiful
429 %%%%%%%%%%%%%%%%%%%%
431 %size(U)
432 %U(1:13)
433 %U(n_mod3-12:end)
434 %p1
436 % résolution de K*X = B }}}
438 % résolution vectorielle en différences finies }}}
440 % {{{ résultats et calculs finaux
442 % récupération des résultats en w et p
443 u = U(1:3:end);
444 w = U(2:3:end);
445 p = U(3:3:end);
447 % calcul des kz explicites
448 % FIXME coefficient à remplacer par la formule analytique, non ?
449 kz2 = (k*k).*((-g/omega/omega).*dlgrho-1) - 0.25/6.4e7;
451 % calcul des rho1 explicites
452 rho1 = (-i.*dlgrho.*w)./BigOmega;
454 % résultats et calculs finaux }}}
456 % {{{ post-processing pour le pas courant
458 wmax(ik, iw) = max(abs(real(w)));
459 umax(ik, iw) = max(abs(real(u)));
460 pmax(ik, iw) = max(abs(real(p)));
461 rhomax(ik, iw) = max(abs(real(rho1)));
462 KZ(ik, iw, :) = sqrt(kz2);
464 if max(abs(real(w))) >= 2.
465 ctrw(ik, iw) = NaN;
466 else
467 ctrw(ik, iw) = 1.;
468 end
470 if sqrt(kz2(j)) >= 2*pi/50.e3
471 A(ik, iw) = 1;
472 else
473 A(ik, iw) = NaN;
474 end
476 figure(100);
477 set(100, 'position', [10 1 750 600])
478 set(gcf, 'Color', 'w');
479 fnt = 20;
480 fnt2 = 18;
481 Zmax = max(Z);
482 Zmin = min(Z);
483 subplot(1, 4, 1);
484 plot(real(w), Z, 'r', 'linewidth', 2);
485 hold on;
486 plot(imag(w), Z, 'g', 'linewidth', 2);
487 hold off;
488 %axis([-2 2 0 max(Z)])
489 v = axis;
490 axis([v(1) v(2) Zmin Zmax])
491 xlabel('Vertical V (m/s)', 'fontsize', 12);
492 ylabel('Altitude (km)', 'fontsize', 16);
493 title(['omega = ', num2str(omega), ' rad/s'], 'fontsize', 12);
494 subplot(1, 4, 2);
495 plot(real(u), Z, 'r', 'linewidth', 2);
496 hold on;
497 plot(imag(u), Z, 'g', 'linewidth', 2);
498 hold off;
499 %axis([-20 20 0 max(Z)])
500 v = axis;
501 axis([v(1) v(2) Zmin Zmax]);
502 xlabel('Horizontal V (m/s)', 'fontsize', 12);
503 subplot(1, 4, 3);
504 plot(real(p), Z, 'r', 'linewidth', 2);
505 hold on;
506 plot(imag(p), Z, 'g', 'linewidth', 2);
507 hold off;
508 %axis([-3000 3000 0 max(Z)])
509 v = axis;
510 axis([v(1) v(2) Zmin Zmax]);
511 xlabel('Pressure (Pa)', 'fontsize', 12);
512 subplot(1, 4, 4);
513 plot(real(rho1), Z, 'r', 'linewidth', 2);
514 hold on;
515 plot(imag(rho1), Z, 'g', 'linewidth', 2);
516 hold off;
517 %axis([-0.1 0.1 0 max(Z)])
518 v = axis;
519 axis([v(1) v(2) Zmin Zmax]);
520 xlabel('Density (kg/m3)', 'fontsize', 12);
521 title(['kx = ', num2str(k), ' rad/m'], 'fontsize', 12);
523 figure(100);
524 F = getframe(gcf);
525 [X, Map] = frame2im(F);
526 file = ['Earth_Mode_L' num2str(round(1/k*1.e-3)) '_T_' num2str(2*pi/omega/60) '.tif'];
527 imwrite(X, file, 'tif', 'Compression', 'none');
529 disp('PRESS ENTER TO CONTINUE...');
530 pause
531 close
533 figure(300);
534 subplot(1, 3, 1);
535 plot(real(sqrt(kz2)), Z, 'r');
536 hold on;
537 plot(imag(sqrt(kz2)), Z, 'g');
538 xlabel('kz (rad/m)', 'fontsize', 12);
539 ylabel('Altitude (km)', 'fontsize', 16);
540 hold off;
541 subplot(1, 3, 2);
542 plot(sqrt(N2), Z, 'r');
543 hold on;
544 v = axis;
545 plot([omega omega], [v(3) v(4)], 'g');
546 hold off;
547 xlabel('N and omega (rad/s)', 'fontsize', 12);
548 subplot(1, 3, 3);
549 plot(2*pi./real(sqrt(kz2))*1.e-3, Z, 'r');
550 hold on;
551 plot(2*pi./imag(sqrt(kz2))*1.e-3, Z, 'g');
552 hold off;
553 xlabel('lambda Z (km)', 'fontsize', 12);
555 pause
556 close
558 figure(2*ik);
559 plot3(I*2*pi/omega/60, real(w), Z, 'k');
560 plot3(I*2*pi/omega/60, imag(w), Z, 'r');
562 figure(2*ik+1);
563 plot3(I*2*pi/omega/60, imag(w), Z, 'k');
565 clear w u p rho1 kz2
567 % post-processing pour le pas courant }}}
569 end
570 end % }}}
572 % boucle sur les trois fréquences }}}
574 % {{{ post-processing
576 Nsol = sqrt(N2(1));
577 Nmax = max(sqrt(N2));
578 Nmin = min(sqrt(N2));
580 njump = 0;
581 njump2 = size(wmax);
583 n1 = 1;
584 n2 = 1;
586 h = figure;
587 set(h, 'position', [2 700 800 600]);
588 set(gcf, 'Color', 'w');
589 font_size = 22;
590 font_sizeN = 18;
592 flag_plot = input('-Plot omega-k or T-lambda ? (0 or 1) ');
594 subplot(2, 2, 1);
595 if flag_plot == 0
596 pcolor(wn(n1:nt), kx(n2:nlon), log10(wmax(:,1+njump:njump2(2))));
597 else
598 pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(wmax(:,1+njump:njump2(2))));
599 end
600 shading flat;
601 v = axis;
602 hold on;
603 plot([Nmax Nmax], [v(3) v(4)], 'k');
604 plot([Nsol Nsol], [v(3) v(4)], 'k');
605 plot([Nmin Nmin], [v(3) v(4)], 'k');
606 colorbar;
607 title('W_r_e_a_l', 'fontsize', font_size);
608 if flag_plot == 0
609 ylabel('kx (rad/m)', 'fontsize', font_size);
610 else
611 ylabel('lambda (km)', 'fontsize', font_size);
612 end
613 set(gca, 'LineWidth', [2])
614 set(gca, 'fontsize', font_sizeN);
616 subplot(2, 2, 2);
617 if flag_plot == 0
618 pcolor(wn(n1:nt), kx(n2:nlon), log10(umax(:,1+njump:njump2(2))));
619 else
620 pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(umax(:,1+njump:njump2(2))));
621 end
622 shading flat;
623 v = axis;
624 hold on;
625 plot([Nmax Nmax], [v(3) v(4)], 'k');
626 plot([Nsol Nsol], [v(3) v(4)], 'k');
627 plot([Nmin Nmin], [v(3) v(4)], 'k');
628 colorbar;
629 title('U_r_e_a_l', 'fontsize', font_size);
630 set(gca, 'LineWidth', [2])
631 set(gca, 'fontsize', font_sizeN);
633 subplot(2, 2, 3);
634 if flag_plot == 0
635 pcolor(wn(n1:nt), kx(n2:nlon), log10(pmax(:,1+njump:njump2(2))));
636 else
637 pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(pmax(:,1+njump:njump2(2))));
638 end
639 shading flat;
640 v = axis;
641 hold on;
642 plot([Nmax Nmax], [v(3) v(4)], 'k');
643 plot([Nsol Nsol], [v(3) v(4)], 'k');
644 plot([Nmin Nmin], [v(3) v(4)], 'k');
645 colorbar;
646 if flag_plot == 0
647 ylabel('kx (rad/m)', 'fontsize', font_size);
648 xlabel('omega (rad/s)', 'fontsize', font_size);
649 else
650 ylabel('lambda (km)', 'fontsize', font_size);
651 xlabel('period T (min)', 'fontsize', font_size);
652 end
653 title('P_r_e_a_l', 'fontsize', font_size);
654 set(gca, 'LineWidth', [2])
655 set(gca, 'fontsize', font_sizeN);
657 subplot(2, 2, 4);
658 if flag_plot == 0
659 pcolor(wn(n1:nt), kx(n2:nlon), log10(rhomax(:,1+njump:njump2(2))));
660 else
661 pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(rhomax(:,1+njump:njump2(2))));
662 end
663 shading flat;
664 v = axis;
665 hold on;
666 plot([Nmax Nmax], [v(3) v(4)], 'k');
667 plot([Nsol Nsol], [v(3) v(4)], 'k');
668 plot([Nmin Nmin], [v(3) v(4)], 'k');
669 colorbar;
670 if flag_plot == 0
671 xlabel('omega (rad/s)', 'fontsize', font_size);
672 else
673 xlabel('period T (min)', 'fontsize', font_size);
674 end
675 title('Rho_r_e_a_l', 'fontsize', font_size);
676 set(gca, 'LineWidth', [2])
677 set(gca, 'fontsize', font_sizeN);
679 % post-processing }}}