matlab_ncl 腳本

%Script to plot out 6-panels of WRFOUT data...
  % Wind magnitude and direction, Heat Flux, 2m Air-Skin Temp (Delta-T)
  % Precipitation (RAINNC + RAINC), Moisture Flux, 2m Mixing Ratio
%
% Joseph B. Zambon
% 26-January-2010

clear all
close all

ncl_file = '/he_data/he/jbzambon/projects/jill/ncl_out.nc';
  %contains: lat/lon, u10, v10, times, slp
wrfout_file = '/he_data/he/jbzambon/projects/jill/wrfout_for_domain1';
  %contains: HFX, QFX, TSK, T2, Q2, RAINC, RAINNC
coastlines = '/he_data/he/jbzambon/matlab/coast.mat';
  %coastline data for world
sp = 20;
  %quiver spacing for wind magnitude + direction plots

ncl = netcdf(ncl_file);
lat = ncl{'lat'}(1,:,:);
lon = ncl{'lon'}(1,:,:);
[X,Y] = meshgrid(lon(1,:),lat(:,1));
U10 = ncl{'u10'}(:);
V10 = ncl{'v10'}(:);
mag10 = sqrt(U10.^2 + V10.^2);
times = ncl{'Times'}(:);

wrf = netcdf(wrfout_file);
hfx = wrf{'HFX'}(:);
qfx = wrf{'QFX'}(:);
tsk = wrf{'TSK'}(:);
T2 = wrf{'T2'}(:);
delT = T2-tsk;
mslp = ncl{'slp'}(:);
Q2 = wrf{'Q2'}(:);
RAINC = wrf{'RAINC'}(:);
RAINNC = wrf{'RAINNC'}(:);
precip = RAINC + RAINNC;

load(coastlines);

set(figure(1),'Position',[0 30 1400 900]);
for i=1:129
  %Squeeze variables
  imag = squeeze(mag10(i,:,:));
  iu = squeeze(U10(i,:,:));
  iv = squeeze(V10(i,:,:));
  ihfx = squeeze(hfx(i,:,:));
  idelT = squeeze(delT(i,:,:));
  ihfx = squeeze(hfx(i,:,:));
  iqfx = squeeze(qfx(i,:,:));
  iq2 = squeeze(Q2(i,:,:));
  %Calculate precip
  if (i<=1)
    iprecip = zeros(235,224);
  else
    iprecip = squeeze(precip(i,:,:))-squeeze(precip((i-1),:,:));
  end
  %Subplot 1,1 - Wind Mag and Vectors
     subplot('position',[0.025 0.57 0.3 0.35])
     pcolor(X,Y,imag); shading flat; caxis([0 20]); hold on;
     quiver(X(1:sp:end),Y(1:sp:end),iu(1:sp:end),iv(1:sp:end),'k')
     plot(ncst(:,1),ncst(:,2),'k')
     pclegend([0 20],[0.025 0.5 0.3 0.02]); hold on;
     title(['Wind Magnitude (ms^-^1) and Direction'])
  %Subplot 1,2 - HFX
     subplot('position',[0.35 0.57 0.3 0.35])
     pcolor(X,Y,ihfx); shading flat; caxis([-200 1000]); hold on;
     plot(ncst(:,1),ncst(:,2),'k')
     pclegend([-200 1000],[0.35 0.5 0.3 0.02]); hold on;
     title(['Heat Flux in Wm^2'])
  %Subplot 1,3 - Delta-T
     subplot('position',[0.675 0.57 0.3 0.35])
     pcolor(X,Y,idelT); shading flat; caxis([-10 2]); hold on;
     plot(ncst(:,1),ncst(:,2),'k')
     pclegend([-10 2],[0.675 0.5 0.3 0.02]); hold on;
     title(['Delta-T (Air - Skin) in \circC'])
  %Subplot 2,1 - Precipitation
     subplot('position',[0.025 0.1 0.3 0.35])
     %change colormap white-to-blue
        load('~/matlab/colormaps/green.mat')
     pcolor(X,Y,iprecip); shading flat; caxis([0 20]); hold on;
     colormap(green); freezeColors;
     plot(ncst(:,1),ncst(:,2),'k')
     B = pclegend([0 20],[0.025 0.03 0.3 0.02]); %colorbar('SouthOutside');
     freezeColors(B);
     title(['Hourly Precip in mm'])
  %Subplot 2,2 - QFX
     subplot('position',[0.35 0.1 0.3 0.35])
     colormap(jet)
     pcolor(X,Y,iqfx); shading flat; caxis([-0.0001 0.0006]); hold on;
     plot(ncst(:,1),ncst(:,2),'k')
     pclegend([-0.0001 0.0006],[0.35 0.03 0.3 0.02]);
     title(['Moisture Flux in kgm^2 * 10^-^3'])
  %Subplot 2,3 - Q2
     subplot('position',[0.675 0.1 0.3 0.35])
     pcolor(X,Y,iq2); shading flat; caxis([0 0.025]); hold on;
     plot(ncst(:,1),ncst(:,2),'k')
     pclegend([0 20],[0.675 0.03 0.3 0.02]);
     title(['2m Mixing Ratio in kgkg^-^1 * 10^-^2'])
  mtit(['Forecast Hour: ' num2str((i-1)*3) '     Valid: ' times(i,1:10) ' ' times(i,12:end)],'yoff',+0.025,'FontSize',18)
  screen2png('test.png')
  eval(['! mv test.png ' num2str(i) '.png'])
  clf
end