[Dart-dev] [4891] DART/trunk/utilities: Added the option to restrict the 'close' computation to a list of particular
nancy at ucar.edu
nancy at ucar.edu
Wed May 4 13:54:22 MDT 2011
Revision: 4891
Author: nancy
Date: 2011-05-04 13:54:22 -0600 (Wed, 04 May 2011)
Log Message:
-----------
Added the option to restrict the 'close' computation to a list of particular
state vector KINDS (e.g. KIND_TEMPERATURE, KIND_U_WIND_COMPONENT). Use
mask arrays and the 'where' construct to do array operations in the code.
Modified Paths:
--------------
DART/trunk/utilities/closest_member_tool.f90
DART/trunk/utilities/closest_member_tool.html
DART/trunk/utilities/closest_member_tool.nml
-------------- next part --------------
Modified: DART/trunk/utilities/closest_member_tool.f90
===================================================================
--- DART/trunk/utilities/closest_member_tool.f90 2011-05-04 19:50:29 UTC (rev 4890)
+++ DART/trunk/utilities/closest_member_tool.f90 2011-05-04 19:54:22 UTC (rev 4891)
@@ -12,29 +12,32 @@
! Program to overwrite the time on each ensemble in a restart file.
-use types_mod, only : r8
-use time_manager_mod, only : time_type, set_time_missing, &
- operator(==), print_time, print_date, &
- set_calendar_type, GREGORIAN, NO_CALENDAR, &
- get_calendar_type
-
-use utilities_mod, only : register_module, do_output, &
- error_handler, nmlfileunit, E_MSG, E_ERR, &
- timestamp, find_namelist_in_file, &
- check_namelist_read, logfileunit, &
- do_nml_file, do_nml_term, open_file, close_file
+use types_mod, only : r8
+use time_manager_mod, only : time_type, set_time_missing, &
+ operator(/=), print_time
+
+use utilities_mod, only : register_module, do_output, &
+ error_handler, nmlfileunit, E_MSG, E_ERR, &
+ timestamp, find_namelist_in_file, &
+ check_namelist_read, logfileunit, &
+ do_nml_file, do_nml_term, open_file, close_file
-use sort_mod, only : slow_index_sort
+use location_mod, only : location_type
-use assim_model_mod, only : static_init_assim_model, get_model_size, &
- open_restart_read, open_restart_write, &
- awrite_state_restart, aread_state_restart, &
- close_restart
+use obs_kind_mod, only : get_num_raw_obs_kinds, get_raw_obs_kind_index, &
+ paramname_length, get_raw_obs_kind_name
-use mpi_utilities_mod, only : initialize_mpi_utilities, task_count, &
- finalize_mpi_utilities
+use sort_mod, only : slow_index_sort
+use assim_model_mod, only : static_init_assim_model, get_model_size, &
+ open_restart_read, open_restart_write, &
+ awrite_state_restart, aread_state_restart, &
+ close_restart, get_state_meta_data
+use mpi_utilities_mod, only : initialize_mpi_utilities, task_count, &
+ finalize_mpi_utilities
+
+
implicit none
! version controlled file description for error handling, do not edit
@@ -43,17 +46,22 @@
revision = "$Revision$", &
revdate = "$Date$"
-integer :: iunit, model_size, io, ens
+integer :: iunit, model_size, io, ens, i, j, kindindex
integer, allocatable :: index_list(:)
-character(len = 128) :: ifile, msgstring
+character(len = 128) :: ifile, msgstring, msgstring1
real(r8), allocatable :: mean(:), member(:), diffs(:)
+logical :: allkinds, done
+logical, allocatable :: usekind(:), useindex(:)
+type(location_type) :: loc
+integer :: num_kinds, stype
type(time_type) :: mean_time, member_time, mean_advance_time, advance_time
+integer, parameter :: max_list_len = 500
character(len=64) :: method_name(4) = (/ &
"Simple Difference ", &
"Normalized Difference", &
- "Simple RMSE ", &
- "Normalized RMSE " /)
+ "Simple RMS Diff ", &
+ "Normalized RMS Diff " /)
!----------------------------------------------------------------
! These variables are namelist-controllable.
@@ -63,13 +71,14 @@
integer :: ens_size = 1
logical :: single_restart_file_in = .true.
integer :: difference_method = 4
+character(len = paramname_length) :: use_only_kinds(max_list_len) = ''
!----------------------------------------------------------------
! different methods to compute 'distance' from mean:
! 1 = simple absolute difference
! 2 = normalized absolute difference
-! 3 = simple rmse difference
-! 4 = normalized rmse difference
+! 3 = simple rms difference
+! 4 = normalized rms difference
!
! (suggest more...)
!----------------------------------------------------------------
@@ -79,7 +88,8 @@
output_file_name, &
ens_size, &
single_restart_file_in, &
- difference_method
+ difference_method, &
+ use_only_kinds
logical :: input_is_model_advance_file = .false.
@@ -91,6 +101,7 @@
!----------------------------------------------------------------
+! program start
!----------------------------------------------------------------
! This program should only be run with a single process
@@ -102,6 +113,7 @@
call register_module(source,revision,revdate)
+
! Read the namelist entry and print it
call find_namelist_in_file("input.nml", "closest_member_tool_nml", iunit)
read(iunit, nml = closest_member_tool_nml, iostat = io)
@@ -122,10 +134,53 @@
write(msgstring, *) 'Computing difference using method: '//trim(method_name(difference_method))
call error_handler(E_MSG,'',msgstring)
-
! make space for the mean and a single member, plus place to sort list for output
-allocate(mean(model_size), member(model_size), index_list(ens_size), diffs(ens_size))
+allocate(mean(model_size), member(model_size))
+allocate(index_list(ens_size), diffs(ens_size))
+! are we adding up only differences from particular kinds, or the entire vector?
+if (use_only_kinds(1) /= '') then
+ allkinds = .false.
+
+ num_kinds = get_num_raw_obs_kinds()
+ allocate(usekind(num_kinds))
+ usekind = .false.
+
+ done = .false.
+ KindList:do i=1, max_list_len
+ if (use_only_kinds(i) == '') then
+ done = .true.
+ exit KindList
+ endif
+ kindindex = get_raw_obs_kind_index(use_only_kinds(i))
+ if (kindindex < 0) then
+ write(msgstring, *) 'unrecognized KIND string: '//trim(use_only_kinds(i))
+ call error_handler(E_ERR,'closest_member_tool', msgstring, &
+ source,revision,revdate)
+ endif
+ usekind(kindindex) = .true.
+
+ enddo KindList
+
+ if (.not. done) then
+ write(msgstring, *) 'cannot have more than ', max_list_len, ' kinds'
+ call error_handler(E_ERR,'closest_member_tool', msgstring, &
+ source,revision,revdate)
+ endif
+
+ write(msgstring, *) 'Computing difference based only on items in state vector items of kind:'
+ call error_handler(E_MSG,'',msgstring)
+ do i=1, num_kinds
+ if (usekind(i)) then
+ write(msgstring, *) ' ', trim(get_raw_obs_kind_name(i))
+ call error_handler(E_MSG,'',msgstring)
+ endif
+ enddo
+
+else
+ allkinds = .true.
+endif
+
mean_time = set_time_missing()
member_time = set_time_missing()
mean_advance_time = set_time_missing()
@@ -141,7 +196,39 @@
endif
call close_restart(iunit)
+! if we are not processing all kinds of state vector items, set up a mask
+! for the kinds we are. do this once at the start so we don't replicate
+! work. the usekind(number_of_different_kinds) array says whether we are
+! going to add differences for this type. the useindex(state_vector_len)
+! array is being precomputed here so it's fast to loop over the ensemble
+! members and only add up differences for the kinds of interest.
+allocate(useindex(model_size))
+if (.not. allkinds) then
+ useindex(:) = .false.
+ j = 0
+ do i=1, model_size
+ call get_state_meta_data(-1 * i, loc, stype)
+ if (stype < 1 .or. stype > num_kinds) then
+ write(msgstring, *) 'bad KIND from get_state_meta_data, ', stype, ' for index ', i
+ write(msgstring1, *) 'must be between 1 and ', num_kinds
+ call error_handler(E_ERR,'closest_member_tool', msgstring, &
+ source,revision,revdate, text2=msgstring1)
+
+ endif
+ if (usekind(stype)) then
+ useindex(i) = .true.
+ j = j + 1
+ endif
+ enddo
+
+ write(msgstring, *) 'using ', j, ' of ', model_size, ' items in the state vector'
+ call error_handler(E_MSG,'closest_member_tool', msgstring)
+else
+ ! use everything.
+ useindex(:) = .true.
+endif
+
! either loop over individual files or open a single file and read a member
! at a time. same functionality; where the file open/close happens differs.
if (single_restart_file_in) then
@@ -158,14 +245,18 @@
call aread_state_restart(member_time, member, iunit)
endif
- ! FIXME: should put in tests for time matching mean to avoid comparing
- ! against the wrong file.
+ if (mean_time /= member_time) then
+ call print_time(mean_time, "time of ensemble mean data")
+ call print_time(member_time, "time of ensemble member data")
+ write(msgstring, *) 'member ', i, ' has a different timestamp than mean'
+ call error_handler(E_ERR,'closest_member_tool', msgstring)
+ endif
- !------------------- Compute scaled RMSE -----------------------
+ !------------------- Compute difference -----------------------
diffs(ens) = compute_diff(mean, member, model_size)
- !------------------- Compute scaled RMSE -----------------------
+ !------------------- Compute difference -----------------------
enddo
call close_restart(iunit)
@@ -187,11 +278,18 @@
call close_restart(iunit)
!------------------- Read restart from file ----------------------
- !------------------- Compute scaled RMSE -----------------------
+ if (mean_time /= member_time) then
+ call print_time(mean_time, "time of ensemble mean data")
+ call print_time(member_time, "time of ensemble member data")
+ write(msgstring, *) 'member ', i, ' has a different timestamp than mean'
+ call error_handler(E_ERR,'closest_member_tool', msgstring)
+ endif
+ !------------------- Compute difference -----------------------
+
diffs(ens) = compute_diff(mean, member, model_size)
- !------------------- Compute scaled RMSE -----------------------
+ !------------------- Compute difference -----------------------
enddo
@@ -206,7 +304,7 @@
call error_handler(E_MSG, '', ' ')
do ens=1, ens_size
- write(msgstring, "(A,I5,A,F18.6)") "Member ", index_list(ens), " difference ", diffs(index_list(ens))
+ write(msgstring, "(A,I5,A,G18.6)") "Member ", index_list(ens), " difference ", diffs(index_list(ens))
call error_handler(E_MSG, '', msgstring)
enddo
@@ -233,7 +331,8 @@
!------------------- Write results to file -----------------------
-deallocate(mean, member, index_list, diffs)
+deallocate(mean, member, index_list, diffs, useindex)
+if (.not. allkinds) deallocate(usekind)
call finalize_mpi_utilities() ! now closes log file, too
@@ -250,77 +349,63 @@
real(r8) :: compute_diff
integer :: i
-real(r8) :: val, r, diff, biggest
+real(r8) :: val, r, diff, biggest, val2
+real(r8), allocatable :: adiff(:)
+integer :: hist(5)
+character(len=72) :: tbuf
+! new strategy: compute an array of differences and sum them at the end.
+! try to use array operations when possible. useindex() is a logical array
+! that can be used as a mask if only some kinds are going to be used.
+! it is set to all .true. if no kinds were set, so it can be used in
+! either case.
+
+allocate(adiff(arraysize))
+adiff = 0.0_r8
+
select case (difference_method)
! simple absolute difference
case (1)
+ where (useindex) adiff = abs(target - candidate)
+
+ compute_diff = sum(adiff)
- val = 0.0
- do i = 1, arraysize
- val = val + abs(target(i) - candidate(i))
- enddo
-
- compute_diff = val
-
! normalized absolute difference
case (2)
- val = 0.0
- do i = 1, arraysize
- if (target(i) == 0.0) then
- if (candidate(i) == 0.0) then
- diff = 0.0
- else
- diff = candidate(i)
- endif
- else
- diff = abs((target(i) - candidate(i)) / target(i))
- endif
- val = val + diff
- enddo
-
- compute_diff = val
+ where (useindex)
+ where (target /= 0.0_r8)
+ adiff = abs((target - candidate) / target)
+ elsewhere
+ adiff = abs(candidate)
+ endwhere
+ endwhere
- ! simple rmse difference
+ compute_diff = sum(adiff)
+
+ ! simple rms difference
case (3)
- val = 0.0
- do i = 1, arraysize
- diff = target(i) - candidate(i)
- val = val + ( diff * diff )
- enddo
-
- if (val > 0) then
- compute_diff = sqrt(val)
- else
- compute_diff = val
- endif
+ where (useindex) adiff = target - candidate
- ! normalized rmse difference
+ compute_diff = sum(adiff * adiff)
+ if (compute_diff > 0.0_r8) compute_diff = sqrt(compute_diff)
+
+ ! normalized rms difference
case (4)
- val = 0.0
- do i = 1, arraysize
- if (target(i) == 0.0) then
- if (candidate(i) == 0.0) then
- diff = 0.0
- else
- diff = candidate(i)
- endif
- else
- diff = (target(i) - candidate(i)) / target(i)
- endif
- val = val + (diff*diff)
- enddo
-
- if (val > 0) then
- compute_diff = sqrt(val)
- else
- compute_diff = val
- endif
+ where (useindex)
+ where (target /= 0.0_r8)
+ adiff = (target - candidate) / target
+ elsewhere
+ adiff = candidate
+ endwhere
+ endwhere
+ compute_diff = sum(adiff * adiff)
+ if (compute_diff > 0.0_r8) compute_diff = sqrt(compute_diff)
+
case default
write(msgstring, *) 'Valid values for difference_method are 1-4, value is', difference_method
call error_handler(E_ERR,'closest_member_tool','Bad value for difference_method', &
@@ -328,6 +413,9 @@
end select
+deallocate(adiff)
+
end function compute_diff
+
end program closest_member_tool
Modified: DART/trunk/utilities/closest_member_tool.html
===================================================================
--- DART/trunk/utilities/closest_member_tool.html 2011-05-04 19:50:29 UTC (rev 4890)
+++ DART/trunk/utilities/closest_member_tool.html 2011-05-04 19:54:22 UTC (rev 4891)
@@ -46,9 +46,11 @@
The ensemble mean is optionally written out from filter by setting a namelist
option in the <a href="../filter/filter.html#Namelist">filter namelist</a>
to output the ensemble mean in restart file format at the end of a run.
-Restart files contain only the bare state vector, so no information
-about the location or variable type is available. The difference is computed
-point by point across the ensemble members.
+The difference is computed point by point across the ensemble members.
+There is an option to restrict the computation to just a subset of the
+entire state vector by listing one or more generic kinds.
+In this case, only state vector items matching one of these kinds
+will contribute to the total difference value.
</P><P>
Available methods are:
<dl>
@@ -65,13 +67,13 @@
normalized by the mean value,
accumulated over the entire state vector.
</dd>
-<dt>3 - simple RMSE difference:</dt>
+<dt>3 - simple RMS difference:</dt>
<dd>
The square root of the accumulated sum of the
square of the difference between each item in the mean vector
and the corresponding item in each ensemble member.
</dd>
-<dt>4 - normalized RMSE difference:</dt>
+<dt>4 - normalized RMS difference:</dt>
<dd>
The square root of the accumulated sum of the
square of the normalized difference between each item in the mean
@@ -119,7 +121,8 @@
<pre>
<em class=call>namelist / closest_member_tool_nml / </em>
input_file_name, output_file_name,
- ens_size, single_restart_file_in, difference_method
+ ens_size, single_restart_file_in, difference_method,
+ use_only_kinds
</pre>
</div>
@@ -168,11 +171,22 @@
<table>
<tr><td>1 </td><td>simple absolute difference</td></tr>
<tr><td>2</td><td>normalized absolute difference</td></tr>
-<tr><td>3</td><td>simple RMSE difference</td></tr>
-<tr><td>4</td><td>normalized RMSE difference</td></tr>
+<tr><td>3</td><td>simple RMS difference</td></tr>
+<tr><td>4</td><td>normalized RMS difference</td></tr>
</table>
Default 4</TD></TR>
+<TR><!--contents--><TD valign=top>use_only_kinds</TD>
+ <!-- type --><TD valign=top>character(len=32)</TD>
+ <!--descript--><TD valign=top>If unspecified, all items in the state vector
+ contribute to the total difference. If one or more kinds
+ are listed here, only items in the state vector of these kinds
+ contribute to the total difference. These are the generic kinds,
+ such as KIND_TEMPERATURE, KIND_U_WIND_COMPONENT, KIND_DENSITY, etc.
+ and not specific types like RADIOSONDE_TEMPERATURE. Consult the
+ model interface code to determine which possible kinds are returned by the
+ <a href="../models/model_mod.html#get_state_meta_data">get_state_meta_data()</a> routine.
+ Default ''.</TD></TR>
</TABLE>
</div>
Modified: DART/trunk/utilities/closest_member_tool.nml
===================================================================
--- DART/trunk/utilities/closest_member_tool.nml 2011-05-04 19:50:29 UTC (rev 4890)
+++ DART/trunk/utilities/closest_member_tool.nml 2011-05-04 19:54:22 UTC (rev 4891)
@@ -5,12 +5,12 @@
# 3 = simple rmse difference
# 4 = normalized rmse difference
-namelist /closest_member_tool_nml/
+&closest_member_tool_nml
input_file_name = 'filter_restart',
output_file_name = 'closest_restart',
ens_size = 20,
single_restart_file_in = .true.,
difference_method = 4,
+ use_only_kinds = '',
/
-
More information about the Dart-dev
mailing list