[Dart-dev] [4550] DART/trunk: Adding first cut at a program to subset an obs sequence file
nancy at ucar.edu
nancy at ucar.edu
Tue Nov 16 16:36:18 MST 2010
Revision: 4550
Author: thoar
Date: 2010-11-16 16:36:18 -0700 (Tue, 16 Nov 2010)
Log Message:
-----------
Adding first cut at a program to subset an obs sequence file
for unique locations and times for a particular observation type.
It is expected to be used to derive a list of 'stations' that
have good temporal coverage to be used for forecast verification.
Modified Paths:
--------------
DART/trunk/models/wrf/work/input.nml
Added Paths:
-----------
DART/trunk/diagnostics/matlab/plot_coverage.m
DART/trunk/models/wrf/work/mkmf_obs_seq_coverage
DART/trunk/models/wrf/work/path_names_obs_seq_coverage
DART/trunk/obs_sequence/obs_seq_coverage.f90
-------------- next part --------------
Added: DART/trunk/diagnostics/matlab/plot_coverage.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_coverage.m (rev 0)
+++ DART/trunk/diagnostics/matlab/plot_coverage.m 2010-11-16 23:36:18 UTC (rev 4550)
@@ -0,0 +1,257 @@
+function obs = plot_coverage(fname)
+%% plot_coverage examines the spatial and temporal coverage for an observation type.
+% obs_seq_coverage must create the netCDF file read by this routine.
+%
+% fname = 'obs_coverage.nc';
+%
+% plot_coverage(fname)
+
+%% 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+if (exist(fname,'file') ~= 2)
+ error('%s does not exist.',fname)
+end
+
+locs = nc_varget(fname,'location');
+
+obs.fname = fname;
+obs.ObsTypeString = nc_attget(fname,nc_global,'obs_of_interest');
+obs.stations = nc_varget(fname,'stations');
+obs.ntimes = nc_varget(fname,'ntimes');
+obs.times = nc_varget(fname,'time');
+obs.timeunits = nc_attget(fname,'time','units');
+
+obs.lons = locs(:,1);
+obs.lats = locs(:,2);
+obs.vert = locs(:,3);
+obs.region = [ min(obs.lons) max(obs.lons) min(obs.lats) max(obs.lats) ];
+timebase = sscanf(obs.timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+obs.timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+obs.times = obs.times + obs.timeorigin;
+
+allt = obs.times(:);
+obs.tmin = datestr(min(allt));
+obs.tmax = datestr(max(allt));
+clear allt;
+
+inds = find(obs.ntimes == 8);
+fprintf('There are %d stations with precisely 8 observation times.\n',length(inds))
+
+%% Now pack the data in the same fashion as the cell array of column labels.
+
+obs.lonindex = 1;
+obs.latindex = 2;
+obs.zindex = 3;
+obs.ntindex = 4;
+obs.Iindex = 5;
+
+obs.colnames{obs.lonindex} = 'longitude';
+obs.colnames{obs.latindex} = 'latitude';
+obs.colnames{obs.zindex} = 'vertical';
+obs.colnames{obs.ntindex} = 'num times';
+obs.colnames{obs.Iindex} = 'stationID';
+
+global obsmat
+obsmat = zeros(length(obs.lons),5);
+obsmat(:,obs.lonindex ) = obs.lons;
+obsmat(:,obs.latindex ) = obs.lats;
+obsmat(:,obs.zindex ) = obs.vert;
+obsmat(:,obs.ntindex ) = obs.ntimes;
+obsmat(:,obs.Iindex ) = obs.stations;
+
+global timemat
+timemat = obs.times;
+
+%% create the linked plots
+
+figure1 = figure(1); clf(figure1); orient tall; wysiwyg
+
+%----------------------------------------------------------------------
+% Top half is a 2D scatterplot
+%----------------------------------------------------------------------
+
+axes0 = axes('Parent',figure1,'OuterPosition',[0 0.4 1 0.57],'FontSize',16);
+view(axes0,[0 90]);
+grid(axes0,'on');
+hold(axes0,'all');
+
+xstring = sprintf('obsmat(:,%d)',obs.lonindex);
+ystring = sprintf('obsmat(:,%d)',obs.latindex);
+
+scatter(obsmat(:,obs.lonindex), obsmat(:,obs.latindex), 18, obsmat(:,obs.ntindex), 'filled', ...
+ 'Parent',axes0,'DisplayName','observation locations', ...
+ 'XDataSource',xstring, ...
+ 'YDataSource',ystring);
+% myworldmap(obs);
+
+xlabel(obs.colnames{obs.lonindex});
+ylabel(obs.colnames{obs.latindex});
+
+hc = colorbar;
+set(get(hc,'YLabel'),'String','# of obs times','Interpreter','none')
+
+ntmax = length(obs.times);
+h = title({sprintf('%s locations',obs.ObsTypeString), ...
+ sprintf('%s ---> %s',obs.tmin,obs.tmax) });
+set(h,'Interpreter','none')
+
+%----------------------------------------------------------------------
+% Bottom half is a "histogram"
+%----------------------------------------------------------------------
+
+axes1 = axes('Parent',figure1,'OuterPosition',[0 0.0 1 0.38],'FontSize',14);
+box(axes1,'on');
+hold(axes1,'all');
+
+xstring = sprintf('obsmat(:,%d)',obs.Iindex);
+ystring = sprintf('obsmat(:,%d)',obs.ntindex);
+
+scatter(obsmat(:,obs.Iindex), obsmat(:,obs.ntindex), ...
+ 'XDataSource',xstring, ...
+ 'YDataSource',ystring);
+xlabel('''station''')
+ylabel('# of obs times')
+
+refreshdata
+linkdata on
+
+%% Create figure for ancillary plots
+% ====================================================================
+
+figure2 = figure(2); clf(figure2); orient tall; wysiwyg
+
+%% Create axes for time VS. station
+axes2 = axes('Parent',figure2,'OuterPosition',[0 0 1 1],'FontSize',14);
+set(axes2,'XAxisLocation','bottom')
+box(axes2,'on');
+hold(axes2,'all');
+
+hist(obsmat(:,obs.ntindex),single([1:max(obsmat(:,obs.ntindex))]))
+h = title(obs.ObsTypeString);
+set(h,'Interpreter','none')
+xlabel('number of observations at a location')
+ylabel('count')
+
+%% Create figure for ancillary plots
+% ====================================================================
+
+figure3 = figure(3); clf(figure3); orient tall; wysiwyg
+
+%% Create axes for time VS. station
+axes3 = axes('Parent',figure3,'OuterPosition',[0 0 1 1],'FontSize',14);
+set(axes3,'XAxisLocation','bottom')
+box(axes3,'on');
+hold(axes3,'all');
+
+ystring = sprintf('obsmat(:,%d)',obs.Iindex);
+
+h = plot(timemat, obsmat(:,obs.Iindex),'kd','Parent',axes3, ...
+ 'DisplayName','time vs station', ...
+ 'YDataSource',ystring, ...
+ 'XDataSource','timemat');
+
+ax = axis;
+datetick(axes3,'x',20);
+axis([min(timemat(:))-0.25 max(timemat(:))+0.25 ax(3) ax(4)])
+h = title(obs.ObsTypeString);
+set(h,'Interpreter','none')
+xlabel('Date - DD/MM/YY')
+ylabel('station ID')
+
+%% thats it folks
+
+% refreshdata
+% linkdata on
+
+
+
+function h = myworldmap(obs)
+
+%%--------------------------------------------------------------------------
+% GET THE ELEVATION DATA AND SET UP THE ASSOCIATED COORDINATE DATA
+%---------------------------------------------------------------------------
+
+load topo; % GET Matlab-native [180x360] ELEVATION DATASET
+lats = -89.5:89.5; % CREATE LAT ARRAY FOR TOPO MATRIX
+lons = 0.5:359.5; % CREATE LON ARRAY FOR TOPO MATRIX
+nlon = length(lons);
+nlat = length(lats);
+
+%%--------------------------------------------------------------------------
+% IF WE NEED TO SWAP HEMISPHERES, DO SO NOW.
+% If we didn't explicitly tell it, make a guess.
+%---------------------------------------------------------------------------
+
+ax = axis;
+
+if (ax(1) < -2)
+ lons = lons - 180.0;
+ topo = [ topo(:,nlon/2+1:nlon) topo(:,1:nlon/2) ];
+end
+
+%%--------------------------------------------------------------------------
+% We need to determine the geographic subset of the elevation matrix.
+%---------------------------------------------------------------------------
+
+lon_ind1 = find(ax(1) <= lons, 1);
+lon_ind2 = find(ax(2) <= lons, 1);
+lat_ind1 = find(ax(3) <= lats, 1);
+lat_ind2 = find(ax(4) <= lats, 1);
+
+if (isempty(lon_ind1)), lon_ind1 = 1; end;
+if (isempty(lon_ind2)), lon_ind2 = nlon; end;
+if (isempty(lat_ind1)), lat_ind1 = 1; end;
+if (isempty(lat_ind2)), lat_ind2 = nlat; end;
+
+elev = topo(lat_ind1:lat_ind2,lon_ind1:lon_ind2);
+x = lons(lon_ind1:lon_ind2);
+y = lats(lat_ind1:lat_ind2);
+
+%%--------------------------------------------------------------------------
+% Contour the "subset" - and give the whole thing an appropriate zlevel
+% so the continents are either at the top of the plot (for depth), or
+% the bottom (for just about everything else.
+% There are differences between 6.5 and 7.0 that make changing the colors
+% of the filled contours a real pain.
+%---------------------------------------------------------------------------
+
+orgholdstate = ishold;
+hold on;
+
+switch lower(obs.Zunits)
+ case 'depth'
+ set(gca,'Zdir','reverse')
+ zlevel = min(ax(5:6));
+ case 'pressure'
+ zlevel = max(ax(5:6));
+ set(gca,'Zdir','reverse')
+ otherwise
+ set(gca,'Zdir','normal')
+ zlevel = min(ax(5:6));
+end
+
+fcolor = [0.7 0.7 0.7]; % light grey
+
+[c,h] = contourf(x,y,elev,[0.0 0.0],'k-');
+
+h_patch = get(h, 'Children');
+
+for i = 1:numel(h_patch)
+ y = get(h_patch(i), 'YData');
+ s = size(y);
+ set(h_patch(i), 'ZData', zlevel*ones(s),'FaceColor',fcolor);
+ set(h_patch(i),'AlphaDataMapping','none','FaceVertexAlphaData',0.3)
+ set(h_patch(i),'FaceAlpha',0.3)
+end
+
+if (orgholdstate == 0), hold off; end;
+
+
Property changes on: DART/trunk/diagnostics/matlab/plot_coverage.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author URL Id
Added: svn:eol-style
+ native
Modified: DART/trunk/models/wrf/work/input.nml
===================================================================
--- DART/trunk/models/wrf/work/input.nml 2010-11-15 23:25:34 UTC (rev 4549)
+++ DART/trunk/models/wrf/work/input.nml 2010-11-16 23:36:18 UTC (rev 4550)
@@ -286,7 +286,18 @@
verbose = .false.,
/
-&restart_file_tool_nml
+&obs_seq_coverage_nml
+ obs_sequence_list = 'obs_coverage_list.txt',
+ obs_sequence_name = '',
+ obs_of_interest = 'METAR_U_10_METER_WIND',
+ lonlim1 = 0.0,
+ lonlim2 = 360.0,
+ latlim1 = -90.0,
+ latlim2 = 90.0,
+ verbose = .true.,
+ /
+
+&restart_file_tool_nml
input_file_name = "filter_restart",
output_file_name = "filter_updated_restart",
ens_size = 1,
Added: DART/trunk/models/wrf/work/mkmf_obs_seq_coverage
===================================================================
--- DART/trunk/models/wrf/work/mkmf_obs_seq_coverage (rev 0)
+++ DART/trunk/models/wrf/work/mkmf_obs_seq_coverage 2010-11-16 23:36:18 UTC (rev 4550)
@@ -0,0 +1,19 @@
+#!/bin/csh
+#
+# 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
+# http://www.image.ucar.edu/DAReS/DART/DART_download
+#
+# $Id$
+
+../../../mkmf/mkmf -p obs_seq_coverage -t ../../../mkmf/mkmf.template -c"-Duse_netCDF" \
+ -a "../../.." $* path_names_obs_seq_coverage
+
+exit $status
+
+# <next few lines under version control, do not edit>
+# $URL:
+# http://subversion.ucar.edu/DAReS/DART/trunk/models/wrf/work/mkmf_obs_seq_coverage $
+# $Revision$
+# $Date$
+
Property changes on: DART/trunk/models/wrf/work/mkmf_obs_seq_coverage
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author URL Id
Added: svn:eol-style
+ native
Added: DART/trunk/models/wrf/work/path_names_obs_seq_coverage
===================================================================
--- DART/trunk/models/wrf/work/path_names_obs_seq_coverage (rev 0)
+++ DART/trunk/models/wrf/work/path_names_obs_seq_coverage 2010-11-16 23:36:18 UTC (rev 4550)
@@ -0,0 +1,14 @@
+assim_model/assim_model_mod.f90
+common/types_mod.f90
+location/threed_sphere/location_mod.f90
+models/wrf/model_mod.f90
+models/wrf/module_map_utils.f90
+mpi_utilities/null_mpi_utilities_mod.f90
+obs_def/obs_def_mod.f90
+obs_kind/obs_kind_mod.f90
+obs_sequence/obs_sequence_mod.f90
+obs_sequence/obs_seq_coverage.f90
+random_nr/random_nr_mod.f90
+random_seq/random_seq_mod.f90
+time_manager/time_manager_mod.f90
+utilities/utilities_mod.f90
Added: DART/trunk/obs_sequence/obs_seq_coverage.f90
===================================================================
--- DART/trunk/obs_sequence/obs_seq_coverage.f90 (rev 0)
+++ DART/trunk/obs_sequence/obs_seq_coverage.f90 2010-11-16 23:36:18 UTC (rev 4550)
@@ -0,0 +1,1128 @@
+! 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
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+program obs_seq_coverage
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+!-----------------------------------------------------------------------
+! This program queries a bunch of obs_seq.xxxx files and tries to
+! figure out 'station coverage' ... what locations are consistently
+! reported through time. Absolutely a 'reverse-engineering exercise'.
+!
+! The observation sequence file only contains lat/lon/level/which_vert,
+! so this is all we have to work with.
+!
+! From: Soyoung Ha <syha at ucar.edu>
+! Date: November 5, 2010 12:00:22 PM MDT
+! To: Tim Hoar <thoar at ucar.edu>, Nancy Collins <nancy at ucar.edu>
+! Subject: obs_seq.out on bluefire
+!
+! Hi Tim and Nancy,
+!
+! I just found that my current DART codes were compiled with r8, not r4.
+! But the MADIS obs netcdf files have lat/lon/elevation information in "float"
+! and only time information in "double".
+! So, although I see the "double" accuracy in my obs_seq.out,
+! I guess those were artificially generated by DART.
+! My obs_seq.out files which were preprocessed using wrf_dart_obs_preprocess
+! every 3 hrs are all available in be:/ptmp/syha/SFCDA/OBS_SEQ/Jun08/Preproc_3hrly/
+! for a month-long period of Jun 2008.
+! Just for your information, I put an ascii header file named "20080601_0000.head.txt"
+! in the same place, which is what I've got from doing "ncdump -h" for the original
+! MADIS metar netcdf file valid at the initial cycle.
+
+!-----------------------------------------------------------------------
+
+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_def, get_copy_meta_data, &
+ get_next_obs, init_obs, init_obs_sequence, &
+ assignment(=), get_num_copies, static_init_obs_sequence, &
+ get_qc, destroy_obs_sequence, read_obs_seq_header, &
+ destroy_obs, get_qc_meta_data
+use obs_def_mod, only : obs_def_type, get_obs_def_time, get_obs_kind, write_obs_def, &
+ get_obs_def_location, set_obs_def_time, &
+ set_obs_def_location, set_obs_def_kind, set_obs_def_error_variance
+use obs_kind_mod, only : max_obs_kinds, get_obs_kind_name, get_obs_kind_index, &
+ write_obs_kind
+use location_mod, only : location_type, get_location, set_location_missing, &
+ write_location, operator(/=), operator(==), &
+ set_location, is_location_in_region, query_location, &
+ nc_write_location_atts, nc_get_location_varids, &
+ nc_write_location, LocationDims
+use time_manager_mod, only : time_type, set_date, set_time, get_time, print_time, &
+ set_time_missing, operator(>), operator(<), operator(==), &
+ operator(<=), operator(-), operator(+), operator(/=)
+use utilities_mod, only : get_unit, 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, get_next_filename, find_textfile_dims, &
+ file_to_text, do_nml_file, do_nml_term
+
+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$'
+
+!---------------------------------------------------------------------
+!---------------------------------------------------------------------
+
+type station
+ integer :: obs_type
+ type(location_type) :: location
+ type(time_type) :: first_time
+ type(time_type) :: last_time
+ integer :: ntimes
+ type(time_type), pointer :: times(:)
+end type station
+
+type(station), allocatable, dimension(:) :: stations
+integer :: num_stations ! This is the current number of unique locations
+integer :: max_stations ! This is the largest possible number of uniq locs
+integer :: station_id ! the index (into stations) of an existing location
+integer :: timeindex ! the index (into the time array of a station)
+
+integer, parameter :: stringlength = 32
+
+!---------------------------------------------------------------------
+! variables associated with the observation
+!---------------------------------------------------------------------
+
+type(obs_sequence_type) :: seq
+type(obs_type) :: obs1, obs2
+type(obs_def_type) :: obs_def
+type(location_type) :: obs_loc, minl, maxl
+real(r8), dimension(LocationDims) :: locarray
+
+character(len = 129) :: obs_seq_in_file_name
+character(len = 129), allocatable, dimension(:) :: obs_seq_filenames
+
+integer :: flavor, flavor_of_interest
+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
+logical :: last_ob_flag
+
+!-----------------------------------------------------------------------
+! Namelist with (some scalar) default values
+!-----------------------------------------------------------------------
+
+character(len = 129) :: obs_sequence_name = 'obs_seq.final'
+character(len = 129) :: obs_sequence_list = ''
+character(len = 129) :: obs_of_interest = 'all'
+
+real(r8) :: lonlim1= MISSING_R8, lonlim2= MISSING_R8
+real(r8) :: latlim1= MISSING_R8, latlim2= MISSING_R8
+
+logical :: debug = .false. ! undocumented ... on purpose
+logical :: verbose = .false.
+
+namelist /obs_seq_coverage_nml/ obs_sequence_name, obs_sequence_list, &
+ lonlim1, lonlim2, latlim1, latlim2, &
+ obs_of_interest, verbose, debug
+
+!-----------------------------------------------------------------------
+! Quantities of interest
+!-----------------------------------------------------------------------
+
+integer, parameter :: Ncopies = 1
+integer :: allNcopies
+character(len=stringlength), dimension(Ncopies) :: copy_names = &
+ (/ 'observation error variance' /)
+
+character(len=stringlength), allocatable, dimension(:) :: module_obs_copy_names
+character(len=stringlength), allocatable, dimension(:) :: obs_copy_names
+character(len=stringlength), allocatable, dimension(:) :: module_qc_copy_names
+character(len=stringlength), allocatable, dimension(:) :: qc_copy_names
+character(len=stringlength), dimension(max_obs_kinds) :: obs_kind_names
+
+real(r8), allocatable, dimension(:) :: qc
+integer, dimension(max_obs_kinds) :: obs_kinds_inds = 0
+
+!-----------------------------------------------------------------------
+! General purpose variables
+!-----------------------------------------------------------------------
+
+integer :: ifile, nread, ngood
+integer :: i, io, ncunit
+
+type(time_type) :: obs_time
+
+character(len = 129) :: ncName, msgstring
+
+!=======================================================================
+! Get the party started
+!=======================================================================
+
+call initialize_utilities('obs_seq_coverage')
+call register_module(source,revision,revdate)
+call static_init_obs_sequence() ! Initialize the obs sequence module
+
+call init_obs(obs1, 0, 0)
+call init_obs(obs2, 0, 0)
+call init_obs_sequence(seq,0,0,0)
+
+! Allocate a hunk of stations. If we fill this up, we will
+! have to create temporary storage, copy, deallocate, reallocate ...
+
+num_stations = 0
+max_stations = 4000
+call initialize_stations(max_stations, stations)
+
+!----------------------------------------------------------------------
+! FIXME : at some point, I'd like to restrict this to just the ones
+! in the obs_seq file ... not the ones supported. Problem is
+! that you have to read all the obs sequence files and process
+! the region information before you know if you have any to
+! output. Maybe allocate more space in the netCDF header and
+! re-enter DEFINE mode after all is said and done ...
+!----------------------------------------------------------------------
+
+num_obs_kinds = max_obs_kinds
+
+do i = 1,max_obs_kinds
+ obs_kind_names(i) = get_obs_kind_name(i)
+enddo
+
+!----------------------------------------------------------------------
+! Read the namelist
+!----------------------------------------------------------------------
+
+call find_namelist_in_file('input.nml', 'obs_seq_coverage_nml', ncunit)
+read(ncunit, nml = obs_seq_coverage_nml, iostat = io)
+call check_namelist_read(ncunit, io, 'obs_seq_coverage_nml')
+
+! Record the namelist values used for the run ...
+if (do_nml_file()) write(nmlfileunit, nml=obs_seq_coverage_nml)
+if (do_nml_term()) write( * , nml=obs_seq_coverage_nml)
+
+if ((obs_sequence_name /= '') .and. (obs_sequence_list /= '')) then
+ write(msgstring,*)'specify "obs_sequence_name" or "obs_sequence_list"'
+ call error_handler(E_MSG, 'obs_seq_coverage', msgstring, source, revision, revdate)
+ write(msgstring,*)'set other to an empty string ... i.e. ""'
+ call error_handler(E_ERR, 'obs_seq_coverage', msgstring, source, revision, revdate)
+endif
+
+flavor_of_interest = get_obs_kind_index(obs_of_interest)
+obs_kinds_inds(flavor_of_interest) = 1
+
+if (flavor_of_interest < 0) then
+ write(msgstring,*)trim(obs_of_interest),' is not a known observation type.'
+ call error_handler(E_ERR, 'obs_seq_coverage', msgstring, source, revision, revdate)
+endif
+
+if (verbose) write(*,*)trim(adjustl(obs_of_interest)),' is type ',flavor_of_interest
+
+!====================================================================
+!====================================================================
+
+minl = set_location( (/ lonlim1, latlim1, 0.0_r8, 1.0_r8 /)) ! vertical unimportant
+maxl = set_location( (/ lonlim2, latlim2, 0.0_r8, 1.0_r8 /)) ! vertical unimportant
+
+!----------------------------------------------------------------------
+! Prepare the variables
+!----------------------------------------------------------------------
+
+allocate(obs_seq_filenames(1000))
+obs_seq_filenames = 'null'
+
+ObsFileLoop : do ifile=1, size(obs_seq_filenames)
+!-----------------------------------------------------------------------
+
+ ! Because of the ability to 'cycle' the ObsFileLoop, we need to
+ ! destroy and deallocate at the top of the loop.
+
+ call destroy_obs(obs1)
+ call destroy_obs(obs2)
+ call destroy_obs_sequence(seq)
+
+ if (allocated(qc)) deallocate(qc)
+ if (allocated(qc_copy_names)) deallocate(qc_copy_names)
+ if (allocated(obs_copy_names)) deallocate(obs_copy_names)
+
+ ! Determine the next input filename ...
+
+ if (obs_sequence_list == '') then
+ obs_seq_in_file_name = next_file(obs_sequence_name,ifile)
+ else
+ obs_seq_in_file_name = get_next_filename(obs_sequence_list,ifile)
+ if (obs_seq_in_file_name == '') exit ObsFileLoop
+ endif
+
+ 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_seq_coverage',msgstring,source,revision,revdate)
+ else
+ write(msgstring,*)trim(obs_seq_in_file_name),&
+ ' does not exist. Finishing up.'
+ call error_handler(E_MSG,'obs_seq_coverage',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(obs2, num_copies, num_qc)
+
+ ! I am taking the observational error variance and making it one of the copies
+
+ allNcopies = num_copies + Ncopies
+
+ if ((num_qc <= 0) .or. (num_copies <=0)) then
+ write(msgstring,*)'need at least 1 qc and 1 observation copy'
+ call error_handler(E_ERR,'obs_seq_coverage',msgstring,source,revision,revdate)
+ endif
+
+ allocate( obs_copy_names(allNcopies), qc_copy_names(num_qc), qc(num_qc))
+
+ 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(*,*)
+ endif
+
+ if (num_obs <= 1) then
+ last_ob_flag = .TRUE.
+ else
+ last_ob_flag = .FALSE.
+ endif
+
+ if ( debug ) write(*,*)'num obs/last_ob_flag are ',num_obs, last_ob_flag
+
+ !--------------------------------------------------------------------
+ ! Read the entire observation sequence - allocates 'seq' internally
+ !--------------------------------------------------------------------
+
+ call read_obs_seq(obs_seq_in_file_name, 0, 0, 0, seq)
+
+ do i=1, num_copies
+ msgstring = trim(get_copy_meta_data(seq,i))//' '
+ obs_copy_names(i) = msgstring(1:stringlength)
+ enddo
+ do i=1, Ncopies
+ obs_copy_names(num_copies+i) = trim(copy_names(i))
+ enddo
+ do i=1, num_qc
+ msgstring = trim(get_qc_meta_data(seq,i))//' '
+ qc_copy_names(i) = msgstring(1:stringlength)
+ enddo
+
+ if ( ifile == 1 ) then ! record the metadata for comparison
+
+ allocate(module_obs_copy_names(allNcopies), &
+ module_qc_copy_names(num_qc) )
+
+ do i=1, num_copies
+ msgstring = trim(get_copy_meta_data(seq,i))//' '
+ module_obs_copy_names(i) = msgstring(1:stringlength)
+ enddo
+ do i=1, Ncopies
+ module_obs_copy_names(num_copies+i) = trim(copy_names(i))
+ enddo
+ do i=1, num_qc
+ msgstring = trim(get_qc_meta_data(seq,i))//' '
+ module_qc_copy_names(i) = msgstring(1:stringlength)
+ enddo
+
+ else ! Compare all subsequent files' metadata to the first one
+
+ do i = 1,allNcopies
+ if (trim(obs_copy_names(i)) /= trim(module_obs_copy_names(i))) then
+ write(msgstring,'(''obs copy '',i3,'' from '',a)') i,trim(obs_seq_in_file_name)
+ call error_handler(E_MSG,'obs_seq_coverage',msgstring,source,revision,revdate)
+ msgstring = 'does not match the same observation copy from the first file.'
+ call error_handler(E_MSG,'obs_seq_coverage',msgstring,source,revision,revdate)
+ write(msgstring,'(''obs copy >'',a,''<'')') trim(obs_copy_names(i))
+ call error_handler(E_MSG,'obs_seq_coverage',msgstring,source,revision,revdate)
+ write(msgstring,'(''expected >'',a,''<'')') trim(module_obs_copy_names(i))
+ call error_handler(E_ERR,'obs_seq_coverage',msgstring,source,revision,revdate)
+ endif
+ enddo
+
+ do i = 1,num_qc
+ if (trim(qc_copy_names(i)) /= trim(module_qc_copy_names(i))) then
+ write(msgstring,'(''qc copy '',i3,'' from '',a)') i,trim(obs_seq_in_file_name)
+ call error_handler(E_MSG,'obs_seq_coverage',msgstring,source,revision,revdate)
+ msgstring = 'does not match the same qc copy from the first file.'
+ call error_handler(E_MSG,'obs_seq_coverage',msgstring,source,revision,revdate)
+ write(msgstring,'(''qc copy '',a)') trim(qc_copy_names(i))
+ call error_handler(E_MSG,'obs_seq_coverage',msgstring,source,revision,revdate)
+ write(msgstring,'(''expected '',a)') trim(module_qc_copy_names(i))
+ call error_handler(E_ERR,'obs_seq_coverage',msgstring,source,revision,revdate)
+ endif
+ enddo
+
+ endif
+
+ !--------------------------------------------------------------------
+ ! Read the first observation in the sequence
+ !--------------------------------------------------------------------
+
+ ngood = 0
+
+ if ( .not. get_first_obs(seq, obs1) ) &
+ call error_handler(E_ERR,'obs_seq_coverage', &
+ 'No first observation in sequence.', &
+ source,revision,revdate)
+
+ !--------------------------------------------------------------------
+ ObservationLoop : do nread = 1,num_obs
+ !--------------------------------------------------------------------
+
+ if ( verbose .and. (mod(nread,1000) == 0) ) &
+ write(*,*)'Processing obs ',nread,' of ',num_obs
+
+ call get_obs_def(obs1, obs_def)
+ call get_qc( obs1, qc)
+
+ flavor = get_obs_kind(obs_def)
+ obs_time = get_obs_def_time(obs_def)
+ obs_loc = get_obs_def_location(obs_def)
+
+ if (verbose .and. (nread == 1)) call print_time(obs_time,'First observation time')
+
+ !-----------------------------------------------------------------
+ ! determine if obs is a new location or time at an existing loc
+ !-----------------------------------------------------------------
+
+ if ( is_location_in_region(obs_loc,minl,maxl) .and. &
+ (flavor == flavor_of_interest) ) then
+
+ ngood = ngood + 1
+ station_id = find_station_location(flavor, obs_loc, stations)
+
+ if ( station_id < 1 ) then
+ station_id = add_new_station(flavor, obs_loc, stations)
+ else
+ if (verbose) write(*,*)'obs(',nread,') matches station ',station_id
+ endif
+
+ if ( is_time_new( obs_time, station_id, stations, timeindex) ) &
+ call update_time( obs_time, station_id, stations, timeindex)
+
+ else
+ locarray = get_location(obs_loc)
+ if (debug) write(*,*)'obs(',nread,') type ',flavor,'is not wanted',locarray
+ endif
+
+ call get_next_obs(seq, obs1, obs2, last_ob_flag)
+ if (.not. last_ob_flag) obs1 = obs2
+
+ if (debug) write(*,*)'obs(',nread,') last_ob_flag ',last_ob_flag
+
+ !--------------------------------------------------------------------
+ enddo ObservationLoop
+ !--------------------------------------------------------------------
+
+ if (verbose) write(*,*)'End of Loop for ',trim(obs_seq_in_file_name)
+
+enddo ObsFileLoop
+
+if (1 == 1) call print_summary
+
+write(ncName,'(''obs_coverage.nc'')')
+
+ncunit = InitNetCDF(ncName)
+call WriteNetCDF(ncunit, ncName, stations)
+call CloseNetCDF(ncunit, ncName)
+
+!-----------------------------------------------------------------------
+! Really, really, done.
+!-----------------------------------------------------------------------
+
+call destroy_obs(obs1)
+call destroy_obs(obs2)
+call destroy_obs_sequence(seq)
+call destroy_stations(stations)
+
+if (allocated(qc)) deallocate(qc)
+if (allocated(qc_copy_names)) deallocate(qc_copy_names)
+if (allocated(obs_copy_names)) deallocate(obs_copy_names)
+if (allocated(module_obs_copy_names)) deallocate(module_obs_copy_names)
+if (allocated(module_qc_copy_names )) deallocate(module_qc_copy_names )
+if (allocated(obs_seq_filenames)) deallocate(obs_seq_filenames)
+
+call timestamp(source,revision,revdate,'end') ! That closes the log file, too.
+
+!======================================================================
+CONTAINS
+!======================================================================
+
+function find_station_location(ObsType, ObsLocation, stationlist) result(station_id)
+! Simply try to find a matching lat/lon for an observation type
+integer, intent(in) :: ObsType
+type(location_type), intent(in) :: ObsLocation
+type(station), dimension(:), intent(inout) :: stationlist
+integer :: station_id
+
+integer :: i
+real(r8), dimension(3) :: obslocarray, stnlocarray
+real(r8) :: londiff, latdiff
+
+station_id = 0;
+
+FindLoop : do i = 1,num_stations
+
+ obslocarray = get_location(ObsLocation)
+ stnlocarray = get_location(stationlist(i)%location)
+
+ londiff = abs(obslocarray(1) - stnlocarray(1))
+ latdiff = abs(obslocarray(2) - stnlocarray(2))
+
+ if ( (londiff <= epsilon(londiff)) .and. &
+ (latdiff <= epsilon(latdiff)) .and. &
+ (ObsType == stationlist(i)%obs_type) ) then
+
+ ! if ( (ObsLocation == stationlist(i)%location) .and. &
+ ! (ObsType == stationlist(i)%obs_type) ) then
+
+ station_id = i
+ exit FindLoop
+ endif
+
+enddo FindLoop
+
+end function find_station_location
+
+
+
+function add_new_station(ObsType, ObsLocation, stationlist) result(station_id)
+
+! Ugh ... if a new location is found, add it. If the stationlist does not have
+! enough space, must copy the info to a temporary list, deallocate/reallocate
+! copy the info back, and deallocate the temporary list. Ugh.
+
+integer, intent(in) :: ObsType
+type(location_type), intent(in) :: ObsLocation
+type(station), allocatable, dimension(:), intent(inout) :: stationlist
+integer :: station_id
+
+integer :: i
+
+type(station), allocatable, dimension(:) :: templist
+
+if ( num_stations >= max_stations ) then ! need to make room
+
+ if (verbose) write(*,*)'Doubling number of possible stations from ', &
+ & num_stations,' to ',2*max_stations
+
+ ! Allocate temporary space; Copy.
+ ! Deallocate/nullify existing space.
+ ! Double the size of the existing space.
+ ! Copy the information back.
+ ! Deallocate/nullify the temporary space.
+ ! Actually add the new station information
+
+ ! Allocate the temporary space.
+ ! We'll worry about the number of time steps later.
+ call initialize_stations(num_stations, templist)
+
+ ! Copy the information to the temporary space.
+ DupLoop1 : do i = 1,num_stations
+
+ templist(i)%obs_type = stationlist(i)%obs_type
+ templist(i)%location = stationlist(i)%location
+ templist(i)%first_time = stationlist(i)%first_time
+ templist(i)%last_time = stationlist(i)%last_time
+ templist(i)%ntimes = stationlist(i)%ntimes
+
+ ! Make sure the time array is the right size, then copy.
+ if (associated(templist(i)%times)) then
+ deallocate( templist(i)%times )
+ nullify( templist(i)%times )
+ endif
+ allocate( templist(i)%times( size(stationlist(i)%times) ) )
+ templist(i)%times = stationlist(i)%times
+
+ enddo DupLoop1
+
+ ! Deallocate the stations, double the array length,
+ ! allocate the new space.
+ call destroy_stations(stationlist)
+ max_stations = 2 * max_stations
+ call initialize_stations(max_stations, stationlist)
+
+ ! Copy the information BACK to the new space.
+ DupLoop2 : do i = 1,num_stations
+
+ stationlist(i)%obs_type = templist(i)%obs_type
+ stationlist(i)%location = templist(i)%location
+ stationlist(i)%first_time = templist(i)%first_time
+ stationlist(i)%last_time = templist(i)%last_time
+ stationlist(i)%ntimes = templist(i)%ntimes
+
+ if (associated(stationlist(i)%times)) then
+ deallocate( stationlist(i)%times )
+ nullify( stationlist(i)%times )
+ endif
+ allocate( stationlist(i)%times( size(templist(i)%times) ) )
+ stationlist(i)%times = templist(i)%times
+
+ enddo DupLoop2
+
+ ! Remove the temporary space.
+ call destroy_stations(templist)
+
+endif
+
+! Add the new station information
+
+num_stations = num_stations + 1
+station_id = num_stations
+
+if (debug) write(*,*)'Adding station ',station_id,' for type ',ObsType
+
+stationlist(station_id)%obs_type = ObsType
+stationlist(station_id)%location = ObsLocation
+
+end function add_new_station
+
+
+
+function is_time_new(ObsTime, stationid, stationlist, timeindex)
+
+! Determine if the observation time is not already in the registry
+! of the times for the particular station.
+
+type(time_type), intent(in) :: ObsTime
+integer, intent(in) :: stationid
+type(station), dimension(:), intent(in) :: stationlist
+integer, intent(out) :: timeindex
+logical :: is_time_new
+
+type(time_type) :: stnhour, obhour, stndelta, obdelta
+integer :: i
+logical :: have_this_hour
+
+have_this_hour = .FALSE.
+timeindex = 0
+
+if ( stationlist(stationid)%ntimes == 0 ) then
+ is_time_new = .TRUE.
+ timeindex = 1
+ return
+endif
+
+TimeLoop : do i = 1,stationlist(stationid)%ntimes
+
+ stnhour = nearest_hour(stationlist(stationid)%times(i))
+ obhour = nearest_hour(ObsTime)
+
+ ! Make sure we only compare observations to the same hour
+ if (stnhour /= obhour) cycle TimeLoop
+
+ have_this_hour = .TRUE.
+
+ ! the time_minus function always returns a positive difference
+ stndelta = stationlist(stationid)%times(i) - stnhour
+ obdelta = ObsTime - obhour
+
+ ! Check to see if the increment is smaller
+ ! if it is, then we know which one to overwrite
+ if (obdelta < stndelta) then
+ if (verbose) call print_time(stationlist(stationid)%times(i),'replacing ')
+ if (verbose) call print_time(ObsTime,'with this observation time')
+ timeindex = i
+ exit TimeLoop
+ endif
+
+enddo TimeLoop
+
+if (have_this_hour .and. (timeindex == 0) ) then
+ ! the time we already have is closer than the candidate
+ is_time_new = .FALSE.
+elseif ( have_this_hour ) then ! candidate is closer
+ is_time_new = .TRUE.
+else ! must be a new observation hour
+ is_time_new = .TRUE.
+ timeindex = stationlist(stationid)%ntimes + 1
+endif
+
+end function is_time_new
+
+
+
+function nearest_hour(sometime)
+! Return the hour nearest to the input time
+type(time_type), intent(in) :: sometime
+type(time_type) :: nearest_hour
+
+type(time_type) :: thirty, tplus30
+integer :: days, hours, secs
+
+thirty = set_time(60*30-1,0) ! almost thirty minutes
+tplus30 = sometime + thirty
+call get_time(tplus30, secs, days)
+hours = secs/(60*60)
+
+nearest_hour = set_time(hours*60*60, days)
+
+if (debug) call print_time(sometime, 'input time')
+if (debug) call print_time(nearest_hour,'top of the hour')
+
+end function nearest_hour
+
+
+
+subroutine update_time(ObsTime, stationid, stationlist, timeindex)
+
+! Add a new time to the station registry.
+! If there is no additional space, must take action.
+
+type(time_type), intent(in) :: ObsTime
+integer, intent(in) :: stationid
+type(station), dimension(:), intent(inout) :: stationlist
+integer, intent(in) :: timeindex
+
+type(time_type), allocatable, dimension(:) :: temptimes
+integer :: ntimes
+
+! Update stuff that seems like a good idea,
+! but I don't really know if I'll use it ...
+if ( stationlist(stationid)%ntimes == 0 ) then
+ stationlist(stationid)%first_time = ObsTime
+ stationlist(stationid)%last_time = ObsTime
+endif
+
+if ( stationlist(stationid)%first_time > ObsTime ) &
+ stationlist(stationid)%first_time = ObsTime
+if ( stationlist(stationid)%last_time < ObsTime ) &
+ stationlist(stationid)%last_time = ObsTime
+
+! Do we need to make room for the new time
+
+ntimes = size(stationlist(stationid)%times)
+
+if ( (stationlist(stationid)%ntimes >= ntimes) .and. &
+ (timeindex > stationlist(stationid)%ntimes) ) then
+
+ allocate(temptimes(ntimes))
+ temptimes = stationlist(stationid)%times
+
+ if (associated(stationlist(stationid)%times)) then
+ deallocate( stationlist(stationid)%times )
+ nullify( stationlist(stationid)%times )
+ endif
+
+ allocate( stationlist(stationid)%times(2*ntimes) )
+ stationlist(stationid)%times(1:ntimes) = temptimes(1:ntimes)
+
+ deallocate(temptimes)
+
+endif
+
+! If the time is new one, increment counter
+if (timeindex > stationlist(stationid)%ntimes) then
+ stationlist(stationid)%ntimes = stationlist(stationid)%ntimes + 1
+ if (debug) write(*,*)'Adding a new time to ',stationid, &
+ ' count ', stationlist(stationid)%ntimes, &
+ ' of ', size(stationlist(stationid)%times)
+endif
+
+! Stuff the time in the appropriate slot ... finally.
+stationlist(stationid)%times(timeindex) = ObsTime
+
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list