[Dart-dev] [5701] DART/branches/development: copied the index sort routine from the sort module and made
nancy at ucar.edu
nancy at ucar.edu
Mon Apr 16 16:03:13 MDT 2012
Revision: 5701
Author: nancy
Date: 2012-04-16 16:03:12 -0600 (Mon, 16 Apr 2012)
Log Message:
-----------
copied the index sort routine from the sort module and made
a time sort routine in the time_manager module. major revamp
of the obs_selection tool - was doing a naive iteration of the
entire selection file which was too slow for real world problems.
i now sort the selection list by time, and then can start at
the last offset and stop searching when i've passed the time
of interest. this has speeded up the tool to where i believe
it's useful again. also fixed some message strings that still
said 'obs_sequence_tool' instead of 'obs_selection', and fixed
a place where i overwrote 'source' instead of using a different
variable name.
Modified Paths:
--------------
DART/branches/development/obs_sequence/obs_selection.f90
DART/branches/development/time_manager/time_manager_mod.f90
-------------- next part --------------
Modified: DART/branches/development/obs_sequence/obs_selection.f90
===================================================================
--- DART/branches/development/obs_sequence/obs_selection.f90 2012-04-16 16:03:56 UTC (rev 5700)
+++ DART/branches/development/obs_sequence/obs_selection.f90 2012-04-16 22:03:12 UTC (rev 5701)
@@ -2,6 +2,19 @@
! provided by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
+! nsc 12apr2012 -
+! was too slow for large lists and large obs_seq files.
+! sorted the obs_def list by time and then started the search
+! at the time of the next obs, and quit looping when past that time.
+! really speeded up - may not have to do anything more complicated
+! at this time.
+! if more speed is needed, next likely place to pick up speed is
+! by sorting all obs at the same time by locations (maybe just
+! by the x coord for starters), or binning in spatial bins if
+! still too slow. but the current fixes should go a long way
+! to making the performance acceptable.
+!
+
program obs_selection
! <next few lines under version control, do not edit>
@@ -27,7 +40,8 @@
read_obs_kind
use time_manager_mod, only : time_type, operator(>), print_time, set_time, &
print_date, set_calendar_type, GREGORIAN, &
- operator(/=), NO_CALENDAR, get_calendar_type
+ operator(/=), operator(<=), NO_CALENDAR, &
+ get_calendar_type, time_index_sort, operator(==)
use obs_sequence_mod, only : obs_sequence_type, obs_type, write_obs_seq, &
init_obs, assignment(=), get_obs_def, &
init_obs_sequence, static_init_obs_sequence, &
@@ -54,18 +68,18 @@
type(obs_sequence_type) :: seq_in, seq_out
type(obs_type) :: obs_in, next_obs_in
type(obs_type) :: obs_out, prev_obs_out
+type(time_type) :: t1, t2
logical :: is_there_one, is_this_last
integer :: size_seq_in, num_copies_in, num_qc_in
integer :: size_seq_out, num_copies_out, num_qc_out
-integer :: num_inserted, iunit, io, i, j, total_num_inserted
+integer :: num_inserted, iunit, io, i, j
+integer :: total_num_inserted, base_index, next_base_index
integer :: max_num_obs, file_id
integer :: first_seq
character(len = 129) :: read_format, meta_data
logical :: pre_I_format, cal
-character(len = 129) :: msgstring
+character(len = 255) :: msgstring, msgstring1, msgstring2, msgstring3
-! could go into namelist if you wanted more control
-integer, parameter :: print_every = 20
!----------------------------------------------------------------
! Namelist input with default values
@@ -77,6 +91,7 @@
integer :: num_input_files = 0
type(obs_def_type), allocatable :: obs_def_list(:)
integer :: obs_def_count
+integer :: print_every_nth_obs = 100
character(len = 129) :: filename_seq(max_num_input_files) = ''
@@ -216,7 +231,7 @@
! Read obs seq to be added, and insert obs from it to the output seq
first_seq = -1
-do i = 1, num_input_files
+FILES: do i = 1, num_input_files
if (.not. process_file(i)) cycle
@@ -268,7 +283,21 @@
if ( is_there_one ) then
- if (good_selection(obs_in, obs_def_list, obs_def_count)) then
+ ! figure out the time of the first obs in this file, and set the
+ ! offset for the first item in the selection file that's at this time.
+ t1 = get_time_from_obs(obs_in)
+ base_index = set_base(t1, obs_def_list, obs_def_count)
+
+ if (base_index < 0) then
+ write(msgstring2, *) 'skipping all obs in ', trim(filename_seq(i))
+ call error_handler(E_MSG, 'obs_selection', &
+ 'first time in obs_sequence file is after all times in selection file', &
+ source, revision, revdate, text2=msgstring2)
+ cycle FILES
+ endif
+ next_base_index = base_index
+
+ if (good_selection(obs_in, obs_def_list, obs_def_count, next_base_index)) then
obs_out = obs_in
call insert_obs_in_seq(seq_out, obs_out) ! new_obs linked list info changes
@@ -280,9 +309,27 @@
call get_next_obs(seq_in, obs_in, next_obs_in, is_this_last)
ObsLoop : do while ( .not. is_this_last)
- obs_in = next_obs_in ! essentially records position in seq_out
+ ! if the time of the next obs is different from this one, bump
+ ! up the start of the search index offset.
+ t1 = get_time_from_obs(obs_in)
+ t2 = get_time_from_obs(next_obs_in)
+ if (t1 /= t2) then
+ base_index = next_base_index
+ next_base_index = set_base(t2, obs_def_list, obs_def_count, base_index)
+
+ if (next_base_index < 0) then
+ ! next obs in selection file is after all the rest of the obs in this
+ ! input obs_seq file, so we can move on now.
+ write(msgstring, *) 'done with ', trim(filename_seq(i))
+ call error_handler(E_MSG, 'obs_selection', msgstring, &
+ source, revision, revdate)
+ cycle FILES
+ endif
+ endif
- if (good_selection(obs_in, obs_def_list, obs_def_count)) then
+ obs_in = next_obs_in ! next obs is the current one now
+
+ if (good_selection(obs_in, obs_def_list, obs_def_count, next_base_index)) then
obs_out = obs_in
! Since the stride through the observation sequence file is always
@@ -299,9 +346,11 @@
prev_obs_out = obs_out ! update position in seq_in for next insert
num_inserted = num_inserted + 1
- if (print_every > 0) then
- if (mod(num_inserted,print_every) == 0) then
- print*, 'inserted number ',num_inserted,' of ',size_seq_in
+ if (print_every_nth_obs > 0) then
+ if (mod(num_inserted,print_every_nth_obs) == 0) then
+ write(msgstring,*) 'inserted number ',num_inserted,' of ',size_seq_in, ' possible obs'
+ call error_handler(E_MSG, 'obs_selection', msgstring, &
+ source, revision, revdate)
endif
endif
@@ -309,6 +358,7 @@
call get_next_obs(seq_in, obs_in, next_obs_in, is_this_last)
+
enddo ObsLoop
total_num_inserted = total_num_inserted + num_inserted
@@ -328,10 +378,10 @@
call destroy_obs_sequence(seq_in)
-enddo
+enddo FILES
write(msgstring,*) 'Starting to process output sequence file ', trim(filename_out)
-call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+call error_handler(E_MSG,'obs_selection',msgstring)
if (.not. print_only) then
print*, 'Total number of obs inserted : ', total_num_inserted
@@ -365,7 +415,7 @@
subroutine obs_seq_modules_used()
! Initialize modules used that require it
-call initialize_utilities('obs_sequence_tool')
+call initialize_utilities('obs_selection')
call register_module(source,revision,revdate)
call static_init_obs_sequence()
@@ -381,7 +431,7 @@
integer :: index
logical :: from_file
-character(len=32) :: source
+character(len=32) :: fsource
! ok, here's the new logic:
! if the user specifies neither filename_seq nor filename_seq_list, we
@@ -399,7 +449,7 @@
if (filename_seq(1) == '' .and. filename_seq_list == '') then
if (num_input_files /= 0 .and. num_input_files /= 1) then
- call error_handler(E_ERR,'obs_sequence_tool', &
+ call error_handler(E_ERR,'obs_selection', &
'if no filenames specified, num_input_files must be 0 or 1', &
source,revision,revdate)
endif
@@ -411,7 +461,7 @@
! make sure the namelist specifies one or the other but not both
if (filename_seq(1) /= '' .and. filename_seq_list /= '') then
- call error_handler(E_ERR,'obs_sequence_tool', &
+ call error_handler(E_ERR,'obs_selection', &
'cannot specify both filename_seq and filename_seq_list', &
source,revision,revdate)
endif
@@ -419,10 +469,10 @@
! if they have specified a file which contains a list, read it into
! the filename_seq array and set the count.
if (filename_seq_list /= '') then
- source = 'filename_seq_list'
+ fsource = 'filename_seq_list'
from_file = .true.
else
- source = 'filename_seq'
+ fsource = 'filename_seq'
from_file = .false.
endif
@@ -432,8 +482,8 @@
if (filename_seq(index) == '') then
if (index == 1) then
- call error_handler(E_ERR,'obs_sequence_tool', &
- trim(source)//' contains no filenames', &
+ call error_handler(E_ERR,'obs_selection', &
+ 'namelist item ', trim(fsource)//' contains no filenames', &
source,revision,revdate)
endif
! leaving num_input_files unspecified (or set to 0) means use
@@ -445,21 +495,21 @@
! if they do give a count, make it match.
if (num_input_files == (index - 1)) return
- write(msgstring, *) 'if num_input_files is 0, the number of files will be automatically computed'
- call error_handler(E_MSG,'obs_sequence_tool', msgstring)
- write(msgstring, *) 'if num_input_files is not 0, it must match the number of filenames specified'
- call error_handler(E_MSG,'obs_sequence_tool', msgstring)
write(msgstring, *) 'num_input_files is ', num_input_files, &
- ' but '//trim(source)//' has filecount ', index - 1
- call error_handler(E_ERR,'obs_sequence_tool', msgstring, &
- source,revision,revdate)
+ ' but namelist item '//trim(fsource)//' has filecount ', index - 1
+
+ write(msgstring2, *) 'if num_input_files is 0, the number of files will be automatically computed'
+ write(msgstring3, *) 'if num_input_files is not 0, it must match the number of filenames specified'
+
+ call error_handler(E_ERR,'obs_selection', msgstring, &
+ source,revision,revdate, text2=msgstring2, text3=msgstring3)
endif
endif
enddo
write(msgstring, *) 'cannot specify more than ',max_num_input_files,' files'
-call error_handler(E_ERR,'obs_sequence_tool', msgstring, &
+call error_handler(E_ERR,'obs_selection', msgstring, &
source,revision,revdate)
end subroutine handle_filenames
@@ -488,7 +538,6 @@
integer :: num_copies2, num_qc2
integer :: num_copies , num_qc, i
character(len=129) :: str1, str2
-character(len=255) :: msgstring1, msgstring2
num_copies1 = get_num_copies(seq1)
num_qc1 = get_num_qc( seq1)
@@ -499,7 +548,7 @@
num_copies = num_copies2
num_qc = num_qc2
-! get this ready in case we have to use it
+! get this ready in case we have to use it below. do not overwrite it!
if (present(fname1) .and. present(fname2)) then
write(msgstring1,*)'Sequence files ', trim(fname1), ' and ', trim(fname2), &
' are not compatible'
@@ -510,42 +559,38 @@
if ( num_copies1 /= num_copies2 ) then
write(msgstring2,*)'Different numbers of data copies found: ', &
num_copies1, ' vs ', num_copies2
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- num_copies = -1
+ call error_handler(E_ERR, 'obs_selection', msgstring1, &
+ source, revision, revdate, text2=msgstring2)
endif
if ( num_qc1 /= num_qc2 ) then
write(msgstring2,*)'Different different numbers of QCs found: ', &
num_qc1, ' vs ', num_qc2
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- num_qc = -1
+ call error_handler(E_MSG, 'obs_selection', msgstring1, &
+ source, revision, revdate, text2=msgstring2)
endif
-if ( num_copies < 0 .or. num_qc < 0 ) then
- call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
-endif
! watch the code flow in this loop and the one below it.
-! the smoothest code order is to determine what the strings are,
-! and then try different things to match them. if a match is found,
-! cycle. if you get to the bottom of the loop, there is no match
-! and a single set of (fatal) error messages is called.
+! if a match is found, cycle. if there is no match
+! a single set of (fatal) error messages is called.
CopyMetaData : do i=1, num_copies
str1 = get_copy_meta_data(seq1,i)
str2 = get_copy_meta_data(seq2,i)
- ! easy case - they match. cycle to next copy.
+ ! they match. cycle to next copy.
if( str1 == str2 ) then
- write(msgstring2,*)'metadata ',trim(str1), ' in both.'
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ ! for now, don't print out if things are ok. this could become
+ ! part of a verbose option.
+ !write(msgstring2,*)'metadata ',trim(str1), ' in both.'
+ !call error_handler(E_MSG, 'obs_selection', msgstring2)
cycle CopyMetaData
endif
! if you get here, the metadata is not the same and the user has not
! given us strings that are ok to match. fail.
write(msgstring2,*)'metadata value mismatch. seq1: ', trim(str1)
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- write(msgstring2,*)'metadata value mismatch. seq2: ', trim(str2)
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
+ write(msgstring3,*)'metadata value mismatch. seq2: ', trim(str2)
+ call error_handler(E_ERR, 'obs_selection', msgstring1, &
+ source, revision, revdate, text2=msgstring2, text3=msgstring3)
enddo CopyMetaData
@@ -553,20 +598,20 @@
str1 = get_qc_meta_data(seq1,i)
str2 = get_qc_meta_data(seq2,i)
- ! easy case - they match. cycle to next copy.
+ ! they match. cycle to next copy.
if( str1 == str2 ) then
- write(msgstring2,*)'metadata ',trim(str1), ' in both.'
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ ! see comment in copy section above about a verbose option.
+ !write(msgstring2,*)'metadata ',trim(str1), ' in both.'
+ !call error_handler(E_MSG, 'obs_selection', msgstring2)
cycle QCMetaData
endif
! if you get here, the metadata is not the same and the user has not
! given us strings that are ok to match. fail.
write(msgstring2,*)'qc metadata value mismatch. seq1: ', trim(str1)
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- write(msgstring2,*)'qc metadata value mismatch. seq2: ', trim(str2)
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
+ write(msgstring3,*)'qc metadata value mismatch. seq2: ', trim(str2)
+ call error_handler(E_ERR, 'obs_selection', msgstring1, &
+ source, revision, revdate, text2=msgstring2, text3=msgstring3)
enddo QCMetaData
@@ -607,7 +652,7 @@
size_seq_in = get_num_obs(seq_in)
if (size_seq_in == 0) then
msgstring = 'Obs_seq file '//trim(filename)//' is empty.'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ call error_handler(E_MSG,'obs_selection',msgstring)
return
endif
@@ -630,7 +675,7 @@
if ( .not. is_there_one ) then
write(msgstring,*)'no first observation in ',trim(filename)
- call error_handler(E_MSG,'obs_sequence_tool', msgstring)
+ call error_handler(E_MSG,'obs_selection', msgstring)
endif
! process it here
@@ -649,8 +694,6 @@
else
type_count(this_obs_kind) = type_count(this_obs_kind) + 1
endif
-! print *, 'obs kind index = ', this_obs_kind
-! if(this_obs_kind > 0)print *, 'obs name = ', get_obs_kind_name(this_obs_kind)
call get_next_obs(seq_in, obs, next_obs, is_this_last)
if (.not. is_this_last) then
@@ -706,7 +749,6 @@
character(len=*), intent(in) :: filename
type(obs_type) :: obs, next_obs
-type(obs_def_type) :: this_obs_def
logical :: is_there_one, is_this_last
integer :: size_seq
integer :: key
@@ -717,7 +759,7 @@
size_seq = get_num_obs(seq)
if (size_seq == 0) then
msgstring = 'Obs_seq file '//trim(filename)//' is empty.'
- call error_handler(E_MSG,'obs_sequence_tool:validate',msgstring)
+ call error_handler(E_MSG,'obs_selection:validate',msgstring)
return
endif
@@ -732,18 +774,15 @@
if ( .not. is_there_one ) then
write(msgstring,*)'no first observation in sequence ' // trim(filename)
- call error_handler(E_MSG,'obs_sequence_tool:validate', msgstring, source, revision, revdate)
+ call error_handler(E_MSG,'obs_selection:validate', msgstring, source, revision, revdate)
endif
-call get_obs_def(obs, this_obs_def)
-last_time = get_obs_def_time(this_obs_def)
+last_time = get_time_from_obs(obs)
-
is_this_last = .false.
ObsLoop : do while ( .not. is_this_last)
- call get_obs_def(obs, this_obs_def)
- this_time = get_obs_def_time(this_obs_def)
+ this_time = get_time_from_obs(obs)
if (last_time > this_time) then
! bad time order of observations in linked list
@@ -753,10 +792,10 @@
if (cal) call print_date(this_time, ' which is date: ')
key = get_obs_key(obs)
- write(msgstring,*)'obs number ', key, ' has earlier time than previous obs'
- call error_handler(E_MSG,'obs_sequence_tool:validate', msgstring)
+ write(msgstring2,*)'obs number ', key, ' has earlier time than previous obs'
write(msgstring,*)'observations must be in increasing time order, file ' // trim(filename)
- call error_handler(E_ERR,'obs_sequence_tool:validate', msgstring, source, revision, revdate)
+ call error_handler(E_ERR,'obs_selection:validate', msgstring, &
+ source, revision, revdate, text2=msgstring2)
endif
last_time = this_time
@@ -784,7 +823,6 @@
integer :: num_copies , num_qc, i
character(len=129) :: str1
-character(len=255) :: msgstring1
num_copies = get_num_copies(seq1)
num_qc = get_num_qc( seq1)
@@ -832,6 +870,9 @@
type(obs_sequence_type) :: seq_in
logical :: is_this_last
real(r8) :: dummy
+ type(obs_def_type), allocatable :: temp_sel_list(:)
+ type(time_type), allocatable :: temp_time(:)
+ integer, allocatable :: sort_index(:)
! if the list of which obs to select comes from the coverage tool,
! it's a list of obs_defs. if it's a full obs_seq file, then
@@ -846,15 +887,36 @@
source,revision,revdate)
endif
- allocate(selection_list(count))
-
! set up the mapping table for the kinds here
call read_obs_kind(iunit, .false.)
+ ! this one stays around and is a return from this subroutine
+ allocate(selection_list(count))
+
+ ! these are temporaries for sorting into time order
+ allocate(temp_sel_list(count), temp_time(count), sort_index(count))
+
+ ! read into array in whatever order these are in the file
do i = 1, count
- call read_obs_def(iunit, selection_list(i), 0, dummy)
+ call read_obs_def(iunit, temp_sel_list(i), 0, dummy)
enddo
+ ! extract just the times into an array
+ do i = 1, count
+ temp_time(i) = get_obs_def_time(temp_sel_list(i))
+ enddo
+
+ ! sort into time order
+ call time_index_sort(temp_time, sort_index, count)
+
+ ! copy into output array using that order
+ do i = 1, count
+ selection_list(i) = temp_sel_list(sort_index(i))
+ enddo
+
+
+ deallocate(temp_sel_list, temp_time, sort_index)
+
call close_file(iunit)
else
@@ -875,6 +937,9 @@
source,revision,revdate)
endif
+ ! in an obs_seq file, using get_next_obs you will
+ ! be guarenteed to get these in increasing time order.
+
is_this_last = .false.
do i = 1, count
if (is_this_last) exit
@@ -891,15 +956,68 @@
selection_count = count
+write(msgstring, *) 'selection file contains ', count, ' entries'
+call error_handler(E_MSG, 'obs_selection', msgstring, source, revision, revdate)
+ call print_time(get_obs_def_time(selection_list(1)), 'time of first selection:')
+if (cal) call print_date(get_obs_def_time(selection_list(1)), ' which is date:')
+ call print_time(get_obs_def_time(selection_list(count)), 'time of last selection:')
+if (cal) call print_date(get_obs_def_time(selection_list(count)), ' which is date:')
+
end subroutine read_selection_list
+!---------------------------------------------------------------------
+function set_base(obs_time, selection_list, selection_count, startindex)
+! find offset of first obs_def entry that has a time >= to the given one
+! if optional startindex is given, start looking there.
+
+ type(time_type), intent(in) :: obs_time
+ type(obs_def_type), intent(in) :: selection_list(:)
+ integer, intent(in) :: selection_count
+ integer, optional, intent(in) :: startindex
+ integer :: set_base
+
+ integer :: i, s
+ type(time_type) :: def_time
+
+ ! if we know an offset to start from, use it. otherwise, start at 1.
+ if (present(startindex)) then
+ s = startindex
+ else
+ s = 1
+ endif
+
+ ! find the index of the first item in this list which is >= given one
+ ! we will start subsequent searches at this offset.
+ do i = s, selection_count
+
+ def_time = get_obs_def_time(selection_list(i))
+
+ ! if we have looped through the obs_def list far enough so
+ ! the next obs_def entry has a time more than or equal to the
+ ! one of the next observation, stop and return this index.
+ if (obs_time <= def_time) then
+ set_base = i
+ return
+ endif
+
+ enddo
+
+! we got all the way through the file and the last obs_def entry
+! was still before the observation time we are looking for.
+set_base = -1
+
+end function set_base
+
+
+
! compare horiz location, time, type - ignores vertical
!---------------------------------------------------------------------
-function good_selection(obs_in, selection_list, selection_count)
+function good_selection(obs_in, selection_list, selection_count, startindex)
type(obs_type), intent(in) :: obs_in
type(obs_def_type), intent(in) :: selection_list(:)
integer, intent(in) :: selection_count
+ integer, intent(in) :: startindex
logical :: good_selection
! first pass, iterate list.
@@ -917,17 +1035,39 @@
base_obs_time = get_obs_def_time(base_obs_def)
base_obs_type = get_obs_kind(base_obs_def)
+ ! this program now time-sorts the selection list first, so we
+ ! are guarenteed the obs_defs will be encountered in time order.
+ ! now optimize in two ways: first, the caller will pass in an
+ ! offset that is less than or equal to the first obs for this time,
+ ! and second this routine will return when the obs_def time is
+ ! larger than the time to match. that should save a lot of looping.
+
+ if (startindex < 1 .or. startindex > selection_count) then
+ write(msgstring2, *) 'startindex ', startindex, ' is not between 1 and ', selection_count
+ call error_handler(E_ERR, 'good_selection', &
+ 'invalid startindex, internal error should not happen', &
+ source, revision, revdate, text2=msgstring2)
+ endif
+
+ ! the obs_def list is time-sorted now. if you get to a larger
+ ! time without a match you can bail early.
+ good_selection = .false.
+
! first select on time - it is an integer comparison
! and quicker than location test. then type, and
! finally location.
- do i = 1, selection_count
+ do i = startindex, selection_count
+ test_obs_time = get_obs_def_time(selection_list(i))
+
+ ! if past the possible times, return now.
+ if (base_obs_time > test_obs_time) return
+
+ if (base_obs_time /= test_obs_time) cycle
+
test_obs_type = get_obs_kind(selection_list(i))
if (base_obs_type /= test_obs_type) cycle
- test_obs_time = get_obs_def_time(selection_list(i))
- if (base_obs_time /= test_obs_time) cycle
-
test_obs_loc = get_obs_def_location(selection_list(i))
if ( .not. horiz_location_equal(base_obs_loc, test_obs_loc)) cycle
@@ -936,8 +1076,7 @@
return
enddo
- ! got to end of selection list without a match, cycle.
- good_selection = .false.
+ ! if you get here, no match, and return value is already false.
end function good_selection
@@ -975,6 +1114,18 @@
end function horiz_location_equal
+!---------------------------------------------------------------------------
+function get_time_from_obs(this_obs)
+ type(obs_type), intent(in) :: this_obs
+ type(time_type) :: get_time_from_obs
+
+type(obs_def_type) :: this_obs_def
+
+call get_obs_def(this_obs, this_obs_def)
+get_time_from_obs = get_obs_def_time(this_obs_def)
+
+end function get_time_from_obs
+
!---------------------------------------------------------------------
end program obs_selection
Modified: DART/branches/development/time_manager/time_manager_mod.f90
===================================================================
--- DART/branches/development/time_manager/time_manager_mod.f90 2012-04-16 16:03:56 UTC (rev 5700)
+++ DART/branches/development/time_manager/time_manager_mod.f90 2012-04-16 22:03:12 UTC (rev 5701)
@@ -50,7 +50,7 @@
! Subroutines and functions operating on time_type
public :: set_time, set_time_missing, increment_time, decrement_time, get_time
-public :: interval_alarm, repeat_alarm, generate_seed
+public :: interval_alarm, repeat_alarm, generate_seed, time_index_sort
! List of available calendar types
public :: THIRTY_DAY_MONTHS, JULIAN, GREGORIAN, NOLEAP, NO_CALENDAR, &
@@ -3158,5 +3158,78 @@
!-------------------------------------------------------------------------
+subroutine time_index_sort(t, index, num)
+
+! Uses a heap sort alogrithm on t, returns array of sorted indices
+! Sorts from earliest (smallest value) time to latest (largest value);
+! uses time > and < operators to do the compare. When index = 1
+! that's the earliest time.
+
+integer, intent(in) :: num
+type(time_type), intent(in) :: t(num)
+integer, intent(out) :: index(num)
+
+integer :: ind, i, j, l_val_index, level
+type(time_type) :: l_val
+
+
+if ( .not. module_initialized ) call time_manager_init
+
+! INITIALIZE THE INDEX ARRAY TO INPUT ORDER
+do i = 1, num
+ index(i) = i
+end do
+
+! Only one element, just send it back
+if(num <= 1) return
+
+level = num / 2 + 1
+ind = num
+
+! Keep looping until finished
+do
+ ! Keep going down levels until bottom
+ if(level > 1) then
+ level = level - 1
+ l_val = t(index(level))
+ l_val_index = index(level)
+ else
+ l_val = t(index(ind))
+ l_val_index = index(ind)
+
+
+ index(ind) = index(1)
+ ind = ind - 1
+ if(ind == 1) then
+ index(1) = l_val_index
+ return
+ endif
+ endif
+
+ i = level
+ j = 2 * level
+
+ do while(j <= ind)
+ if(j < ind) then
+ if(t(index(j)) < t(index(j + 1))) j = j + 1
+ endif
+ if(l_val < t(index(j))) then
+ index(i) = index(j)
+ i = j
+ j = 2 * j
+ else
+ j = ind + 1
+ endif
+
+ end do
+ index(i) = l_val_index
+
+end do
+
+end subroutine time_index_sort
+
+
+!-------------------------------------------------------------------------
+
end module time_manager_mod
More information about the Dart-dev
mailing list