ninto/MATLAB/viscosity/modesT_visc.m

   1
   2                 %%%%%%%%%%%%%%%%%%%%%%
   3                 % TEST GRAVITY WAVES %
   4                 %    cas visqueux    %
   5                 %%%%%%%%%%%%%%%%%%%%%%
   6
   7 % pour le moment, il n'y a pas distinction entre les deux axes horizontaux,
   8 % c'est du 3D planaire disons :)
   9
  10 % {{{ preprocessing
  11
  12 clear all;
  13 close all;
  14
  15 % enable viscosity?
  16 has_viscosity = 0;
  17
  18 % modele1:
  19 %   altitude (0 - 500km)
  20 %   densité
  21 %   dérivée de la densité
  22 %   dérivée du logarithme de la densité
  23 load modele1
  24
  25
  26 %Z = RHOout(:,1);
  27 %Z       = modele1(:,1);
  28 %nnn     = max(size(Z));     % 9999 (1D)
  29 %rho     = modele1(:,2);
  30 %drho    = modele1(:,3);
  31 %dlgrho  = modele1(:,4);
  32
  33 % tests à plus basse résolution
  34 res     = 10;
  35 binf    = 1;
  36 bsup    = 9999;
  37 Z       = modele1(binf:res:bsup,1);
  38 rho     = modele1(binf:res:bsup,2);
  39 drho    = modele1(binf:res:bsup,3);
  40 dlgrho  = modele1(binf:res:bsup,4);
  41
  42 xx         = size(Z);       % 9999
  43 n_mod      = xx(1);         % 9999 le nombre de modes (autant que de tranches
  44 % d'altitudes ?)
  45 I(1:n_mod) = 1;
  46
  47 % calcul de la dérivée spatiale verticale ?
  48 for j = 1:n_mod-1
  49     dz(j) = (Z(j+1)-Z(j))*1.e3;     % pas spatial en mètres, utilisé dans les
  50     % dérivées en dz j'imagine
  51 end
  52 dz(n_mod) = dz(n_mod-1);          % FIXME ajouté pour la version vectorisée, bof
  53
  54 dz = dz';
  55
  56 % conditions limites sur la vitesse horizontale : pour le calcul du gradient de vent
  57 u0b = 0;      % vitesse en bas  (bottom)
  58 u0t = 0;      % vitesse en haut (top)
  59 % coefficient d'un gradient linéaire des vents égal à d(u0)/dz
  60 alpha = (u0t-u0b)/((Z(n_mod)-Z(1))*1.e3);
  61 % FIXME gruik pour le gradient linéaire
  62 alpha_local = alpha.*ones(n_mod,1);
  63 % profil des vitesses moyennes
  64 u0 = Z.*(1.e3*alpha) + u0b;
  65
  66 M = 5.9768e+24;     % masse de la Terre
  67 G = 6.67e-11;       % constante gravitationnelle
  68 R = 6378.e3;        % rayon de la Terre
  69
  70 %----------------------------------
  71 %   omega, kx & ky à la main
  72 %----------------------------------
  73
  74 T  = [50 100 150]*60    % trois périodes de temps (s) : 3000, 6000, 9000 secondes
  75 %T  = [10 20 30]*60
  76 wn = 2*pi./T;           % les trois fréquences associées...
  77
  78 h  = 2500;              % hauteur de l'océan
  79 g  = 9.8;
  80 c  = sqrt(g*h);         % vitesse de phase des ondes acoustiques simples
  81 kx = wn/c;              % nombres d'ondes horizontaux
  82 ky = kx;                % .
  83
  84 %----------------------------------
  85 % itérateurs spatio-temporels
  86
  87 nlat = max(size(ky));   % en fait, on a fait en sorte que nlon = nlat...
  88 nlon = max(size(kx));
  89 nt   = max(size(wn));
  90
  91 %----------------------------------
  92 if has_viscosity
  93     % viscosité cinématique
  94     % FIXME à récupérer comme data du modèle d'atmosphère
  95     visc = 1.3e-5; % valeur moyenne pour commencer
  96 end
  97
  98 %----------------------------------
  99 % bibliothèques
 100
 101 load BestView2
 102
 103 % }}}
 104
 105 % {{{ résolution en différences finies
 106
 107 % {{{ boucle sur les nombres d'onde
 108 % (sur les longeurs d'onde, ici trois)
 109 % une figure par mode ik pour le Hwater(kx)
 110 for ik = 1:nlon % (nlon-1)/2 + 2:nlon
 111
 112     ik
 113     k    = kx(ik);          % nombre d'onde horizontal courant
 114
 115     %drawnow;        % DRAWNOW causes figure windows and their children to update and
 116                     %% flushes the system event queue. Any callbacks generated by incoming
 117                     %% events - e.g. mouse or key events - will be dispatched before
 118                     %% DRAWNOW returns.
 119
 120     %figure(2*ik);   % FIGURE(H) makes H the current figure,  forces it to become visible,
 121                     %% and raises it above all other figures on the screen.  If Figure H
 122                     %% does not exist,  and H is an integer,  a new figure is created with
 123                     %% handle H.
 124
 125     %set(2*ik, 'position', [44   250   500   750]);
 126     %set(gcf, 'Color', 'w');
 127     %set(gca, 'LineWidth', [2])
 128     %set(gca, 'fontsize', 18);
 129     %view(VV);
 130     %title(['Hwater = ', num2str(h), 'm and kx = ',  num2str(k), ' rad/m'], 'fontsize', 18);
 131     %%     xlabel('omega w (rad/s)', 'fontsize', 18);
 132     %xlabel('period T (min)', 'fontsize', 18);
 133     %ylabel('Vertical V_r_e_a_l (m/s)', 'fontsize', 18);
 134     %zlabel('altitude (km)', 'fontsize', 18);
 135     %grid on;
 136     %hold on;
 137
 138     %figure(2*ik+1);
 139     %set(2*ik+1, 'position', [44   250   500   750]);
 140     %set(gcf, 'Color', 'w');
 141     %set(gca, 'LineWidth', [2])
 142     %set(gca, 'fontsize', 18);
 143     %view(VV);
 144     %title(['Hwater = ', num2str(h), 'm   and   kx = ',  num2str(k), ' rad/m'], 'fontsize', 18);
 145     %%    xlabel('omega (rad/s)', 'fontsize', 18);
 146     %xlabel('period T (min)', 'fontsize', 18);
 147     %ylabel('Vertical V_i_m_a_g (m/s)', 'fontsize', 18);
 148     %zlabel('altitude (km)', 'fontsize', 18);
 149     %grid on;
 150     %hold on;
 151
 152     % }}}
 153
 154     % {{{ boucle sur les trois fréquences par mode
 155     %     sachant qu'un mode est une altitude ici
 156     for iw = 1:nt %(nt-1)/2 + 2:nt
 157         % avec nt = max(size(wn)) == 3 ici
 158
 159         % {{{ constantes et utilitaires
 160
 161         N2         = -g.*dlgrho;
 162         N2over2g   = (N2./(2*g)).^2;
 163         omega      = wn(iw)
 164         BigOmega   = omega - k.*u0;                % FIXME k à vectoriser quand je passerai en full 3D
 165         kxy2       = k^2;                          % FIXME no full 3D
 166         phase      = k - omega;                    % (pas unitaires) FIXME à modifier pour la full 3D
 167         expniphase = exp(-i*phase);                % FIXME vérifier cette formule
 168         kz2        = k.*k.*(N2./omega./omega - 1) - N2over2g;
 169
 170         % }}}
 171
 172         % {{{ conditions initiales
 173
 174         % TODO pour les conditions initiales :
 175         % il est compliqué de propager directement en visqueux car il faut deux conditions
 176         % initiales. Pour palier ce problème, sous l'hypothèse que la viscosité joue un très
 177         % faible rôle sur les premiers kilomètres de propagation verticale, on commencera
 178         % la propagation en utilisant le propagateur non visqueux, puis la version visqueuse
 179         % prendra le relai. Il faudra déterminer à quelle altitude on fait cette transition.
 180         % On peut aussi faire tout en visqueux sauf les deux premiers pas, mais dans l'idée,
 181         % il est intéressant de ne pas considérer la viscosité aussi longtemps que possible,
 182         % car les calculs sont plus coûteux.
 183
 184         % {{{ première condition initiale ; détermination du type d'onde
 185         for j = 1:1
 186
 187             % C.I.1 perturbations verticale et de pression
 188             % on se donne une perturbation unitaire, dans l'optique de l'étude des modes
 189             w1   = complex(1,0);
 190             w(j) = w1;
 191
 192             % C.I.1 densité
 193             % FIXME je trouve un facteur 1/rho0 pour rho1
 194             rho1(j) = -i/BigOmega(j)*dlgrho(j)*w(j);
 195
 196             % calcul de la pression initiale sous une contrainte de continuité à l'interface
 197             if N2(j) >= omega*omega
 198                 % cas d'une onde propagative
 199                 p(j) = i/(k*k)*(alpha*k - i*sqrt(kz2(j))*(omega - k*u0(1)) - 0.5*(omega - k*u0(j))*dlgrho(j))*w1;
 200             else
 201                 % cas d'une onde non propagative
 202                 p(j) = i/(k*k)*(alpha*k -   sqrt(kz2(j))*(omega - k*u0(j)) - 0.5*(omega - k*u0(j))*dlgrho(j))*w1;
 203             end
 204
 205             % C.I.1 perturbation horizontale
 206             % FIXME je trouve un facteur 1/rho0 pour u(j)
 207             u(j)   = 1/BigOmega(j) * (k*p(j) - i*alpha*w(j));
 208             %v(j)   = u(j);
 209
 210             if has_viscosity
 211                 w_diff(j) = 0; u_diff(j) = 0; p_diff(j) = 0; rho1_diff(j) = 0;
 212             end
 213
 214         end % première condition intiale }}}
 215
 216         % {{{ seconde condition initiale pour l'ordre 2
 217         for j = 2:2
 218
 219             % idée : prendre les dérivées premières
 220             % pour les perturbations, on peut considérer que leurs accélérations sont nulles
 221             % i.e. utiliser des conditions de Dirichlet homogènes
 222             %u(j) = 0;
 223             %v(j) = 0;    % FIXME no full 3D
 224
 225             % on utilise le propagateur non visqueux au premier ordre pour initier le calcul
 226             % FIXME dans l'idée, on l'utilisera aussi après, tant que la viscosité
 227             % est négligeable : ça diminue le coût en calcul :)
 228             w(j)    = -i*k*u(j-1)*dz(j-1) + 0.5*dlgrho(j-1)*dz(j-1)*w(j-1) + w(j-1);
 229
 230             p(j)    = i*BigOmega(j-1)*dz(j-1)*w(j-1) - 0.5*dlgrho(j-1)*dz(j-1)*p(j-1) - rho1(j-1)*dz(j-1)*g + p(j-1);
 231
 232             u(j)    =  1/BigOmega(j)*(-i*w(j)*alpha + k*p(j));
 233
 234             rho1(j) = -i/BigOmega(j)*dlgrho(j)*w(j);
 235
 236             if has_viscosity
 237                 w_diff(j) = 0; u_diff(j) = 0; p_diff(j) = 0; rho1_diff(j) = 0;
 238             end
 239
 240         end
 241         % seconde condition initiale }}}
 242
 243         % conditions initiales }}}
 244
 245         % {{{ propagation verticale
 246         for j = 3:n_mod
 247
 248             % {{{
 249             % cas visqueux, en plusieurs étapes :
 250             % 1. calcul des perturbations horizontales u(j) et v(j)
 251             % 2. calcul du nouveau w(j) en fonction de u(j) et v(j)
 252             % 3. calcul du nouveau p(j) en fonction de w(j)
 253
 254             % calcul des perturbations horizontales
 255             % avec (1) et (2)
 256             %u(j) = 2*u(j-1) ...
 257             %     - u(j-2) ...
 258             %     - dz(j-1)^2 * ...
 259             %       ( (1/visc) * (i*BigOmega(j) * u(j-1) - w(j-1) * ((u0(j) - u0(j-2))/(2*dz(j-1))) - (1/rho(j-1))*(i*k*p(j-1))) ...
 260             %         - kxy2 * u(j-1) ...
 261             %         + expniphase * (u0(j) - 2*u0(j-1) + u0(j-2))/(dz(j-1)^2) );
 262
 263             % FIXME commenté pour le moment, 3D == 2D avec une seule dimension horizontale pour le moment
 264             % si réintégré, check le facteur alpha qui est 1D pour le moment
 265             %v(j) = 2*v(j-1) ...
 266             %    - v(j-2) ...
 267             %    - dz(j-1)^2 * ...
 268             %      ( (1/visc) * (i*BigOmega(j) * v(j-1) - w(j-1) * ((v0(j) - v0(j-2))/(2*dz(j-1))) - (1/rho(j-1))*(i*kx*p(j-1))) ...
 269             %        - kxy2 * v(j-1) ...
 270             %        + expniphase * (v0(j) - 2*v0(j-1) + v0(j-2))/(dz(j-1)^2) );
 271
 272             % calcul de la perturbation verticale
 273             % avec (1), (2) et (4)
 274             %BigOmegaVisc = (BigOmega(j-1) + i*visc*kxy2);
 275             %w(j) = w(j-1)*(dz(j-1)*dlgrho(j-1) - (2*dz(j-1)*k*alpha)/BigOmegaVisc) ...
 276             %     - w(j-2) ...
 277             %     - 2*dz(j-1)*((1/BigOmegaVisc)*(i*kxy2 + visc*sqrt(rho(j-1))*( ...
 278             %       k*( (u(j) - 2*u(j-1) + u(j-2))/(dz(j-1)^2) + expniphase*((u0(j) - 2*u0(j-1) + u0(j-2))/(dz(j-1)^2)) ))));
 279                    %+ ky*( (v(j) - 2*v(j-1) + v(j-2))/(dz(j-1)^2) + expniphase*((v0(j) - 2*v0(j-1) + v0(j-2))/(dz(j-1)^2))
 280
 281             % calcul de la perturbation de pression
 282             % avec (3) et (5)
 283             %p(j) = p(j-1) - dz(j-1)*dlgrho(j-1)*p(j-2) ...
 284             %     + 2*dz(j-1)*(i*BigOmega(j-1) - visc*kxy2 + (i*g*dlgrho(j-1))/(BigOmega(j-1)))*w(j-1) ...
 285             %     + 2*visc*sqrt(rho(j-1))*((w(j) - 2*w(j-1) + w(j-2))/(dz(j-1)));
 286
 287             % calcul de la perturbation de densité
 288             %rho1(j) = -i/BigOmega(j)*dlgrho(j)*w(j);
 289             % }}}
 290
 291             % FIXME retour à la version non visqueuse pour le moment
 292             % TODO  if pour intégrer le cas visqueux, sur un critère à déterminer
 293
 294             w(j)    = -i*k*u(j-1)*dz(j-1) + 0.5*dlgrho(j-1)*dz(j-1)*w(j-1) + w(j-1);
 295
 296             p(j)    = i*BigOmega(j-1)*dz(j-1)*w(j-1) - 0.5*dlgrho(j-1)*dz(j-1)*p(j-1) - rho1(j-1)*dz(j-1)*g + p(j-1);
 297
 298             u(j)    =  1/BigOmega(j)*(-i*w(j)*alpha + k*p(j));
 299
 300             rho1(j) = -i/BigOmega(j)*dlgrho(j)*w(j);
 301
 302             if has_viscosity
 303
 304                 w_old = w(j); p_old = p(j); u_old = u(j); rho1_old = rho1(j);
 305
 306                 % FIXME cas visqueux
 307
 308                 % précalculs
 309                 % TODO à faire une fois les formules fixées
 310
 311                 % calculs
 312
 313                 % même formule, mais u(visc)
 314                 w(j)      = -i*k*u(j-1)*2*dz(j-1) ...
 315                           + dlgrho(j-1)*dz(j-1)*w(j-1) ...
 316                           + w(j-2);
 317
 318                 w_diff(j) = w(j) - w_old;
 319
 320                 % nouveaux termes
 321                 p_visc(j) = -2*dz(j-1)*w(j-1)*visc*kxy2 ...
 322                           + (2*visc)/(dz(j-1))*(w(j) - 2*w(j-1) + w(j-2));
 323
 324                 p(j)      = -dlgrho(j-1)*dz(j-1)*p(j-1) ...
 325                           + 2*dz(j-1)*w(j-1)*i*BigOmega(j-1) ...
 326                           - rho1(j-1)*2*dz(j-1)*g ...
 327                           + p(j-2) ...
 328                           + p_visc(j);
 329
 330                 p_diff(j) = p(j) - p_old;
 331
 332                 % nouveaux termes
 333                 u(j)      = ((i*BigOmega(j))/visc - kxy2 + (dz(j-1))^(-2))^(-1) ...
 334                             * ( alpha/visc*w(j) ...
 335                                 + (i*k)/(visc)*p(j) ...
 336                                 - sqrt(rho(j)*expniphase*((alpha_local(j)-alpha_local(j-2))/(2*dz(j-1)))) ...
 337                                 - (dz(j-1))^(-2)*(u(j-2) - 2*u(j-1)));
 338
 339                 u_diff(j) = u(j) - u_old;
 340
 341                 % même formule, mais w(u(visc))
 342                 rho1(j)   = -i/BigOmega(j)*dlgrho(j)*w(j);
 343
 344                 rho1_diff(j) = rho1(j) - rho1_old;
 345
 346                 % TODO
 347                 % calculer l'effet visqueux à part et le plotter
 348                 % vérifier son ordre de grandeur et sa dimension analytique
 349                 % changer les fréquences pour voir la dépendance
 350
 351             end
 352
 353         end % }}}
 354
 355         % {{{ cas finaux
 356         % }}}
 357
 358         % {{{ figures récapitulatives
 359
 360         % bornes supérieures des quantités propagées
 361         wmax(ik, iw)   = max(abs(real(w)));
 362         umax(ik, iw)   = max(abs(real(u)));
 363         pmax(ik, iw)   = max(abs(real(p)));
 364         rhomax(ik, iw) = max(abs(real(rho1)));
 365         % nombre d'onde vertical
 366         KZ(ik, iw, :)  = sqrt(kz2);
 367
 368         % FIXME pas compris
 369         if max(abs(real(w))) >=  2.
 370             ctrw(ik, iw) = NaN;
 371         else
 372             ctrw(ik, iw) = 1.;
 373         end
 374
 375         if sqrt(kz2(j)) >=  2*pi/50.e3
 376             A(ik, iw) = 1;
 377         else
 378             A(ik, iw) = NaN;
 379         end
 380
 381         % figure principale, dessin des modes dans des subfigures
 382         figure(100);
 383         set(100, 'position', [10 1 750 600])
 384         set(gcf, 'Color', 'w');
 385         fnt   =  20;
 386         fnt2  =  18;
 387         Zmax  =  max(Z);
 388         Zmin  =  min(Z);
 389
 390         % perturbation de la vitesse verticale, parties réelles et imaginaires
 391         subplot(1, 4, 1);
 392         plot(real(w), Z, 'r', 'linewidth', 2);
 393         hold on;
 394         plot(imag(w), Z, 'g', 'linewidth', 2);
 395         hold off;
 396         v = axis;
 397         axis([v(1) v(2) Zmin Zmax])
 398         xlabel('Vertical V (m/s)', 'fontsize', 12);
 399         ylabel('Altitude (km)', 'fontsize', 16);
 400         title(['Period = ', num2str(2*pi/omega/60), ' mn' ' (omega = ', num2str(omega), ' rad/s)'], 'fontsize', 12);
 401
 402         % perturbation de la vitesse horizontale, parties réelles et imaginaires
 403         subplot(1, 4, 2);
 404         plot(real(u), Z, 'r', 'linewidth', 2);
 405         hold on;
 406         plot(imag(u), Z, 'g', 'linewidth', 2);
 407         hold off;
 408         v = axis;
 409         axis([v(1) v(2) Zmin Zmax]);
 410         xlabel('Horizontal V (m/s)', 'fontsize', 12);
 411
 412         % perturbation de pression, parties réelles et imaginaires
 413         subplot(1, 4, 3);
 414         plot(real(p), Z, 'r', 'linewidth', 2);
 415         hold on;
 416         plot(imag(p), Z, 'g', 'linewidth', 2);
 417         hold off;
 418         v = axis;
 419         axis([v(1) v(2) Zmin Zmax]);
 420         xlabel('Pressure (Pa)', 'fontsize', 12);
 421
 422         % perturbation de la densité, parties réelles et imaginaires
 423         subplot(1, 4, 4);
 424         plot(real(rho1), Z, 'r', 'linewidth', 2);
 425         hold on;
 426         plot(imag(rho1), Z, 'g', 'linewidth', 2);
 427         hold off;
 428         v = axis;
 429         axis([v(1) v(2) Zmin Zmax]);
 430         xlabel('Density (kg/m3)', 'fontsize', 12);
 431         title(['Lambda = ', num2str(round(1/k*1e-3)), ' km'  ' (kx = ',  num2str(k), ' rad/m)'], 'fontsize', 12);
 432
 433         % export TIF
 434         figure(100);
 435         F = getframe(gcf);
 436         [X, Map]  =  frame2im(F);
 437         file  =  ['mode_L' num2str(round(1/k*1.e-3)) '_T' num2str(2*pi/omega/60) '.tif'];
 438         imwrite(X, file, 'tif', 'Compression', 'none');
 439
 440         %% paramètres
 441         %disp('paramètres !');
 442         %figure(300);
 443         %% nombres d'onde
 444         %subplot(1, 3, 1);
 445         %plot(real(sqrt(kz2)), Z, 'r');
 446         %hold on;
 447         %plot(imag(sqrt(kz2)), Z, 'g');
 448         %xlabel('kz (rad/m)', 'fontsize', 12);
 449         %ylabel('Altitude (km)', 'fontsize', 16);
 450         %hold off;
 451         %% fréquences
 452         %subplot(1, 3, 2);
 453         %plot(sqrt(N2), Z, 'r');
 454         %hold on;
 455         %v = axis;
 456         %plot([omega omega], [v(3) v(4)], 'g');
 457         %hold off;
 458         %xlabel('N and omega (rad/s)', 'fontsize', 12);
 459         %% longueurs d'onde
 460         %subplot(1, 3, 3);
 461         %plot(2*pi./real(sqrt(kz2))*1.e-3, Z, 'r');
 462         %hold on;
 463         %plot(2*pi./imag(sqrt(kz2))*1.e-3, Z, 'g');
 464         %hold off;
 465         %xlabel('lambda Z (km)', 'fontsize', 12);
 466
 467         %pause
 468         %close
 469
 470         if has_viscosity
 471
 472             % perturbations, parties visqueuses
 473             figure(500);
 474             subplot(1,4,1);
 475             plot(real(w_diff), Z, 'r', 'linewidth', 2);
 476             hold on;
 477             plot(imag(w_diff), Z, 'g', 'linewidth', 2);
 478             hold off;
 479             v = axis;
 480             axis([v(1) v(2) Zmin Zmax])
 481             xlabel('Vertical V, viscosity related (m/s)', 'fontsize', 12);
 482             ylabel('Altitude (km)', 'fontsize', 16);
 483             title(['omega =  ', num2str(omega), ' rad/s'], 'fontsize', 12);
 484
 485             subplot(1,4,2);
 486             plot(real(u_diff), Z, 'r', 'linewidth', 2);
 487             hold on;
 488             plot(imag(u_diff), Z, 'g', 'linewidth', 2);
 489             hold off;
 490             v = axis;
 491             axis([v(1) v(2) Zmin Zmax])
 492             xlabel('Horizontal V (m/s), viscosity related', 'fontsize', 12);
 493
 494             subplot(1,4,3);
 495             plot(real(p_diff), Z, 'r', 'linewidth', 2);
 496             hold on;
 497             plot(imag(p_diff), Z, 'g', 'linewidth', 2);
 498             hold off;
 499             v = axis;
 500             axis([v(1) v(2) Zmin Zmax])
 501             xlabel('Pressure (Pa), viscosity related', 'fontsize', 12);
 502
 503             subplot(1, 4, 4);
 504             plot(real(rho1_diff), Z, 'r', 'linewidth', 2);
 505             hold on;
 506             plot(imag(rho1_diff), Z, 'g', 'linewidth', 2);
 507             hold off;
 508             v = axis;
 509             axis([v(1) v(2) Zmin Zmax]);
 510             xlabel('Density (kg/m3), viscosity related', 'fontsize', 12);
 511             title(['kx = ',  num2str(k), ' rad/m'], 'fontsize', 12);
 512
 513         end
 514
 515         % graphs 3D
 516         % perturbations de la vitesse verticale
 517         %figure(2*ik);
 518         %plot3(I*2*pi/omega/60, real(w), Z, 'k');
 519         %plot3(I*2*pi/omega/60, imag(w), Z, 'r');
 520
 521         %pause
 522         %close
 523
 524         disp('PRESS ENTER TO CONTINUE...');
 525         pause
 526         close
 527
 528         % next!
 529         clear w u p rho1 kz2
 530
 531         % }}}
 532
 533     end
 534     % }}}
 535 end
 536
 537 % }}}
 538
 539 % {{{ post-processing
 540
 541 Nsol = sqrt(N2(1));
 542 Nmax = max(sqrt(N2));
 543 Nmin = min(sqrt(N2));
 544
 545 njump  = 0;
 546 njump2 = size(wmax);
 547
 548 n1 = 1;
 549 n2 = 1;
 550
 551 h = figure;
 552 set(h, 'position', [2 700 800 600]);
 553 set(gcf, 'Color', 'w');
 554 font_size  =  22;
 555 font_sizeN =  18;
 556
 557 flag_plot  =  input('-Plot omega-k or T-lambda ? (0 or 1) ');
 558
 559 subplot(2, 2, 1);
 560 if flag_plot == 0
 561     pcolor(wn(n1:nt), kx(n2:nlon), log10(wmax(:,1+njump:njump2(2))));
 562 else
 563     pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(wmax(:,1+njump:njump2(2))));
 564 end
 565 shading flat;
 566 v = axis;
 567 hold on;
 568 plot([Nmax Nmax], [v(3) v(4)], 'k');
 569 plot([Nsol Nsol], [v(3) v(4)], 'k');
 570 plot([Nmin Nmin], [v(3) v(4)], 'k');
 571 colorbar;
 572 title('W_r_e_a_l', 'fontsize', font_size);
 573 if flag_plot == 0
 574     ylabel('kx (rad/m)', 'fontsize', font_size);
 575 else
 576     ylabel('lambda (km)', 'fontsize', font_size);
 577 end
 578 set(gca, 'LineWidth', [2])
 579 set(gca, 'fontsize', font_sizeN);
 580
 581 subplot(2, 2, 2);
 582 if flag_plot == 0
 583     pcolor(wn(n1:nt), kx(n2:nlon), log10(umax(:,1+njump:njump2(2))));
 584 else
 585     pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(umax(:,1+njump:njump2(2))));
 586 end
 587 shading flat;
 588 v = axis;
 589 hold on;
 590 plot([Nmax Nmax], [v(3) v(4)], 'k');
 591 plot([Nsol Nsol], [v(3) v(4)], 'k');
 592 plot([Nmin Nmin], [v(3) v(4)], 'k');
 593 colorbar;
 594 title('U_r_e_a_l', 'fontsize', font_size);
 595 set(gca, 'LineWidth', [2])
 596 set(gca, 'fontsize', font_sizeN);
 597
 598 subplot(2, 2, 3);
 599 if flag_plot == 0
 600     pcolor(wn(n1:nt), kx(n2:nlon), log10(pmax(:,1+njump:njump2(2))));
 601 else
 602     pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(pmax(:,1+njump:njump2(2))));
 603 end
 604 shading flat;
 605 v = axis;
 606 hold on;
 607 plot([Nmax Nmax], [v(3) v(4)], 'k');
 608 plot([Nsol Nsol], [v(3) v(4)], 'k');
 609 plot([Nmin Nmin], [v(3) v(4)], 'k');
 610 colorbar;
 611 if flag_plot == 0
 612     ylabel('kx (rad/m)', 'fontsize', font_size);
 613     xlabel('omega (rad/s)', 'fontsize', font_size);
 614 else
 615     ylabel('lambda (km)', 'fontsize', font_size);
 616     xlabel('period T (min)', 'fontsize', font_size);
 617 end
 618 title('P_r_e_a_l', 'fontsize', font_size);
 619 set(gca, 'LineWidth', [2])
 620 set(gca, 'fontsize', font_sizeN);
 621
 622 subplot(2, 2, 4);
 623 if flag_plot == 0
 624     pcolor(wn(n1:nt), kx(n2:nlon), log10(rhomax(:,1+njump:njump2(2))));
 625 else
 626     pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(rhomax(:,1+njump:njump2(2))));
 627 end
 628 shading flat;
 629 v = axis;
 630 hold on;
 631 plot([Nmax Nmax], [v(3) v(4)], 'k');
 632 plot([Nsol Nsol], [v(3) v(4)], 'k');
 633 plot([Nmin Nmin], [v(3) v(4)], 'k');
 634 colorbar;
 635 if flag_plot == 0
 636     xlabel('omega (rad/s)', 'fontsize', font_size);
 637 else
 638     xlabel('period T (min)', 'fontsize', font_size);
 639 end
 640 title('Rho_r_e_a_l', 'fontsize', font_size);
 641 set(gca, 'LineWidth', [2])
 642 set(gca, 'fontsize', font_sizeN);
 643 % }}}
 644