[Dart-dev] [6567] DART/trunk/obs_sequence: major update to this tool.
nancy at ucar.edu
nancy at ucar.edu
Mon Nov 4 15:48:15 MST 2013
Revision: 6567
Author: nancy
Date: 2013-11-04 15:48:15 -0700 (Mon, 04 Nov 2013)
Log Message:
-----------
major update to this tool. it enforces the restriction that it
should be used on completely compatible sets of obs_seq.final files,
and it can now compare more than 4 files at a single time.
it will take a list of groups of files and compare them N at a time.
previous versions of this tool supported only a 2-way and then
a 4-way compare.
Modified Paths:
--------------
DART/trunk/obs_sequence/obs_common_subset.f90
DART/trunk/obs_sequence/obs_common_subset.html
DART/trunk/obs_sequence/obs_common_subset.nml
-------------- next part --------------
Modified: DART/trunk/obs_sequence/obs_common_subset.f90
===================================================================
--- DART/trunk/obs_sequence/obs_common_subset.f90 2013-11-04 22:45:21 UTC (rev 6566)
+++ DART/trunk/obs_sequence/obs_common_subset.f90 2013-11-04 22:48:15 UTC (rev 6567)
@@ -6,25 +6,42 @@
program obs_common_subset
-! this program expects to be given pairs, or lists of pairs, of obs_seq
-! files as input and does an obs-by-obs comparison. any obs which does not
-! match in type, time, location, or (if present) DART QC is discarded.
-! in addition to matching, the DART QCs have to be below a given threshold.
-! the intent is to run two experiments with identical input observations
-! but some difference in namelist settings or other filter differences.
-! then take the two output obs_seq.final files and remove any obs which
-! were not assimilated in both experiments for any reason. then comparing
-! the output obs_seq files will be comparing exactly the same obs, and only
-! those which were assimilated. differences in the output diagnostics will
-! be because of the experiment differences, not differences in the number of
-! obs assimilated.
+! This program takes as input one or more sets of N obs_seq filenames.
+! N is usually 2, but can be any number of files to compare together.
+! It opens the obs_seq files N at a time and does an obs-by-obs comparison.
+! It is intended to pull out a common subset of obs which were successfully
+! assimilated in all experiments so the subsequent diagnostics can give
+! an apples-to-apples comparison of the relative error or fit to obs
+! between different experiments.
+!
+! This tool does not search in the input obs_seq files to try to match obs
+! in different orders or to skip over obs present in only one input file,
+! so all experiments should be done with identical input obs_seq files.
+! The obs to be assimilated or evaluated can be selected by namelist
+! at run time if they are different, and different experiment settings
+! can be selected by other namelist values.
+!
+! If the next input obs don't match exactly (identical types, times, locations,
+! and DART QC) they are discarded. If they match, and if the DART QC indicates
+! they have been successfully assimilated in all experiments, then each of the
+! N input obs are copied through to the N output obs_seq files. At the end of
+! running this tool it creates N output obs_seq files which contain a consistent
+! set of assimilated obs across the set.
-! running this program creates 2 new output files for each pair of input files;
-! the names of the input obs_seq files with a suffix appended. it can take
-! list of obs_seq files, either explicitly listed or in a separate file with
-! one input filename per line.
+! The outputs are suitable for input to the obs_diag or obs_seq_to_netcdf tools
+! and there are matlab scripts for comparing the results of N experiments.
+! Differences in the output diagnostics will now be because of the experiment
+! differences, not differences in the number of obs assimilated.
+! Running this program creates sets of N new output files corresponding to each
+! of N input files. The output names are the input filenames with a suffix appended.
+! It can take a list of obs_seq files, either given directly in the namelist or
+! from a separate ascii file with one input filename per line. If more than N files
+! are specified it will process them N at a time and loop until the end of the list.
+! The filenames should be listed so all files from experiment 1 are first, then all
+! files from experiment 2, up to N experiments.
+
use types_mod, only : r8, missing_r8, metadatalength, obstypelength
use utilities_mod, only : register_module, initialize_utilities, &
find_namelist_in_file, check_namelist_read, &
@@ -57,75 +74,74 @@
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate = "$Date$"
-type(obs_sequence_type) :: seq_in1, seq_in2, seq_out1, seq_out2
-type(obs_sequence_type) :: seq_in3, seq_in4, seq_out3, seq_out4
-type(obs_type) :: obs_in1, next_obs_in1, obs_in2, next_obs_in2
-type(obs_type) :: obs_in3, next_obs_in3, obs_in4, next_obs_in4
-type(obs_type) :: obs_out1, prev_obs_out1, obs_out2, prev_obs_out2
-type(obs_type) :: obs_out3, prev_obs_out3, obs_out4, prev_obs_out4
-logical :: is_this_last1, is_this_last2
-logical :: is_this_last3, is_this_last4
-logical :: comparing3 = .true.
-logical :: comparing4 = .true.
-logical :: wanted = .false.
-integer :: size_seq_in1, num_copies_in1, num_qc_in1
-integer :: size_seq_in2, num_copies_in2, num_qc_in2
-integer :: size_seq_in3, num_copies_in3, num_qc_in3
-integer :: size_seq_in4, num_copies_in4, num_qc_in4
-integer :: num_copies_in, num_qc_in, size_seq_out
-integer :: num_inserted, iunit, io, i, j
-integer :: max_num_obs1, max_num_obs2, file_id
-integer :: max_num_obs3, max_num_obs4
-integer :: num_rejected_badqc, num_rejected_diffqc
-integer :: num_rejected_other
+! max number of files to compare against each other. this program will loop
+! over sets of files for as long as there are input files in the list. if you
+! need to compare more than 50 files together, let us know what you are doing
+! because i'd be very interested to hear. but you can make the maxseq larger
+! in that case.
+
+integer, parameter :: maxseq = 50 ! max number of files compared at once
+
+type(obs_sequence_type) :: seq_in(maxseq)
+type(obs_sequence_type) :: seq_out(maxseq)
+type(obs_type) :: obs_in(maxseq), next_obs_in(maxseq)
+type(obs_type) :: obs_out(maxseq), prev_obs_out(maxseq)
+logical :: is_this_last(maxseq)
+logical :: wanted
+integer :: size_seq_in(maxseq), num_copies_in(maxseq), num_qc_in(maxseq)
+integer :: size_seq_out, num_inserted, iunit, io, i, j, k
+integer :: max_num_obs(maxseq), file_id, atonce, nsets, other, offset
+integer :: num_rejected_badqc, num_rejected_diffqc, num_rejected_other
+integer :: num_mismatch_loc, num_mismatch_time, num_mismatch_type
character(len = 129) :: read_format
logical :: pre_I_format, cal
-character(len = 256) :: msgstring, msgstring1, msgstring2
-character(len = 164) :: filename_out1, filename_out2
-character(len = 164) :: filename_out3, filename_out4
+character(len = 256) :: msgstring, msgstring1, msgstring2, msgstring3
+character(len = 256) :: filename_in(maxseq)
+character(len = 300) :: filename_out(maxseq) ! filename_in + . + suffix
-character(len = metadatalength) :: dart_qc_meta_data = 'DART quality control'
+character(len = metadatalength), parameter :: dart_qc_meta_data = 'DART quality control'
character(len = metadatalength) :: meta_data
+
integer :: qc_index
-integer :: num_input_files = 0
+integer :: num_input_sets = 0
-! could go into namelist if you wanted more control
-integer, parameter :: print_every = 200
-integer, parameter :: qc_threshold = 1 ! ok if <= this
-
!----------------------------------------------------------------
! Namelist input with default values
! max_num_input_files : maximum number of input sequence files to be processed
-! lazy, pick big numbers. make them bigger if too small.
-integer, parameter :: max_num_input_files = 1000
-integer, parameter :: max_obs_input_types = 500
+! if filenames are specified in the namelist. if filenames are specified in a
+! separate ascii input file, the list length can be infinite. for names in the
+! namelist, the len must be > num_to_compare_at_once * number of sets of files
+! to do. e.g. if you are comparing 3 files at a time, and want to compare a
+! series of 100 sets of 3 files, that would be 300 filenames. make these
+! numbers bigger if these values are too small, and recompile.
-character(len = 129) :: filename_seq1(max_num_input_files) = ''
-character(len = 129) :: filename_seq2(max_num_input_files) = ''
-character(len = 129) :: filename_seq3(max_num_input_files) = ''
-character(len = 129) :: filename_seq4(max_num_input_files) = ''
-character(len = 129) :: filename_seq_list1 = ''
-character(len = 129) :: filename_seq_list2 = ''
-character(len = 129) :: filename_seq_list3 = ''
-character(len = 129) :: filename_seq_list4 = ''
+integer, parameter :: max_num_input_files = 5000
+
+character(len = 256) :: filename_seq(max_num_input_files) = ''
+character(len = 256) :: filename_seq_list = ''
character(len = 32) :: filename_out_suffix = '.common'
-logical :: print_only = .false.
-character(len=32) :: calendar = 'Gregorian'
+integer :: num_to_compare_at_once = 2
+logical :: print_only = .false.
+! not expected to be changed often, but its here if needed.
+integer :: print_every = 1000
+integer :: dart_qc_threshold = 3 ! ok if DART QC <= this
+
+character(len=32) :: calendar = 'Gregorian'
+
namelist /obs_common_subset_nml/ &
- filename_seq1, filename_seq_list1, &
- filename_seq2, filename_seq_list2, &
- filename_seq3, filename_seq_list3, &
- filename_seq4, filename_seq_list4, &
- filename_out_suffix, print_only, calendar
+ num_to_compare_at_once, &
+ filename_seq, filename_seq_list, &
+ filename_out_suffix, print_only, &
+ print_every, dart_qc_threshold, calendar
!----------------------------------------------------------------
! Start of the program:
!
-! Process each input observation sequence file in turn, optionally
-! selecting observations to insert into the output sequence file.
+! Process each N input observation sequence file in turn,
+! selecting observations to insert into the output sequence files.
!----------------------------------------------------------------
call setup()
@@ -140,18 +156,36 @@
if (do_nml_term()) write( * , nml=obs_common_subset_nml)
! the logic here is slightly different than the obs_sequence_tool.
-! the user gives us 4 lists of obs_seq files; either in the namelist
+! the user gives us a list of obs_seq files; either in the namelist
! or as the name of a file which contains the list, one per line.
-! either way, the lists must be the same length.
+! the files are compared N at a time, until the list is exhausted.
! as in the other tools, it is an error to specify both an explicit
! list and the name of file for input.
-call handle_filenames(filename_seq1, filename_seq_list1, &
- filename_seq2, filename_seq_list2, &
- filename_seq3, filename_seq_list3, &
- filename_seq4, filename_seq_list4, &
- num_input_files)
+write(msgstring, '(A,I3,A)') "Comparing ", num_to_compare_at_once, " obs_seq files at a time"
+call error_handler(E_MSG, "obs_common_subset", msgstring, source,revision,revdate)
+if ((num_to_compare_at_once < 2) .or. (num_to_compare_at_once > maxseq)) then
+ write(msgstring, *) "num_to_compare_at_once must be >= 2 and <= ", maxseq
+ call error_handler(E_ERR, "obs_common_subset", msgstring, source,revision,revdate)
+endif
+
+! if there is only a single file from each experiment to compare, the order
+! is obvious. if each experiment has a list of input obs_seq.final files from
+! a multi-step assimilation, list all the files from experiment 1, then all
+! files from experiment 2, etc.
+call handle_filenames(num_to_compare_at_once, filename_seq, filename_seq_list, &
+ num_input_sets)
+
+if (num_input_sets > 1) then
+ write(msgstring, '(A,I3,A)') "Will loop ", num_input_sets, " times to process the entire input filelist"
+ call error_handler(E_MSG, "obs_common_subset", msgstring, source,revision,revdate)
+endif
+
+! because i am a lazy typist
+atonce = num_to_compare_at_once
+nsets = num_input_sets
+
! the default is a gregorian calendar. if you are using a different type
! set it in the namelist. this only controls how it prints out the first
! and last timestamps in the obs_seq files.
@@ -163,365 +197,258 @@
! end of namelist processing and setup
-! single pass algorithm (unlike other obs tools). process the files in
-! pairs, making sure the metadata matches. the number of obs does NOT
-! have to match. iterate them pairwise and keep any matching obs which
-! share the same QC values, as long as the QCs are <= threshold.
-! (QC tests are skipped if DART QC isn't in the metadata.)
+! process the files in sets, making sure the metadata matches. keep matching
+! obs which share the same QC values, as long as the QCs are <= threshold.
+! since this tool has very little utility without the DART QC values, fail
+! if they aren't present.
-! count of input files was set in the handle_filenames() routine above.
-do i = 1, num_input_files
+! count of input sets was set in the handle_filenames() routine above.
+NUMSETS: do j = 1, nsets
- if ((len(filename_seq1(i)) == 0) .or. (filename_seq1(i) == "")) then
- write(msgstring, *) 'entry ', i, 'for sequence1 is empty or null'
- call error_handler(E_ERR,'obs_common_subset', msgstring, &
- source,revision,revdate) ! shouldn't happen
- endif
- if ((len(filename_seq2(i)) == 0) .or. (filename_seq2(i) == "")) then
- write(msgstring, *) 'entry ', i, 'for sequence2 is empty or null'
- call error_handler(E_ERR,'obs_common_subset', msgstring, &
- source,revision,revdate) ! shouldn't happen
- endif
+ ! if num_to_compare_at_once = 3 and num_input_sets = 100, then this
+ ! assumes the input file list is organized in this order:
+ ! a1,b1,c1, a2,b2,c2, ... a100,b100,c100
+ ! where a, b, c are the 3 separate experiments to compare,
+ ! and 1, 2, ... 100 are the obs_seq.final output files from 100 steps
+ ! of a multi-step assimilation.
- if ((len(filename_seq3(i)) == 0) .or. (filename_seq3(i) == "")) then
- comparing3 = .false.
- write(msgstring, *) 'entry ', i, 'for sequence3 is empty or null'
- call error_handler(E_MSG,'obs_common_subset', msgstring, &
- source,revision,revdate) ! may happen
- endif
+ do i = 1, atonce
+
+ offset = (i-1)*num_input_sets + j
- if ((len(filename_seq4(i)) == 0) .or. (filename_seq4(i) == "")) then
- comparing4 = .false.
- write(msgstring, *) 'entry ', i, 'for sequence4 is empty or null'
- call error_handler(E_MSG,'obs_common_subset', msgstring, &
- source,revision,revdate) ! may happen
- endif
+ ! fill in the names of the current set of obs_seq files
+ filename_in(i) = filename_seq(offset)
- ! read in the next pair of files.
+ enddo
- call read_obs_seq_header(filename_seq1(i), num_copies_in1, num_qc_in1, &
- size_seq_in1, max_num_obs1, file_id, read_format, pre_I_format, &
- close_the_file = .true.)
+ ! now filename_in(:) has all the names for this set
+ do i = 1, atonce
+
+ if ((len(filename_in(i)) == 0) .or. (filename_in(i) == "")) then
+ write(msgstring, *) 'sequence name ', i, ' in set ', j, ' is empty or null'
+ call error_handler(E_ERR,'obs_common_subset', msgstring, &
+ source,revision,revdate) ! shouldn't happen because already tested for this
+ endif
- if (max_num_obs1 == 0) then
- write(msgstring,*) 'No obs in input sequence file ', trim(filename_seq1(i))
- call error_handler(E_MSG,'obs_common_subset',msgstring)
- cycle
- endif
+ ! read in the next set of files.
- call read_obs_seq_header(filename_seq2(i), num_copies_in2, num_qc_in2, &
- size_seq_in2, max_num_obs2, file_id, read_format, pre_I_format, &
- close_the_file = .true.)
-
- if (max_num_obs2 == 0) then
- write(msgstring,*) 'No obs in input sequence file ', trim(filename_seq2(i))
- call error_handler(E_MSG,'obs_common_subset',msgstring)
- cycle
- endif
-
- if (comparing3) then
- call read_obs_seq_header(filename_seq3(i), num_copies_in3, num_qc_in3, &
- size_seq_in3, max_num_obs3, file_id, read_format, pre_I_format, &
+ call read_obs_seq_header(filename_in(i), num_copies_in(i), num_qc_in(i), &
+ size_seq_in(i), max_num_obs(i), file_id, read_format, pre_I_format, &
close_the_file = .true.)
- if (max_num_obs3 == 0) then
- write(msgstring,*) 'No obs in input sequence file ', trim(filename_seq3(i))
- call error_handler(E_MSG,'obs_common_subset',msgstring)
- cycle
+ if (size_seq_in(i) == 0) then
+ write(msgstring1, *) 'This tool cannot process empty obs_seq files'
+ write(msgstring,*) 'Obs in input sequence file ', trim(filename_in(i))
+ call error_handler(E_ERR,'obs_common_subset',msgstring, &
+ source,revision,revdate, text2=msgstring1)
endif
- endif
- if (comparing4) then
- call read_obs_seq_header(filename_seq4(i), num_copies_in4, num_qc_in4, &
- size_seq_in4, max_num_obs4, file_id, read_format, pre_I_format, &
- close_the_file = .true.)
+ enddo
- if (max_num_obs4 == 0) then
- write(msgstring,*) 'No obs in input sequence file ', trim(filename_seq4(i))
- call error_handler(E_MSG,'obs_common_subset',msgstring)
- cycle
+ do i = 2, atonce
+ if (size_seq_in(1) /= size_seq_in(i)) then
+ write(msgstring,*) 'Input obs_seq files cannot have different observation counts'
+ write(msgstring2,*) 'Counts for file 1 and file ', i,' are ', size_seq_in(1), &
+ ' obs and ', size_seq_in(i), ' obs'
+ call error_handler(E_ERR,'obs_common_subset',msgstring, &
+ source,revision,revdate, text2=msgstring2)
endif
- endif
+ enddo
- write(*, *) 'Starting to process input sequence files: '
- write(*,*) trim(filename_seq1(i))
- write(*,*) trim(filename_seq2(i))
- if (comparing3) write(*,*) trim(filename_seq3(i))
- if (comparing4) write(*,*) trim(filename_seq4(i))
-
+ write( *, *) 'Starting to process input sequence files: '
write(logfileunit, *) 'Starting to process input sequence files: '
- write(logfileunit,*) trim(filename_seq1(i))
- write(logfileunit,*) trim(filename_seq2(i))
- if (comparing3) write(logfileunit,*) trim(filename_seq3(i))
- if (comparing4) write(logfileunit,*) trim(filename_seq4(i))
- call read_obs_seq(filename_seq1(i), 0, 0, 0, seq_in1)
- call read_obs_seq(filename_seq2(i), 0, 0, 0, seq_in2)
- if (comparing3) call read_obs_seq(filename_seq3(i), 0, 0, 0, seq_in3)
- if (comparing4) call read_obs_seq(filename_seq4(i), 0, 0, 0, seq_in4)
+ do i = 1, atonce
+ write( * , *) trim(filename_in(i))
+ write(logfileunit, *) trim(filename_in(i))
+ call read_obs_seq(filename_in(i), 0, 0, 0, seq_in(i))
+ enddo
+
! make sure the files have compatible metadata. this errors out here if not.
- call compare_metadata(seq_in1, seq_in2, filename_seq1(i), filename_seq2(i))
+ call compare_metadata(atonce, seq_in, filename_in)
- ! if we get here, the number of copies and qcs must match. use a single
- ! variable for both.
- num_copies_in = num_copies_in1
- num_qc_in = num_qc_in1
-
! sanity check - ensure the linked list times are in increasing time order
- call validate_obs_seq_time(seq_in1, filename_seq1(i))
- call validate_obs_seq_time(seq_in2, filename_seq2(i))
- if (comparing3) call validate_obs_seq_time(seq_in3, filename_seq3(i))
- if (comparing4) call validate_obs_seq_time(seq_in4, filename_seq4(i))
+ do i = 1, atonce
+ call validate_obs_seq_time(seq_in(i), filename_in(i))
+ enddo
- ! the output can have no more than the minimum number of input obs.
- size_seq_out = min(size_seq_in1, size_seq_in2)
- if (comparing3) size_seq_out = min(size_seq_out, size_seq_in3)
- if (comparing4) size_seq_out = min(size_seq_out, size_seq_in4)
+ ! the output can have no more than the number of input obs, and we already
+ ! tested that all sequences have the same obs count.
+ size_seq_out = size_seq_in(1)
- ! find which DART qc copy has the DART quality control. set to -1 if not
- ! present, which will skip the threshold test. and at this point we know
- ! the two sequences are the same so just pick one and query it.
+ ! find which DART QC copy has the DART quality control. this version of
+ ! the tool fails if one is not found. there seems very little utility in
+ ! this tool if it tries to process files without a DART QC.
qc_index = -1
- do j=1, get_num_qc(seq_in1)
- if(index(get_qc_meta_data(seq_in1,j), dart_qc_meta_data) > 0) then
- qc_index = j
- exit
+ QCLOOP: do i=1, get_num_qc(seq_in(1))
+ if(index(get_qc_meta_data(seq_in(1),i), dart_qc_meta_data) > 0) then
+ qc_index = i
+ exit QCLOOP
endif
- enddo
+ enddo QCLOOP
+ ! if we do decide there is utility in running these on files without a
+ ! DART QC, remove this error exit here. nsc 24 oct 2013
+ if (qc_index < 0) then
+ write(msgstring,*) 'DART QC was not found in the input obs_seq files.'
+ write(msgstring2,*) 'This tool only works on output files from an assimilation.'
+ call error_handler(E_ERR,'obs_common_subset',msgstring, &
+ source,revision,revdate, text2=msgstring2)
+ endif
+
+
! blank line, start of actually creating output file
call error_handler(E_MSG,' ',' ')
- ! Initialize individual observation variables
- call init_obs( obs_in1, num_copies_in, num_qc_in)
- call init_obs( next_obs_in1, num_copies_in, num_qc_in)
- call init_obs( obs_out1, num_copies_in, num_qc_in)
- call init_obs(prev_obs_out1, num_copies_in, num_qc_in)
+ do i = 1, atonce
+ ! Initialize individual observation variables
+ call init_obs( obs_in(i), num_copies_in(1), num_qc_in(1))
+ call init_obs( next_obs_in(i), num_copies_in(1), num_qc_in(1))
+ call init_obs( obs_out(i), num_copies_in(1), num_qc_in(1))
+ call init_obs(prev_obs_out(i), num_copies_in(1), num_qc_in(1))
- call init_obs( obs_in2, num_copies_in, num_qc_in)
- call init_obs( next_obs_in2, num_copies_in, num_qc_in)
- call init_obs( obs_out2, num_copies_in, num_qc_in)
- call init_obs(prev_obs_out2, num_copies_in, num_qc_in)
+ ! create the output sequences
+ call init_obs_sequence(seq_out(i), num_copies_in(1), num_qc_in(1), size_seq_out)
- if (comparing3) then
- call init_obs( obs_in3, num_copies_in, num_qc_in)
- call init_obs( next_obs_in3, num_copies_in, num_qc_in)
- call init_obs( obs_out3, num_copies_in, num_qc_in)
- call init_obs(prev_obs_out3, num_copies_in, num_qc_in)
- endif
+ do k=1, num_copies_in(1)
+ meta_data = get_copy_meta_data(seq_in(1), k)
+ call set_copy_meta_data(seq_out(i), k, meta_data)
+ enddo
+ do k=1, num_qc_in(1)
+ meta_data = get_qc_meta_data(seq_in(1), k)
+ call set_qc_meta_data(seq_out(i), k, meta_data)
+ enddo
- if (comparing4) then
- call init_obs( obs_in4, num_copies_in, num_qc_in)
- call init_obs( next_obs_in4, num_copies_in, num_qc_in)
- call init_obs( obs_out4, num_copies_in, num_qc_in)
- call init_obs(prev_obs_out4, num_copies_in, num_qc_in)
- endif
-
- ! create the output sequences here
- call init_obs_sequence(seq_out1, num_copies_in, num_qc_in, size_seq_out)
- call init_obs_sequence(seq_out2, num_copies_in, num_qc_in, size_seq_out)
- if (comparing3) call init_obs_sequence(seq_out3, num_copies_in, num_qc_in, size_seq_out)
- if (comparing4) call init_obs_sequence(seq_out4, num_copies_in, num_qc_in, size_seq_out)
-
- do j=1, num_copies_in
- meta_data = get_copy_meta_data(seq_in1, j)
- call set_copy_meta_data(seq_out1, j, meta_data)
- call set_copy_meta_data(seq_out2, j, meta_data)
- if (comparing3) call set_copy_meta_data(seq_out3, j, meta_data)
- if (comparing4) call set_copy_meta_data(seq_out4, j, meta_data)
enddo
- do j=1, num_qc_in
- meta_data = get_qc_meta_data(seq_in1, j)
- call set_qc_meta_data(seq_out1, j, meta_data)
- call set_qc_meta_data(seq_out2, j, meta_data)
- if (comparing3) call set_qc_meta_data(seq_out3, j, meta_data)
- if (comparing4) call set_qc_meta_data(seq_out4, j, meta_data)
- enddo
- ! not sure what this should do - print after it has generated the common set?
- if (print_only) call print_obs_seq(seq_in1, filename_seq1(i))
- if (print_only) call print_obs_seq(seq_in2, filename_seq2(i))
- if (print_only .and. comparing3) call print_obs_seq(seq_in3, filename_seq3(i))
- if (print_only .and. comparing4) call print_obs_seq(seq_in4, filename_seq4(i))
-
- !-------------------------------------------------------------
! Start to insert obs from sequence_in into sequence_out
- !
! NOTE: insert_obs_in_seq CHANGES the obs passed in.
! Must pass a copy of incoming obs to insert_obs_in_seq.
- !--------------------------------------------------------------
+
num_inserted = 0
num_rejected_badqc = 0
num_rejected_diffqc = 0
num_rejected_other = 0
+ num_mismatch_time = 0
+ num_mismatch_type = 0
+ num_mismatch_loc = 0
- if ( get_first_obs(seq_in1, obs_in1) .and. get_first_obs(seq_in2, obs_in2) ) then
-
- if (comparing3) then
- if (.not. get_first_obs(seq_in3, obs_in3)) then
- write(msgstring, *)'no first observation in ',trim(filename_seq3(i))
- call error_handler(E_ERR,'obs_common_subset',msgstring,source,revision,revdate)
- endif
+ ! get the first obs from each file. since we have already checked
+ ! for empty sequence files, this is not expected to fail at this point.
+ do i = 1, atonce
+ if ( .not. get_first_obs(seq_in(i), obs_in(i))) then
+ write(msgstring, *)'no first observation in ',trim(filename_in(i))
+ call error_handler(E_ERR,'obs_common_subset',msgstring,source,revision,revdate)
endif
+ is_this_last(i) = .false.
+ next_obs_in(i) = obs_in(i)
+ enddo
- if (comparing4) then
- if (.not. get_first_obs(seq_in4, obs_in4)) then
- write(msgstring, *)'no first observation in ',trim(filename_seq4(i))
- call error_handler(E_ERR,'obs_common_subset',msgstring,source,revision,revdate)
- endif
- endif
- is_this_last1 = .false.
- is_this_last2 = .false.
- is_this_last3 = .false.
- is_this_last4 = .false.
- next_obs_in1 = obs_in1
- next_obs_in2 = obs_in2
- if (comparing3) next_obs_in3 = obs_in3
- if (comparing4) next_obs_in4 = obs_in4
+ ! since we've checked to be sure that all the input obs_seq files have
+ ! the same obs counts, i can test the result from the first seq only.
+ ObsLoop : do while ( .not. is_this_last(1) )
- ObsLoop : do while ( .not. is_this_last1 .and. .not. is_this_last2)
+ ! move the next obs into position to be worked on
+ do i = 1, atonce
+ obs_in(i) = next_obs_in(i)
+ enddo
- ! essentially record position in seq_out
- obs_in1 = next_obs_in1
- obs_in2 = next_obs_in2
- if (comparing3) obs_in3 = next_obs_in3
- if (comparing4) obs_in4 = next_obs_in4
+ ! see if we want to keep it or not
+ wanted = all_good(obs_in, atonce, qc_index, dart_qc_threshold)
- if (comparing4) then
- wanted = all_good(obs_in1, obs_in2, qc_index, qc_threshold, obs_in3, obs_in4)
- elseif (comparing3) then
- wanted = all_good(obs_in1, obs_in2, qc_index, qc_threshold, obs_in3)
+ if ( wanted ) then
+ do i = 1, atonce
+ obs_out(i) = obs_in(i)
+ enddo
+
+ ! Since the stride through the observation sequence file is always
+ ! guaranteed to be in temporally-ascending order, we can use the
+ ! 'previous' observation as the starting point to search for the
+ ! correct insertion point. This speeds up the insert code a lot.
+
+ if (num_inserted == 0) then
+ do i = 1, atonce
+ call insert_obs_in_seq(seq_out(i), obs_out(i))
+ enddo
else
- wanted = all_good(obs_in1, obs_in2, qc_index, qc_threshold)
+ do i = 1, atonce
+ call insert_obs_in_seq(seq_out(i), obs_out(i), prev_obs_out(i))
+ enddo
endif
- if ( wanted ) then
- obs_out1 = obs_in1
- obs_out2 = obs_in2
- if (comparing3) obs_out3 = obs_in3
- if (comparing4) obs_out4 = obs_in4
+ ! update position in seq for next insert
+ do i = 1, atonce
+ prev_obs_out(i) = obs_out(i)
+ enddo
- ! Since the stride through the observation sequence file is always
- ! guaranteed to be in temporally-ascending order, we can use the
- ! 'previous' observation as the starting point to search for the
- ! correct insertion point. This speeds up the insert code a lot.
+ num_inserted = num_inserted + 1
- if (num_inserted > 0) then
- call insert_obs_in_seq(seq_out1, obs_out1, prev_obs_out1)
- call insert_obs_in_seq(seq_out2, obs_out2, prev_obs_out2)
- if (comparing3) call insert_obs_in_seq(seq_out3, obs_out3, prev_obs_out3)
- if (comparing4) call insert_obs_in_seq(seq_out4, obs_out4, prev_obs_out4)
- else
- call insert_obs_in_seq(seq_out1, obs_out1)
- call insert_obs_in_seq(seq_out2, obs_out2)
- if (comparing3) call insert_obs_in_seq(seq_out3, obs_out3)
- if (comparing4) call insert_obs_in_seq(seq_out4, obs_out4)
+ if (print_every > 0) then
+ if (mod(num_inserted,print_every) == 0) then
+ print*, 'inserted obs number ',num_inserted,' of ',size_seq_out
endif
+ endif
- ! update position in seq for next insert
- prev_obs_out1 = obs_out1
- prev_obs_out2 = obs_out2
- if (comparing3) prev_obs_out3 = obs_out3
- if (comparing4) prev_obs_out4 = obs_out4
- num_inserted = num_inserted + 1
+ endif
- if (print_every > 0) then
- if (mod(num_inserted,print_every) == 0) then
- print*, 'inserted number ',num_inserted,' of ',size_seq_out
- endif
- endif
+ do i = 1, atonce
+ call get_next_obs(seq_in(i), obs_in(i), next_obs_in(i), is_this_last(i))
+ enddo
- endif
+ enddo ObsLoop
- call get_next_obs(seq_in1, obs_in1, next_obs_in1, is_this_last1)
- call get_next_obs(seq_in2, obs_in2, next_obs_in2, is_this_last2)
- if (comparing3) call get_next_obs(seq_in3, obs_in3, next_obs_in3, is_this_last3)
- if (comparing4) call get_next_obs(seq_in4, obs_in4, next_obs_in4, is_this_last4)
- enddo ObsLoop
+ ! finished creating output sequences. print info and write out files
+ write(msgstring, *) 'Number of obs in each output seq file :', get_num_key_range(seq_out(1))
+ call error_handler(E_MSG,'obs_common_subset',msgstring)
+
+ do i = 1, atonce
+ filename_out(i) = trim(filename_in(i))//trim(filename_out_suffix)
+ call print_obs_seq(seq_out(i), filename_out(i))
+ enddo
+
+ if (print_only) then
+ write(msgstring,*) 'Output sequence files not created; print_only in namelist is .true.'
+ call error_handler(E_MSG,'', msgstring)
else
- write(msgstring, *)'no first observation in ',trim(filename_seq1(i))
- write(msgstring1,*)' or in ',trim(filename_seq2(i))
- call error_handler(E_MSG,'obs_common_subset', msgstring, text2=msgstring1)
- endif
- if (.not. print_only) then
- print*, '--------- Obs seq file pair # : ', i
- print*, 'Number of obs in sequence 1 : ', size_seq_in1
- print*, 'Number of obs in sequence 2 : ', size_seq_in2
+ print*, '--------- Obs seq file set # : ', j
+ print*, 'Number of obs in each input : ', size_seq_in(1)
print*, 'Number of obs rejected diff qc : ', num_rejected_diffqc
print*, 'Number of obs rejected bad qc : ', num_rejected_badqc
print*, 'Number of obs rejected other : ', num_rejected_other
+ if (num_rejected_other > 0) then
+ print*, ' Number of mismatches in time : ', num_mismatch_time
+ print*, ' Number of mismatches in type : ', num_mismatch_type
+ print*, ' Number of mismatches in loc : ', num_mismatch_loc
+ endif
print*, 'Number of obs copied to output : ', num_inserted
print*, '---------------------------------------------------------'
- endif
- call destroy_obs_sequence(seq_in1)
- call destroy_obs_sequence(seq_in2)
- if (comparing3) call destroy_obs_sequence(seq_in3)
- if (comparing4) call destroy_obs_sequence(seq_in4)
+ write(msgstring, *) 'Starting to write output sequence files'
+ call error_handler(E_MSG,'obs_common_subset',msgstring)
- filename_out1 = trim(filename_seq1(i))//trim(filename_out_suffix)
- filename_out2 = trim(filename_seq2(i))//trim(filename_out_suffix)
- if (comparing3) filename_out3 = trim(filename_seq3(i))//trim(filename_out_suffix)
- if (comparing4) filename_out4 = trim(filename_seq4(i))//trim(filename_out_suffix)
-
- write(msgstring, *) 'Starting to write output sequence files'
- call error_handler(E_MSG,'obs_common_subset',msgstring)
-
- print*, 'Number of obs in the output seq file :', get_num_key_range(seq_out1)
- print*, 'and :', get_num_key_range(seq_out2)
- if (comparing3) print*, 'and :', get_num_key_range(seq_out3)
- if (comparing4) print*, 'and :', get_num_key_range(seq_out4)
-
- call print_obs_seq(seq_out1, filename_out1)
- call print_obs_seq(seq_out2, filename_out2)
- if (comparing3) call print_obs_seq(seq_out3, filename_out3)
- if (comparing4) call print_obs_seq(seq_out4, filename_out4)
- if (.not. print_only) then
- call write_obs_seq(seq_out1, filename_out1)
- call write_obs_seq(seq_out2, filename_out2)
- if (comparing3) call write_obs_seq(seq_out3, filename_out3)
- if (comparing4) call write_obs_seq(seq_out4, filename_out4)
- else
- write(msgstring,*) 'Output sequence files not created; print_only in namelist is .true.'
- call error_handler(E_MSG,'', msgstring)
+ do i = 1, atonce
+ call write_obs_seq(seq_out(i), filename_out(i))
+ enddo
endif
- ! clean up
+ ! clean up - prev_obs_out is a copy of obs_out, so don't call destroy on it
+ ! because it's already been deleted.
- call destroy_obs_sequence(seq_out1)
- call destroy_obs( obs_out1)
- call destroy_obs( obs_in1)
- call destroy_obs( next_obs_in1)
+ do i = 1, atonce
+ call destroy_obs_sequence(seq_in(i))
+ call destroy_obs_sequence(seq_out(i))
+ call destroy_obs(obs_in(i))
+ call destroy_obs(obs_out(i))
+ call destroy_obs(next_obs_in(i))
+ enddo
- call destroy_obs_sequence(seq_out2)
- call destroy_obs( obs_out2)
- call destroy_obs( obs_in2)
- call destroy_obs( next_obs_in2)
+enddo NUMSETS
- if (comparing3) then
- call destroy_obs_sequence(seq_out3)
- call destroy_obs( obs_out3)
- call destroy_obs( obs_in3)
- call destroy_obs( next_obs_in3)
- endif
-
- if (comparing4) then
- call destroy_obs_sequence(seq_out4)
- call destroy_obs( obs_out4)
- call destroy_obs( obs_in4)
- call destroy_obs( next_obs_in4)
- endif
-
- !call destroy_obs(prev_obs_out1) ! these are copies of something already
- !call destroy_obs(prev_obs_out2) ! deleted; duplicate dels if called.
-
-enddo
-
call shutdown()
!---------------------------------------------------------------------
@@ -533,6 +460,7 @@
!---------------------------------------------------------------------
+
subroutine setup()
! Initialize modules used that require it
@@ -543,6 +471,7 @@
end subroutine setup
!---------------------------------------------------------------------
+
subroutine shutdown()
call finalize_utilities('obs_common_subset')
@@ -550,257 +479,180 @@
end subroutine shutdown
!---------------------------------------------------------------------
-subroutine handle_filenames(filename_seq1, filename_seq_list1, &
- filename_seq2, filename_seq_list2, &
- filename_seq3, filename_seq_list3, &
- filename_seq4, filename_seq_list4, &
- num_input_files)
-! sort out the input lists, set the length as a return in num_input_files,
-! and make sure what's specified is consistent.
+subroutine handle_filenames(setsize, filename_seq, filename_seq_list, &
+ num_input_sets)
-character(len=*), intent(inout) :: filename_seq1(:), filename_seq2(:)
-character(len=*), intent(inout) :: filename_seq3(:), filename_seq4(:)
-character(len=*), intent(in) :: filename_seq_list1, filename_seq_list2
-character(len=*), intent(in) :: filename_seq_list3, filename_seq_list4
-integer, intent(out) :: num_input_files
+! sort out the input lists, set the length as a return in num_input_sets,
+! fill in filename_seq if list was given, and make sure what's specified is consistent.
-integer :: indx
-logical :: from_file1, from_file2, from_file3, from_file4
-character(len=32) :: source1, source2, source3, source4
+integer, intent(in) :: setsize
+character(len=*), intent(inout) :: filename_seq(:)
+character(len=*), intent(in) :: filename_seq_list
+integer, intent(out) :: num_input_sets
+integer :: indx, num_input_files
+logical :: from_file
+character(len=32) :: source
+
! if the user specifies neither filename_seq nor filename_seq_list, error
! if the user specifies both, error.
-! if the two list lengths aren't equal, error.
-! if the user gives one or more filelist names, we make sure the length(s) are
-! not more than maxfiles and read it/them into the explicit list and continue.
-! set num_input_files to the number of pairs in the lists
+! if list length isn't a whole multiple of num_to_compare_at_once, error.
+! if user gives a filelist name, make sure the length is not more than
+! maxfiles and read it into the explicit list and continue.
+! set num_input_sets to the count of sets in the list
-if (filename_seq1(1) == '' .and. filename_seq_list1 == '') then
+if (filename_seq(1) == '' .and. filename_seq_list == '') then
+ msgstring2='One of filename_seq or filename_seq_list must be specified in namelist'
call error_handler(E_ERR,'handle_filenames', &
- 'no filenames specified as input 1', &
- source,revision,revdate)
+ 'no filenames specified as input', &
+ source,revision,revdate, text2=msgstring2)
endif
-if (filename_seq2(1) == '' .and. filename_seq_list2 == '') then
- call error_handler(E_ERR,'handle_filenames', &
- 'no filenames specified as input 2', &
- source,revision,revdate)
-endif
-if (filename_seq3(1) == '' .and. filename_seq_list3 == '') then
- comparing3 = .false.
-endif
-if (filename_seq4(1) == '' .and. filename_seq_list4 == '') then
- comparing4 = .false.
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list