[Dart-dev] [3285] group/posters/trunk/elements/DampedInflationFigure.m: figure showing stable (for one month) adaptive inflations - image plot with associated time series above

[thoar at subversion.ucar.edu]

An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080331/5e467f00/attachment.html
-------------- next part --------------
Added: group/posters/trunk/elements/DampedInflationFigure.m
===================================================================
--- group/posters/trunk/elements/DampedInflationFigure.m	                        (rev 0)
+++ group/posters/trunk/elements/DampedInflationFigure.m	2008-03-31 20:07:36 UTC (rev 3285)
@@ -0,0 +1,135 @@
+function h = DampedInflationFigure(level)
+%
+% h = DampedInflationFigure(15);
+%
+% /ptmp/raeder/J/T85_2/01_*/P*.nc      no inflation
+% .............CSL/Damp.9.6/Prior_extract_1-31.nc   damped adaptive prior state space inflation 
+% ..................../62s2_500_all has the obs_space files in both T85_N directories.
+% 01_62 is the last obs_seq file in both cases.
+% Odd #s have 6Z and 12Z, even have 18Z and "24Z", aka 00Z.
+
+if ( exist('get_copy_index') ~=2 )
+   addpath /fs/image/home/thoar/DART/matlab
+end
+
+basedir = '/ptmp/raeder/CSL';
+
+%----------------------------------------------------------------------
+
+varlist = {'P0','PS','T','U','V','Q','CLDLIQ','CLDICE'};
+varlist = {'P0','T','U','V'};
+varlist = {'P0','PS','Q','T','US','VS'};
+
+%----------------------------------------------------------------------
+% Get initial inflation values
+%----------------------------------------------------------------------
+
+ddir = 'Damp.9.6';
+
+inffname = sprintf('%s/%s/Prior_extract_1-31.nc',basedir,ddir);
+
+lon     = getnc(inffname,'lon');
+lat     = getnc(inffname,'lat');
+levels  = getnc(inffname,'lev');
+timearr = getnc(inffname,'time');
+nT      = length(timearr);
+
+% (145,-30) is about Australia
+% (255, 40) is about Boulder
+lon1ind = find(abs(lon-145) == min(abs(lon-145)));
+lat1ind = find(abs(lat +30) == min(abs(lat+ 30)));
+lon2ind = find(abs(lon-255) == min(abs(lon-255)));
+lat2ind = find(abs(lat -40) == min(abs(lat- 40)));
+
+disp(sprintf('lon,lat is %.3f %.3f',lon(lon1ind),lat(lat1ind)))
+disp(sprintf('lon,lat is %.3f %.3f',lon(lon2ind),lat(lat2ind)))
+disp(sprintf('level %d is %f',level,levels(level)))
+
+inf_mean_index = get_copy_index(inffname,'inflation mean');
+corner = [nT, inf_mean_index, -1, -1, level];
+endpnt = [nT, inf_mean_index, -1, -1, level];
+
+for varnum=5:5
+
+   figure(varnum)
+   orient landscape
+   wysiwyg
+
+   % HARDWIRING for the moment ... 'US' is staggered zonal wind ... aka 'U'
+
+   varname = 'U'
+
+   %-------------------------------------------------------------------
+   % Get final inflation values from last timestep
+   %-------------------------------------------------------------------
+
+   finval = getnc(inffname, varlist{varnum}, corner, endpnt);
+
+   maxval = max(finval(:));
+
+   subplot('position',[0.05,0.05,0.85,0.60]);
+   set(gca,'FontSize',18,'FontWeight','normal')
+
+   imagesc(lon,lat,finval,[1,maxval]); hold on;
+   set(gca,'YDir','normal');
+   pos = get(gca,'Position');
+
+   h = colorbar([0.91,0.05,0.02,0.6]);
+   set(h,'FontSize',18);
+   set(get(h,'YLabel'),'String','inflation factor','FontSize',18);
+
+   h1 = worldmap;
+   set(h1,'Color','w','LineWidth',2.0);
+ % text(20,110,sprintf('%s at level %d ',varlist{varnum},levels(level)))
+   text(20,110,sprintf('%s at %.0f hPa, t = 08/01',varname,levels(level)),'FontSize',18,'FontWeight','normal')
+
+   h1 = plot(lon(lon1ind),lat(lat1ind),'kx');
+   set(h1,'MarkerSize',18,'LineWidth',4.0);
+
+   h1 = plot(lon(lon2ind),lat(lat2ind),'wx');
+   set(h1,'MarkerSize',18,'LineWidth',4.0);
+
+   hold off;
+
+   %-------------------------------------------------------------------
+   % Get initial inflation values from first timestep
+   %-------------------------------------------------------------------
+
+   inf_mean_index = get_copy_index(inffname,'inflation mean');
+
+   corner = [1, inf_mean_index, -1, -1, level];
+   endpnt = [1, inf_mean_index, -1, -1, level];
+
+   startval = getnc(inffname, varlist{varnum}, corner, endpnt);
+
+   subplot('position',[0.10,0.75,0.2,0.125])
+   imagesc(lon,lat,startval,[1,maxval]);
+   h  = worldmap;
+   set(h ,'Color',[0.95 0.95 0.95],'LineWidth',1.5);
+   set(gca,'YDir','normal');
+   axis off;
+
+   %-------------------------------------------------------------------
+   % Get Time series for two locations
+   %-------------------------------------------------------------------
+
+   corner = [-1, inf_mean_index, lat1ind, lon1ind, level];
+   endpnt = [-1, inf_mean_index, lat1ind, lon1ind, level];
+   australia = getnc(inffname, varlist{varnum}, corner, endpnt);
+
+   corner = [-1, inf_mean_index, lat2ind, lon2ind, level];
+   endpnt = [-1, inf_mean_index, lat2ind, lon2ind, level];
+   boulder = getnc(inffname, varlist{varnum}, corner, endpnt);
+   
+   subplot('position',[0.40,0.725,0.50,0.175])
+      set(gca,'FontSize',18,'FontWeight','normal')
+      plot(timearr, australia, 'k--', ...
+           timearr, boulder, 'r-', 'LineWidth', 2.0)
+      axis([-Inf Inf 0 maxval])
+      datetick('x',6,'keeplimits','keepticks');
+      h1 = legend('New South Wales','Colorado');
+      set(h1,'Box','off','Location','NorthWest');
+
+end
+
+colormap(hotres)