[Dart-dev] [4899] DART/trunk/diagnostics: Extended obs_diag to finalize the netCDF file even when using
nancy at ucar.edu
nancy at ucar.edu
Fri May 6 16:35:46 MDT 2011
Revision: 4899
Author: thoar
Date: 2011-05-06 16:35:46 -0600 (Fri, 06 May 2011)
Log Message:
-----------
Extended obs_diag to finalize the netCDF file even when using
obs_seq.out files as input (the tables relating observation types
and integers is useful in decoding the locations file, if requested).
The output is restricted to only put the observation types USED as
netcdf global attributes.
The wind_obs_* stuff is completely superceded by the obs_seq_to_netcdf
and this version of plot_wind_vectors.
Modified Paths:
--------------
DART/trunk/diagnostics/matlab/plot_wind_vectors.m
DART/trunk/diagnostics/threed_sphere/obs_diag.f90
Removed Paths:
-------------
DART/trunk/diagnostics/threed_sphere/wind_obs_to_table.f90
DART/trunk/diagnostics/threed_sphere/wind_obs_to_table.html
DART/trunk/diagnostics/threed_sphere/wind_obs_to_table.nml
-------------- next part --------------
Modified: DART/trunk/diagnostics/matlab/plot_wind_vectors.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_wind_vectors.m 2011-05-05 16:23:45 UTC (rev 4898)
+++ DART/trunk/diagnostics/matlab/plot_wind_vectors.m 2011-05-06 22:35:46 UTC (rev 4899)
@@ -1,33 +1,43 @@
-function obs = plot_wind_vectors( fname, ncname, platform, varargin )
+function data = plot_wind_vectors( fname, platform, CopyString, QCString, varargin )
%% plot_wind_vectors creates maps with wind vector images overlain.
%
% The required arguments are:
-% fname ... the file name containing the wind observation pairs
-% ncname ... the output netCDF file from obs_diag.f90 (for some metadata)
+% fname ... the output netCDF file from obs_seq_to_netcdf
% platform ... a string to represent the observation platform.
+% usually 'RADIOSONDE', 'SAT', 'METAR', ... get hints from:
+% ncdump -v ObsTypesMetaData obs_epoch_xxx.nc | grep _U_
+% CopyString ... which observation copy is of interest?
+% ncdump -v CopyMetaData obs_epoch_xxx.nc
+% QCString ... which QC copy is of interest?
+% ncdump -v QCMetaData obs_epoch_xxx.nc
%
% The optional arguments are:
-% levels ... specifies the vertical area of interest. If not present,
-% all vertical levels are used.
-% region ... specifies that horizontal area of interest. If not present,
+% region ... specifies that horizontal & vertical area of interest. If not present,
% all available observations are used.
% scalefactor ... provides control over the plotted size of the
% wind vectors. A smaller number results in a
% bigger wind vector. If not present, a value of 10.0 is
% used.
%
-% fname = 'wind_vectors.006.dat';
-% ncname = 'obs_to_table_output.nc';
-% platform = 'RADIOSONDE';
-% levels = [1020 500];
-% region = [0 360 0 90]; %
+% EXAMPLE 1:
+% fname = 'obs_epoch_001.nc';
+% platform = 'SAT'; % usually 'RADIOSONDE', 'SAT', 'METAR', ...
+% CopyString = 'NCEP BUFR observation';
+% QCString = 'DART quality control';
+% region = [0 360 0 90 1020 500]; %
% scalefactor = 5; % reference arrow magnitude
%
-% obs = plot_wind_vectors(fname, ncname, platform, ...
-% 'levels', levels, 'region', region, 'scalefactor', scalefactor);
+% obs = plot_wind_vectors(fname, platform, CopyString, QCString, ...
+% 'region', region, 'scalefactor', scalefactor);
%
+%
+% EXAMPLE 2 (CONUS domain):
+%
+% region = [210 310 12 65 -Inf Inf];
+% obs = plot_wind_vectors('obs_epoch_001.nc', 'SAT', ...
+% 'NCEP BUFR observation', 'DART quality control','region',region);
-%% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+%% DART software - Copyright � 2004 - 2010 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
% http://www.image.ucar.edu/DAReS/DART/DART_download
%
@@ -37,37 +47,56 @@
% $Revision$
% $Date$
-narg = nargin;
+% Set sensible defaults
-if narg == 3
+region = [0 360 -90 90 -Inf Inf];
+scalefactor = 10;
- levels = [];
- region = [];
- scalefactor = 10;
+% Harvest input
-else
+if nargin ~= 4
+ [region, scalefactor] = parseinput(varargin{:});
+end
- [levels, region, scalefactor] = parseinput(varargin{:});
+%% Start the ball rolling
+if (exist(fname,'file') ~= 2)
+ error('%s does not exist',fname)
end
data.filename = fname;
-data.ncname = ncname;
data.platform = platform;
-data.levels = levels;
+data.copystring = CopyString;
+data.qcstring = QCString;
data.region = region;
data.scalefactor = scalefactor;
-f = netcdf(ncname,'nowrite');
-data.platforms = f{'ObservationTypes'}(:);
-data.timeunits = f{'time_bounds'}.units(:);
-close(f);
+%% Read the observation sequence
-timebase = sscanf(data.timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
-data.timeorigin = datenum(timebase(1),timebase(2),timebase(3));
-obs = Load_Parse(data);
+[UtypeString, VtypeString] = FindObsType( fname, platform );
+verbose = 0;
+Uobs = read_obs_netcdf(fname, UtypeString, region, CopyString, QCString, verbose);
+Vobs = read_obs_netcdf(fname, VtypeString, region, CopyString, QCString, verbose);
-if isempty(obs)
+if (length(Uobs.obs) ~= length(Vobs.obs))
+ error('Houston, we have a problem.')
+end
+
+lonmismatch = Uobs.lons ~= Vobs.lons;
+latmismatch = Uobs.lats ~= Vobs.lats;
+ zmismatch = Uobs.z ~= Vobs.z;
+
+if ( sum(lonmismatch) ~= 0),
+ warning('DART:UVcollocaton','There are %d mismatched (in longitude) observations',sum(lonmismatch))
+end
+if ( sum(latmismatch) ~= 0),
+ warning('DART:UVcollocaton','There are %d mismatched (in latitude) observations',sum(latmismatch))
+end
+if ( sum(zmismatch) ~= 0),
+ warning('DART:UVcollocaton','There are %d mismatched (in vertical) observations',sum(zmismatch))
+end
+
+if (sum(lonmismatch) == length(lonmismatch))
clf;
axis([0 1 0 1]); axis image
h = text(0.5,0.5,sprintf('%s has no %s data',data.filename,data.platform));
@@ -75,58 +104,86 @@
return
end
-obs.times = obs.times + data.timeorigin;
-t1 = datestr(min(obs.times));
-t2 = datestr(max(obs.times));
+%% must only use the observations that are co-located to
+% set up the plotting structure.
-clf;
-axlims = DrawBackground( obs );
+inds = ((Uobs.lons == Vobs.lons) & ...
+ (Uobs.lats == Vobs.lats) & ...
+ (Uobs.z == Vobs.z));
+data.time = Uobs.time(inds);
+data.lon = Uobs.lons(inds);
+data.lat = Uobs.lats(inds);
+data.z = Uobs.z(inds);
+data.Uqc = Uobs.qc(inds);
+data.Vqc = Vobs.qc(inds);
+data.U = Uobs.obs(inds);
+data.V = Vobs.obs(inds);
+data.level1 = min(Uobs.z);
+data.levelN = max(Uobs.z);
+data.levstring = sprintf('%.2f to %.2f',data.level1,data.levelN);
-goodUV = find( (obs.Uqc < 1) & (obs.Vqc < 1));
-baadUV = find( (obs.Uqc > 1) & (obs.Vqc > 1));
-goodU = find( (obs.Uqc < 1) & (obs.Vqc > 0));
-goodV = find( (obs.Vqc < 1) & (obs.Uqc > 0));
+%% Start the Plotting
+clf;
+axish = gca;
+axlims = DrawBackground( data );
+set(axish,'Layer','top')
+
+if ( isfinite(strfind(lower(QCString),'dart')) )
+ % We know how to interpret QC codes
+ goodUV = find( (data.Uqc < 2) & (data.Vqc < 2));
+ baadUV = find( (data.Uqc > 1) & (data.Vqc > 1));
+ goodU = find( (data.Uqc < 2) & (data.Vqc > 1));
+ goodV = find( (data.Uqc > 1) & (data.Vqc < 2));
+else
+ % We do not know how to interpret QC codes, so they
+ % are all 'good'
+ baadUV = [];
+ goodU = [];
+ goodV = [];
+end
+
legh = [];
legstr = {};
if ~ isempty(goodUV)
- hgood = obs2plot(obs, goodUV, [0 0 0] );
+ hgood = obs2plot(data, goodUV, [0 0 0] );
legh = [legh; hgood];
- legstr{length(legstr)+1} = sprintf('%d good',length(goodUV));
+ legstr{length(legstr)+1} = sprintf('%d ''good''',length(goodUV));
end
if ~ isempty(baadUV)
- hbaadUV = obs2plot(obs, baadUV, [1 0 0] );
+ hbaadUV = obs2plot(data, baadUV, [1 0 0] );
legh = [legh; hbaadUV];
- legstr{length(legstr)+1} = sprintf('%d badUbadV',length(baadUV));
+ legstr{length(legstr)+1} = sprintf('%d ''badU badV''',length(baadUV));
end
if ~ isempty(goodU)
- hgoodU = obs2plot(obs, goodU, [0 1 0] );
+ hgoodU = obs2plot(data, goodU, [0 1 0] );
legh = [legh; hgoodU];
- legstr{length(legstr)+1} = sprintf('%d goodUbadV',length(goodU));
+ legstr{length(legstr)+1} = sprintf('%d ''goodU badV''',length(goodU));
end
if ~ isempty(goodV)
- hgoodV = obs2plot(obs, goodV, [0 0 1] );
+ hgoodV = obs2plot(data, goodV, [0 0 1] );
legh = [legh; hgoodV];
- legstr{length(legstr)+1} = sprintf('%d badUgoodV',length(goodV));
+ legstr{length(legstr)+1} = sprintf('%d ''badU goodV''',length(goodV));
end
+t1 = datestr(min(data.time),'yyyy-mm-dd HH:MM:SS');
+t2 = datestr(max(data.time),'yyyy-mm-dd HH:MM:SS');
h = title({sprintf('%s %s %s',t1,platform,t2), ...
- sprintf('levels %s ',obs.levelstring)});
+ sprintf('levels %s ',data.levstring)});
set(h,'FontSize',18)
h = xlabel(data.filename); set(h,'Interpreter','none');
-legend(legh,legstr,'Location','NorthWestOutside')
+legend(legh,legstr,'Location','Best')
hold off;
-
function axlims = DrawBackground( obs )
%======================================================================
@@ -134,14 +191,14 @@
% Figure out bounds of the data
axlims = [ min(obs.lon) max(obs.lon) min(obs.lat) max(obs.lat) ] ;
else
- axlims = obs.region;
+ axlims = obs.region(1:4);
end
% It is nice to have a little padding around the perimeter
dx = 0.05 * (axlims(2) - axlims(1));
dy = 0.05 * (axlims(4) - axlims(3));
-
-axis(axlims + [-dx dx -dy dy])
+axlims(1:4) = axlims(1:4) + [-dx dx -dy dy];
+axis(axlims)
axis image
% It is nice to know where the land is
@@ -181,237 +238,75 @@
-function obs = Load_Parse( data )
+function [ustring, vstring] = FindObsType( ncname, platform )
%======================================================================
% Makes no attempt to find/replace/identify MISSING values
%
% data.filename incoming filename with the data
% data.ncname incoming filename with the metadata
% data.platform the observation platform of interest
-% data.levels the top/bottom levels of interest
% data.region the region of interest [xmin xmax ymin ymax]
% data.scalefactor the reference wind vector magnitude
% data.platforms the observation platforms in the incoming file
% data.timeunits the units string e.g. 'days since yyyy-mm-dd'
% data.timeorigin
-obsmat = load(data.filename);
+ObsTypeStrings = nc_varget(ncname,'ObsTypesMetaData');
% Find the types of data in the obs_sequence file for this epoch.
+% Turns out all the winds are either xxxx_U_WIND_COMPONENT or xxxx_U_10_METER_WIND
+% so either way, they start out xxxx_U_
-platformIDs = unique(obsmat(:,1));
-uid = floor(platformIDs/100);
-vid = platformIDs - uid*100;
-Ustrings = data.platforms(uid,:);
-Vstrings = data.platforms(vid,:);
+utarget = sprintf('%s_U_',strtrim(platform));
+vtarget = sprintf('%s_V_',strtrim(platform));
+n = length(utarget);
+ustring = [];
+vstring = [];
-nplatforms = length(uid);
-pid = [];
-obs = [];
-levelstring = [];
-regionstring = [];
+for i = 1:size(ObsTypeStrings,1)
-% This block divines the platform string and companion obs_kinds
-% from the netCDF file metadata - the only reason we need the netCDF file.
+ utf = strncmpi(ObsTypeStrings(i,:),utarget,n);
+ vtf = strncmpi(ObsTypeStrings(i,:),vtarget,n);
-for i = 1:nplatforms
- uindex = findstr(Ustrings(i,:),'_U_WIND_COMPONENT') - 1;
- vindex = findstr(Vstrings(i,:),'_V_WIND_COMPONENT') - 1;
-
- if (isempty(uindex) | isempty(vindex))
- uindex = findstr(Ustrings(i,:),'_U_10_METER_WIND') - 1;
- vindex = findstr(Vstrings(i,:),'_V_10_METER_WIND') - 1;
+ if ( utf )
+ ustring = ObsTypeStrings(i,:);
+ uindex = i;
end
+ if ( vtf )
+ vstring = ObsTypeStrings(i,:);
+ vindex = i;
+ end
- Ubase = Ustrings(i,1:uindex);
- Vbase = Vstrings(i,1:vindex);
-
- if ~strcmp(Ubase,Vbase)
- error('U and V wind component platforms do not match')
- end
-
- % Determine what numeric pid corresponds to the desired platform string
- if strcmp(lower(Ubase),lower(data.platform))
- pid = platformIDs(i);
- end
-
- % echo a little informational statement about how many obs of
- % this type there are in this assimilation period.
- inds = find(obsmat(:,1) == platformIDs(i));
- nobs = length(inds);
-
- disp(sprintf('%6d %14s observations in %s (%4d)',nobs,Ubase,data.filename,platformIDs(i)))
-
end
-% This block extracts just the desired observations based on platform.
-
-if isempty(pid)
- disp(sprintf('no %s observations in %s', data.platform, data.filename))
- return
+if ( isempty(ustring) || isempty(vstring) )
+ error('no %s winds in %s',platform,ncname)
end
-inds = find(obsmat(:,1) == pid);
+% echo a little informational statement about the number of obs
-if isempty(inds)
- disp(sprintf('no %s observations (type %d) in %s',data.platform, pid, data.filename))
- return
-end
+obs_type = nc_varget(ncname,'obs_type');
+numU = sum(obs_type == uindex);
+numV = sum(obs_type == vindex);
+fprintf('%8d %s observations in %s\n', numU, ustring, ncname)
+fprintf('%8d %s observations in %s\n', numV, vstring, ncname)
-platform = obsmat(inds, 1);
-day = obsmat(inds, 2);
-seconds = obsmat(inds, 3);
-lon = obsmat(inds, 4);
-lat = obsmat(inds, 5);
-level = obsmat(inds, 6);
-Uqc = obsmat(inds, 7);
-Vqc = obsmat(inds, 8);
-U = obsmat(inds, 9);
-V = obsmat(inds,10);
-if ( size(obsmat,2) > 10 )
- Upr = obsmat(inds,11);
- Vpr = obsmat(inds,12);
- Upo = obsmat(inds,13);
- Vpo = obsmat(inds,14);
+if (numU ~= numV)
+ error('Different number of U,V observations. Dying ...')
end
-times = day + seconds/86400;
-%--------------------------------------------------
-% Subset the levels of interest
-%--------------------------------------------------
-if ( isempty(data.levels) )
- % then we want all levels, do nothing ...
- level1 = min(level);
- levelN = max(level);
- levelstring = sprintf('all (%.2f to %.2f)',level1,levelN);
-else
- level1 = min(data.levels);
- levelN = max(data.levels);
-
- levelstring = sprintf('%.2f to %.2f',level1,levelN);
-
- inds = find ((level >= level1) & (level <= levelN));
-
- if (length(inds) == 0)
- disp(sprintf('no %s observations in %s', data.platform, levelstring))
- return
- end
-
- platform = platform(inds);
- day = day(inds);
- seconds = seconds(inds);
- lon = lon(inds);
- lat = lat(inds);
- level = level(inds);
- Uqc = Uqc(inds);
- Vqc = Vqc(inds);
- U = U(inds);
- V = V(inds);
- if ( size(obsmat,2) > 10 )
- Upr = Upr(inds);
- Vpr = Vpr(inds);
- Upo = Upo(inds);
- Vpo = Vpo(inds);
- end
- times = times(inds);
-end
-
-%--------------------------------------------------
-% Subset the region of interest
-% for the moment, we are not supporting wrapping at Greenwich.
-%--------------------------------------------------
-
-if isempty(data.region)
- % then we want the entire dataset, do nothing ...
- regionstring = 'global';
-else
- lon1 = data.region(1);
- lonN = data.region(2);
- lat1 = data.region(3);
- latN = data.region(4);
-
- regionstring = sprintf('(%.2f -> %.2f, %.2f -> %.2f)',lon1,lonN,lat1,latN);
-
- inds = find ((lon >= lon1) & (lon <= lonN) & ...
- (lat >= lat1) & (lat <= latN));
-
- if (length(inds) == 0)
- disp(sprintf('no %s observations in %s', data.platform, regionstring))
- return
- end
-
- platform = platform(inds);
- day = day(inds);
- seconds = seconds(inds);
- lon = lon(inds);
- lat = lat(inds);
- level = level(inds);
- Uqc = Uqc(inds);
- Vqc = Vqc(inds);
- U = U(inds);
- V = V(inds);
- if ( size(obsmat,2) > 10 )
- Upr = Upr(inds);
- Vpr = Vpr(inds);
- Upo = Upo(inds);
- Vpo = Vpo(inds);
- end
- times = times(inds);
-end
-
-%--------------------------------------------------
-% Insert NaN's for missing values
-%--------------------------------------------------
-
-U( U < -900) = NaN;
-V( V < -900) = NaN;
-if ( size(obsmat,2) > 10 )
- Upr( Upr < -900) = NaN;
- Vpr( Vpr < -900) = NaN;
- Upo( Upo < -900) = NaN;
- Vpo( Vpo < -900) = NaN;
-end
-
-%--------------------------------------------------
-% Create the output structure.
-%--------------------------------------------------
-
-obs.platform = platform;
-obs.day = day;
-obs.seconds = seconds;
-obs.lon = lon;
-obs.lat = lat;
-obs.level = level;
-obs.Uqc = Uqc;
-obs.Vqc = Vqc;
-obs.U = U;
-obs.V = V;
-if ( size(obsmat,2) > 10 )
- obs.Upr = Upr;
- obs.Vpr = Vpr;
- obs.Upo = Upo;
- obs.Vpo = Vpo;
-end
-obs.times = times;
-obs.levels = data.levels;
-obs.region = data.region;
-obs.scalefactor = data.scalefactor;
-obs.levelstring = levelstring;
-
-
-function [levels, region, scalefactor] = parseinput(varargin)
+function [region, scalefactor] = parseinput(varargin)
%======================================================================
% try to parse the input pairs ... which must be pairs
if (mod(length(varargin),2) ~= 0)
- error('Wrong number (%d) of optional arguments. Must be parameter/value pairs: ''levels'',[1000 500]',length(varargin))
+ error('Wrong number (%d) of optional arguments. Must be parameter/value pairs: ''region'',[0 360 -90 90 1020 500]',length(varargin))
end
npairs = length(varargin)/2;
-levels = [];
region = [];
scalefactor = 10.0;
@@ -428,15 +323,12 @@
end
end
-% Make sure the levels array has a top/bottom
-if ~ isempty(levels)
+if (length(region) < 4)
+ region = [0 360 -90 90 -Inf Inf];
+elseif (length(region) == 4)
+ region = [region -Inf Inf];
+elseif (length(region) ~= 6)
+ warning('DART:region input','region must be length 4 or 6')
+ error('Unable to interpret region - %s',num2str(region))
end
-% Make sure the geographic array makes sense
-if ~ isempty(region)
-end
-
-% Make sure the scalefactor makes sense
-if ~ isempty(scalefactor)
-end
-
Modified: DART/trunk/diagnostics/threed_sphere/obs_diag.f90
===================================================================
--- DART/trunk/diagnostics/threed_sphere/obs_diag.f90 2011-05-05 16:23:45 UTC (rev 4898)
+++ DART/trunk/diagnostics/threed_sphere/obs_diag.f90 2011-05-06 22:35:46 UTC (rev 4899)
@@ -113,7 +113,7 @@
integer :: prior_spread_index, posterior_spread_index
integer :: flavor, wflavor ! THIS IS THE (global) 'KIND' in the obs_def_mod list.
integer :: num_copies, num_qc, num_obs, max_num_obs, obs_seq_file_id
-integer :: num_obs_kinds
+integer :: num_obs_types
! variables used primarily/exclusively for the rank histogram
integer :: ens_size, rank_histogram_bin
@@ -233,7 +233,7 @@
real(r8), dimension(:,:,:,:), pointer :: observation, ens_mean
integer, dimension(:,:,:,:), pointer :: NDartQC_0, NDartQC_1, NDartQC_2, NDartQC_3
integer, dimension(:,:,:,:), pointer :: NDartQC_4, NDartQC_5, NDartQC_6, NDartQC_7
- integer, dimension(:,:,:,:,:), pointer :: hist_bin
+ integer, dimension(:,:,:,:,:), pointer :: hist_bin => NULL()
end type TLRV_type
type LRV_type
@@ -289,7 +289,9 @@
real(r8), allocatable, dimension(:) :: scale_factor ! to convert to plotting units
integer, allocatable, dimension(:) :: ob_defining_vert ! obs index defining vert coord type
-character(len = stringlength), pointer, dimension(:) :: my_obs_kind_names
+! List of observations types augmented with 'WIND' types if need be.
+! Replace calls to 'get_obs_kind_name' ---> index into 'obs_type_strings'
+character(len = stringlength), pointer, dimension(:) :: obs_type_strings
! These pairs of variables are used when we diagnose which observations
! are far from the background.
@@ -323,14 +325,15 @@
! Define/Append the 'horizontal wind' obs_kinds to supplant the list declared
! in obs_kind_mod.f90 i.e. if there is a RADIOSONDE_U_WIND_COMPONENT
! and a RADIOSONDE_V_WIND_COMPONENT, there must be a RADIOSONDE_HORIZONTAL_WIND
-! Replace calls to 'get_obs_kind_name' with variable 'my_obs_kind_names'
+! Replace calls to 'get_obs_kind_name' with variable 'obs_type_strings'
!----------------------------------------------------------------------
-num_obs_kinds = grok_observation_names(my_obs_kind_names)
+num_obs_types = grok_observation_names(obs_type_strings)
-allocate( which_vert(num_obs_kinds), &
- scale_factor(num_obs_kinds), &
- ob_defining_vert(num_obs_kinds))
+allocate( which_vert(num_obs_types), &
+ scale_factor(num_obs_types), &
+ ob_defining_vert(num_obs_types) )
+
which_vert = VERTISUNDEF
scale_factor = 1.0_r8
ob_defining_vert = -1
@@ -408,7 +411,7 @@
! SetTime rectifies user input and the final binning sequence.
!----------------------------------------------------------------------
-call SetTime(beg_time, end_time, binsep, binwidth, halfbinwidth, &
+call SetTime(beg_time, end_time, binsep, halfbinwidth, &
TimeMin, TimeMax, Nepochs, bincenter, binedges, epoch_center, epoch_edges, &
obs_used_in_epoch)
@@ -430,27 +433,27 @@
allocate(obs_seq_filenames(Nepochs*400))
obs_seq_filenames = 'null'
-allocate(guess%rmse( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%bias( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%spread( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%totspread( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%observation(Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%ens_mean( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%Nposs( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%Nused( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NbigQC( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NbadIZ( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NbadUV( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NbadLV( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NbadDartQC( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NDartQC_0( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NDartQC_1( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NDartQC_2( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NDartQC_3( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NDartQC_4( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NDartQC_5( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NDartQC_6( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- guess%NDartQC_7( Nepochs, Nlevels, Nregions, num_obs_kinds) )
+allocate(guess%rmse( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%bias( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%spread( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%totspread( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%observation(Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%ens_mean( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%Nposs( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%Nused( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NbigQC( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NbadIZ( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NbadUV( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NbadLV( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NbadDartQC( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NDartQC_0( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NDartQC_1( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NDartQC_2( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NDartQC_3( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NDartQC_4( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NDartQC_5( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NDartQC_6( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%NDartQC_7( Nepochs, Nlevels, Nregions, num_obs_types) )
guess%rmse = 0.0_r8
guess%bias = 0.0_r8
@@ -478,29 +481,29 @@
guess%num_times = Nepochs
guess%num_levels = Nlevels
guess%num_regions = Nregions
-guess%num_variables = num_obs_kinds
+guess%num_variables = num_obs_types
-allocate(analy%rmse( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%bias( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%spread( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%totspread( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%observation(Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%ens_mean( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%Nposs( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%Nused( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NbigQC( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NbadIZ( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NbadUV( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NbadLV( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NbadDartQC( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NDartQC_0( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NDartQC_1( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NDartQC_2( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NDartQC_3( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NDartQC_4( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NDartQC_5( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NDartQC_6( Nepochs, Nlevels, Nregions, num_obs_kinds), &
- analy%NDartQC_7( Nepochs, Nlevels, Nregions, num_obs_kinds) )
+allocate(analy%rmse( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%bias( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%spread( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%totspread( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%observation(Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%ens_mean( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%Nposs( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%Nused( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NbigQC( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NbadIZ( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NbadUV( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NbadLV( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NbadDartQC( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NDartQC_0( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NDartQC_1( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NDartQC_2( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NDartQC_3( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NDartQC_4( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NDartQC_5( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NDartQC_6( Nepochs, Nlevels, Nregions, num_obs_types), &
+ analy%NDartQC_7( Nepochs, Nlevels, Nregions, num_obs_types) )
analy%rmse = 0.0_r8
analy%bias = 0.0_r8
@@ -528,29 +531,29 @@
analy%num_times = Nepochs
analy%num_levels = Nlevels
analy%num_regions = Nregions
-analy%num_variables = num_obs_kinds
+analy%num_variables = num_obs_types
-allocate(guessAVG%rmse( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%bias( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%spread( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%totspread( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%observation(Nlevels, Nregions, num_obs_kinds), &
- guessAVG%ens_mean( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%Nposs( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%Nused( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NbigQC( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NbadIZ( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NbadUV( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NbadLV( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NbadDartQC( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NDartQC_0( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NDartQC_1( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NDartQC_2( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NDartQC_3( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NDartQC_4( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NDartQC_5( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NDartQC_6( Nlevels, Nregions, num_obs_kinds), &
- guessAVG%NDartQC_7( Nlevels, Nregions, num_obs_kinds) )
+allocate(guessAVG%rmse( Nlevels, Nregions, num_obs_types), &
+ guessAVG%bias( Nlevels, Nregions, num_obs_types), &
+ guessAVG%spread( Nlevels, Nregions, num_obs_types), &
+ guessAVG%totspread( Nlevels, Nregions, num_obs_types), &
+ guessAVG%observation(Nlevels, Nregions, num_obs_types), &
+ guessAVG%ens_mean( Nlevels, Nregions, num_obs_types), &
+ guessAVG%Nposs( Nlevels, Nregions, num_obs_types), &
+ guessAVG%Nused( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NbigQC( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NbadIZ( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NbadUV( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NbadLV( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NbadDartQC( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NDartQC_0( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NDartQC_1( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NDartQC_2( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NDartQC_3( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NDartQC_4( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NDartQC_5( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NDartQC_6( Nlevels, Nregions, num_obs_types), &
+ guessAVG%NDartQC_7( Nlevels, Nregions, num_obs_types) )
guessAVG%rmse = 0.0_r8
guessAVG%bias = 0.0_r8
@@ -577,29 +580,29 @@
guessAVG%string = 'VPguess'
guessAVG%num_levels = Nlevels
guessAVG%num_regions = Nregions
-guessAVG%num_variables = num_obs_kinds
+guessAVG%num_variables = num_obs_types
-allocate(analyAVG%rmse( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%bias( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%spread( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%totspread( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%observation(Nlevels, Nregions, num_obs_kinds), &
- analyAVG%ens_mean( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%Nposs( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%Nused( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NbigQC( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NbadIZ( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NbadUV( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NbadLV( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NbadDartQC( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NDartQC_0( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NDartQC_1( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NDartQC_2( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NDartQC_3( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NDartQC_4( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NDartQC_5( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NDartQC_6( Nlevels, Nregions, num_obs_kinds), &
- analyAVG%NDartQC_7( Nlevels, Nregions, num_obs_kinds) )
+allocate(analyAVG%rmse( Nlevels, Nregions, num_obs_types), &
+ analyAVG%bias( Nlevels, Nregions, num_obs_types), &
+ analyAVG%spread( Nlevels, Nregions, num_obs_types), &
+ analyAVG%totspread( Nlevels, Nregions, num_obs_types), &
+ analyAVG%observation(Nlevels, Nregions, num_obs_types), &
+ analyAVG%ens_mean( Nlevels, Nregions, num_obs_types), &
+ analyAVG%Nposs( Nlevels, Nregions, num_obs_types), &
+ analyAVG%Nused( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NbigQC( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NbadIZ( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NbadUV( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NbadLV( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NbadDartQC( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NDartQC_0( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NDartQC_1( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NDartQC_2( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NDartQC_3( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NDartQC_4( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NDartQC_5( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NDartQC_6( Nlevels, Nregions, num_obs_types), &
+ analyAVG%NDartQC_7( Nlevels, Nregions, num_obs_types) )
analyAVG%rmse = 0.0_r8
analyAVG%bias = 0.0_r8
@@ -626,7 +629,7 @@
analyAVG%string = 'VPanaly'
analyAVG%num_levels = Nlevels
analyAVG%num_regions = Nregions
-analyAVG%num_variables = num_obs_kinds
+analyAVG%num_variables = num_obs_types
!----------------------------------------------------------------------
! Open file for histogram of innovations, as a function of standard deviation.
@@ -677,8 +680,7 @@
! Read in information about observation sequence so we can allocate
! observations. We need info about how many copies, qc values, etc.
- obs_seq_in_file_name = trim(adjustl(obs_seq_in_file_name)) ! Lahey requirement
- obs_seq_filenames(ifile) = trim(adjustl(obs_seq_in_file_name))
+ obs_seq_filenames(ifile) = trim( obs_seq_in_file_name)
call read_obs_seq_header(obs_seq_in_file_name, &
num_copies, num_qc, num_obs, max_num_obs, &
@@ -701,14 +703,14 @@
write(logfileunit,*)'num_qc is ',num_qc
write(logfileunit,*)'num_obs is ',num_obs
write(logfileunit,*)'max_num_obs is ',max_num_obs
- write(logfileunit,*)'obs_seq_read_format is ',trim(adjustl(obs_seq_read_format))
+ write(logfileunit,*)'obs_seq_read_format is ',trim(obs_seq_read_format)
write(logfileunit,*)'pre_I_format is ',pre_I_format
write( * ,*)
write( * ,*)'num_copies is ',num_copies
write( * ,*)'num_qc is ',num_qc
write( * ,*)'num_obs is ',num_obs
write( * ,*)'max_num_obs is ',max_num_obs
- write( * ,*)'obs_seq_read_format is ',trim(adjustl(obs_seq_read_format))
+ write( * ,*)'obs_seq_read_format is ',trim(obs_seq_read_format)
write( * ,*)'pre_I_format is ',pre_I_format
endif
@@ -777,9 +779,9 @@
else
if (verbose) then
write(logfileunit,*)'seqTN > TimeMin ... using ', &
- trim(adjustl(obs_seq_in_file_name))
+ trim(obs_seq_in_file_name)
write( * ,*)'seqTN > TimeMin ... using ', &
- trim(adjustl(obs_seq_in_file_name))
+ trim(obs_seq_in_file_name)
endif
endif
@@ -798,7 +800,7 @@
if (allocated(copyvals)) deallocate( copyvals )
exit ObsFileLoop
else
- if (verbose) write(*,*)'seqT1 < TimeMax ... using ',trim(adjustl(obs_seq_in_file_name))
+ if (verbose) write(*,*)'seqT1 < TimeMax ... using ',trim(obs_seq_in_file_name)
endif
!--------------------------------------------------------------------
@@ -823,12 +825,13 @@
endif
else
! This should happen exactly once, if at all.
- allocate(guess%hist_bin( Nepochs, Nlevels, Nregions, num_obs_kinds, ens_size+1))
+ allocate(guess%hist_bin( Nepochs, Nlevels, Nregions, num_obs_types, ens_size+1))
allocate(ens_copy_index(ens_size))
guess%hist_bin = 0
call init_random_seq(ran_seq, seed=23)
endif
create_rank_histogram = .true.
+ if (verbose) write(*,*) 'Creating rank histogram with ',ens_size+1,' bins.'
else
write(*,*) 'Cannot create rank histogram.'
create_rank_histogram = .false.
@@ -876,9 +879,9 @@
! Append epoch number to name
write(locName,'(a,i3.3,a)') 'observation_locations.', iepoch, '.dat'
if (file_exist(locName)) then
- lunit = open_file(trim(adjustl(locName)),form='formatted',action='append')
+ lunit = open_file(locName, form='formatted', action='append')
else
- lunit = open_file(trim(adjustl(locName)),form='formatted',action='rewind')
+ lunit = open_file(locName, form='formatted', action='rewind')
write(lunit, '(a)') ' lon lat lev kind key QCval'
endif
endif
@@ -1120,7 +1123,7 @@
if ( create_rank_histogram ) then
call get_obs_values(observation, copyvals)
rank_histogram_bin = Rank_Histogram(copyvals, obs_index, &
- obs_error_variance, ens_size, ens_copy_index)
+ obs_error_variance, ens_copy_index)
endif
!--------------------------------------------------------------
@@ -1380,8 +1383,8 @@
enddo EpochLoop
if (verbose) then
- write(logfileunit,*)'End of EpochLoop for ',trim(adjustl(obs_seq_in_file_name))
- write( * ,*)'End of EpochLoop for ',trim(adjustl(obs_seq_in_file_name))
+ write(logfileunit,*)'End of EpochLoop for ',trim(obs_seq_in_file_name)
+ write( * ,*)'End of EpochLoop for ',trim(obs_seq_in_file_name)
endif
call destroy_obs(obs1)
@@ -1401,8 +1404,8 @@
if (verbose) then
do ivar = 1,SIZE(which_vert)
- write(logfileunit,*)'which_vert(',ivar,' of ',num_obs_kinds,') = ',which_vert(ivar)
- write( * ,*)'which_vert(',ivar,' of ',num_obs_kinds,') = ',which_vert(ivar)
+ write(logfileunit,*)'which_vert(',ivar,' of ',num_obs_types,') = ',which_vert(ivar)
+ write( * ,*)'which_vert(',ivar,' of ',num_obs_types,') = ',which_vert(ivar)
enddo
endif
@@ -1411,7 +1414,7 @@
write( * ,*)'Normalizing time-level-region-variable quantities.'
endif
-do ivar = 1,num_obs_kinds
+do ivar = 1,num_obs_types
do iregion= 1,Nregions
do ilev = 1,Nlevels
do iepoch = 1,Nepochs
@@ -1509,7 +1512,7 @@
write( * ,*)'Normalize quantities for all levels.'
endif
-do ivar=1,num_obs_kinds
+do ivar=1,num_obs_types
do iregion=1, Nregions
do ilev=1, Nlevels
@@ -1669,14 +1672,19 @@
if ( sum(obs_used_in_epoch) == 0 ) then
- call print_date(AllseqT1,'First observation date')
- call print_date( TimeMin,'First requested date')
- call print_date(AllseqTN,'Last observation date')
- call print_date( TimeMax,'Last requested date')
+ call print_date(AllseqT1,' First observation date')
+ call print_date(AllseqTN,' Last observation date')
+ call print_date( TimeMin,' First requested date')
+ call print_date( TimeMax,' Last requested date')
- write(msgstring,*)'No observations in requested time bins.'
- call error_handler(E_ERR,'obs_diag',msgstring,source,revision,revdate)
+ write( * ,*)
+ write(logfileunit,*)
+ write(msgstring,*)'WARNING: NO OBSERVATIONS in requested time bins.'
+ call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+
endif
!----------------------------------------------------------------------
@@ -1699,9 +1707,10 @@
guess%NDartQC_0, guess%NDartQC_1, guess%NDartQC_2, guess%NDartQC_3, &
guess%NDartQC_4, guess%NDartQC_5, guess%NDartQC_6, guess%NDartQC_7)
-if ( create_rank_histogram ) &
-deallocate(guess%hist_bin, ens_copy_index)
+if (associated(guess%hist_bin)) deallocate(guess%hist_bin)
+if (allocated(ens_copy_index)) deallocate(ens_copy_index)
+
deallocate(analy%rmse, analy%bias, analy%spread, analy%totspread, &
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list