[Dart-dev] [4326] DART/trunk/models/wrf/matlab: These routines are ANCIENT - they use netCDF variables with a 'state'
nancy at ucar.edu
nancy at ucar.edu
Thu Mar 25 16:01:44 MDT 2010
Revision: 4326
Author: thoar
Date: 2010-03-25 16:01:44 -0600 (Thu, 25 Mar 2010)
Log Message:
-----------
These routines are ANCIENT - they use netCDF variables with a 'state'
variable (i.e. it has not been parsed into prognostic variables). All I
did was replace the now-unsupported getnc() and netcdf_toolbox() functions
with the equivalent snctools() functions. There is no telling if these
ever worked, or if they work now ... I felt compelled to update them in
the face of the new release.
Modified Paths:
--------------
DART/trunk/models/wrf/matlab/correl.m
DART/trunk/models/wrf/matlab/map_wrf_diff.m
DART/trunk/models/wrf/matlab/map_wrf_diff_time_vect.m
DART/trunk/models/wrf/matlab/map_wrf_vect.m
DART/trunk/models/wrf/matlab/rms_cross_time.m
DART/trunk/models/wrf/matlab/stats_wrf_prof_mem.m
DART/trunk/models/wrf/matlab/stats_wrf_prof_vect.m
DART/trunk/models/wrf/matlab/stats_wrf_time_bd.m
DART/trunk/models/wrf/matlab/stats_wrf_time_vect.m
DART/trunk/models/wrf/matlab/stats_wrf_time_vect_point.m
DART/trunk/models/wrf/matlab/w2_wrf.m
-------------- next part --------------
Modified: DART/trunk/models/wrf/matlab/correl.m
===================================================================
--- DART/trunk/models/wrf/matlab/correl.m 2010-03-23 17:14:07 UTC (rev 4325)
+++ DART/trunk/models/wrf/matlab/correl.m 2010-03-25 22:01:44 UTC (rev 4326)
@@ -12,75 +12,86 @@
clear
-fname = 'Prior';
-tlon = getnc(fname, 'XLON');
-we = size(tlon, 2)
-tlat = getnc(fname, 'XLAT');
-sn = size(tlat, 1)
-level = getnc(fname, 'level');
-bt = size(level, 1)
-copy = getnc(fname, 'copy');
-ens_size = size(copy, 1) - 2;
+fname = 'Prior_Diag.nc';
-true_times = getnc(fname, 'time');
-num_true_times = size(true_times, 1)
+if (exist(fname,'file') ~= 2)
+ error('%s does not exist.',fname)
+end
+if (~ nc_isvar(fname,'state'))
+ disp('This is an old-school routine ...')
+ error('%s does not have a "state" variable',fname)
+end
+
+tlon = nc_varget(fname, 'XLON_d01'); we = size( tlon, 2);
+tlat = nc_varget(fname, 'XLAT_d01'); sn = size( tlat, 1);
+level = nc_varget(fname, 'level_d01'); bt = size(level, 1);
+
+ens_size = get_ens_size(fname);
+
vartype = input('Input variable type for correlation, 1=U, 2=V, 3=W, 4=GZ, 5=T, 6=MU, 7=QV, 8=QC, 9=QR: ');
+vartime = input('Input variable time (index): ');
+vari = input('Input variable x location : ');
+varj = input('Input variable y location : ');
-vartime = input('Input variable time : ');
-
-vari = input('Input variable x location : ');
-varj = input('Input variable y location : ');
if vartype == 6
-vark = 1;
+ vark = 1;
else
-vark = input('Input variable z location : ');
+ vark = input('Input variable z location : ');
end
start_var = 1;
-nx = we + 1;
-ny = sn;
+nx = we + 1;
+ny = sn;
if vartype > 1
-start_var = start_var + bt*(we + 1)*sn;
-nx = we;
-ny = sn + 1;
+ start_var = start_var + bt*(we + 1)*sn;
+ nx = we;
+ ny = sn + 1;
end
if vartype > 2
-start_var = start_var + bt*we*(sn + 1);
-nx = we;
-ny = sn;
+ start_var = start_var + bt*we*(sn + 1);
+ nx = we;
+ ny = sn;
end
if vartype > 3
-start_var = start_var + (bt + 1)*we*sn;
+ start_var = start_var + (bt + 1)*we*sn;
end
if vartype > 4
-start_var = start_var + (bt + 1)*we*sn;
+ start_var = start_var + (bt + 1)*we*sn;
end
if vartype > 5
-start_var = start_var + bt*we*sn;
+ start_var = start_var + bt*we*sn;
end
if vartype > 6
-start_var = start_var + we*sn;
+ start_var = start_var + we*sn;
end
if vartype > 7
-start_var = start_var + bt*we*sn*(vartype-7);
+ start_var = start_var + bt*we*sn*(vartype-7);
end
-start_var = start_var + vari + varj*nx + (vark-1)*nx*ny
+start_var = start_var + vari + varj*nx + (vark-1)*nx*ny;
+count = [1 1 1];
+ens_var = zeros(1,ens_size);
+
for imem = 1:ens_size
-ens_var(imem) = getnc(fname, 'state',[vartime imem start_var],[vartime imem start_var],[1 1 1]);
+
+ memstring = sprintf('ensemble member %d',imem);
+ memindex = get_copy_index(fname,memstring);
+
+ start = [vartime memindex start_var] - 1 ;
+ ens_var(imem) = nc_varget(fname, 'state', start, count);
end
-iso = [0.1:0.1:1];
+iso = 0.1:0.1:1;
-% Select field to plot (U, V, W, GZ, T, MU, QV, QC, QR)
+%% Select field to plot (U, V, W, GZ, T, MU, QV, QC, QR)
field_num = input('Input field type, 1=U, 2=V, 3=W, 4=GZ, 5=T, 6=MU, 7=QV, 8=QC, 9=QR: ');
stime = input('Initial time : ');
ftime = input('End time : ');
-% Get level for free atmosphere fields
+%% Get level for free atmosphere fields
if field_num == 6
field_level = 1;
else
@@ -88,113 +99,113 @@
end
start_var = 1;
-nx = we + 1;
-ny = sn;
+nx = we + 1;
+ny = sn;
var_units = 'U (m/s)';
if field_num > 1
-start_var = start_var + bt*(we + 1)*sn;
-nx = we;
-ny = sn + 1;
-var_units = 'V (m/s)';
+ start_var = start_var + bt*(we + 1)*sn;
+ nx = we;
+ ny = sn + 1;
+ var_units = 'V (m/s)';
end
if field_num > 2
-start_var = start_var + bt*we*(sn + 1);
-nx = we;
-ny = sn;
-var_units = 'W (m/s)';
+ start_var = start_var + bt*we*(sn + 1);
+ nx = we;
+ ny = sn;
+ var_units = 'W (m/s)';
end
if field_num > 3
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'GZ (m^2/s^2)';
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'GZ (m^2/s^2)';
end
if field_num > 4
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'T (K)';
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'T (K)';
end
if field_num > 5
-start_var = start_var + bt*we*sn;
-var_units = 'MU (Pa)';
+ start_var = start_var + bt*we*sn;
+ var_units = 'MU (Pa)';
end
if field_num > 6
-start_var = start_var + we*sn;
-var_units = 'QV (kg/kg)';
+ start_var = start_var + we*sn;
+ var_units = 'QV (kg/kg)';
end
if field_num > 7
-start_var = start_var + bt*we*sn*(field_num-7);
-var_units = 'QC (kg/kg)';
+ start_var = start_var + bt*we*sn*(field_num-7);
+ var_units = 'QC (kg/kg)';
end
if field_num > 8
-var_units = 'QR (kg/kg)';
+ var_units = 'QR (kg/kg)';
end
scrsz = get(0,'ScreenSize');
figure('Position',[1 scrsz(4)/2 0.9*scrsz(4) 0.9*scrsz(4)])
- m = ceil(sqrt(ftime-stime+1));
-
+m = ceil(sqrt(ftime-stime+1));
start_var = start_var + nx*ny*(field_level - 1);
-end_var = start_var + nx*ny - 1;
+end_var = start_var + nx*ny - 1;
+pane = 1;
- pane = 1;
-
for itime = stime:ftime
-plot_title = [var_units ' Level: ' num2str(field_level) ' Time: ' num2str(itime)];
+ plot_title = [var_units ' Level: ' num2str(field_level) ' Time: ' num2str(itime)];
-% Extract field
+ %% Extract field
-field_vec = zeros(nx*ny,ens_size);
+ field_vec = zeros(nx*ny,ens_size);
+ correlation = zeros(nx*ny,1);
+ count = [1 1 1];
-for imem = 1:ens_size
+ for imem = 1:ens_size
+ memstring = sprintf('ensemble member %d',imem);
+ memindex = get_copy_index(fname,memstring);
+ start = [itime memindex start_var] - 1;
+ state_vec_prior = nc_varget(fname, 'state', start, count);
+ field_vec(:,imem) = state_vec_prior;
+ end
-state_vec_prior = getnc(fname, 'state',[itime imem start_var],[itime imem end_var],[1 1 1]);
+ iplot = 1;
-field_vec(:,imem) = state_vec_prior;
+ for i=1:nx*ny
+ if ( (cov(ens_var) ~= 0.0) && (cov(field_vec(i,:)) ~= 0.0) )
+ corrmat = corrcoef(ens_var,field_vec(i,:));
+ correlation(i) = corrmat(1,2);
+ else
+ iplot=0;
+ end
+ end
-end
+ if iplot
-iplot = 1;
+ field = reshape(correlation, [nx, ny]);
-for i=1:nx*ny
-if ( (cov(ens_var) ~= 0.0) & (cov(field_vec(i,:)) ~= 0.0) )
- corrmat = corrcoef(ens_var,field_vec(i,:));
- correlation(i) = corrmat(1,2);
-else
- iplot=0;
-end
-end
+ % Plot field
-if iplot
+ subplot(m,m,pane);
-field = reshape(correlation, [nx, ny]);
+ %nc=5
-% Plot field
+ %colormap = (prism(nc))
-subplot(m,m,pane);
+ if field_num > 2
+ [C,h] = contour(tlon,tlat, field', iso );
+ hold on
+ [Cm,hm] = contour (tlon,tlat, field', -iso, ':');
+ plot(tlon(vari+1),tlat(varj+1),'kx','LineWidth',2,'MarkerSize',15)
+ else
+ %[C, h] = contourf(field');
+ [C,h] = contour (field', iso);
+ hold on
+ [Cm,hm] = contour (field', -iso, '--');
+ plot(vari+1,varj+1,'kx','LineWidth',2,'MarkerSize',15)
+ end
+ title(plot_title)
+ %colorbar('vert')
+ clabel(C, h);
+ clabel(Cm, hm);
-%nc=5
+ end
-%colormap = (prism(nc))
+ pane = pane + 1;
-if field_num > 2
- [C,h] = contour(tlon,tlat, field', iso );
- hold on
- [Cm,hm] = contour (tlon,tlat, field', -iso, ':');
- plot(tlon(vari+1),tlat(varj+1),'kx','LineWidth',2,'MarkerSize',15)
-else
- %[C, h] = contourf(field');
- [C,h] = contour (field', iso);
- hold on
- [Cm,hm] = contour (field', -iso, '--');
- plot(vari+1,varj+1,'kx','LineWidth',2,'MarkerSize',15)
end
-title(plot_title)
-%colorbar('vert')
-clabel(C, h);
-clabel(Cm, hm);
-
-end
-
-pane = pane + 1;
-
-end
Modified: DART/trunk/models/wrf/matlab/map_wrf_diff.m
===================================================================
--- DART/trunk/models/wrf/matlab/map_wrf_diff.m 2010-03-23 17:14:07 UTC (rev 4325)
+++ DART/trunk/models/wrf/matlab/map_wrf_diff.m 2010-03-25 22:01:44 UTC (rev 4326)
@@ -10,25 +10,42 @@
% $Revision$
% $Date$
-% Select field to plot (U, V, W, GZ, T, MU, QV, QC, QR)
+%% Select field to plot (U, V, W, GZ, T, MU, QV, QC, QR)
+map_proj = {'lambert', 'ups', 'mercator'};
+
field_num = input('Input field type, 1=U, 2=V, 3=W, 4=GZ, 5=T, 6=MU, 7=QV, 8=QC, 9=QR: ');
+member = input('Input ensemble member index: ');
+itime = input('Time index : ');
-map_proj = {'lambert', 'ups', 'mercator'};
+trfname = 'True_State.nc';
+prfname = 'Prior_Diag.nc';
+pofname = 'Posterior_Diag.nc';
-member = input('Input ensemble member: ');
+if (exist(trfname,'file') ~= 2)
+ error('%s does not exist.',trfname)
+end
+if (exist(prfname,'file') ~= 2)
+ error('%s does not exist.',prfname)
+end
+if (exist(pofname,'file') ~= 2)
+ error('%s does not exist.',pofname)
+end
-itime = input('Time: ');
+if (exist(prfname,'file') ~= 2)
+ error('%s does not exist',prfname)
+end
-fname = 'Prior_Diag';
-tlon = getnc(fname, 'XLON');
-we = size(tlon, 2);
-tlat = getnc(fname, 'XLAT');
-sn = size(tlat, 1);
-level = getnc(fname, 'level');
-bt = size(level, 1);
+if (~ nc_isvar(prfname,'state'))
+ disp('This is an old-school routine ...')
+ error('%s does not have a "state" variable',prfname)
+end
-% Get level for free atmosphere fields
+tlon = nc_varget(prfname, 'XLON_d01'); we = size( tlon, 2);
+tlat = nc_varget(prfname, 'XLAT_d01'); sn = size( tlat, 1);
+level = nc_varget(prfname, 'level_d01'); bt = size(level, 1);
+
+%% Get level for free atmosphere fields
if field_num == 6
field_level = 1;
else
@@ -36,80 +53,85 @@
end
start_var = 1;
-nx = we + 1;
-ny = sn;
+nx = we + 1;
+ny = sn;
var_units = 'U (m/s)';
+
if field_num > 1
-start_var = start_var + bt*(we + 1)*sn;
-nx = we;
-ny = sn + 1;
-var_units = 'V (m/s)';
+ start_var = start_var + bt*(we + 1)*sn;
+ nx = we;
+ ny = sn + 1;
+ var_units = 'V (m/s)';
end
if field_num > 2
-start_var = start_var + bt*we*(sn + 1);
-nx = we;
-ny = sn;
-var_units = 'W (m/s)';
+ start_var = start_var + bt*we*(sn + 1);
+ nx = we;
+ ny = sn;
+ var_units = 'W (m/s)';
end
if field_num > 3
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'GZ (m^2/s^2)';
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'GZ (m^2/s^2)';
end
if field_num > 4
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'T (K)';
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'T (K)';
end
if field_num > 5
-start_var = start_var + bt*we*sn;
-var_units = 'MU (Pa)';
+ start_var = start_var + bt*we*sn;
+ var_units = 'MU (Pa)';
end
if field_num > 6
-start_var = start_var + we*sn;
-var_units = 'QV (kg/kg)';
+ start_var = start_var + we*sn;
+ var_units = 'QV (kg/kg)';
end
if field_num > 7
-start_var = start_var + bt*we*sn*(field_num-7);
-var_units = 'QC (kg/kg)';
+ start_var = start_var + bt*we*sn*(field_num-7);
+ var_units = 'QC (kg/kg)';
end
if field_num > 8
-var_units = 'QR (kg/kg)';
+ var_units = 'QR (kg/kg)';
end
plot_title = [var_units ' Level: ' num2str(field_level) ' Time: ' num2str(itime)];
start_var = start_var + nx*ny*(field_level - 1);
-end_var = start_var + nx*ny - 1;
+end_var = start_var + nx*ny - 1;
-% Extract field
+memstring = sprintf('ensemble member %d',member);
+mmbr_ind = get_copy_index(prfname,memstring);
+true_ind = get_copy_index(trfname,'true state');
-fname = 'True_State';
-state_vec_truth = getnc(fname, 'state',[itime -1 start_var],[itime -1 end_var],[1 1 1]);
+%% Extract field
-fname = 'Prior_Diag';
-state_vec_prior = getnc(fname, 'state',[itime member start_var],[itime member end_var],[1 1 1]);
+count = [itime 1 (end_var-start_var+1)];
-fname = 'Posterior_Diag';
-state_vec_posterior = getnc(fname, 'state',[itime member start_var],[itime member end_var],[1 1 1]);
+start = [itime true_ind start_var ] - 1;
+state_vec_truth = nc_varget(trfname, 'state', start, count);
+start = [itime mmbr_ind start_var ] - 1;
+state_vec_prior = nc_varget(prfname, 'state', start, count);
+state_vec_poste = nc_varget(pofname, 'state', start, count);
-field_vec = state_vec_posterior - state_vec_truth;
-%field_vec = state_vec_posterior - state_vec_prior;
-%field_vec = state_vec_posterior;
+field_vec = state_vec_poste - state_vec_truth;
+%field_vec = state_vec_poste - state_vec_prior;
+%field_vec = state_vec_poste;
field = reshape(field_vec, [nx, ny]);
-% Plot field
+%% Plot field
%nc=5
-%colormap = (prism(nc))
+% colormap = (prism(nc))
if field_num > 2
-[C, h] = contourf(tlon,tlat,field');
+ [C, h] = contourf(tlon, tlat, field');
else
-%[C, h] = contourf(field');
-[C,h] = contour ( field' , [0.5:1:5] );
-hold on
-[Cm,hm] = contour ( field' ,- [0.5:1:5] , 'k:');
+ %[C, h] = contourf(field');
+ [C,h] = contour(field', 0.5:1:5 );
+ hold on
+ [Cm,hm] = contour(field', -0.5:1:5, 'k:');
end
+
title(plot_title)
colorbar('vert')
clabel(C, h);
Modified: DART/trunk/models/wrf/matlab/map_wrf_diff_time_vect.m
===================================================================
--- DART/trunk/models/wrf/matlab/map_wrf_diff_time_vect.m 2010-03-23 17:14:07 UTC (rev 4325)
+++ DART/trunk/models/wrf/matlab/map_wrf_diff_time_vect.m 2010-03-25 22:01:44 UTC (rev 4326)
@@ -10,24 +10,37 @@
% $Revision$
% $Date$
+trfname = 'True_State.nc';
+prfname = 'Prior_Diag.nc';
+pofname = 'Posterior_Diag.nc';
+
+if (exist(trfname,'file') ~= 2)
+ error('%s does not exist.',trfname)
+end
+if (exist(prfname,'file') ~= 2)
+ error('%s does not exist.',prfname)
+end
+if (exist(pofname,'file') ~= 2)
+ error('%s does not exist.',pofname)
+end
+
+if (~ nc_isvar(prfname,'state'))
+ disp('This is an old-school routine ...')
+ error('%s does not have a "state" variable',prfname)
+end
+
% Select field to plot (U, V, W, GZ, T, MU, QV, QC, QR)
field_num = input('Input field type, 1=U, 2=V, 3=W, 4=GZ, 5=T, 6=MU, 7=QV, 8=QC, 9=QR: ');
-member = input('Input ensemble member: ');
+member = input('Input ensemble member index: ');
-fname = 'Prior_Diag';
-tlon = getnc(fname, 'XLON');
-we = size(tlon, 2);
-tlat = getnc(fname, 'XLAT');
-sn = size(tlat, 1);
-level = getnc(fname, 'level');
-bt = size(level, 1);
-true_times = getnc(fname, 'time');
-num_true_times = size(true_times, 1)
+tlon = nc_varget(prfname, 'XLON_d01'); we = size( tlon, 2);
+tlat = nc_varget(prfname, 'XLAT_d01'); sn = size( tlat, 1);
+level = nc_varget(prfname, 'level_d01'); bt = size(level, 1);
- stime = input('Initial time : ');
- ftime = input('End time : ');
+stime = input('Initial time (index): ');
+ftime = input('End time (index): ');
% Get level for free atmosphere fields
if field_num == 6
@@ -37,110 +50,108 @@
end
start_var = 1;
-nx = we + 1;
-ny = sn;
+nx = we + 1;
+ny = sn;
var_units = 'U (m/s)';
-iso = [0.5:1:5];
+iso = 0.5:1:5;
+
if field_num > 1
-start_var = start_var + bt*(we + 1)*sn;
-nx = we;
-ny = sn + 1;
-var_units = 'V (m/s)';
+ start_var = start_var + bt*(we + 1)*sn;
+ nx = we;
+ ny = sn + 1;
+ var_units = 'V (m/s)';
end
if field_num > 2
-start_var = start_var + bt*we*(sn + 1);
-nx = we;
-ny = sn;
-var_units = 'W (m/s)';
-iso = [0.01:0.01:0.1];
+ start_var = start_var + bt*we*(sn + 1);
+ nx = we;
+ ny = sn;
+ var_units = 'W (m/s)';
+ iso = 0.01:0.01:0.1;
end
if field_num > 3
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'GZ (m^2/s^2)';
-iso = [50:50:300];
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'GZ (m^2/s^2)';
+ iso = 50:50:300;
end
if field_num > 4
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'T (K)';
-iso = [0.5:0.5:5];
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'T (K)';
+ iso = 0.5:0.5:5;
end
if field_num > 5
-start_var = start_var + bt*we*sn;
-var_units = 'MU (Pa)';
-iso = [100:100:600];
+ start_var = start_var + bt*we*sn;
+ var_units = 'MU (Pa)';
+ iso = 100:100:600;
end
if field_num > 6
-start_var = start_var + we*sn;
-var_units = 'QV (kg/kg)';
-iso = [0.0001:0.0001:0.001];
+ start_var = start_var + we*sn;
+ var_units = 'QV (kg/kg)';
+ iso = 0.0001:0.0001:0.001;
end
if field_num > 7
-start_var = start_var + bt*we*sn*(field_num-7);
-var_units = 'QC (kg/kg)';
-iso = [0.00001:0.00001:0.0001];
+ start_var = start_var + bt*we*sn*(field_num-7);
+ var_units = 'QC (kg/kg)';
+ iso = 0.00001:0.00001:0.0001;
end
if field_num > 8
-var_units = 'QR (kg/kg)';
-iso = [0.00001:0.00001:0.0001];
+ var_units = 'QR (kg/kg)';
+ iso = 0.00001:0.00001:0.0001;
end
scrsz = get(0,'ScreenSize');
figure('Position',[1 scrsz(4)/2 0.9*scrsz(4) 0.9*scrsz(4)])
- m = ceil(sqrt(ftime-stime+1));
+m = ceil(sqrt(ftime-stime+1));
start_var = start_var + nx*ny*(field_level - 1);
-end_var = start_var + nx*ny - 1;
+end_var = start_var + nx*ny - 1;
- pane = 1;
+pane = 1;
for itime = stime:ftime
-plot_title = [var_units ' Level: ' num2str(field_level) ' Time: ' num2str(itime)];
+ plot_title = [var_units ' Level: ' num2str(field_level) ' Time: ' num2str(itime)];
-% Extract field
+ %% Extract field
+ count = [ 1 1 (end_var-start_var+1)];
+ tstart = [itime true_ind start_var] - 1;
+ pstart = [itime member start_var] - 1;
-fname = 'True_State';
-state_vec_truth = getnc(fname, 'state',[itime -1 start_var],[itime -1 end_var],[1 1 1]);
+ state_vec_truth = nc_varget(trfname, 'state', tstart, count);
+ state_vec_prior = nc_varget(prfname, 'state', pstart, count);
+ state_vec_poste = nc_varget(pofname, 'state', pstart, count);
-fname = 'Prior_Diag';
-state_vec_prior = getnc(fname, 'state',[itime member start_var],[itime member end_var],[1 1 1]);
+ %field_vec = state_vec_prior - state_vec_truth;
+ field_vec = state_vec_poste - state_vec_prior;
+ %field_vec = state_vec_poste;
-fname = 'Posterior_Diag';
-%fname = 'True_Bdy';
-state_vec_posterior = getnc(fname, 'state',[itime member start_var],[itime member end_var],[1 1 1]);
+ field = reshape(field_vec, [nx, ny]);
-%field_vec = state_vec_prior - state_vec_truth;
-field_vec = state_vec_posterior - state_vec_prior;
-%field_vec = state_vec_posterior;
+ % Plot field
-field = reshape(field_vec, [nx, ny]);
+ subplot(m,m,pane);
-% Plot field
+ %nc=5
-subplot(m,m,pane);
+ %colormap = (prism(nc))
+ if field_num > 2
+ [C,h] = contour(tlon,tlat, field', iso );
+ hold on
+ [Cm,hm] = contour (tlon,tlat, field', -iso, ':');
+ else
+ %[C, h] = contourf(field');
+ [C,h] = contour (field', iso);
+ hold on
+ [Cm,hm] = contour (field', -iso, '--');
+ end
+ title(plot_title)
+ %colorbar('vert')
+ clabel(C, h);
+ clabel(Cm, hm);
-%nc=5
+ pane = pane + 1;
-%colormap = (prism(nc))
-if field_num > 2
- [C,h] = contour(tlon,tlat, field', iso );
- hold on
- [Cm,hm] = contour (tlon,tlat, field', -iso, ':');
-else
- %[C, h] = contourf(field');
- [C,h] = contour (field', iso);
- hold on
- [Cm,hm] = contour (field', -iso, '--');
end
-title(plot_title)
-%colorbar('vert')
-clabel(C, h);
-clabel(Cm, hm);
-pane = pane + 1;
-
-end
-
% Loop for another try
%map_wrf;
Modified: DART/trunk/models/wrf/matlab/map_wrf_vect.m
===================================================================
--- DART/trunk/models/wrf/matlab/map_wrf_vect.m 2010-03-23 17:14:07 UTC (rev 4325)
+++ DART/trunk/models/wrf/matlab/map_wrf_vect.m 2010-03-25 22:01:44 UTC (rev 4326)
@@ -10,84 +10,92 @@
% $Revision$
% $Date$
-% Select field to plot (U, V, W, GZ, T, MU, QV, QC, QR)
+%% Select field to plot (U, V, W, GZ, T, MU, QV, QC, QR)
field_num = input('Input field type, 1=U, 2=V, 3=W, 4=GZ, 5=T, 6=MU, 7=QV, 8=QC, 9=QR: ');
% Get file name of true state file
-fname = 'True_State';
-tlon = getnc(fname, 'XLON');
-we = size(tlon, 2);
-tlat = getnc(fname, 'XLAT');
-sn = size(tlat, 1);
-level = getnc(fname, 'level');
-bt = size(level, 1);
+fname = 'True_State';
-state_vec = getnc(fname, 'state');
+if (exist(fname,'file') ~= 2)
+ error('%s does not exist.',fname)
+end
-% Get a time level from the user
+if (~ nc_isvar(fname,'state'))
+ disp('This is an old-school routine ...')
+ error('%s does not have a "state" variable',fname)
+end
+
+tlon = nc_varget(fname, 'XLON_d01'); we = size( tlon, 2);
+tlat = nc_varget(fname, 'XLAT_d01'); sn = size( tlat, 1);
+level = nc_varget(fname, 'level_d01'); bt = size(level, 1);
+
+state_vec = nc_varget(fname, 'state');
+
+%% Get a time level from the user
itime = input('Input time level: ');
single_state = state_vec(itime, :);
-% Get level for free atmosphere fields
+%% Get level for free atmosphere fields
if field_num == 6
field_level = 1;
else
field_level = input('Input level: ');
end
-start_var = 1
-nx = we + 1
-ny = sn
-var_units = 'U (m/s)'
+start_var = 1;
+nx = we + 1;
+ny = sn;
+var_units = 'U (m/s)';
+
if field_num > 1
- start_var = start_var + bt*(we + 1)*sn
- nx = we
- ny = sn + 1
-var_units = 'V (m/s)'
+ start_var = start_var + bt*(we + 1)*sn;
+ nx = we;
+ ny = sn + 1;
+ var_units = 'V (m/s)';
end
if field_num > 2
- start_var = start_var + bt*we*(sn + 1)
- nx = we
- ny = sn
-var_units = 'W (m/s)'
+ start_var = start_var + bt*we*(sn + 1);
+ nx = we;
+ ny = sn;
+ var_units = 'W (m/s)';
end
if field_num > 3
- start_var = start_var + (bt + 1)*we*sn
-var_units = 'GZ (m^2/s^2)'
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'GZ (m^2/s^2)';
end
if field_num > 4
- start_var = start_var + (bt + 1)*we*sn
-var_units = 'T (K)'
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'T (K)';
end
if field_num > 5
- start_var = start_var + bt*we*sn
-var_units = 'MU (Pa)'
+ start_var = start_var + bt*we*sn;
+ var_units = 'MU (Pa)';
end
if field_num > 6
- start_var = start_var + we*sn
-var_units = 'QV (kg/kg)'
+ start_var = start_var + we*sn;
+ var_units = 'QV (kg/kg)';
end
if field_num > 7
- start_var = start_var + bt*we*sn*(field_num-7)
-var_units = 'QC (kg/kg)'
+ start_var = start_var + bt*we*sn*(field_num-7);
+ var_units = 'QC (kg/kg)';
end
if field_num > 8
-var_units = 'QR (kg/kg)'
+ var_units = 'QR (kg/kg)';
end
-plot_title = ['True state ' var_units ' Level: ' num2str(field_level) ' Time: ' num2str(itime)]
+plot_title = ['True state ' var_units ' Level: ' num2str(field_level) ' Time: ' num2str(itime)];
-start_var = start_var + nx*ny*(field_level - 1)
+start_var = start_var + nx*ny*(field_level - 1);
-% Extract field
+%% Extract field
field_vec = single_state(start_var : start_var + nx*ny - 1);
field = reshape(field_vec, [nx, ny]);
-% Plot field
+%% Plot field
%nc=5
Modified: DART/trunk/models/wrf/matlab/rms_cross_time.m
===================================================================
--- DART/trunk/models/wrf/matlab/rms_cross_time.m 2010-03-23 17:14:07 UTC (rev 4325)
+++ DART/trunk/models/wrf/matlab/rms_cross_time.m 2010-03-25 22:01:44 UTC (rev 4326)
@@ -10,132 +10,136 @@
% $Revision$
% $Date$
-fname = 'Prior_Diag';
-tlon = getnc(fname, 'XLON');
-we = size(tlon, 2);
-tlat = getnc(fname, 'XLAT');
-sn = size(tlat, 1);
-level = getnc(fname, 'level');
-bt = size(level, 1);
-copy = getnc(fname, 'copy');
-ens_size = size(copy, 1) - 2;
+prfname = 'Prior_Diag.nc';
+pofname = 'Posterior_Diag.nc';
+trfname = 'True_State.nc';
-true_times = getnc(fname, 'time');
-num_true_times = size(true_times, 1);
-%num_true_times = 27
-
-dx = 200.;
-
-mean_ind = ens_size + 1;
-std_ind = mean_ind + 1;
-
-scrsz = get(0,'ScreenSize');
-figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)])
-
-for field_num = [1:9]
-
-start_var = 1;
-nx = we + 1;
-ny = sn;
-var_units = 'U (m/s)';
-maxlev = bt;
-maxval = 8.0;
-if field_num > 1
-start_var = start_var + bt*(we + 1)*sn ;
-nx = we;
-ny = sn + 1;
-var_units = 'V (m/s)';
+if (exist(trfname,'file') ~= 2)
+ error('%s does not exist.',trfname)
end
-if field_num > 2
-start_var = start_var + bt*we*(sn + 1);
-nx = we;
-ny = sn;
-var_units = 'W (m/s)';
-maxlev = bt + 1;
-maxval = 0.03;
+if (exist(pofname,'file') ~= 2)
+ error('%s does not exist.',pofname)
end
-if field_num > 3
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'GZ (m^2/s^2)';
-maxval = 700.0;
+if (exist(prfname,'file') ~= 2)
+ error('%s does not exist.',prfname)
end
-if field_num > 4
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'T (K)';
-maxlev = bt;
-maxval = 4.0;
+
+if (~ nc_isvar(pofname,'state'))
+ disp('This is an old-school routine ...')
+ error('%s does not have a "state" variable',pofname)
end
-if field_num > 5
-start_var = start_var + bt*we*sn;
-var_units = 'MU (Pa)';
-maxlev = 1;
-maxval = 700.0;
-end
-if field_num > 6
-start_var = start_var + we*sn;
-var_units = 'QV (kg/kg)';
-maxlev = bt;
-maxval = 0.001;
-end
-if field_num > 7
-start_var = start_var + bt*we*sn*(field_num-7);
-var_units = 'QC (kg/kg)';
-maxval = 0.00003;
-end
-if field_num > 8
-var_units = 'QR (kg/kg)';
-maxval = 0.00007;
-end
-x = [1:nx];
-rmse = zeros(nx,num_true_times);
-field1d = zeros(maxlev,1);
+tlon = nc_varget(pofname, 'XLON_d01');
+tlat = nc_varget(pofname, 'XLAT_d01');
+level = nc_varget(pofname, 'level_d01');
+copy = nc_varget(pofname, 'copy');
+true_times = nc_varget(pofname, 'time');
-f_size = nx*ny*maxlev;
+we = size( tlon, 2);
+sn = size( tlat, 1);
+bt = size( level, 1);
+ens_size = size( copy, 1) - 2;
+num_true_times = size(true_times, 1);
-end_var = start_var + f_size - 1;
+dx = 200.0;
-for itime = num_true_times-1:num_true_times
+mean_ind = get_copy_index(pofname,'ensemble mean');
-% Extract field
+scrsz = get(0,'ScreenSize');
+figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)])
-fname = 'True_State';
-field_vec_truth = getnc(fname, 'state',[itime -1 start_var],[itime -1 end_var],[1 1 1]);
+for field_num = 1:9
-fname = 'Prior_Diag';
-field_vec_prior = getnc(fname, 'state',[itime mean_ind start_var],[itime mean_ind end_var],[1 1 1]);
+ start_var = 1;
+ nx = we + 1;
+ ny = sn;
+ var_units = 'U (m/s)';
+ maxlev = bt;
+ maxval = 8.0;
-fname = 'Posterior_Diag';
-field_vec_posterior = getnc(fname, 'state',[itime mean_ind start_var],[itime mean_ind end_var],[1 1 1]);
+ if field_num > 1
+ start_var = start_var + bt*(we + 1)*sn ;
+ nx = we;
+ ny = sn + 1;
+ var_units = 'V (m/s)';
+ end
+ if field_num > 2
+ start_var = start_var + bt*we*(sn + 1);
+ nx = we;
+ ny = sn;
+ var_units = 'W (m/s)';
+ maxlev = bt + 1;
+ maxval = 0.03;
+ end
+ if field_num > 3
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'GZ (m^2/s^2)';
+ maxval = 700.0;
+ end
+ if field_num > 4
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'T (K)';
+ maxlev = bt;
+ maxval = 4.0;
+ end
+ if field_num > 5
+ start_var = start_var + bt*we*sn;
+ var_units = 'MU (Pa)';
+ maxlev = 1;
+ maxval = 700.0;
+ end
+ if field_num > 6
+ start_var = start_var + we*sn;
+ var_units = 'QV (kg/kg)';
+ maxlev = bt;
+ maxval = 0.001;
+ end
+ if field_num > 7
+ start_var = start_var + bt*we*sn*(field_num-7);
+ var_units = 'QC (kg/kg)';
+ maxval = 0.00003;
+ end
+ if field_num > 8
+ var_units = 'QR (kg/kg)';
+ maxval = 0.00007;
+ end
-field_vec = field_vec_posterior - field_vec_truth;
+ x = 1:nx;
+ rmse = zeros(nx,num_true_times);
+ field1d = zeros(maxlev,1);
+ f_size = nx*ny*maxlev;
+ end_var = start_var + f_size - 1;
-field3d = reshape(field_vec, [nx, ny, maxlev]);
+ for itime = num_true_times-1:num_true_times
-for ix = 1:nx
-field1d(:) = field3d(ix,sn/2,:);
-rmse(ix,itime) = sqrt((field1d'*field1d)/maxlev);
-x(ix) = (ix-1)*dx;
-end
+ % Extract field
-end
+ count = [ 1 1 (end_var - start_var + 1)];
+ start = [itime 1 start_var] -1;
+ field_vec_truth = nc_varget(trfname, 'state', start, count);
-subplot(3,3,field_num);
+ start = [itime mean_ind start_var];
+ field_vec_prior = nc_varget(prfname, 'state', start, count);
+ field_vec_posterior = nc_varget(pofname, 'state', start, count);
-plot(x,rmse)
+ field_vec = field_vec_posterior - field_vec_truth;
+ field3d = reshape(field_vec, [nx, ny, maxlev]);
- plot_title = [var_units];
+ for ix = 1:nx
+ field1d(:) = field3d(ix,sn/2,:);
+ rmse(ix,itime) = sqrt((field1d'*field1d)/maxlev);
+ x(ix) = (ix-1)*dx;
+ end
- title(plot_title)
+ end
-% room = (x(2*num_true_times)-x(1))/10;
+ subplot(3,3,field_num);
+ plot(x,rmse)
+ title(var_units)
-% axis ([(x(1)-room) (x(2*num_true_times)+room) 0.0 maxval])
-
end
xlabel('distance (km)')
% Loop for another try
%map_wrf;
-
Modified: DART/trunk/models/wrf/matlab/stats_wrf_prof_mem.m
===================================================================
--- DART/trunk/models/wrf/matlab/stats_wrf_prof_mem.m 2010-03-23 17:14:07 UTC (rev 4325)
+++ DART/trunk/models/wrf/matlab/stats_wrf_prof_mem.m 2010-03-25 22:01:44 UTC (rev 4326)
@@ -10,119 +10,122 @@
% $Revision$
% $Date$
-fname = 'Prior_Diag';
-tlon = getnc(fname, 'XLON');
-we = size(tlon, 2);
-tlat = getnc(fname, 'XLAT');
-sn = size(tlat, 1);
-level = getnc(fname, 'level');
-bt = size(level, 1);
-copy = getnc(fname, 'copy');
-ens_size = size(copy, 1) - 2;
+trfname = 'True_State.nc';
+prfname = 'Prior_Diag.nc';
+pofname = 'Posterior_Diag.nc';
-state_vec_prior = getnc(fname, 'state');
+if (exist(trfname,'file') ~= 2)
+ error('%s does not exist.',trfname)
+end
+if (exist(prfname,'file') ~= 2)
+ error('%s does not exist.',prfname)
+end
+if (exist(pofname,'file') ~= 2)
+ error('%s does not exist.',pofname)
+end
-fname = 'Posterior_Diag';
+if (~ nc_isvar(prfname,'state'))
+ disp('This is an old-school routine ...')
+ error('%s does not have a "state" variable',prfname)
+end
-state_vec_posterior = getnc(fname, 'state');
+tlon = nc_varget(prfname, 'XLON_d01'); we = size( tlon, 2);
+tlat = nc_varget(prfname, 'XLAT_d01'); sn = size( tlat, 1);
+level = nc_varget(prfname, 'level_d01'); bt = size(level, 1);
+times = nc_varget(prfname, 'time'); num_times = size(times, 1);
-fname = 'True_State';
+state_vec_truth = nc_varget(trfname, 'state');
+state_vec_prior = nc_varget(prfname, 'state');
+state_vec_poste = nc_varget(pofname, 'state');
-state_vec_truth = getnc(fname, 'state');
+dim = size(state_vec_truth);
-dim = size (state_vec_truth);
-
len = dim(2);
-true_times = getnc(fname, 'time');
-num_true_times = size(true_times, 1);
+state_vec_truth = reshape(state_vec_truth, [num_times, 1, len]);
-state_vec_truth = reshape(state_vec_truth, [num_true_times, 1, len]);
-
member = input('Input ensemble member: ');
+itime = input('Time (index): ');
-itime = input('Time: ');
+poste_tru = state_vec_poste(itime,member,:) - state_vec_truth(itime,1,:);
+prior_tru = state_vec_prior(itime,member,:) - state_vec_truth(itime,1,:);
-poste_tru = state_vec_posterior(itime,member,:) - state_vec_truth(itime,1,:);
-prior_tru = state_vec_prior(itime,member, :) - state_vec_truth(itime,1,:);
-
scrsz = get(0,'ScreenSize');
figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)])
for field_num = [1:5 7:9]
-start_var = 1
-nx = we + 1
-ny = sn
-var_units = 'U (m/s)'
-maxlev = bt
-if field_num > 1
- start_var = start_var + bt*(we + 1)*sn
- nx = we
- ny = sn + 1
-var_units = 'V (m/s)'
-end
-if field_num > 2
- start_var = start_var + bt*we*(sn + 1)
- nx = we
- ny = sn
-var_units = 'W (m/s)'
-maxlev = bt + 1
-end
-if field_num > 3
- start_var = start_var + (bt + 1)*we*sn
-var_units = 'GZ (m^2/s^2)'
-end
-if field_num > 4
- start_var = start_var + (bt + 1)*we*sn
-var_units = 'T (K)'
-maxlev = bt
-end
-if field_num > 5
- start_var = start_var + bt*we*sn
-var_units = 'MU (Pa)'
-maxlev = 1
-end
-if field_num > 6
- start_var = start_var + we*sn
-var_units = 'QV (kg/kg)'
-maxlev = bt
-end
-if field_num > 7
- start_var = start_var + bt*we*sn*(field_num-7)
-var_units = 'QC (kg/kg)'
-end
-if field_num > 8
-var_units = 'QR (kg/kg)'
-end
+ start_var = 1;
+ nx = we + 1;
+ ny = sn;
+ var_units = 'U (m/s)';
+ maxlev = bt;
-y = [1:maxlev];
-rmse_prior = y;
-rmse_poste = y;
+ if field_num > 1
+ start_var = start_var + bt*(we + 1)*sn;
+ nx = we;
+ ny = sn + 1;
+ var_units = 'V (m/s)';
+ end
+ if field_num > 2
+ start_var = start_var + bt*we*(sn + 1);
+ nx = we;
+ ny = sn;
+ var_units = 'W (m/s)';
+ maxlev = bt + 1;
+ end
+ if field_num > 3
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'GZ (m^2/s^2)';
+ end
+ if field_num > 4
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'T (K)';
+ maxlev = bt;
+ end
+ if field_num > 5
+ start_var = start_var + bt*we*sn;
+ var_units = 'MU (Pa)';
+ maxlev = 1;
+ end
+ if field_num > 6
+ start_var = start_var + we*sn;
+ var_units = 'QV (kg/kg)';
+ maxlev = bt;
+ end
+ if field_num > 7
+ start_var = start_var + bt*we*sn*(field_num-7);
+ var_units = 'QC (kg/kg)';
+ end
+ if field_num > 8
+ var_units = 'QR (kg/kg)';
+ end
-for field_level = 1:maxlev
+ y = 1:maxlev;
+ rmse_prior = y;
+ rmse_poste = y;
-% Extract field
+ for field_level = 1:maxlev
-field_vec = prior_tru(start_var : start_var + nx*ny - 1);
-rmse_prior(field_level) = sqrt((field_vec*field_vec')/(nx*ny));
+ % Extract field
-field_vec = poste_tru(start_var : start_var + nx*ny - 1);
-rmse_poste(field_level) = sqrt((field_vec*field_vec')/(nx*ny));
+ field_vec = prior_tru(start_var : start_var + nx*ny - 1);
+ rmse_prior(field_level) = sqrt((field_vec*field_vec')/(nx*ny));
-start_var = start_var + nx*ny;
+ field_vec = poste_tru(start_var : start_var + nx*ny - 1);
+ rmse_poste(field_level) = sqrt((field_vec*field_vec')/(nx*ny));
-end
+ start_var = start_var + nx*ny;
-subplot(3,3,field_num); plot(rmse_prior,y,rmse_poste,y)
+ end
-plot_title = [var_units]
+ subplot(3,3,field_num);
+ plot(rmse_prior,y,rmse_poste,y)
+ title(var_units)
-title(plot_title)
-
end
-legend('Prior','Posterior')
+legend('Prior','Posterior');
% Loop for another try
%map_wrf;
Modified: DART/trunk/models/wrf/matlab/stats_wrf_prof_vect.m
===================================================================
--- DART/trunk/models/wrf/matlab/stats_wrf_prof_vect.m 2010-03-23 17:14:07 UTC (rev 4325)
+++ DART/trunk/models/wrf/matlab/stats_wrf_prof_vect.m 2010-03-25 22:01:44 UTC (rev 4326)
@@ -10,131 +10,138 @@
% $Revision$
% $Date$
-fname = 'Prior_Diag';
-tlon = getnc(fname, 'XLON');
-we = size(tlon, 2);
-tlat = getnc(fname, 'XLAT');
-sn = size(tlat, 1);
-level = getnc(fname, 'level');
-bt = size(level, 1);
-copy = getnc(fname, 'copy');
-ens_size = size(copy, 1) - 2;
+trfname = 'True_State.nc';
+prfname = 'Prior_Diag.nc';
+pofname = 'Posterior_Diag.nc';
-mean_ind = ens_size + 1;
-std_ind = mean_ind + 1;
-
-itime = input('Time: ');
-
-figure_title = ['Time: ',num2str(itime)];
-
-scrsz = get(0,'ScreenSize');
-figure('Position',[1 scrsz(4)/2 scrsz(3)/2 0.9*scrsz(4)],'Name',figure_title)
-
-for field_num = [1:5 7:9]
-
-start_var = 1;
-nx = we + 1;
-ny = sn;
-var_units = 'U (m/s)';
-maxlev = bt;
-maxval = 8.0;
-if field_num > 1
-start_var = start_var + bt*(we + 1)*sn;
-nx = we;
-ny = sn + 1;
-var_units = 'V (m/s)';
+if (exist(trfname,'file') ~= 2)
+ error('%s does not exist.',trfname)
end
-if field_num > 2
-start_var = start_var + bt*we*(sn + 1);
-nx = we;
-ny = sn;
-var_units = 'W (m/s)';
-maxlev = bt + 1;
-maxval = 0.03;
+if (exist(prfname,'file') ~= 2)
+ error('%s does not exist.',prfname)
end
-if field_num > 3
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'GZ (m^2/s^2)';
-maxval = 700.0;
+if (exist(pofname,'file') ~= 2)
+ error('%s does not exist.',pofname)
end
-if field_num > 4
-start_var = start_var + (bt + 1)*we*sn;
-var_units = 'T (K)';
-maxlev = bt;
-maxval = 4.0;
+
+if (~ nc_isvar(prfname,'state'))
+ disp('This is an old-school routine ...')
+ error('%s does not have a "state" variable',prfname)
end
-if field_num > 5
-start_var = start_var + bt*we*sn;
-var_units = 'MU (Pa)';
-maxlev = 1;
-maxval = 700.0;
-end
-if field_num > 6
-start_var = start_var + we*sn;
-var_units = 'QV (kg/kg)';
-maxlev = bt;
-maxval = 0.001;
-end
-if field_num > 7
-start_var = start_var + bt*we*sn*(field_num-7);
-var_units = 'QC (kg/kg)';
-maxval = 0.00004;
-end
-if field_num > 8
-var_units = 'QR (kg/kg)';
-maxval = 0.00007;
-end
-f_size = nx*ny;
+tlon = nc_varget(prfname, 'XLON_d01'); we = size( tlon, 2);
+tlat = nc_varget(prfname, 'XLAT_d01'); sn = size( tlat, 1);
+level = nc_varget(prfname, 'level_d01'); bt = size(level, 1);
-y = [1:maxlev];
-rmse_prior = y;
-rmse_poste = y;
-spread_prior = y;
-spread_poste = y;
+true_ind = get_copy_index(trfname,'true state');
+mean_ind = get_copy_index(prfname,'ensemble mean');
+sprd_ind = get_copy_index(prfname,'ensemble spread');
-for field_level = 1:maxlev
+itime = input('Time (index): ');
-% Extract field
+figure_title = ['Time: ',num2str(itime)];
-end_var = start_var + nx*ny - 1;
+scrsz = get(0,'ScreenSize');
+figure('Position',[1 scrsz(4)/2 scrsz(3)/2 0.9*scrsz(4)],'Name',figure_title)
-fname = 'True_State';
-state_vec_truth = getnc(fname, 'state',[itime -1 start_var],[itime -1 end_var],[1 1 1]);
+for field_num = [1:5 7:9]
-fname = 'Prior_Diag';
-field_vec_prior = getnc(fname, 'state',[itime mean_ind start_var],[itime mean_ind end_var],[1 1 1]);
+ start_var = 1;
+ nx = we + 1;
+ ny = sn;
+ var_units = 'U (m/s)';
+ maxlev = bt;
+ maxval = 8.0;
-field_vec = getnc(fname, 'state',[itime std_ind start_var],[itime std_ind end_var],[1 1 1]);
+ if field_num > 1
+ start_var = start_var + bt*(we + 1)*sn;
+ nx = we;
+ ny = sn + 1;
+ var_units = 'V (m/s)';
+ end
+ if field_num > 2
+ start_var = start_var + bt*we*(sn + 1);
+ nx = we;
+ ny = sn;
+ var_units = 'W (m/s)';
+ maxlev = bt + 1;
+ maxval = 0.03;
+ end
+ if field_num > 3
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'GZ (m^2/s^2)';
+ maxval = 700.0;
+ end
+ if field_num > 4
+ start_var = start_var + (bt + 1)*we*sn;
+ var_units = 'T (K)';
+ maxlev = bt;
+ maxval = 4.0;
+ end
+ if field_num > 5
+ start_var = start_var + bt*we*sn;
+ var_units = 'MU (Pa)';
+ maxlev = 1;
+ maxval = 700.0;
+ end
+ if field_num > 6
+ start_var = start_var + we*sn;
+ var_units = 'QV (kg/kg)';
+ maxlev = bt;
+ maxval = 0.001;
+ end
+ if field_num > 7
+ start_var = start_var + bt*we*sn*(field_num-7);
+ var_units = 'QC (kg/kg)';
+ maxval = 0.00004;
+ end
+ if field_num > 8
+ var_units = 'QR (kg/kg)';
+ maxval = 0.00007;
+ end
-spread_prior(field_level) = sqrt((field_vec'*field_vec)/(f_size));
+ f_size = nx*ny;
+ y = 1:maxlev;
+ rmse_prior = y;
+ rmse_poste = y;
+ sprd_prior = y;
+ sprd_poste = y;
-fname = 'Posterior_Diag';
-field_vec_posterior = getnc(fname, 'state',[itime mean_ind start_var],[itime mean_ind end_var],[1 1 1]);
+ for field_level = 1:maxlev
-field_vec = getnc(fname, 'state',[itime std_ind start_var],[itime std_ind end_var],[1 1 1]);
+ %% Extract field
-spread_poste(field_level) = sqrt((field_vec'*field_vec)/(f_size));
+ start = [itime true_ind start_var] - 1;
+ count = [ 1 1 nx*ny];
-field_vec = field_vec_prior - state_vec_truth;
-rmse_prior(field_level) = sqrt((field_vec'*field_vec)/(f_size));
+ state_vec_truth = nc_varget(trfname, 'state', start, count);
-field_vec = field_vec_posterior - state_vec_truth;
-rmse_poste(field_level) = sqrt((field_vec'*field_vec)/(f_size));
+ start = [itime mean_ind start_var] - 1;
-start_var = start_var + nx*ny;
+ field_vec_prior = nc_varget(prfname, 'state', start, count);
+ field_vec_poste = nc_varget(pofname, 'state', start, count);
-end
+ start = [itime sprd_ind start_var] - 1;
-subplot(3,3,field_num);
-plot(rmse_prior,y,'b',rmse_poste,y,'g',spread_prior,y,'--b',spread_poste,y,'--g');
+ field_vec = nc_varget(prfname, 'state',start, count);
+ sprd_prior(field_level) = sqrt((field_vec'*field_vec)/(f_size));
+ field_vec = nc_varget(pofname, 'state', start, count);
+ sprd_poste(field_level) = sqrt((field_vec'*field_vec)/(f_size));
-plot_title = [var_units];
+ field_vec = field_vec_prior - state_vec_truth;
+ rmse_prior(field_level) = sqrt((field_vec'*field_vec)/(f_size));
+ field_vec = field_vec_poste - state_vec_truth;
+ rmse_poste(field_level) = sqrt((field_vec'*field_vec)/(f_size));
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list