[Dart-dev] [5513] DART/branches/development/obs_sequence/obs_common_subset.f90: Untested version that will compare UP TO 4 files and output JUST those
nancy at ucar.edu
nancy at ucar.edu
Tue Jan 17 15:16:12 MST 2012
Revision: 5513
Author: thoar
Date: 2012-01-17 15:16:11 -0700 (Tue, 17 Jan 2012)
Log Message:
-----------
Untested version that will compare UP TO 4 files and output JUST those
observations that were successfully assimilated in all N files.
The operative term here is 'untested'.
Modified Paths:
--------------
DART/branches/development/obs_sequence/obs_common_subset.f90
-------------- next part --------------
Modified: DART/branches/development/obs_sequence/obs_common_subset.f90
===================================================================
--- DART/branches/development/obs_sequence/obs_common_subset.f90 2012-01-17 18:26:09 UTC (rev 5512)
+++ DART/branches/development/obs_sequence/obs_common_subset.f90 2012-01-17 22:16:11 UTC (rev 5513)
@@ -28,7 +28,8 @@
find_namelist_in_file, check_namelist_read, &
error_handler, E_ERR, E_MSG, nmlfileunit, &
do_nml_file, do_nml_term, get_next_filename, &
- open_file, close_file, finalize_utilities
+ open_file, close_file, finalize_utilities, &
+ logfileunit
use location_mod, only : location_type, get_location, set_location, &
LocationName, read_location, operator(/=), &
write_location
@@ -61,20 +62,31 @@
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
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 = metadatalength) :: dart_qc_meta_data = 'DART quality control'
character(len = metadatalength) :: meta_data
@@ -93,20 +105,25 @@
integer, parameter :: max_num_input_files = 1000
integer, parameter :: max_obs_input_types = 500
-
character(len = 129) :: filename_seq1(max_num_input_files) = ''
-character(len = 129) :: filename_seq_list1 = ''
character(len = 129) :: filename_seq2(max_num_input_files) = ''
-character(len = 129) :: filename_seq_list2 = ''
+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 = ''
character(len = 32) :: filename_out_suffix = '.common'
logical :: print_only = .false.
character(len=32) :: calendar = 'Gregorian'
-
namelist /obs_common_subset_nml/ &
- filename_seq1, filename_seq_list1, filename_out_suffix, &
- filename_seq2, filename_seq_list2, print_only, calendar
+ 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
!----------------------------------------------------------------
! Start of the program:
@@ -127,7 +144,7 @@
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 2 lists of obs_seq files; either in the namelist
+! the user gives us 4 lists 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.
! as in the other tools, it is an error to specify both an explicit
@@ -135,6 +152,8 @@
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)
! the default is a gregorian calendar. if you are using a different type
@@ -157,17 +176,33 @@
! count of input files was set in the handle_filenames() routine above.
do i = 1, num_input_files
- if (len(filename_seq1(i)) .eq. 0 .or. filename_seq1(i) .eq. "" ) then
+ if ((len(filename_seq1(i)) .eq. 0) .or. (filename_seq1(i) .eq. "")) 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)) .eq. 0 .or. filename_seq2(i) .eq. "" ) then
+ if ((len(filename_seq2(i)) .eq. 0) .or. (filename_seq2(i) .eq. "")) 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
+ write(*,*)'filename_seq3(i) is ',filename_seq3(i), len(filename_seq3(i)), (filename_seq3(i) .eq. "")
+
+ if ((len(filename_seq3(i)) .eq. 0) .or. (filename_seq3(i) .eq. "")) 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
+
+ if ((len(filename_seq4(i)) .eq. 0) .or. (filename_seq4(i) .eq. "")) 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
+
! read in the next pair of files.
call read_obs_seq_header(filename_seq1(i), num_copies_in1, num_qc_in1, &
@@ -190,15 +225,47 @@
cycle
endif
- write(msgstring, *) 'Starting to process input sequence files: '
- write(msgstring1,*) trim(filename_seq1(i))
- write(msgstring2,*) trim(filename_seq2(i))
- call error_handler(E_MSG,'obs_common_subset',msgstring, &
- text2=msgstring1, text3=msgstring2)
+ 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, &
+ close_the_file = .true.)
- 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 (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
+ 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.)
+
+ 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
+ endif
+ endif
+
+ 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(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)
+
! 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))
@@ -208,14 +275,16 @@
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))
+ 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))
- ! the output can have no more than the smaller number of input obs
+ ! the output can have no more than the max number of input obs
! these can be different sizes - we are going to copy only the matching obs
size_seq_in1 = max_num_obs1
size_seq_in2 = max_num_obs2
- size_seq_out = min(size_seq_in1, size_seq_in2)
+ size_seq_out = max(size_seq_in1, size_seq_in2)
! 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
@@ -232,27 +301,49 @@
call error_handler(E_MSG,' ',' ')
! Initialize individual observation variables
- call init_obs( obs_in1, num_copies_in, num_qc_in)
- call init_obs( obs_in2, num_copies_in, num_qc_in)
- call init_obs(next_obs_in1, num_copies_in, num_qc_in)
- call init_obs(next_obs_in2, num_copies_in, num_qc_in)
+ 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)
+
+ 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_out1, num_copies_in, num_qc_in)
call init_obs(prev_obs_out2, num_copies_in, num_qc_in)
+
+ 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
+ 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)
+ 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)
+ 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)
+ 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 why these would be needed?
@@ -263,8 +354,10 @@
!size_seq_in = get_num_key_range(seq_in) !current size of seq_in
! 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) 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
@@ -279,35 +372,73 @@
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
+ endif
+
+ 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.
- next_obs_in1 = obs_in1
- next_obs_in2 = obs_in2
+ 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
ObsLoop : do while ( .not. is_this_last1 .and. .not. is_this_last2)
- obs_in1 = next_obs_in1 ! essentially records position in seq_out
- obs_in2 = next_obs_in2
+ ! 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
- if (good_match(obs_in1, obs_in2, qc_index, qc_threshold)) then
- obs_out1 = obs_in1
- obs_out2 = obs_in2
+ 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)
+ else
+ wanted = all_good(obs_in1, obs_in2, qc_index, qc_threshold)
+ endif
+ if ( wanted ) then
+ obs_out1 = obs_in1
+ obs_out2 = obs_in2
+ if (comparing3) obs_out3 = obs_in3
+ if (comparing4) obs_out4 = obs_in4
+
! 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
- call insert_obs_in_seq(seq_out1, obs_out1, prev_obs_out1)
- call insert_obs_in_seq(seq_out2, obs_out2, prev_obs_out2)
+ 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)
+ 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)
endif
- prev_obs_out1 = obs_out1 ! update position in seq for next insert
- prev_obs_out2 = obs_out2
+ ! 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
if (print_every > 0) then
@@ -318,8 +449,10 @@
endif
- 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)
+ 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
@@ -340,26 +473,33 @@
print*, '---------------------------------------------------------'
endif
- call destroy_obs_sequence(seq_in1)
- call destroy_obs_sequence(seq_in2)
+ 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)
-
- filename_out1 = trim(filename_seq1(i))//trim(filename_out_suffix)
- filename_out2 = trim(filename_seq2(i))//trim(filename_out_suffix)
+ 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 process output sequence files ', &
- trim(filename_out1)
- write(msgstring1,*) 'and ', trim(filename_out2)
- call error_handler(E_MSG,'obs_common_subset',msgstring, text2=msgstring1)
+ 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)
+ 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)
+ 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)
+ 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)
@@ -368,13 +508,29 @@
! clean up
call destroy_obs_sequence(seq_out1)
+ call destroy_obs( obs_out1)
+ call destroy_obs( obs_in1)
+ call destroy_obs( next_obs_in1)
+
call destroy_obs_sequence(seq_out2)
- call destroy_obs( obs_in1 )
- call destroy_obs( obs_in2 )
- call destroy_obs(next_obs_in1 )
- call destroy_obs(next_obs_in2 )
- call destroy_obs( obs_out1)
- call destroy_obs( obs_out2)
+ call destroy_obs( obs_out2)
+ call destroy_obs( obs_in2)
+ call destroy_obs( next_obs_in2)
+
+ 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.
@@ -409,18 +565,23 @@
!---------------------------------------------------------------------
subroutine handle_filenames(filename_seq1, filename_seq_list1, &
- filename_seq2, filename_seq_list2, num_input_files)
+ 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.
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
integer :: indx
-logical :: from_file1, from_file2
-character(len=32) :: source1, source2
+logical :: from_file1, from_file2, from_file3, from_file4
+character(len=32) :: source1, source2, source3, source4
! if the user specifies neither filename_seq nor filename_seq_list, error
! if the user specifies both, error.
@@ -439,10 +600,22 @@
'no filenames specified as input 2', &
source,revision,revdate)
endif
+if (filename_seq3(1) == '' .and. filename_seq_list3 == '') then
+ call error_handler(E_ERR,'obs_common_subset', &
+ 'no filenames specified as input 3', &
+ source,revision,revdate)
+endif
+if (filename_seq4(1) == '' .and. filename_seq_list4 == '') then
+ call error_handler(E_ERR,'obs_common_subset', &
+ 'no filenames specified as input 4', &
+ source,revision,revdate)
+endif
! make sure the namelist specifies one or the other but not both
if ((filename_seq1(1) /= '' .and. filename_seq_list1 /= '') .or. &
- (filename_seq2(1) /= '' .and. filename_seq_list2 /= '')) then
+ (filename_seq2(1) /= '' .and. filename_seq_list2 /= '') .or. &
+ (filename_seq3(1) /= '' .and. filename_seq_list3 /= '') .or. &
+ (filename_seq4(1) /= '' .and. filename_seq_list4 /= '')) then
call error_handler(E_ERR,'obs_common_subset', &
'cannot specify both filename_seq and filename_seq_list', &
source,revision,revdate)
@@ -457,6 +630,7 @@
source1 = 'filename_seq1'
from_file1 = .false.
endif
+
if (filename_seq_list2 /= '') then
source2 = 'filename_seq_list2'
from_file2 = .true.
@@ -465,12 +639,37 @@
from_file2 = .false.
endif
+if (filename_seq_list3 /= '') then
+ source3 = 'filename_seq_list3'
+ from_file3 = .true.
+else
+ source3 = 'filename_seq3'
+ from_file3 = .false.
+endif
+
+if (filename_seq_list4 /= '') then
+ source4 = 'filename_seq_list4'
+ from_file4 = .true.
+else
+ source4 = 'filename_seq4'
+ from_file4 = .false.
+endif
+
+write(*,*)'filename_seq1(1) ',trim(filename_seq1(1)),' ', trim(source1), ' ',from_file1
+write(*,*)'filename_seq2(1) ',trim(filename_seq2(1)),' ', trim(source2), ' ',from_file2
+write(*,*)'filename_seq3(1) ',trim(filename_seq3(1)),' ', trim(source3), ' ',from_file3
+write(*,*)'filename_seq4(1) ',trim(filename_seq4(1)),' ', trim(source4), ' ',from_file4
+
! the point of this loop is to count up how many pairs of input seq files we have.
do indx = 1, max_num_input_files
if (from_file1) &
filename_seq1(indx) = get_next_filename(filename_seq_list1, indx)
if (from_file2) &
filename_seq2(indx) = get_next_filename(filename_seq_list2, indx)
+ if (from_file3) &
+ filename_seq3(indx) = get_next_filename(filename_seq_list3, indx)
+ if (from_file4) &
+ filename_seq4(indx) = get_next_filename(filename_seq_list4, indx)
! a pair of empty names ends the list and we return with the count.
! (unless both lists are empty and then we're unhappy)
@@ -873,72 +1072,132 @@
! compare location, time, type, QC, plus make sure QC <= threshold
!---------------------------------------------------------------------
-function good_match(obs_in1, obs_in2, qc_index, qc_threshold)
- type(obs_type), intent(in) :: obs_in1, obs_in2
- integer, intent(in) :: qc_index, qc_threshold
- logical :: good_match
+function all_good(obs_in1, obs_in2, qc_index, qc_threshold, obs_in3, obs_in4)
+type(obs_type), intent(in) :: obs_in1, obs_in2
+integer, intent(in) :: qc_index, qc_threshold
+type(obs_type), optional, intent(in) :: obs_in3, obs_in4
+logical :: all_good
-type(obs_def_type) :: base_obs_def, test_obs_def
-integer :: base_obs_type, test_obs_type
-type(time_type) :: base_obs_time, test_obs_time
-type(location_type) :: base_obs_loc, test_obs_loc
-integer :: base_qc, test_qc
+type(obs_def_type) :: test1_obs_def, test2_obs_def, test3_obs_def, test4_obs_def
+integer :: test1_obs_type, test2_obs_type, test3_obs_type, test4_obs_type
+type(time_type) :: test1_obs_time, test2_obs_time, test3_obs_time, test4_obs_time
+type(location_type) :: test1_obs_loc, test2_obs_loc, test3_obs_loc, test4_obs_loc
+integer :: test1_qc, test2_qc, test3_qc, test4_qc
real(r8) :: temp(1)
-
! assume failure so we can return as soon as we know they don't match.
-good_match = .false.
+all_good = .false.
-call get_obs_def(obs_in1, base_obs_def)
-base_obs_loc = get_obs_def_location(base_obs_def)
-base_obs_time = get_obs_def_time(base_obs_def)
-base_obs_type = get_obs_kind(base_obs_def)
+call get_obs_def(obs_in1, test1_obs_def)
+test1_obs_loc = get_obs_def_location(test1_obs_def)
+test1_obs_time = get_obs_def_time(test1_obs_def)
+test1_obs_type = get_obs_kind(test1_obs_def)
if (qc_index < 0) then
- base_qc = 0
+ test1_qc = 0
else
call get_qc(obs_in1, temp, qc_index)
- base_qc = nint(temp(1))
+ test1_qc = nint(temp(1))
endif
-call get_obs_def(obs_in2, test_obs_def)
-test_obs_loc = get_obs_def_location(test_obs_def)
-test_obs_time = get_obs_def_time(test_obs_def)
-test_obs_type = get_obs_kind(test_obs_def)
+call get_obs_def(obs_in2, test2_obs_def)
+test2_obs_loc = get_obs_def_location(test2_obs_def)
+test2_obs_time = get_obs_def_time(test2_obs_def)
+test2_obs_type = get_obs_kind(test2_obs_def)
if (qc_index < 0) then
- test_qc = 0
+ test2_qc = 0
else
call get_qc(obs_in2, temp, qc_index)
- test_qc = nint(temp(1))
+ test2_qc = nint(temp(1))
endif
+if (present(obs_in3)) then
+ call get_obs_def(obs_in3, test3_obs_def)
+ test3_obs_loc = get_obs_def_location(test3_obs_def)
+ test3_obs_time = get_obs_def_time(test3_obs_def)
+ test3_obs_type = get_obs_kind(test3_obs_def)
+ if (qc_index < 0) then
+ test3_qc = 0
+ else
+ call get_qc(obs_in3, temp, qc_index)
+ test3_qc = nint(temp(1))
+ endif
+else
+ test3_qc = test1_qc
+ test3_obs_type = test1_obs_type
+ test3_obs_time = test1_obs_time
+ test3_obs_loc = test1_obs_loc
+endif
+
+if (present(obs_in4)) then
+ call get_obs_def(obs_in4, test4_obs_def)
+ test4_obs_loc = get_obs_def_location(test4_obs_def)
+ test4_obs_time = get_obs_def_time(test4_obs_def)
+ test4_obs_type = get_obs_kind(test4_obs_def)
+ if (qc_index < 0) then
+ test4_qc = 0
+ else
+ call get_qc(obs_in4, temp, qc_index)
+ test4_qc = nint(temp(1))
+ endif
+else
+ test4_qc = test1_qc
+ test4_obs_type = test1_obs_type
+ test4_obs_time = test1_obs_time
+ test4_obs_loc = test1_obs_loc
+endif
+
! first compare the integer values; it's quicker than location.
! in this case leave the QC until last so we can get statistics on why
! an obs is rejected. if a lot of obs fail the other tests, then the input
! files aren't as identical as they're supposed to be.
-if ((base_obs_type /= test_obs_type) .or. &
- (base_obs_time /= test_obs_time) .or. &
- (base_obs_loc /= test_obs_loc)) then
+if ((test1_obs_type /= test2_obs_type) .or. &
+ (test1_obs_type /= test3_obs_type) .or. &
+ (test1_obs_type /= test4_obs_type) .or. &
+ (test2_obs_type /= test3_obs_type) .or. &
+ (test2_obs_type /= test4_obs_type) .or. &
+ (test3_obs_type /= test4_obs_type)) then
num_rejected_other = num_rejected_other + 1
return
endif
-if (base_qc /= test_qc) then
+if ((test1_obs_time /= test2_obs_time) .or. &
+ (test1_obs_time /= test3_obs_time) .or. &
+ (test1_obs_time /= test4_obs_time) .or. &
+ (test2_obs_time /= test3_obs_time) .or. &
+ (test2_obs_time /= test4_obs_time) .or. &
+ (test3_obs_time /= test4_obs_time)) then
+ num_rejected_other = num_rejected_other + 1
+ return
+endif
+
+if ((test1_obs_loc /= test2_obs_loc) .or. &
+ (test1_obs_loc /= test3_obs_loc) .or. &
+ (test1_obs_loc /= test4_obs_loc) .or. &
+ (test2_obs_loc /= test3_obs_loc) .or. &
+ (test2_obs_loc /= test4_obs_loc) .or. &
+ (test3_obs_loc /= test4_obs_loc)) then
+ num_rejected_other = num_rejected_other + 1
+ return
+endif
+
+! TJH FIXME left off here ...
+if (test1_qc /= test2_qc) then
num_rejected_diffqc = num_rejected_diffqc + 1
return
endif
! this assumes we have already tested for both qcs being equal, so it can
! just test one and know whether they are both over the threshold or not.
-if (base_qc > qc_threshold) then
+if (test1_qc > qc_threshold) then
num_rejected_badqc = num_rejected_badqc + 1
return
endif
! all match - good return.
-good_match = .true.
+all_good = .true.
-end function good_match
+end function all_good
!---------------------------------------------------------------------
! currently unused:
More information about the Dart-dev
mailing list