[Dart-dev] [3825] DART/trunk: This change supports the creation of ASCII files ( and supporting plotting
nancy at ucar.edu
nancy at ucar.edu
Fri Apr 17 13:13:26 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090417/fd607869/attachment-0001.html
-------------- next part --------------
Copied: DART/trunk/diagnostics/matlab/PlotWindObs.m (from rev 3819, DART/trunk/models/wrf/matlab/PlotWindObs.m)
===================================================================
--- DART/trunk/diagnostics/matlab/PlotWindObs.m (rev 0)
+++ DART/trunk/diagnostics/matlab/PlotWindObs.m 2009-04-17 19:13:25 UTC (rev 3825)
@@ -0,0 +1,27 @@
+%
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2007, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+for periods = 2:37;
+
+ fname = sprintf('wind_vectors.%03d.dat',periods);
+ ncname = 'obs_diag_output.nc';
+ platform = 'SAT';
+ level = -1;
+
+ obs = plot_wind_vectors(fname,ncname,platform,level);
+
+ disp('Pausing, hit any key to continue ...')
+ pause
+
+end
+
Copied: DART/trunk/diagnostics/matlab/plot_wind_vectors.m (from rev 3823, DART/trunk/models/wrf/matlab/plot_wind_vectors.m)
===================================================================
--- DART/trunk/diagnostics/matlab/plot_wind_vectors.m (rev 0)
+++ DART/trunk/diagnostics/matlab/plot_wind_vectors.m 2009-04-17 19:13:25 UTC (rev 3825)
@@ -0,0 +1,441 @@
+function obs = plot_wind_vectors( fname, ncname, platform, 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)
+% platform ... a string to represent the observation platform.
+%
+% 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,
+% 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]; %
+% scalefactor = 5; % reference arrow magnitude
+%
+% obs = plot_wind_vectors(fname, ncname, platform, ...
+% 'levels', levels, 'region', region, 'scalefactor', scalefactor);
+%
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2007, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+narg = nargin;
+
+if narg == 3
+
+ levels = [];
+ region = [];
+ scalefactor = 10;
+
+else
+
+ [levels, region, scalefactor] = parseinput(varargin{:});
+
+end
+
+data.filename = fname;
+data.ncname = ncname;
+data.platform = platform;
+data.levels = levels;
+data.region = region;
+data.scalefactor = scalefactor;
+
+f = netcdf(ncname,'nowrite');
+data.platforms = f{'ObservationTypes'}(:);
+data.timeunits = f{'time_bounds'}.units(:);
+close(f);
+
+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);
+
+if isempty(obs)
+ clf;
+ axis([0 1 0 1]); axis image
+ h = text(0.5,0.5,sprintf('%s has no %s data',data.filename,data.platform));
+ set(h,'Interpreter','none')
+ return
+end
+
+obs.times = obs.times + data.timeorigin;
+t1 = datestr(min(obs.times));
+t2 = datestr(max(obs.times));
+
+clf;
+axlims = DrawBackground( obs );
+
+
+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));
+
+legh = [];
+legstr = {};
+
+if ~ isempty(goodUV)
+ hgood = obs2plot(obs, goodUV, [0 0 0] );
+ legh = [legh; hgood];
+ legstr{length(legstr)+1} = sprintf('%d good',length(goodUV));
+end
+
+if ~ isempty(baadUV)
+ hbaadUV = obs2plot(obs, baadUV, [1 0 0] );
+ legh = [legh; hbaadUV];
+ legstr{length(legstr)+1} = sprintf('%d badUbadV',length(baadUV));
+end
+
+if ~ isempty(goodU)
+ hgoodU = obs2plot(obs, goodU, [0 1 0] );
+ legh = [legh; hgoodU];
+ legstr{length(legstr)+1} = sprintf('%d goodUbadV',length(goodU));
+end
+
+if ~ isempty(goodV)
+ hgoodV = obs2plot(obs, goodV, [0 0 1] );
+ legh = [legh; hgoodV];
+ legstr{length(legstr)+1} = sprintf('%d badUgoodV',length(goodV));
+end
+
+h = title({sprintf('%s %s %s',t1,platform,t2), ...
+ sprintf('levels %s ',obs.levelstring)});
+set(h,'FontSize',18)
+
+h = xlabel(data.filename); set(h,'Interpreter','none');
+
+legend(legh,legstr,'Location','NorthWestOutside')
+
+hold off;
+
+
+
+function axlims = DrawBackground( obs )
+%======================================================================
+
+if isempty(obs.region)
+ % Figure out bounds of the data
+ axlims = [ min(obs.lon) max(obs.lon) min(obs.lat) max(obs.lat) ] ;
+else
+ axlims = obs.region;
+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])
+axis image
+
+% It is nice to know where the land is
+worldmap('light');
+hold on;
+
+% Plot a wind vector pair of known length for reference
+tx = axlims(1)+dx;
+ty = axlims(4)-dy;
+
+U = 10/obs.scalefactor;
+V = 0/obs.scalefactor;
+h = quiver(tx, ty, U, V, 0.0 ,'LineWidth',4.0,'Color','k');
+set(h,'LineWidth',3.0)
+
+U = 0/obs.scalefactor; V = 10/obs.scalefactor;
+h = quiver(tx, ty, U, V, 0.0 ,'LineWidth',4.0,'Color','k');
+set(h,'LineWidth',3.0)
+h = text(tx, ty-0.1*dy,'10 m/s','VerticalAlignment','top');
+
+
+
+function h1 = obs2plot(obs, mask, colspec )
+%======================================================================
+% Actually plots the wind vectors. to defeat Matlab's automatic
+% scaling, we use a scalefactor of 0.0 and manually scale the
+% data from the user input (or default).
+
+lon = obs.lon(mask);
+lat = obs.lat(mask);
+U = obs.U(mask)/obs.scalefactor;
+V = obs.V(mask)/obs.scalefactor;
+h1 = quiver(lon, lat, U, V, 0.0);
+h2 = plot(lon, lat, '.','MarkerSize',4);
+set(h1,'Color',colspec)
+set(h2,'Color',colspec)
+
+
+
+function obs = Load_Parse( data )
+%======================================================================
+% 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);
+
+% Find the types of data in the obs_sequence file for this epoch.
+
+platformIDs = unique(obsmat(:,1));
+uid = floor(platformIDs/100);
+vid = platformIDs - uid*100;
+Ustrings = data.platforms(uid,:);
+Vstrings = data.platforms(vid,:);
+
+nplatforms = length(uid);
+pid = [];
+obs = [];
+levelstring = [];
+regionstring = [];
+
+% This block divines the platform string and companion obs_kinds
+% from the netCDF file metadata - the only reason we need the netCDF file.
+
+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;
+ 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
+end
+
+inds = find(obsmat(:,1) == pid);
+
+if isempty(inds)
+ disp(sprintf('no %s observations (type %d) in %s',data.platform, pid, data.filename))
+ return
+end
+
+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);
+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);
+ Upr = Upr(inds);
+ Vpr = Vpr(inds);
+ Upo = Upo(inds);
+ Vpo = Vpo(inds);
+ 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)
+%======================================================================
+% 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))
+end
+
+npairs = length(varargin)/2;
+
+levels = [];
+region = [];
+scalefactor = 10.0;
+
+for i = 1:2:length(varargin)
+ switch lower(varargin{i})
+ case 'levels'
+ levels = varargin{i+1};
+ case 'region'
+ region = varargin{i+1};
+ case 'scalefactor'
+ scalefactor = varargin{i+1};
+ otherwise
+ error('Unknown parameter %s',varargin{i})
+ end
+end
+
+% Make sure the levels array has a top/bottom
+if ~ isempty(levels)
+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 2009-04-17 19:08:42 UTC (rev 3824)
+++ DART/trunk/diagnostics/threed_sphere/obs_diag.f90 2009-04-17 19:13:25 UTC (rev 3825)
@@ -46,7 +46,9 @@
use utilities_mod, only : open_file, close_file, register_module, &
file_exist, error_handler, E_ERR, E_WARN, E_MSG, &
initialize_utilities, logfileunit, nmlfileunit, timestamp, &
- find_namelist_in_file, check_namelist_read, nc_check
+ find_namelist_in_file, check_namelist_read, nc_check, &
+ next_file, find_textfile_dims, file_to_text
+
use sort_mod, only : sort
use typeSizes
@@ -503,14 +505,8 @@
ObsFileLoop : do ifile=1, Nepochs*4
!-----------------------------------------------------------------------
- ! TODO: NextFile only handles relative filenames (as of 12/20/2007).
- ! The assumption is that the file are organized obs_0001/obs_seq.final
- ! NextFile finds the first slash instead of the last slash. If everything
- ! were slaved off the last slash I suspect relative and absolute would work.
- ! TJH 2007/12/20
+ obs_seq_in_file_name = next_file(obs_sequence_name,ifile)
- obs_seq_in_file_name = NextFile(obs_sequence_name,ifile)
-
if ( file_exist(trim(obs_seq_in_file_name)) ) then
write(msgstring,*)'opening ', trim(obs_seq_in_file_name)
call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
@@ -3221,108 +3217,6 @@
- Function NextFile(fname,ifile)
- !----------------------------------------------------------------------
- ! The file name can take one of three forms:
- ! /absolute/path/to/nirvana/obs_001/obs_seq.final (absolute path)
- ! obs_0001/obs_seq.final (relative path)
- ! obs_seq.final (no path ... local)
- !
- ! If there is a '/' in the file name, we grab the portion before the
- ! slash and look for an underscore. Anything following the underscore
- ! is presumed to be the portion to increment.
- !
- ! If there is no slash AND ifile is > 1 ... we have already read
- ! the 'one and only' obs_seq.final file and we return 'done'.
- !----------------------------------------------------------------------
-
- character(len=129), intent(in) :: fname
- integer, intent(in) :: ifile
- character(len=129) :: NextFile
-
- integer, SAVE :: filenum = 0
- integer, SAVE :: dir_prec = 0
- character(len=129), SAVE :: dir_base
- character(len=129), SAVE :: filename
- character(len=stringlength), SAVE :: dir_ext
-
- character(len=129) :: dir_name
- integer :: slashindex, splitindex, i, strlen
-
- if (ifile == 1) then ! First time through ... find things.
-
- ! Start looking (right-to-left) for the 'slash'.
- ! Anything to the right of it must be a filename.
- ! Anything to the left must be the part that gets incremented.
-
- filename = adjustl(fname)
- NextFile = trim(filename)
- strlen = len_trim(filename)
- slashindex = 0
-
- SlashLoop : do i = strlen,1,-1
- if ( NextFile(i:i) == '/' ) then
- slashindex = i
- exit SlashLoop
- endif
- enddo SlashLoop
-
- if (slashindex > 0) then ! we have a directory structure
-
- dir_name = trim(fname(1:slashindex-1))
- filename = trim(fname(slashindex+1:129))
- strlen = len_trim(dir_name)
- splitindex = 0
-
- SplitLoop : do i = strlen,1,-1
- if ( dir_name(i:i) == '_' ) then
- splitindex = i
- exit SplitLoop
- endif
- enddo SplitLoop
-
- if (splitindex <= 0) then
- filenum = -1 ! indicates no next file
- else
- dir_base = dir_name(1:splitindex-1)
- dir_ext = dir_name(splitindex+1:slashindex-1)
- dir_prec = slashindex - splitindex - 1
- read(dir_ext,*) filenum
- endif
-
- else ! we have one single file - on the first trip through
-
- filenum = -1 ! indicates no next file
-
- endif
-
- else
-
- if (filenum < 0) then
- NextFile = 'doneDONEdoneDONE'
- else
-
- filenum = filenum + 1
- if (dir_prec == 1) then
- write(NextFile,'(a,''_'',i1.1,''/'',a)') trim(dir_base),filenum,trim(filename)
- elseif (dir_prec == 2) then
- write(NextFile,'(a,''_'',i2.2,''/'',a)') trim(dir_base),filenum,trim(filename)
- elseif (dir_prec == 3) then
- write(NextFile,'(a,''_'',i3.3,''/'',a)') trim(dir_base),filenum,trim(filename)
- elseif (dir_prec == 4) then
- write(NextFile,'(a,''_'',i4.4,''/'',a)') trim(dir_base),filenum,trim(filename)
- else
- write(NextFile,'(a,''_'',i5.5,''/'',a)') trim(dir_base),filenum,trim(filename)
- endif
-
- endif
-
- endif
-
- end Function NextFile
-
-
-
Subroutine ObsLocationsExist( printswitch )
! This routine checks for the existence of observation location files.
Added: DART/trunk/diagnostics/threed_sphere/obs_to_table.f90
===================================================================
--- DART/trunk/diagnostics/threed_sphere/obs_to_table.f90 (rev 0)
+++ DART/trunk/diagnostics/threed_sphere/obs_to_table.f90 2009-04-17 19:13:25 UTC (rev 3825)
@@ -0,0 +1,924 @@
+! Data Assimilation Research Testbed -- DART
+! Copyright 2004-2007, Data Assimilation Research Section
+! University Corporation for Atmospheric Research
+! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+
+program obs_to_table
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+!-----------------------------------------------------------------------
+! The programs defines a series of epochs (periods of time)
+!
+! All 'possible' obs_kinds are treated separately.
+!-----------------------------------------------------------------------
+
+use types_mod, only : r4, r8, digits12, MISSING_R8, MISSING_R4
+use obs_sequence_mod, only : read_obs_seq, obs_type, obs_sequence_type, get_first_obs, &
+ get_obs_from_key, get_obs_def, get_copy_meta_data, &
+ get_obs_time_range, get_time_range_keys, get_num_obs, &
+ get_next_obs, get_num_times, get_obs_values, init_obs, &
+ assignment(=), get_num_copies, static_init_obs_sequence, &
+ get_qc, destroy_obs_sequence, read_obs_seq_header, &
+ get_last_obs, destroy_obs, get_num_qc, get_qc_meta_data
+use obs_def_mod, only : obs_def_type, get_obs_def_error_variance, get_obs_def_time, &
+ get_obs_def_location, get_obs_kind, get_obs_name
+use obs_kind_mod, only : max_obs_kinds, get_obs_kind_var_type, get_obs_kind_name, &
+ do_obs_form_pair, add_wind_names, &
+ KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT
+use location_mod, only : location_type, get_location, set_location_missing, &
+ write_location, operator(/=), operator(==), &
+ set_location, is_location_in_region, VERTISUNDEF
+use time_manager_mod, only : time_type, set_date, set_time, get_time, print_time, &
+ set_calendar_type, get_calendar_string, print_date, &
+ operator(*), operator(+), operator(-), &
+ operator(>), operator(<), operator(/), &
+ operator(/=), operator(<=)
+use schedule_mod, only : schedule_type, set_regular_schedule, get_schedule_length, &
+ get_time_from_schedule
+use utilities_mod, only : open_file, close_file, register_module, &
+ file_exist, error_handler, E_ERR, E_WARN, E_MSG, &
+ initialize_utilities, nmlfileunit, timestamp, &
+ find_namelist_in_file, check_namelist_read, nc_check, &
+ next_file, find_textfile_dims, file_to_text
+
+use typeSizes
+use netcdf
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = '$URL$', &
+ revision = '$Revision$', &
+ revdate = '$Date$'
+
+!---------------------------------------------------------------------
+!---------------------------------------------------------------------
+
+integer, parameter :: stringlength = 32
+logical, parameter :: DEBUG = .false.
+
+!---------------------------------------------------------------------
+! variables associated with the observation
+!---------------------------------------------------------------------
+
+type(obs_sequence_type) :: seq
+type(obs_type) :: observation, next_obs
+type(obs_type) :: obs1, obsN
+type(obs_def_type) :: obs_def
+type(location_type) :: obs_loc, minl, maxl
+
+character(len = 129) :: obs_seq_in_file_name
+character(len = 129), allocatable, dimension(:) :: obs_seq_filenames
+
+! Storage with fixed size for observation space diagnostics
+real(r8), dimension(1) :: prior_mean, posterior_mean, prior_spread, posterior_spread
+real(r8) :: pr_mean, po_mean ! same as above, without useless dimension
+real(r8) :: pr_sprd, po_sprd ! same as above, without useless dimension
+
+! We are treating winds as a vector pair, but we are handling the
+! observations serially. Consequently, we exploit the fact that
+! the U observations are _followed_ by the V observations.
+
+real(r8) :: U_obs = 0.0_r8
+real(r8) :: U_obs_err_var = 0.0_r8
+type(location_type) :: U_obs_loc
+integer :: U_flavor
+integer :: U_type = KIND_V_WIND_COMPONENT ! intentional mismatch
+real(r8) :: U_pr_mean = 0.0_r8
+real(r8) :: U_pr_sprd = 0.0_r8
+real(r8) :: U_po_mean = 0.0_r8
+real(r8) :: U_po_sprd = 0.0_r8
+integer :: U_qc = 0
+
+integer :: obs_index, prior_mean_index, posterior_mean_index
+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
+character(len=129) :: obs_seq_read_format
+logical :: pre_I_format
+
+integer, dimension(2) :: key_bounds
+real(r8), dimension(1) :: obs
+real(r8) :: obs_err_var
+
+integer, allocatable, dimension(:) :: keys
+
+logical :: out_of_range, is_there_one, keeper
+
+!---------------------------------------------------------------------
+! variables associated with quality control
+!
+! qc_index reflects the 'original' QC value of the observation, if any.
+! Most frequently represents the value NCEP assigned to their
+! observations.
+!
+! dart_qc_index
+! 0 observation assimilated
+! 1 observation evaluated only
+! --- everything above this means the prior and posterior are OK
+! 2 assimilated, but the posterior forward operator failed
+! 3 Evaluated only, but the posterior forward operator failed
+! --- everything above this means only the prior is OK
+! 4 prior forward operator failed
+! 5 not used
+! 6 prior QC rejected
+! 7 outlier rejected
+! 8+ reserved for future use
+
+integer :: qc_index, dart_qc_index
+integer :: qc_integer
+integer, parameter :: QC_MAX = 7
+integer, parameter :: QC_MAX_PRIOR = 3
+integer, parameter :: QC_MAX_POSTERIOR = 1
+integer, dimension(0:QC_MAX) :: qc_counter = 0
+real(r8), allocatable, dimension(:) :: qc
+real(r8), allocatable, dimension(:) :: copyvals
+
+!-----------------------------------------------------------------------
+! Namelist with (some scalar) default values
+!-----------------------------------------------------------------------
+
+character(len = 129) :: obs_sequence_name = 'obs_seq.final'
+
+real(r8) :: lonlim1= MISSING_R8, lonlim2= MISSING_R8
+real(r8) :: latlim1= MISSING_R8, latlim2= MISSING_R8
+
+logical :: verbose = .false.
+
+namelist /obs_to_table_nml/ obs_sequence_name, lonlim1, lonlim2, &
+ latlim1, latlim2, verbose
+
+!-----------------------------------------------------------------------
+! Variables used to accumulate the statistics.
+!-----------------------------------------------------------------------
+
+integer, parameter :: Ncopies = 10
+character(len = stringlength), dimension(Ncopies) :: copy_names = &
+ (/ 'Nposs ', 'Nused ', 'NbadQC ', 'NbadIZ ', 'NbadUV ', &
+ 'NbadLV ', 'rmse ', 'bias ', 'spread ', 'totalspread' /)
+
+!-----------------------------------------------------------------------
+! General purpose variables
+!-----------------------------------------------------------------------
+
+integer :: iepoch, ivar, ifile, num_obs_in_epoch
+real(r8) :: obslon, obslat, obslevel, obsloc3(3)
+
+integer :: obsindex, iunit, io
+integer :: Nepochs
+
+character(len = stringlength), pointer, dimension(:) :: my_obs_kind_names
+
+type(schedule_type) :: schedule
+type(time_type) :: TimeMin, TimeMax ! of the entire period of interest
+type(time_type) :: beg_time, end_time ! of the particular bin
+type(time_type) :: seqT1, seqTN ! first,last time in entire observation sequence
+type(time_type) :: obs_time
+
+character(len = 129) :: ncName, windName, msgstring, calendarstring
+
+integer :: Nidentity = 0 ! identity observations
+
+!=======================================================================
+! Get the party started
+!=======================================================================
+
+call initialize_utilities('obs_to_table')
+call register_module(source,revision,revdate)
+call static_init_obs_sequence() ! Initialize the obs sequence module
+
+!----------------------------------------------------------------------
+! 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'
+!----------------------------------------------------------------------
+
+num_obs_kinds = add_wind_names(my_obs_kind_names)
+
+!----------------------------------------------------------------------
+! Read the namelist
+!----------------------------------------------------------------------
+
+call find_namelist_in_file('input.nml', 'obs_to_table_nml', iunit)
+read(iunit, nml = obs_to_table_nml, iostat = io)
+call check_namelist_read(iunit, io, 'obs_to_table_nml')
+
+! Record the namelist values used for the run ...
+write(nmlfileunit, nml=obs_to_table_nml)
+write( * , nml=obs_to_table_nml)
+
+!----------------------------------------------------------------------
+! SetSchedule rectifies user input and the final binning sequence.
+!----------------------------------------------------------------------
+
+call set_regular_schedule(schedule) ! also sets calendar type
+
+Nepochs = get_schedule_length(schedule)
+call get_time_from_schedule(TimeMin, schedule, 1, 1)
+call get_time_from_schedule(TimeMax, schedule, Nepochs, 1)
+call get_calendar_string(calendarstring)
+
+U_obs_loc = set_location_missing()
+minl = set_location(lonlim1, latlim1, 0.0_r8, VERTISUNDEF) ! vertical unimportant
+maxl = set_location(lonlim2, latlim2, 0.0_r8, VERTISUNDEF) ! vertical unimportant
+
+!----------------------------------------------------------------------
+! Prepare the variables
+!----------------------------------------------------------------------
+
+allocate(obs_seq_filenames(Nepochs*4))
+obs_seq_filenames = 'null'
+
+ObsFileLoop : do ifile=1, Nepochs*4
+!-----------------------------------------------------------------------
+
+ obs_seq_in_file_name = next_file(obs_sequence_name,ifile)
+
+ if ( file_exist(trim(obs_seq_in_file_name)) ) then
+ write(msgstring,*)'opening ', trim(obs_seq_in_file_name)
+ call error_handler(E_MSG,'obs_to_table',msgstring,source,revision,revdate)
+ else
+ write(msgstring,*)trim(obs_seq_in_file_name),&
+ ' does not exist. Finishing up.'
+ call error_handler(E_MSG,'obs_to_table',msgstring,source,revision,revdate)
+ exit ObsFileLoop
+ endif
+
+ ! 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(obs_seq_in_file_name) ! Lahey requirement
+ 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, &
+ obs_seq_file_id, obs_seq_read_format, pre_I_format, &
+ close_the_file = .true.)
+
+ ! Initialize some (individual) observation variables
+
+ call init_obs( obs1, num_copies, num_qc) ! First obs in sequence
+ call init_obs( obsN, num_copies, num_qc) ! Last obs in sequence
+ call init_obs(observation, num_copies, num_qc) ! current obs
+ call init_obs( next_obs, num_copies, num_qc) ! duh ...
+
+ if (num_qc > 0) allocate( qc(num_qc) )
+ if (num_copies > 0) allocate( copyvals(num_copies) )
+
+ if ( DEBUG ) then
+ 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(obs_seq_read_format)
+ write(*,*)'pre_I_format is ',pre_I_format
+ write(*,*)
+ endif
+
+ ! Read in the entire observation sequence
+
+ call read_obs_seq(obs_seq_in_file_name, 0, 0, 0, seq)
+
+ !--------------------------------------------------------------------
+ ! Determine the time encompassed in the observation sequence.
+ !--------------------------------------------------------------------
+
+ is_there_one = get_first_obs(seq, obs1)
+ if ( .not. is_there_one ) then
+ call error_handler(E_ERR,'obs_to_table','No first observation in sequence.', &
+ source,revision,revdate)
+ endif
+ call get_obs_def(obs1, obs_def)
+ seqT1 = get_obs_def_time(obs_def)
+
+ is_there_one = get_last_obs(seq, obsN)
+ if ( .not. is_there_one ) then
+ call error_handler(E_ERR,'obs_to_table','No last observation in sequence.', &
+ source,revision,revdate)
+ endif
+ call get_obs_def(obsN, obs_def)
+ seqTN = get_obs_def_time(obs_def)
+
+ if ( verbose ) then
+ call print_time(seqT1,'First observation time')
+ call print_time(seqTN,'Last observation time')
+ call print_time(TimeMin,'TimeMin from input')
+ call print_time(TimeMax,'TimeMax from input')
+ call print_date(seqT1,'First observation date')
+ call print_date(seqTN,'Last observation date')
+ call print_date(TimeMin,'DateMin from input')
+ call print_date(TimeMax,'DateMax from input')
+ endif
+
+ !--------------------------------------------------------------------
+ ! If the last observation is before the period of interest, move on.
+ !--------------------------------------------------------------------
+
+ if ( seqTN < TimeMin ) then
+ if (verbose) write(*,*)'seqTN < TimeMin ... trying next file.'
+ call destroy_obs(obs1)
+ call destroy_obs(obsN)
+ call destroy_obs(observation)
+ call destroy_obs(next_obs)
+ call destroy_obs_sequence(seq)
+ if (allocated(qc)) deallocate( qc )
+ if (allocated(copyvals)) deallocate( copyvals )
+ cycle ObsFileLoop
+ else
+ if (verbose) write(*,*)'seqTN > TimeMin ... using ', trim(obs_seq_in_file_name)
+ endif
+
+ !--------------------------------------------------------------------
+ ! If the first observation is after the period of interest, finish.
+ !--------------------------------------------------------------------
+
+ if ( seqT1 > TimeMax ) then
+ if (verbose) write(*,*)'seqT1 > TimeMax ... stopping.'
+ call destroy_obs(obs1)
+ call destroy_obs(obsN)
+ call destroy_obs(observation)
+ call destroy_obs(next_obs)
+ call destroy_obs_sequence(seq)
+ if (allocated(qc)) deallocate( qc )
+ if (allocated(copyvals)) deallocate( copyvals )
+ exit ObsFileLoop
+ else
+ if (verbose) write(*,*)'seqT1 < TimeMax ... using ',trim(obs_seq_in_file_name)
+ endif
+
+ !--------------------------------------------------------------------
+ ! Find the index of obs, ensemble mean, spread ... etc.
+ !--------------------------------------------------------------------
+ ! Only require obs_index to be present; this allows the program
+ ! to be run on obs_seq.in files which have no means or spreads.
+ ! You can still plot locations, but that's it.
+ !--------------------------------------------------------------------
+
+ call SetIndices( obs_index, qc_index, dart_qc_index )
+
+ !====================================================================
+ EpochLoop : do iepoch = 1, Nepochs
+ !====================================================================
+
+ call get_time_from_schedule(beg_time, schedule, iepoch, 1)
+ call get_time_from_schedule(end_time, schedule, iepoch, 2)
+
+ call get_obs_time_range(seq, beg_time, end_time, key_bounds, &
+ num_obs_in_epoch, out_of_range)
+
+ if( num_obs_in_epoch == 0 ) then
+ if (verbose) write(*,*)' No observations in epoch ',iepoch,' cycling ...'
+ cycle EpochLoop
+ endif
+
+ write(*,*)'num_obs_in_epoch (', iepoch, ') = ', num_obs_in_epoch
+
+ ! Append epoch number to output file name
+ write(windName,'(a,i3.3,a)') 'wind_vectors.', iepoch, '.dat'
+ iunit = open_file(windName, form='formatted', action='rewind')
+
+ allocate(keys(num_obs_in_epoch))
+
+ call get_time_range_keys(seq, key_bounds, num_obs_in_epoch, keys)
+
+ !-----------------------------------------------------------------
+ ObservationLoop : do obsindex = 1, num_obs_in_epoch
+ !-----------------------------------------------------------------
+
+ ! 'flavor' is from the 'master list' in the obs_kind_mod.f90
+ ! each obs_seq.final file has their own private kind - which
+ ! gets mapped to the 'master list', if you will.
+
+ call get_obs_from_key(seq, keys(obsindex), observation)
+ call get_obs_values(observation, obs, obs_index)
+ call get_obs_def(observation, obs_def)
+ call get_qc(observation, qc)
+
+ flavor = get_obs_kind(obs_def)
+ obs_err_var = get_obs_def_error_variance(obs_def)
+ obs_time = get_obs_def_time(obs_def)
+ obs_loc = get_obs_def_location(obs_def)
+ obsloc3 = get_location(obs_loc)
+
+ obslon = obsloc3(1) ! [ 0, 360]
+ obslat = obsloc3(2) ! [-90, 90]
+ obslevel = obsloc3(3) ! variable-dependent
+
+ !--------------------------------------------------------------
+ ! We have one Region of interest
+ !--------------------------------------------------------------
+
+ keeper = is_location_in_region( obs_loc, minl, maxl )
+
+ if ( .not. keeper ) cycle ObservationLoop
+
+ !--------------------------------------------------------------
+ ! Convert the DART QC data
+ !--------------------------------------------------------------
+
+ if (dart_qc_index > 0) then
+ qc_integer = min( nint(qc(dart_qc_index)), QC_MAX )
+ else if (qc_index > 0) then ! If there is no dart_qc in obs_seq ...
+ qc_integer = nint(qc(qc_index))
+ else ! If there is no original QC value ...
+ qc_integer = 999
+ endif
+
+ !--------------------------------------------------------------
+ ! Summary of observation knowledge at this point
+ !--------------------------------------------------------------
+
+ if ( DEBUG ) then
+ write(*,*)'observation # ',obsindex
+ write(*,*)'key ',keys(obsindex)
+ write(*,*)'obs_flavor ',flavor
+ write(*,*)'obs_err_var ',obs_err_var
+ write(*,*)'lon/lat/level ',obslon,obslat,obslevel
+ write(*,*)'obs(1) ',obs(1)
+ write(*,*)'qc ',qc
+ endif
+
+ !--------------------------------------------------------------
+ ! If it is a U wind component, all we need to do is save it.
+ ! It will be matched up with the subsequent V component.
+ ! At some point we have to remove the dependency that the
+ ! U component MUST preceed the V component.
+ !--------------------------------------------------------------
+
+ if ( get_obs_kind_var_type(flavor) == KIND_U_WIND_COMPONENT ) then
+
+ U_obs = obs(1)
+ U_obs_err_var = obs_err_var
+ U_obs_loc = obs_loc
+ U_flavor = flavor
+ U_type = KIND_U_WIND_COMPONENT
+ U_qc = qc_integer
+
+ cycle ObservationLoop
+
+ endif
+
+
+ !-----------------------------------------------------------
+ ! Additional work for horizontal wind (given U,V)
+ !-----------------------------------------------------------
+
+ ObsIsWindCheck: if ( get_obs_kind_var_type(flavor) == KIND_V_WIND_COMPONENT ) then
+
+ ! The big assumption is that the U wind component has
+ ! immediately preceeded the V component and has been saved.
+ !
+ ! We check for observation compatibility and the index for
+ ! this wind 'kind' ... not originally in the max_obs_kind namelist.
+ ! this will be the 'wflavor' (wind) flavor.
+
+ if ((obs_loc == U_obs_loc) .and. &
+ do_obs_form_pair(flavor,U_flavor,keys(obsindex),my_obs_kind_names,wflavor)) then
+
+ call WritePairs( iunit, U_flavor, flavor, &
+ obs_time, obslon, obslat, obslevel, &
+ U_qc, qc_integer, U_obs, obs(1))
+ else
+ write(*,*)'time series : V with no U obs index ', keys(obsindex)
+ wflavor = -99
+ endif
+
+ endif ObsIsWindCheck
+
+ !-----------------------------------------------------------------
+ enddo ObservationLoop
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list