source: Roms_tools/Run_v2.1/TEST_CASES/plot_vortex.m @ 1

Last change on this file since 1 was 1, checked in by cholod, 13 years ago

import Roms_Agrif

File size: 3.7 KB
Line 
1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2%
3%  Make a plot from the results of the VORTEX test case
4%
5%  Further Information: 
6%  http://www.brest.ird.fr/Roms_tools/
7
8%  This file is part of ROMSTOOLS
9%
10%  ROMSTOOLS is free software; you can redistribute it and/or modify
11%  it under the terms of the GNU General Public License as published
12%  by the Free Software Foundation; either version 2 of the License,
13%  or (at your option) any later version.
14%
15%  ROMSTOOLS is distributed in the hope that it will be useful, but
16%  WITHOUT ANY WARRANTY; without even the implied warranty of
17%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18%  GNU General Public License for more details.
19%
20%  You should have received a copy of the GNU General Public License
21%  along with this program; if not, write to the Free Software
22%  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
23%  MA  02111-1307  USA
24%
25%  Copyright (c) 2005-2006 by Pierrick Penven
26%  e-mail:Pierrick.Penven@ird.fr 
27%
28%  Ref: Penven, P., L. Debreu, P. Marchesiello and J.C. McWilliams,
29%       Application of the ROMS embedding procedure for the Central
30%      California Upwelling System,  Ocean Modelling, 2006.
31%
32%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33clear all
34close all
35%
36% User defined parameters
37%
38vname='temp';
39tindex=11;
40%
41% Caxis depending of the variable name
42%
43if length(vname)>1
44  if vname(1:2)=='te'
45    type='r';
46    ddd=1;
47    cmin=17;
48    dc=0.1;
49    cmax=22;
50    cff=1;
51  elseif vname(1:2)=='ze'
52    type='r';
53    ddd=0;
54    cmin=-100;
55    dc=5;
56    cmax=100;
57    cff=100;
58  elseif vname(1:2)=='ub'
59    type='u';
60    ddd=0;
61    cmin=-50;
62    dc=5;
63    cmax=50;
64    cff=100;
65  elseif vname(1:2)=='vb'
66    type='v';
67    ddd=0;
68    cmin=-50;
69    dc=5;
70    cmax=50;
71    cff=100;
72  end
73else
74  if vname(1)=='u'
75    type='u';
76    ddd=1;
77    cmin=-100;
78    dc=20;
79    cmax=100;
80    cff=100;
81  elseif vname(1)=='v'
82    type='v';
83    ddd=1;
84    cmin=-100;
85    dc=20;
86    cmax=100;
87    cff=100;
88  end
89end
90%
91% Parent
92%
93nc=netcdf('vortex_his.nc');
94N=length(nc('s_rho'));
95time=round(nc{'scrum_time'}(tindex)/(24*3600));
96disp(['Day : ',num2str(time)])
97X=1e-3*nc{'x_rho'}(:);
98Y=1e-3*nc{'y_rho'}(:);
99if ddd==1
100  t1=cff*squeeze(nc{vname}(tindex,N,:,:));
101else
102  t1=cff*squeeze(nc{vname}(tindex,:,:));
103end
104close(nc)
105[Xu,Xv,Xrp]=rho2uvp(X);
106[Yu,Yv,Yrp]=rho2uvp(Y);
107if type=='r'
108  X1=X;
109  Y1=Y;
110elseif type=='u'
111  X1=Xu;
112  Y1=Yu;
113elseif type=='v'
114  X1=Xv;
115  Y1=Yv;
116end
117%
118% Child
119%
120nc=netcdf('vortex_his.nc.1');
121nestvortex=0;
122if ~isempty(nc)
123  nestvortex=1;
124  X=1e-3*nc{'x_rho'}(:);
125  Y=1e-3*nc{'y_rho'}(:);
126  if ddd==1
127    t=cff*squeeze(nc{vname}(tindex,N,:,:));
128  else
129    t=cff*squeeze(nc{vname}(tindex,:,:));
130  end
131  close(nc)
132  [Xu,Xv,Xp]=rho2uvp(X);
133  [Yu,Yv,Yp]=rho2uvp(Y);
134  if type=='r'
135    X2=X;
136    Y2=Y;
137  elseif type=='u'
138    X2=Xu;
139    Y2=Yu;
140  elseif type=='v'
141    X2=Xv;
142    Y2=Yv;
143  end
144  Xbox=cat(1,Xp(1:end,1),  ...
145                Xp(end,1:end)' ,...
146                Xp(end:-1:1,end),...
147                Xp(1,end:-1:1)');
148  Ybox=cat(1,Yp(1:end,1),  ...
149                Yp(end,1:end)' ,...
150                Yp(end:-1:1,end),...
151                Yp(1,end:-1:1)');
152end
153%
154% Plots
155%
156contour(X1,Y1,t1,[cmin:dc:cmax],'r')
157xlabel('X [km]')
158ylabel('Y [km]')
159title([vname,' - Day ',num2str(time)])
160if nestvortex==1
161  hold on
162  contour(X2,Y2,t,[cmin:dc:cmax],'k--')
163  plot(Xbox,Ybox,'k');
164  hold off
165  figure
166  tint=interp2(X1,Y1,t1,X2,Y2,'cubic');
167  tdiff=tint-t;
168  disp(['max difference = ',num2str(max(max(abs(tdiff))),2)])
169  disp(['relative difference = ',...
170        num2str(100*max(max(abs(tdiff)))/max(max(abs(tint))),2),' %'])
171  imagesc(flipud(tdiff))
172  colorbar
173  title(['Parent - Child : ',vname,' - Day ',num2str(time)])
174end
Note: See TracBrowser for help on using the repository browser.