[Dart-dev] [4327] DART/trunk/models/wrf/matlab: These are slightly less ancient, but nonetheless fundamentally untested.
nancy at ucar.edu
nancy at ucar.edu
Thu Mar 25 16:10:15 MDT 2010
Revision: 4327
Author: thoar
Date: 2010-03-25 16:10:15 -0600 (Thu, 25 Mar 2010)
Log Message:
-----------
These are slightly less ancient, but nonetheless fundamentally untested.
They are syntactically OK, logically is another matter. Their origin is
lost - these should be considered as a decent starting point _only_
I simply removed dependencies on unsupported toolboxes and replaced them
with supported calls.
Modified Paths:
--------------
DART/trunk/models/wrf/matlab/compute_pressure.m
DART/trunk/models/wrf/matlab/compute_reflectivity.m
DART/trunk/models/wrf/matlab/get_aux_fields_for_p.m
DART/trunk/models/wrf/matlab/get_aux_fields_for_ref.m
DART/trunk/models/wrf/matlab/get_constants.m
DART/trunk/models/wrf/matlab/interp_to_height.m
DART/trunk/models/wrf/matlab/interp_to_pressure.m
DART/trunk/models/wrf/matlab/map_wrf.m
DART/trunk/models/wrf/matlab/map_wrf_diff_time.m
DART/trunk/models/wrf/matlab/map_wrf_spread_time.m
DART/trunk/models/wrf/matlab/psfc_map.m
DART/trunk/models/wrf/matlab/psfc_movie.m
DART/trunk/models/wrf/matlab/psfc_series.m
DART/trunk/models/wrf/matlab/stats_wrf_prof.m
DART/trunk/models/wrf/matlab/stats_wrf_time.m
-------------- next part --------------
Modified: DART/trunk/models/wrf/matlab/compute_pressure.m
===================================================================
--- DART/trunk/models/wrf/matlab/compute_pressure.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/compute_pressure.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -30,9 +30,8 @@
% $Revision$
% $Date$
- phi_eta = ...
- ( phi(2:end,:,:) - phi(1:end-1,:,:) ) ./ repmat( dnw(:), [1 size(mu)] );
- alpha = -phi_eta ./ repmat( reshape( mu, [1 size(mu)] ), [ length(dnw) 1 1 ] );
+phi_eta = phi(2:end,:,:) - phi(1:end-1,:,:) ) ./ repmat( dnw(:), [1 size(mu)];
+alpha = -phi_eta ./ repmat( reshape( mu, [1 size(mu)] ), [ length(dnw) 1 1 ] );
- % Gas law:
- pres = p0 * ( Rd * theta .* (1 + (Rd/Rv)*qv ) ./ ( p0 * alpha ) ).^gamma ;
+%% Gas law:
+pres = p0 * ( Rd * theta .* (1 + (Rd/Rv)*qv ) ./ ( p0 * alpha ) ).^gamma ;
Modified: DART/trunk/models/wrf/matlab/compute_reflectivity.m
===================================================================
--- DART/trunk/models/wrf/matlab/compute_reflectivity.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/compute_reflectivity.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -21,32 +21,45 @@
% $Revision$
% $Date$
-dief = 0.224;
+dief = 0.224;
+n0r = 8.0e6;
+n0g = 4.0e4;
+n0s = 3.0e6;
+rho_r = 1000.0;
+rho_g = 917.0;
+rho_s = 100.0;
-n0r = 8.0e6; n0g = 4.0e4; n0s = 3.0e6;
-rho_r = 1000.0; rho_g = 917.0; rho_s = 100.0;
-
[Nk Nj Ni] = size(temp) ;
ref = zeros(Nk, Nj, Ni) ;
ar = 7.2e20/(((pi*rho_r)^1.75)*(n0r^0.75)) ;
ag_dry = dief*((rho_g/rho_r)^2)*7.2e20 / (((pi*rho_g)^1.75)*(n0g^0.75)) ;
-% This is appropriate for 10-cm radar.
+
+%% This is appropriate for 10-cm radar.
ag_wet = (7.2e20/(((pi*rho_g)^1.75)*(n0g^0.75)))^0.95 ;
as_wet = 7.2e20/(((pi*rho_s)^1.75)*(n0s^0.75)) ;
as_dry = dief*((rho_s/rho_r)^2)*as_wet ;
-% RAIN
+%% RAIN
precip = rho .* qr ;
-for ii = 1:Ni; for jj = 1:Nj; for kk = 1:Nk;
- if( precip(kk,jj,ii) > 0.0 ) ref(kk,jj,ii) = ref(kk,jj,ii) + ar * (precip(kk,jj,ii)^1.75) ; end
-end; end; end
+for ii = 1:Ni
+for jj = 1:Nj
+for kk = 1:Nk
+ if( precip(kk,jj,ii) > 0.0 )
+ ref(kk,jj,ii) = ref(kk,jj,ii) + ar * (precip(kk,jj,ii)^1.75);
+ end
+end
+end
+end
-% HAIL / GRAUPEL
+%% HAIL / GRAUPEL
precip = rho .* qg ;
-for ii = 1:Ni; for jj = 1:Nj; for kk = 1:Nk;
+
+for ii = 1:Ni
+for jj = 1:Nj
+for kk = 1:Nk
if( precip(kk,jj,ii) > 0.0 )
if ( temp(kk,jj,ii) < 273.15 )
ref(kk,jj,ii) = ref(kk,jj,ii) + ag_dry * (precip(kk,jj,ii)^1.75) ;
@@ -54,11 +67,16 @@
ref(kk,jj,ii) = ref(kk,jj,ii) + ag_wet * (precip(kk,jj,ii)^1.6625) ;
end
end
-end; end; end
+end
+end
+end
-% SNOW
+%% SNOW
precip = rho .* qs ;
-for ii = 1:Ni; for jj = 1:Nj; for kk = 1:Nk;
+
+for ii = 1:Ni
+for jj = 1:Nj
+for kk = 1:Nk
if( precip(kk,jj,ii) > 0.0 )
if ( temp(kk,jj,ii) < 273.15 )
ref(kk,jj,ii) = ref(kk,jj,ii) + as_dry * (precip(kk,jj,ii)^1.75) ;
@@ -66,12 +84,18 @@
ref(kk,jj,ii) = ref(kk,jj,ii) + as_wet * (precip(kk,jj,ii)^1.75) ;
end
end
-end; end; end
+end
+end
+end
-for ii = 1:Ni; for jj = 1:Nj; for kk = 1:Nk;
+for ii = 1:Ni
+for jj = 1:Nj
+for kk = 1:Nk
if( ref(kk,jj,ii) > 0.0 )
ref(kk,jj,ii) = 10.0 * log10(ref(kk,jj,ii)) ;
else
ref(kk,jj,ii) = NaN ;
end
-end; end; end
+end
+end
+end
Modified: DART/trunk/models/wrf/matlab/get_aux_fields_for_p.m
===================================================================
--- DART/trunk/models/wrf/matlab/get_aux_fields_for_p.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/get_aux_fields_for_p.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -5,13 +5,22 @@
%
% Other inputs: T0 = wrf caries theta as deviation from this value
% varargin = {} (if fields come from wrfinput files)
-% {time_index, member_index}
+% {time_index, copy_index}
% (if fields come from DART diagnostic files)
% Outputs: mu = full mu (2d)
% dnw = intervals between w levels (1d)
% phi = full geopotential (3d)
% theta = full theta (3d)
% qv = water-vapor mixing ratio (3d)
+%
+% Example to read a DART netcdf file.
+% filename = 'Prior_Diag.nc';
+% time_index = 1;
+% copy_index = 1;
+% dom_id = 1;
+% T0 = 270.0; % this is an example - I have no idea what it should be.
+% [ mu, dnw, phi, theta, qv ] = ...
+% get_aux_fields_for_p( filename, T0, time_index, copy_index, dom_id )
%% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -23,43 +32,88 @@
% $Revision$
% $Date$
- % Retrieve required fields
- if isempty( varargin )
- % Read all data
- nc = netcdf( filename , 'nowrite' ) ;
+if (exist(filename,'file') ~= 2)
+ error('%s does not exist.',filename)
+end
- mu = nc{'MU'}(:,:) + nc{'MUB'} ;
- dnw = nc{'DNW'}(:) ;
- phi = nc{'PH'}(:,:,:) + nc{'PHB'} ;
- theta = nc{'T'}(:,:,:) + T0 ;
- qv = nc{'QVAPOR'}(:,:,:) ;
- close(nc);
- elseif length( varargin ) == 2
- time_index = varargin{1} ; mem_index = varargin{2} ;
- % note that DART diagnostic files include multiple times, members
- % Read all data
- nc = netcdf( filename , 'nowrite' ) ;
+%% Retrieve required fields
- mu = squeeze(nc{'MU'}(time_index,mem_index, :,:)) + nc{'MUB'} ;
- dnw = nc{'DNW'}(:) ;
- phi = squeeze(nc{'PH'}(time_index,mem_index,:,:,:)) + nc{'PHB'} ;
- theta = squeeze(nc{'T'}(time_index,mem_index,:,:,:)) + T0 ;
- qv = squeeze(nc{'QVAPOR'}(time_index,mem_index,:,:,:)) ;
- close(nc);
- elseif length( varargin ) == 3
- time_index = varargin{1} ; mem_index = varargin{2} ; id = varargin{3} ;
-% note that DART diagnostic files include multiple times, members, domains
- mu = getnc(filename,['MU_d0',int2str(id)],[time_index mem_index -1 -1], ...
- [time_index mem_index -1 -1],[1 1 1 1]) + ...
- getnc(filename,['MUB_d0',int2str(id)],[-1 -1],[-1 -1],[1 1]);
- dnw = getnc(filename,['DNW_d0',int2str(id)],[-1],[-1],[-1]);
- phi = getnc(filename,['PH_d0',int2str(id)],[time_index mem_index -1 -1 -1], ...
- [time_index mem_index -1 -1 -1],[1 1 1 1 1]) + ...
- getnc(filename,['PHB_d0',int2str(id)],[-1 -1 -1],[-1 -1 -1],[1 1 1]);
- theta = getnc(filename,['T_d0',int2str(id)],[time_index mem_index -1 -1 -1], ...
- [time_index mem_index -1 -1 -1],[1 1 1 1 1]) + T0;
- qv = getnc(filename,['QVAPOR_d0',int2str(id)],[time_index mem_index -1 -1 -1], ...
- [time_index mem_index -1 -1 -1],[1 1 1 1 1]);
- else
+if isempty( varargin )
+
+ %% Read all data from 'wrfinput' style files
+
+ varexist(filename,{'MU','MUB','DNW','PH','PHB','T','QVAPOR'})
+
+ dnw = nc_varget(filename, 'DNW');
+ qv = nc_varget(filename,'QVAPOR');
+ theta = nc_varget(filename, 'T') + T0;
+ mu = nc_varget(filename, 'MU') + nc_varget(filename,'MUB');
+ phi = nc_varget(filename, 'PH') + nc_varget(filename,'PHB');
+
+elseif length( varargin ) == 2
+
+ %% Read all data from DART diagnostic file - no domain
+ % note that DART diagnostic files include multiple times, copies
+
+ varexist(filename,{'MU','MUB','DNW','PH','PHB','T','QVAPOR'})
+
+ time_index = varargin{1};
+ copy_index = varargin{2};
+
+ start = [time_index, copy_index, 1, 1] - 1;
+ count = [ 1, 1, -1, -1];
+
+ dnw = nc_varget(filename, 'DNW');
+ qv = nc_varget(filename,'QVAPOR',start,count);
+ theta = nc_varget(filename, 'T',start,count) + T0;
+ mu = nc_varget(filename, 'MU',start,count) + nc_varget(filename,'MUB');
+ phi = nc_varget(filename, 'PH',start,count) + nc_varget(filename,'PHB');
+
+elseif length( varargin ) == 3
+
+ %% note that DART diagnostic files include multiple times, copies, domains
+
+ time_index = varargin{1};
+ copy_index = varargin{2};
+ id = varargin{3};
+
+ varexist(filename,{[ 'MU_d0',int2str(id)], ['MUB_d0',int2str(id)], ...
+ [ 'PH_d0',int2str(id)], ['PHB_d0',int2str(id)], ...
+ [ 'T_d0',int2str(id)], ...
+ [ 'DNW_d0',int2str(id)], ...
+ ['QVAPOR_d0',int2str(id)]} )
+
+ start = [time_index, copy_index, 1, 1] - 1;
+ count = [ 1, 1, -1, -1];
+
+ dnw = nc_varget(filename,['DNW_d0',int2str(id)]);
+ mu = nc_varget(filename,[ 'MU_d0',int2str(id)], start, count) + ...
+ nc_varget(filename,['MUB_d0',int2str(id)]);
+
+ start = [time_index, copy_index, 1, 1, 1] - 1;
+ count = [ 1, 1, -1, -1, -1];
+
+ qv = nc_varget(filename,['QVAPOR_d0',int2str(id)], start, count);
+ theta = nc_varget(filename,[ 'T_d0',int2str(id)], start, count) + T0;
+ phi = nc_varget(filename,[ 'PH_d0',int2str(id)], start, count) + ...
+ nc_varget(filename,[ 'PHB_d0',int2str(id)]);
+
+else
disp('*** Incorrect number of input arguments in get_aux_fields_for_p')
- end
+end
+
+
+function varexist(filename, varnames)
+%% We already know the file exists by this point.
+% Lets check to make sure that file contains all needed variables.
+
+nvars = length(varnames);
+
+for i = 1:nvars
+
+ if ( ~ nc_isvar(filename,varnames{i}) )
+ error('%s is not a variable in %s',varnames{i},filename)
+ end
+
+end
+
Modified: DART/trunk/models/wrf/matlab/get_aux_fields_for_ref.m
===================================================================
--- DART/trunk/models/wrf/matlab/get_aux_fields_for_ref.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/get_aux_fields_for_ref.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -10,6 +10,13 @@
% Outputs: qr = rain water mixing ratio (3d)
% qg = graupel mixing ratio (3d)
% qs = snow mixing ratio (3d)
+%
+% Example to read a DART netcdf file.
+% filename = 'Prior_Diag.nc';
+% time_index = 1;
+% copy_index = 1;
+% dom_id = 1;
+% [ qr, qg, qs ] = get_aux_fields_for_ref( filename, time_index, copy_index, dom_id );
%% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -21,35 +28,72 @@
% $Revision$
% $Date$
- % Retrieve required fields
- if isempty( varargin )
- % Read all data
- nc = netcdf( filename , 'nowrite' ) ;
+if (exist(filename,'file') ~= 2)
+ error('%s does not exist.',filename)
+end
- qr = nc{'QRAIN'}(:,:,:) ;
- qg = nc{'QGRAUP'}(:,:,:) ;
- qs = nc{'QSNOW'}(:,:,:) ;
- close(nc);
- elseif length( varargin ) == 2
- time_index = varargin{1} ; mem_index = varargin{2} ;
- % note that DART diagnostic files include multiple times, members
- % Read all data
- nc = netcdf( filename , 'nowrite' ) ;
+%% Retrieve required fields
- qr = squeeze(nc{'QRAIN'}(time_index,mem_index,:,:,:)) ;
- qg = squeeze(nc{'QGRAUP'}(time_index,mem_index,:,:,:)) ;
- qs = squeeze(nc{'QSNOW'}(time_index,mem_index,:,:,:)) ;
- close(nc);
- elseif length( varargin ) == 3
- time_index = varargin{1} ; mem_index = varargin{2} ; id = varargin{3} ;
-% note that DART diagnostic files include multiple times, members, domains
+if isempty( varargin )
- qr = getnc(filename,['QRAIN_d0',int2str(id)],[time_index mem_index -1 -1 -1], ...
- [time_index mem_index -1 -1 -1],[1 1 1 1 1]);
- qg = getnc(filename,['QGRAUP_d0',int2str(id)],[time_index mem_index -1 -1 -1], ...
- [time_index mem_index -1 -1 -1],[1 1 1 1 1]);
- qs = getnc(filename,['QSNOW_d0',int2str(id)],[time_index mem_index -1 -1 -1], ...
- [time_index mem_index -1 -1 -1],[1 1 1 1 1]);
- else
+ %% Read all data from 'wrfinput' style files
+ varexist(filename, {'QRAIN','QGRAUP','QSNOW'});
+
+ qr = nc_varget(filename, 'QRAIN');
+ qg = nc_varget(filename,'QGRAUP');
+ qs = nc_varget(filename, 'QSNOW');
+
+elseif length( varargin ) == 2
+
+ %% Read all data from DART diagnostic file - no domain
+ % note that DART diagnostic files include multiple times, copies
+
+ time_index = varargin{1};
+ copy_index = varargin{2};
+
+ varexist(filename, {'QRAIN','QGRAUP','QSNOW'});
+
+ start = [time_index, copy_index, 1, 1, 1] - 1;
+ count = [ 1, 1, -1, -1, -1];
+ qr = nc_varget(filename, 'QRAIN',start,count);
+ qg = nc_varget(filename,'QGRAUP',start,count);
+ qs = nc_varget(filename, 'QSNOW',start,count);
+
+elseif length( varargin ) == 3
+
+ %% note that DART diagnostic files include multiple times, copies, domains
+
+ time_index = varargin{1};
+ copy_index = varargin{2};
+ id = varargin{3};
+
+ start = [time_index, copy_index, 1, 1, 1] - 1;
+ count = [ 1, 1, -1, -1, -1];
+
+ varexist(filename, {[ 'QRAIN_d0',int2str(id)], ...
+ ['QGRAUP_d0',int2str(id)], ...
+ [ 'QSNOW_d0',int2str(id)]} );
+
+ qr = nc_varget(filename,[ 'QRAIN_d0',int2str(id)], start, count);
+ qg = nc_varget(filename,['QGRAUP_d0',int2str(id)], start, count);
+ qs = nc_varget(filename,[ 'QSNOW_d0',int2str(id)], start, count);
+
+else
disp('*** Incorrect number of input arguments in get_aux_fields_for_p')
- end
+end
+
+
+function varexist(filename, varnames)
+%% We already know the file exists by this point.
+% Lets check to make sure that file contains all needed variables.
+
+nvars = length(varnames);
+
+for i = 1:nvars
+
+ if ( ~ nc_isvar(filename,varnames{i}) )
+ error('%s is not a variable in %s',varnames{i},filename)
+ end
+
+end
+
Modified: DART/trunk/models/wrf/matlab/get_constants.m
===================================================================
--- DART/trunk/models/wrf/matlab/get_constants.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/get_constants.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -13,12 +13,13 @@
% $Revision$
% $Date$
- %--Useful constants
- Rd = 287.0;
- Cp = 7.0*Rd/2.0;
- gamma = Cp / (Cp - Rd) ;
- Rv = 461;
- g = 9.81;
- L_c = 2.25e6;
- T0 = 300;
- p0 = 1000.e2;
+%% Useful constants
+
+Rd = 287.0;
+Cp = 7.0*Rd/2.0;
+gamma = Cp / (Cp - Rd) ;
+Rv = 461;
+g = 9.81;
+L_c = 2.25e6;
+T0 = 300;
+p0 = 1000.e2;
Modified: DART/trunk/models/wrf/matlab/interp_to_height.m
===================================================================
--- DART/trunk/models/wrf/matlab/interp_to_height.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/interp_to_height.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -17,26 +17,29 @@
% $Revision$
% $Date$
- [Nk Nj Ni] = size(heights) ;
+[Nk Nj Ni] = size(heights);
+below = zeros(Nj,Ni);
+var_below = zeros(Nj,Ni);
+dz = zeros(Nj,Ni);
+dvar = zeros(Nj,Ni);
- below = zeros(Nj,Ni) ; var_below = zeros(Nj,Ni) ;
- dz = zeros(Nj,Ni) ; dvar = zeros(Nj,Ni) ;
+% % at each horiz. location, find highest height beneath level
+for ii = 1:Ni
+for jj = 1:Nj
- % at each horiz. location, find highest height beneath level
- for ii = 1:Ni; for jj = 1:Nj;
% kk(jj,ii) = max( find( heights(:,jj,ii) - level < 0 ) ) ;
kk = max( find( heights(:,jj,ii) - level < 0 ) ) ;
% if isempty(kk); disp([ kk , jj, ii ] ); end
if isempty(kk); kk = 1; end % level is below surface here
- below(jj,ii) = heights(max(1 ,kk ),jj,ii) ; % height below level
- dz(jj,ii) = heights(min(Nk,kk+1),jj,ii) ...
- - heights(max(1 ,kk ),jj,ii) ;
- var_below(jj,ii) = var_in(max(1 ,kk ),jj,ii) ;
- dvar (jj,ii) = var_in(min(Nk,kk+1),jj,ii) - var_in(max(1 ,kk ),jj,ii) ;
- end; end
+ below(jj,ii) = heights(max(1 ,kk ),jj,ii) ; % height below level
+ dz(jj,ii) = heights(min(Nk,kk+1),jj,ii) - heights(max(1 ,kk ),jj,ii) ;
+ var_below(jj,ii) = var_in( max(1 ,kk ),jj,ii) ;
+ dvar (jj,ii) = var_in( min(Nk,kk+1),jj,ii) - var_in(max(1 ,kk ),jj,ii) ;
+end
+end
- var_interp = ( dvar ./ dz ) .* ( level - below ) + var_below ;
+var_interp = ( dvar ./ dz ) .* ( level - below ) + var_below ;
- var_interp( level < heights(1,:,:) ) = NaN ;
- % level is beneath lowest mass level
+var_interp( level < heights(1,:,:) ) = NaN ;
+% level is beneath lowest mass level
Modified: DART/trunk/models/wrf/matlab/interp_to_pressure.m
===================================================================
--- DART/trunk/models/wrf/matlab/interp_to_pressure.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/interp_to_pressure.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -13,26 +13,28 @@
% $Revision$
% $Date$
- [Nk Nj Ni] = size(pressure) ;
+[Nk Nj Ni] = size(pressure) ;
+below = zeros(Nj,Ni);
+var_below = zeros(Nj,Ni);
+dlogp = zeros(Nj,Ni);
+dvar = zeros(Nj,Ni);
- below = zeros(Nj,Ni) ; var_below = zeros(Nj,Ni) ;
- dlogp = zeros(Nj,Ni) ; dvar = zeros(Nj,Ni) ;
-
- % at each horiz. location, find highest pressure level beneath p_level
- for ii = 1:Ni; for jj = 1:Nj;
+% at each horiz. location, find highest pressure level beneath p_level
+for ii = 1:Ni
+for jj = 1:Nj
% kk(jj,ii) = max( find( pressure(:,jj,ii) - p_level > 0 ) ) ;
kk = max( find( pressure(:,jj,ii) - p_level > 0 ) ) ;
% if isempty(kk); disp([ kk , jj, ii ] ); end
if isempty(kk); kk = 1; end % p_level is below surface here
- below(jj,ii) = log( pressure(max(1 ,kk ),jj,ii) ) ; % pressure level below p_level
- dlogp(jj,ii) = log( pressure(min(Nk,kk+1),jj,ii) ) ...
- - log( pressure(max(1 ,kk ),jj,ii) ) ;
- var_below(jj,ii) = var_in(max(1 ,kk ),jj,ii) ;
- dvar (jj,ii) = var_in(min(Nk,kk+1),jj,ii) - var_in(max(1 ,kk ),jj,ii) ;
- end; end
+ below(jj,ii) = log( pressure(max(1 ,kk ),jj,ii) ) ; % pressure level below p_level
+ dlogp(jj,ii) = log( pressure(min(Nk,kk+1),jj,ii) ) - below(jj,ii);
+ var_below(jj,ii) = var_in(max(1 ,kk ),jj,ii);
+ dvar(jj,ii) = var_in(min(Nk,kk+1),jj,ii) - var_below(jj,ii);
+end
+end
- var_interp = ( dvar ./ dlogp ) .* ( log(p_level) - below ) + var_below ;
+var_interp = ( dvar ./ dlogp ) .* ( log(p_level) - below ) + var_below ;
- var_interp( p_level > pressure(1,:,:) ) = NaN ;
- % p_level is beneath surface
+var_interp( p_level > pressure(1,:,:) ) = NaN ;
+ % p_level is beneath surface
Modified: DART/trunk/models/wrf/matlab/map_wrf.m
===================================================================
--- DART/trunk/models/wrf/matlab/map_wrf.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/map_wrf.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -1,14 +1,15 @@
function field = map_wrf(fname, varname, copystring, levelindx, timeindx )
-%% map_wrf simply creates a map of a WRF field.
+%% map_wrf creates a map of a WRF field - please check what is supposed to
+% happen with the staggering.
%
+% (rather slow) Example:
% fname = 'Prior_Diag.nc';
% varname = 'U_d01';
% copystring = 'ensemble mean';
% levelindx = 10;
% timeindx = 1;
%
-% Example:
-% map_wrf(fname, varname, copystring, levelindx, timeindx )
+% map_wrf(fname, varname, copystring, levelindx, timeindx );
%% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -32,10 +33,10 @@
if (exist(fname,'file') ~= 2)
error('%s does not exist',fname)
end
-if ( ~ nc_isvar(fname,varname) )
- error('%s is not a variable in %s',varname,fname)
-end
+varexist(fname,{varname,'time','XLON_d01','XLAT_d01','level_d01', ...
+ 'TRUELAT1','TRUELAT2', 'CEN_LAT', 'CEN_LON', 'MAP_PROJ'})
+
copyindex = get_copy_index(fname,copystring);
%% Read map metadata
@@ -43,6 +44,7 @@
xlon = nc_varget(fname, 'XLON_d01'); we = size( xlon, 2);
xlat = nc_varget(fname, 'XLAT_d01'); sn = size( xlat, 1);
level = nc_varget(fname,'level_d01'); bt = size(level, 1);
+
stdlat1 = nc_varget(fname, 'TRUELAT1');
stdlat2 = nc_varget(fname, 'TRUELAT2');
cen_lat = nc_varget(fname, 'CEN_LAT');
@@ -127,29 +129,51 @@
% subplot(m,m,pane);
- axesm(map_proj{mp},'Origin',[0 cen_lon 0],'MapParallels',[stdlat1 stdlat2],...
- 'MapLatLimit',[minlat maxlat],'MapLonLimit',[minlon maxlon]);
+ axesm(map_proj{mp}, 'Origin', [0 cen_lon 0], ...
+ 'MapParallels',[stdlat1 stdlat2],...
+ 'MapLatLimit',[minlat maxlat], ...
+ 'MapLonLimit',[minlon maxlon]);
% framem;
- states = shaperead('usastatelo', 'UseGeoCoords', true);
- geoshow([states.Lat],[states.Lon],'Color','black');
- % linem(states.Lon, states.Lat, 'color', [0 0 0]);
+% states = shaperead('usastatelo', 'UseGeoCoords', true);
+% geoshow([states.Lat],[states.Lon],'Color','black');
+
+ landareas = shaperead('landareas', 'UseGeoCoords', true);
+ geoshow([landareas.Lat],[landareas.Lon],'Color','black');
- % axis( [-0.6 0.6 .2 1.2 ]) % This works pretty well for present CONUS domain
-
if min(field(:)) ~= max(field(:))
- disp('contouring ...')
+% disp('contouring (slow) ...')
+% [C h]=contourm(xlat,xlon,field) ;
+% h = clabelm(C,h,'labelspacing',288); set(h,'Fontsize',9);
- [C h]=contourm(xlat,xlon,field) ;
- h = clabelm(C,h,'labelspacing',288); set(h,'Fontsize',9);
+ h1 = surfm(xlat,xlon,field);
- % Can get lat,lon lines via, e.g., contourm(xlat,xlon,xlat)
- % [C h]=contourm(xlat,xlon,xlat,'k') ;
- % h = clabelm(C,h,'labelspacing',288); set(h,'Fontsize',9);
-
title(plot_title,'Interpreter','none')
- %colorbar('vert')
+ colorbar('vert')
+ else
+ error('field is all uniform value')
end
end
+
+
+
+function varexist(filename, varnames)
+%% We already know the file exists by this point.
+% Lets check to make sure that file contains all needed variables.
+
+nvars = length(varnames);
+gotone = ones(1,nvars);
+
+for i = 1:nvars
+ gotone(i) = nc_isvar(filename,varnames{i});
+ if ( ~ gotone(i) )
+ fprintf('\n%s is not a variable in %s\n',varnames{i},filename)
+ end
+end
+
+if ~ all(gotone)
+ error('missing required variable ... exiting')
+end
+
Modified: DART/trunk/models/wrf/matlab/map_wrf_diff_time.m
===================================================================
--- DART/trunk/models/wrf/matlab/map_wrf_diff_time.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/map_wrf_diff_time.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -12,10 +12,10 @@
field_name = input('Input field type (U, V, W, PH, T, MU, QV, QC, QR, XLAND, VECT, HDIV): ');
-map_proj = {'lambert', 'ups', 'mercator'};
+map_proj = {'lambert', 'ups', 'mercator'};
state_name = {'True_State' 'Prior_Diag' 'Posterior_Diag'};
-disp('Plotting horizontal map of a linear combinaison of'); disp (state_name)
+disp('Plotting horizontal map of a linear combination of'); disp (state_name)
fact = input('Input coefficients (between [ ]):');
@@ -29,12 +29,18 @@
error('Nothing to plot.')
end
-dx = getnc(fname, 'DX');
-stdlat1 = getnc(fname, 'TRUELAT1');
-stdlat2 = getnc(fname, 'TRUELAT2');
-cen_lon = getnc(fname, 'CEN_LON');
-mp = getnc(fname, 'MAP_PROJ');
+if (exist(fname,'file') ~= 2)
+ error('%s does not exist.',fname)
+end
+truelat1 = nc_attget(fname, nc_global, 'TRUELAT1');
+truelat2 = nc_attget(fname, nc_global, 'TRUELAT2');
+cen_lat = nc_attget(fname, nc_global, 'CEN_LAT');
+cen_lon = nc_attget(fname, nc_global, 'CEN_LON');
+stand_lon = nc_attget(fname, nc_global, 'STAND_LON');
+mp = nc_attget(fname, nc_global, 'MAP_PROJ');
+dx = nc_attget(fname, nc_global, 'DX' );
+
num_domains = size(dx,1);
if (num_domains > 1)
@@ -48,24 +54,21 @@
end
-xlon = getnc(fname, ['XLON_d0',int2str(id)]);
-we = size(xlon, 2);
-xlat = getnc(fname, ['XLAT_d0',int2str(id)]);
-sn = size(xlat, 1);
-level = getnc(fname, ['level_d0',int2str(id)]);
-bt = size(level, 1);
+xlon = nc_varget(fname, [ 'XLON_d0',int2str(id)]); we = size( xlon, 2);
+xlat = nc_varget(fname, [ 'XLAT_d0',int2str(id)]); sn = size( xlat, 1);
+level = nc_varget(fname, ['level_d0',int2str(id)]); bt = size(level, 1);
cop = zeros(1,3);
cop(1) = -1;
-if (fact(2) ~= 0.0) | (fact(3) ~= 0.0)
+if (fact(2) ~= 0.0) || (fact(3) ~= 0.0)
if fact(2) ~= 0.0
fname = char(state_name(2));
else
fname = char(state_name(3));
end
- ncopy = size(getnc(fname, 'copy'), 1);
- copy_name = deblank(getnc(fname, 'CopyMetaData'));
+ ncopy = size(nc_varget(fname, 'copy'), 1);
+ copy_name = deblank(nc_varget(fname, 'CopyMetaData'));
for icopy = 1:ncopy
disp(['Copy # ' int2str(icopy) ' is ' copy_name(icopy,:)])
end
@@ -75,41 +78,44 @@
string2 = sprintf('True state');
end
+%% FIXME ... this is horrid! use vector _something_ ...
+
for iy = 1:sn
- for ix = 1:we
- if(xlon(iy,ix) > 0.0)
+for ix = 1:we
+ if(xlon(iy,ix) > 0.0)
xlon(iy,ix) = xlon(iy,ix) - 360.0;
- end
end
end
+end
minlat = min(xlat(:)); maxlat = max(xlat(:));
minlon = min(xlon(:)); maxlon = max(xlon(:));
-true_times = getnc(fname, 'time');
-num_true_times = size(true_times, 1)
+true_times = nc_varget(fname, 'time');
+num_true_times = size(true_times, 1);
stime = input('Initial time : ');
ftime = input('End time : ');
-% Get level for free atmosphere fields
-if strcmp(field_name,'MU') | strcmp(field_name,'XLAND')
+%% Get level for free atmosphere fields
+
+if strcmp(field_name,'MU') || strcmp(field_name,'XLAND')
field_level = 1;
vert_coord = 1;
lev_units = '';
else
vert_coord = 0;
- while vert_coord ~= 1 & vert_coord ~= 2 & vert_coord ~= 3
+ while vert_coord ~= 1 && vert_coord ~= 2 && vert_coord ~= 3
vert_coord = input('Input vertical coordinate: model(1), pressure(2), height(3): ');
if vert_coord == 1
field_level = input('Input model level: ');
- lev_units = 'model level';
+ lev_units = 'model level';
elseif vert_coord == 2
field_level = input('Input pressure level (hPa): ');
- lev_units = 'hPa';
+ lev_units = 'hPa';
elseif vert_coord == 3
field_level = input('Input height (meters): ');
- lev_units = 'm';
+ lev_units = 'm';
else
disp('Invalid choice. Please try again')
end
@@ -120,65 +126,67 @@
uname = ['U_d0',int2str(id)];
vname = ['V_d0',int2str(id)];
-if strcmp(field_name,'U') | strcmp(field_name,'V') | strcmp(field_name,'VECT')
+if strcmp(field_name,'U') || strcmp(field_name,'V') || strcmp(field_name,'VECT')
var_units = 'm/s';
- iso = [0.5:1:5];
+ iso = 0.5:1:5;
end
if strcmp(field_name,'W')
var_units = 'm/s';
- iso = [0.01:0.01:0.1];
+ iso = 0.01:0.01:0.1;
end
if strcmp(field_name,'PH')
var_units = 'm^2/s^2';
- iso = [50:50:300];
+ iso = 50:50:300;
end
if strcmp(field_name,'T')
var_units = 'K';
- iso = [250:5:310];
+ iso = 250:5:310;
end
if strcmp(field_name,'MU')
var_units = 'Pa';
- iso = [100:100:600];
+ iso = 100:100:600;
end
if strcmp(field_name,'QV')
var_units = 'kg/kg';
- iso = [0.0001:0.0001:0.001];
+ iso = 0.0001:0.0001:0.001;
end
if strcmp(field_name,'QC')
var_units = 'kg/kg';
- iso = [0.00001:0.00001:0.0001];
+ iso = 0.00001:0.00001:0.0001;
end
if strcmp(field_name,'QR')
var_units = 'kg/kg';
- iso = [0.00001:0.00001:0.0001];
+ iso = 0.00001:0.00001:0.0001;
end
if strcmp(field_name,'XLAND')
var_units = '-';
end
if strcmp(field_name,'HDIV')
var_units = '(m/s)/km';
- iso = [-0.5:0.05:0.5];
+ iso = -0.5:0.05:0.5;
end
if strcmp(field_name,'REF')
var_units = 'dBZ';
- iso = [-10:5:80];
+ iso = -10:5:80;
end
m = ceil(sqrt(ftime-stime+1));
-%iso = [-3:0.25:3];
-%iso = [0.1:0.2:5.25];
-%iso = [0.01:0.02:5.25];
+%iso = -3:0.25:3;
+%iso = 0.1:0.2:5.25;
+%iso = 0.01:0.02:5.25;
pane = 1;
for itime = stime:ftime
- plot_title = sprintf('%s (%s) %i (%s) day = %i sec = %i %s',field_name,var_units,field_level,lev_units,fix(true_times(itime)),fix((true_times(itime)-fix(true_times(itime)))*86400),string2);
+ plot_title = sprintf('%s (%s) %i (%s) day = %i sec = %i %s', ...
+ field_name, var_units, field_level, lev_units, fix(true_times(itime)), ...
+ fix((true_times(itime)-fix(true_times(itime)))*86400), string2);
-% Extract field
+ %% Extract field
- field = zeros(sn,we);
+ field = zeros(sn,we);
fieldu = zeros(sn,we);
fieldv = zeros(sn,we);
@@ -187,58 +195,60 @@
if fact(istate) ~= 0.0
fname = char(state_name(istate));
- %--Set up, compute pressure
+
+ %% Set up, compute pressure
[ Cp, Rd, gamma, Rv, L_c, g, T0, p0 ] = get_constants ;
[ mu, dnw, phi, theta, qv ] = ...
- get_aux_fields_for_p( fname, T0, itime, cop(istate), id ) ;
+ get_aux_fields_for_p( fname, T0, itime, cop(istate), id ) ;
pres = compute_pressure( mu, dnw, phi, theta, qv, Rd,Rv,gamma,p0 ) ;
var_in = zeros(bt,sn,we);
- %--Retrieve specified variable from netcdf file
+ %--Retrieve specified variable from netcdf file
if vert_coord == 1
if strcmp(field_name,'MU')
- corner = [itime cop(istate) -1 -1];
- end_point = [itime cop(istate) -1 -1];
- stride = [1 1 1 1];
+ start = [itime cop(istate) 1 1] - 1;
+ count = [ 1 1 -1 -1];
elseif strcmp(field_name,'XLAND')
- corner = [-1 -1];
- end_point = [-1 -1];
- stride = [1 1];
+ start = [ 1 1];
+ count = [-1 -1];
else
- corner = [itime cop(istate) field_level -1 -1];
- end_point = [itime cop(istate) field_level -1 -1];
- stride = [1 1 1 1 1];
+ start = [itime cop(istate) field_level 1 1] - 1;
+ count = [ 1 1 1 -1 -1];
end
if strcmp(field_name,'VECT')
- stag_field = getnc(fname,uname,corner,end_point,stride);
+
+ stag_field = nc_varget(fname,uname,start,count);
for iy = 1:sn
for ix = 1:we
fieldu(iy,ix) = fieldu(iy,ix) + ...
- fact(istate)*(stag_field(iy,ix) + stag_field(iy,ix+1))/2.0;
+ fact(istate)*(stag_field(iy,ix) + ...
+ stag_field(iy,ix+1))/2.0;
end
end
- stag_field = getnc(fname,vname,corner,end_point,stride);
+ stag_field = nc_varget(fname,vname,start,count);
for iy = 1:sn
for ix = 1:we
fieldv(iy,ix) = fieldv(iy,ix) + ...
fact(istate)*(stag_field(iy,ix) + stag_field(iy+1,ix))/2.0;
end
end
+
elseif strcmp(field_name,'HDIV')
- stag_field = getnc(fname,uname,corner,end_point,stride);
+
+ stag_field = nc_varget(fname,uname,start,count);
for iy = 1:sn
for ix = 1:we
field(iy,ix) = field(iy,ix) - ...
fact(istate)*(stag_field(iy,ix+1) - stag_field(iy,ix))/dx(id);
end
end
- stag_field = getnc(fname,vname,corner,end_point,stride);
+ stag_field = nc_varget(fname,vname,start,count);
for iy = 1:sn
for ix = 1:we
field(iy,ix) = field(iy,ix) - ...
@@ -246,9 +256,10 @@
end
end
field = field*1000.0;
+
else
- stag_field = getnc(fname, var_name,corner,end_point,stride);
+ stag_field = nc_varget(fname, var_name,start,count);
if strcmp(field_name,'U')
for iy = 1:sn
for ix = 1:we
@@ -271,11 +282,14 @@
elseif vert_coord == 2
+ start = [itime cop(istate) 1 1 1] -1;
+ count = [ 1 1 -1 -1 -1];
+
if strcmp(field_name,'VECT')
var_inu = zeros(bt,sn,we);
var_inv = zeros(bt,sn,we);
- stag_field = getnc(fname,uname,[itime cop(istate) -1 -1 -1], ...
- [itime cop(istate) -1 -1 -1],[1 1 1 1 1]);
+
+ stag_field = nc_varget(fname, uname, start, count);
for iz = 1:bt
for iy = 1:sn
for ix = 1:we
@@ -283,8 +297,8 @@
end
end
end
- stag_field = getnc(fname,vname,[itime cop(istate) -1 -1 -1], ...
- [itime cop(istate) -1 -1 -1],[1 1 1 1 1]);
+
+ stag_field = nc_varget(fname, vname, start, count);
for iz = 1:bt
for iy = 1:sn
for ix = 1:we
@@ -292,9 +306,10 @@
end
end
end
+
elseif strcmp(field_name,'HDIV')
- stag_field = getnc(fname,uname,[itime cop(istate) -1 -1 -1], ...
- [itime cop(istate) -1 -1 -1],[1 1 1 1 1]);
+
+ stag_field = nc_varget(fname, uname, start, count);
for iz = 1:bt
for iy = 1:sn
for ix = 1:we
@@ -303,8 +318,8 @@
end
end
end
- stag_field = getnc(fname,vname,[itime cop(istate) -1 -1 -1], ...
- [itime cop(istate) -1 -1 -1],[1 1 1 1 1]);
+
+ stag_field = nc_varget(fname, vname, start, count);
for iz = 1:bt
for iy = 1:sn
for ix = 1:we
@@ -322,8 +337,7 @@
elseif strcmp(field_name,'T') % Same here
stag_field = theta ;
else
- stag_field = getnc(fname,var_name,[itime cop(istate) -1 -1 -1], ...
- [itime cop(istate) -1 -1 -1],[1 1 1 1 1]);
+ stag_field = nc_varget(fname, var_name, start, count);
end
if strcmp(field_name,'U')
@@ -342,7 +356,7 @@
end
end
end
- elseif strcmp(field_name,'W') | strcmp(field_name,'PH')
+ elseif strcmp(field_name,'W') || strcmp(field_name,'PH')
for iz = 1:bt
for iy = 1:sn
for ix = 1:we
@@ -374,15 +388,12 @@
height = compute_height( phi, g ) ;
- if strcmp(field_name,'REF')
+ if strcmp(field_name,'REF')
- [ qr, qg, qs ] = ...
- get_aux_fields_for_ref( fname, itime, cop(istate), id ) ;
+ [ qr, qg, qs ] = get_aux_fields_for_ref( fname, itime, cop(istate), id ) ;
- rho = compute_density( mu, dnw, phi ) ;
-
- temp = compute_temperature( pres, theta, Cp, Rd, p0 ) ;
-
+ rho = compute_density( mu, dnw, phi ) ;
+ temp = compute_temperature( pres, theta, Cp, Rd, p0 ) ;
var_in = compute_reflectivity( qr, qg, qs, rho, temp ) ;
end
@@ -396,12 +407,12 @@
end
-% Plot field
+%% Plot field
subplot(m,m,pane);
axesm(map_proj{mp(id)},'Origin',[0 cen_lon(id) 0],'MapParallels', ...
- [stdlat1(id) stdlat2(id)],...
+ [truelat1(id) truelat2(id)],...
'MapLatLimit',[minlat maxlat],'MapLonLimit',[minlon maxlon]);
framem;
@@ -417,10 +428,10 @@
scale = ceil(we/20);
-% legend
+ % legend
- quiverm(xlat([1:scale:sn],[1:scale:we]),xlon([1:scale:sn],[1:scale:we]), ...
- fieldv([1:scale:sn],[1:scale:we]),fieldu([1:scale:sn],[1:scale:we]),'k')
+ quiverm(xlat(1:scale:sn,1:scale:we), xlon(1:scale:sn,1:scale:we), ...
+ fieldv(1:scale:sn,1:scale:we),fieldu(1:scale:sn,1:scale:we),'k')
else
@@ -438,23 +449,20 @@
% hm = clabelm(Cm,hm,'labelspacing',288); set(hm,'Fontsize',12);
% end
- if strcmp(field_name,'REF')
- eval('dbz_colors')
- h = pcolorm(xlat,xlon,field);
+ if strcmp(field_name,'REF')
+ eval('dbz_colors')
+ h = pcolorm(xlat,xlon,field);
+ else
+ [C h] = contourfm(xlat,xlon,field, iso); caxis([min(iso(:)),max(iso(:))]);
+ end
- else
+ %[C,h] = contour (field, iso);
+ %hold on
+ %[Cm,hm] = contour (field, -iso, '--');
+ cb = colorbar('vert'); set(cb,'Fontsize',12);
+ %clabel(C, h);
+ %clabel(Cm, hm);
- [C h] = contourfm(xlat,xlon,field, iso); caxis([min(iso(:)),max(iso(:))]);
-
- end
-
-%[C,h] = contour (field, iso);
-%hold on
-%[Cm,hm] = contour (field, -iso, '--');
-cb = colorbar('vert'); set(cb,'Fontsize',12);
-%clabel(C, h);
-%clabel(Cm, hm);
-
end
end
Modified: DART/trunk/models/wrf/matlab/map_wrf_spread_time.m
===================================================================
--- DART/trunk/models/wrf/matlab/map_wrf_spread_time.m 2010-03-25 22:01:44 UTC (rev 4326)
+++ DART/trunk/models/wrf/matlab/map_wrf_spread_time.m 2010-03-25 22:10:15 UTC (rev 4327)
@@ -10,163 +10,132 @@
% $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: ');
-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);
-ncopy = getnc(fname, 'copy');
-ens_size = size(ncopy, 1) - 2;
+prfname = 'Prior_Diag.nc';
+pofname = 'Posterior_Diag.nc';
-mean_ind = ens_size + 1;
-std_ind = mean_ind + 1;
+if (exist(prfname,'file') ~= 2)
+ error('%s does not exist.',prfname)
+end
+if (exist(pofname,'file') ~= 2)
+ error('%s does not exist.',pofname)
+end
-icopy = std_ind;
+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);
-true_times = getnc(fname, 'time');
-num_true_times = size(true_times, 1)
+ens_size = get_ens_size(prfname);
+sprd_ind = get_copy_index(prfname,'ensemble spread');
- stime = input('Initial time : ');
- ftime = input('End time : ');
+stime = input('Initial time (index): ');
+ftime = input('End time (index): ');
-% 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
-nx = we + 1;
-ny = sn;
+nx = we + 1;
+ny = sn;
var_units = 'U (m/s)';
-var_name = 'U';
-maxlev = bt;
-iso = [0.5:1:5];
+var_name = 'U';
+maxlev = bt;
+iso = 0.5:1:5;
+
if field_num > 1
-nx = we;
-ny = sn + 1;
-var_units = 'V (m/s)';
-var_name = 'V';
+ nx = we;
+ ny = sn + 1;
+ var_units = 'V (m/s)';
+ var_name = 'V';
end
if field_num > 2
-nx = we;
-ny = sn;
-var_units = 'W (m/s)';
-var_name = 'W';
-maxlev = bt + 1;
-iso = [0.01:0.01:0.1];
+ nx = we;
+ ny = sn;
+ var_units = 'W (m/s)';
+ var_name = 'W';
+ maxlev = bt + 1;
+ iso = 0.01:0.01:0.1;
end
if field_num > 3
-var_units = 'GZ (m^2/s^2)';
-var_name = 'PH';
-iso = [50:50:300];
+ var_units = 'GZ (m^2/s^2)';
+ var_name = 'PH';
+ iso = 50:50:300;
end
if field_num > 4
-var_units = 'T (K)';
-var_name = 'T';
-maxlev = bt;
-iso = [0.5:0.5:5];
+ var_units = 'T (K)';
+ var_name = 'T';
+ maxlev = bt;
+ iso = 0.5:0.5:5;
end
if field_num > 5
-var_units = 'MU (Pa)';
-var_name = 'MU';
-maxlev = 1;
-iso = [100:100:600];
+ var_units = 'MU (Pa)';
+ var_name = 'MU';
+ maxlev = 1;
+ iso = 100:100:600;
end
if field_num > 6
-var_units = 'QV (kg/kg)';
-var_name = 'QVAPOR';
-maxlev = bt;
-iso = [0.0001:0.0001:0.001];
+ var_units = 'QV (kg/kg)';
+ var_name = 'QVAPOR';
+ maxlev = bt;
+ iso = 0.0001:0.0001:0.001;
end
if field_num > 7
-var_units = 'QC (kg/kg)';
-var_name = 'QCLOUD';
-iso = [0.00001:0.00001:0.0001];
+ var_units = 'QC (kg/kg)';
+ var_name = 'QCLOUD';
+ iso = 0.00001:0.00001:0.0001;
end
if field_num > 8
-var_units = 'QR (kg/kg)';
-var_name = 'QRAIN';
-iso = [0.00001:0.00001:0.0001];
+ var_units = 'QR (kg/kg)';
+ var_name = 'QRAIN';
+ 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 = ftime-stime+1;
+m = ftime-stime+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)];
-if maxlev > 1
-corner_m = [itime icopy field_level -1 -1];
-end_point_m = [itime icopy field_level -1 -1];
-stride = [1 1 1 1 1];
-else
-corner_m = [itime icopy -1 -1];
-end_point_m = [itime icopy -1 -1];
-stride = [1 1 1 1];
-end
+ if maxlev > 1
+ start = [itime sprd_ind field_level 1 1] - 1;
+ count = [ 1 1 1 -1 -1];
+ else
+ start = [itime sprd_ind 1 1] -1;
+ count = [ 1 1 -1 -1];
+ end
-% Extract field
+ %% Extract fields
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list