[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 &nbsp;</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