[Dart-dev] [6199] DART/branches/development/adaptive_inflate: fix the log messages to be accurate for inflation types
nancy at ucar.edu
nancy at ucar.edu
Thu May 30 11:57:32 MDT 2013
Revision: 6199
Author: nancy
Date: 2013-05-30 11:57:31 -0600 (Thu, 30 May 2013)
Log Message:
-----------
fix the log messages to be accurate for inflation types
1 and 3; use open_file() and close_file() routines instead
of directly calling open() and close(); remove unused routines
from the 'use' blocks; add error checking to get and set routines
that only apply to obs-space inflation; change one test to
be <= 0 instead of < 0 because it results in the same behavior
with less computation. update the docs to point out that for
inflation type 3 only the first index value will be used for
all items in the state vector.
Modified Paths:
--------------
DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90
DART/branches/development/adaptive_inflate/adaptive_inflate_mod.html
-------------- next part --------------
Modified: DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90
===================================================================
--- DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90 2013-05-30 16:45:12 UTC (rev 6198)
+++ DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90 2013-05-30 17:57:31 UTC (rev 6199)
@@ -12,14 +12,14 @@
!
! Operations and storage required for various adaptive inflation algorithms
-use types_mod, only : r8, PI
+use types_mod, only : r8, PI, missing_r8
use time_manager_mod, only : time_type, get_time, set_time
-use utilities_mod, only : file_exist, get_unit, register_module, &
+use utilities_mod, only : register_module, open_file, close_file, &
error_handler, E_ERR, E_MSG
use random_seq_mod, only : random_seq_type, random_gaussian, init_random_seq
-use ensemble_manager_mod, only : ensemble_type, all_vars_to_all_copies, all_copies_to_all_vars, &
- read_ensemble_restart, write_ensemble_restart, get_copy_owner_index, &
- prepare_to_write_to_vars, prepare_to_read_from_vars, prepare_to_update_vars, map_pe_to_task
+use ensemble_manager_mod, only : ensemble_type, read_ensemble_restart, write_ensemble_restart, &
+ get_copy_owner_index, prepare_to_write_to_vars, &
+ prepare_to_read_from_vars, prepare_to_update_vars, map_pe_to_task
use mpi_utilities_mod, only : my_task_id, send_to, receive_from
implicit none
@@ -58,7 +58,7 @@
end type adaptive_inflate_type
! Module storage for writing error messages
-character(len = 129) :: errstring
+character(len = 255) :: msgstring
! Flag indicating whether module has been initialized
logical :: initialized = .false.
@@ -100,11 +100,13 @@
if(.not. initialized) then
initialized = .true.
call register_module(source, revision, revdate)
+endif
- ! If non-deterministic inflation is being done, need to initialize random sequence.
- ! use the task id number (plus 1 since they start at 0) to set the initial seed.
- ! NOTE: non-deterministic inflation does NOT reproduce as process count is varied!
- if(.not. deterministic) call init_random_seq(inflate_handle%ran_seq, my_task_id()+1)
+! If non-deterministic inflation is being done, need to initialize random sequence.
+! use the task id number (plus 1 since they start at 0) to set the initial seed.
+! NOTE: non-deterministic inflation does NOT reproduce as process count is varied!
+if(.not. deterministic) then
+ call init_random_seq(inflate_handle%ran_seq, my_task_id()+1)
endif
! more information for users to document what they selected in the nml:
@@ -156,14 +158,14 @@
! Cannot support non-determistic inflation and an inf_lower_bound < 1
if(.not. deterministic .and. inf_lower_bound < 1.0_r8) then
- write(errstring, *) 'Cannot have non-deterministic inflation and inf_lower_bound < 1'
- call error_handler(E_ERR, 'adaptive_inflate_init', errstring, source, revision, revdate)
+ write(msgstring, *) 'Cannot have non-deterministic inflation and inf_lower_bound < 1'
+ call error_handler(E_ERR, 'adaptive_inflate_init', msgstring, source, revision, revdate)
endif
-! give this an initial value which defaults to time-constant inflation
-! change it below if the value(s) are > 0
-minmax_mean(:) = 0.0_r8
-minmax_sd(:) = 0.0_r8
+! give these distinctive values; if inflation is being used
+! (e.g. inf_flavor > 0) then they should be set in all cases.
+minmax_mean(:) = missing_r8
+minmax_sd(:) = missing_r8
!------ Block for state space inflation initialization ------
@@ -174,16 +176,12 @@
! Verify that indices are contiguous
if(ss_inflate_sd_index /= ss_inflate_index + 1) then
- write(errstring, *) 'ss_inflate_index = ', ss_inflate_index, &
+ write(msgstring, *) 'ss_inflate_index = ', ss_inflate_index, &
' and ss_inflate_sd_index = ', ss_inflate_sd_index, ' must be continguous'
call error_handler(E_ERR, 'adaptive_inflate_init', &
- errstring, source, revision, revdate)
+ msgstring, source, revision, revdate)
endif
- ! set this, and then below if sd values are from a file, compute them
- ! and send the value to PE0 if not already there.
- minmax_sd = sd_initial
-
! Read in initial values from file OR get from namelist arguments
! If either mean, sd, or both are to be read from the restart file, read them in.
@@ -220,6 +218,25 @@
endif
endif
+ ! Inflation type 3 is spatially-constant. Make sure the entire array is set to that
+ ! value. the computation only uses index 1, but the diagnostics write out the entire
+ ! array and it will be misleading if not constant. the inf values were set above.
+ ! if they were set by namelist, this code changes nothing. but if they were read in
+ ! from a file, then it is possible the values vary across the array. these lines
+ ! ensure the entire array contains a single constant value to match what the code uses.
+ if(inf_flavor == 3) then
+ call get_copy_owner_index(ss_inflate_index, owner, owners_index)
+ if (owner == ens_handle%my_pe) then
+ call prepare_to_update_vars(ens_handle)
+ ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
+ endif
+ call get_copy_owner_index(ss_inflate_sd_index, owner, owners_index)
+ if (owner == ens_handle%my_pe) then
+ call prepare_to_update_vars(ens_handle)
+ ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
+ endif
+ endif
+
! this block figures out what the min/max value of the mean/sd is
! if we are reading in the values from a restart file. it is used
! in diagnostic output so it needs to get to PE0. we also could check it
@@ -271,24 +288,6 @@
endif
endif
- ! Inflation type 3 is spatially-constant. Make sure the entire array is set to that
- ! value. the computation only uses index 1, but the diagnostics write out the entire
- ! array and it will be misleading if not constant. the inf values were set above.
- ! if they were set by namelist, this code changes nothing. but if they were read in
- ! from a file, then it is possible the values vary across the array. these lines
- ! ensure the entire array contains a single constant value to match what the code uses.
- if(inf_flavor == 3) then
- call get_copy_owner_index(ss_inflate_index, owner, owners_index)
- if (owner == ens_handle%my_pe) then
- call prepare_to_update_vars(ens_handle)
- ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
- endif
- call get_copy_owner_index(ss_inflate_sd_index, owner, owners_index)
- if (owner == ens_handle%my_pe) then
- call prepare_to_update_vars(ens_handle)
- ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
- endif
- endif
!------ Block for obs. space inflation initialization ------
@@ -299,28 +298,22 @@
! Only single values for inflation, inflation_sd (not arrays)
if(mean_from_restart .or. sd_from_restart) then
! Open the file
- restart_unit = get_unit()
- open(unit = restart_unit, file = in_file_name, action = 'read', &
- form = 'formatted', iostat = io)
- if (io /= 0) then
- write(errstring, *) 'unable to open inflation restart file ', &
- trim(in_file_name), ' for reading'
- call error_handler(E_ERR, 'adaptive_inflate_init', &
- errstring, source, revision, revdate)
- endif
+ restart_unit = open_file(in_file_name, form='formatted', action='read')
read(restart_unit, *, iostat = io) inflate_handle%inflate, inflate_handle%sd
if (io /= 0) then
- write(errstring, *) 'unable to read inflation restart values from ', &
+ write(msgstring, *) 'unable to read inflation restart values from ', &
trim(in_file_name)
call error_handler(E_ERR, 'adaptive_inflate_init', &
- errstring, source, revision, revdate)
+ msgstring, source, revision, revdate)
endif
- close(restart_unit)
+ call close_file(restart_unit)
endif
! If using the namelist values, set (or overwrite) them here.
if (.not. mean_from_restart) inflate_handle%inflate = inf_initial
if (.not. sd_from_restart) inflate_handle%sd = sd_initial
+ minmax_mean(:) = inflate_handle%inflate
+ minmax_sd(:) = inflate_handle%sd
endif
! Write to log file what kind of inflation is being used.
@@ -359,25 +352,40 @@
sadapt = ' spatially-constant,'
akind = ' state-space'
case default
- write(errstring, *) 'Illegal inflation value for ', label
- call error_handler(E_ERR, 'adaptive_inflate_init', errstring, source, revision, revdate)
+ write(msgstring, *) 'Illegal inflation value for ', label
+ call error_handler(E_ERR, 'adaptive_inflate_init', msgstring, source, revision, revdate)
end select
! say in plain english what kind of inflation was selected.
-write(errstring, '(4A)') &
+write(msgstring, '(4A)') &
trim(det), trim(tadapt), trim(sadapt), trim(akind)
-call error_handler(E_MSG, trim(label) // ' inflation:', errstring, source, revision, revdate)
+call error_handler(E_MSG, trim(label) // ' inflation:', msgstring, source, revision, revdate)
+! if inflation was set from a namelist, that message has already been
+! printed (further up in this routine). for values set from a restart
+! file, if the inflation flavor is 2, the values printed are the min and
+! max from the entire array. for flavors 1 and 3 there is only a single
+! value to print out.
if (inf_flavor > 0) then
if (mean_from_restart) then
- write(errstring, '(A, F8.3, A, F8.3)') &
- 'inf mean from restart file: min value: ', minmax_mean(1), ' max value: ', minmax_mean(2)
- call error_handler(E_MSG, trim(label) // ' inflation:', errstring, source, revision, revdate)
+ if (inf_flavor == 2) then
+ write(msgstring, '(A, F8.3, A, F8.3)') &
+ 'inf mean from restart file: min value: ', minmax_mean(1), ' max value: ', minmax_mean(2)
+ else
+ write(msgstring, '(A, F8.3, A, F8.3)') &
+ 'inf mean from restart file: value: ', minmax_mean(1)
+ endif
+ call error_handler(E_MSG, trim(label) // ' inflation:', msgstring, source, revision, revdate)
endif
if (sd_from_restart) then
- write(errstring, '(A, F8.3, A, F8.3)') &
- 'inf stddev from restart file: min value: ', minmax_sd(1), ' max value: ', minmax_sd(2)
- call error_handler(E_MSG, trim(label) // ' inflation:', errstring, source, revision, revdate)
+ if (inf_flavor == 2) then
+ write(msgstring, '(A, F8.3, A, F8.3)') &
+ 'inf stddev from restart file: min value: ', minmax_sd(1), ' max value: ', minmax_sd(2)
+ else
+ write(msgstring, '(A, F8.3, A, F8.3)') &
+ 'inf stddev from restart file: value: ', minmax_sd(1)
+ endif
+ call error_handler(E_MSG, trim(label) // ' inflation:', msgstring, source, revision, revdate)
endif
endif
@@ -399,10 +407,10 @@
if(do_varying_ss_inflate(inflate_handle) .or. do_single_ss_inflate(inflate_handle)) then
! Verify that indices are contiguous
if(ss_inflate_sd_index /= ss_inflate_index + 1) then
- write(errstring, *) 'ss_inflate_index = ', ss_inflate_index, &
+ write(msgstring, *) 'ss_inflate_index = ', ss_inflate_index, &
' and ss_inflate_sd_index = ', ss_inflate_sd_index, ' must be continguous'
call error_handler(E_ERR, 'adaptive_inflate_end', &
- errstring, source, revision, revdate)
+ msgstring, source, revision, revdate)
endif
! Write the inflate and inflate_sd as two copies for a restart
@@ -412,28 +420,21 @@
! Flavor 1 is observation space, write its restart directly
else if(do_obs_inflate(inflate_handle)) then
! Open the restart file
- restart_unit = get_unit()
- open(unit = restart_unit, file = inflate_handle%out_file_name, &
- action = 'write', form = 'formatted', iostat = io)
- if (io /= 0) then
- write(errstring, *) 'unable to open inflation restart file ', &
- trim(inflate_handle%out_file_name), ' for writing'
- call error_handler(E_ERR, 'adaptive_inflate_end', &
- errstring, source, revision, revdate)
- endif
+ restart_unit = open_file(inflate_handle%out_file_name, &
+ form = 'formatted', action='write')
write(restart_unit, *, iostat = io) inflate_handle%inflate, inflate_handle%sd
if (io /= 0) then
- write(errstring, *) 'unable to write into inflation restart file ', &
+ write(msgstring, *) 'unable to write into inflation restart file ', &
trim(inflate_handle%out_file_name)
call error_handler(E_ERR, 'adaptive_inflate_end', &
- errstring, source, revision, revdate)
+ msgstring, source, revision, revdate)
endif
- close(unit = restart_unit)
+ call close_file(restart_unit)
endif
endif
! Need to close diagnostic files for observation space if in use
-if(inflate_handle%obs_diag_unit > -1) close(inflate_handle%obs_diag_unit)
+if(inflate_handle%obs_diag_unit > -1) call close_file(inflate_handle%obs_diag_unit)
end subroutine adaptive_inflate_end
@@ -501,7 +502,12 @@
real(r8) :: get_inflate
type(adaptive_inflate_type), intent(in) :: inflate_handle
-get_inflate = inflate_handle%inflate
+if(do_obs_inflate(inflate_handle)) then
+ get_inflate = inflate_handle%inflate
+else
+ write(msgstring, *) 'This routine can only be used with obs_space inflation'
+ call error_handler(E_ERR, 'get_inflate', msgstring, source, revision, revdate)
+endif
end function get_inflate
@@ -515,7 +521,12 @@
real(r8) :: get_sd
type(adaptive_inflate_type), intent(in) :: inflate_handle
-get_sd = inflate_handle%sd
+if(do_obs_inflate(inflate_handle)) then
+ get_sd = inflate_handle%sd
+else
+ write(msgstring, *) 'This routine can only be used with obs_space inflation'
+ call error_handler(E_ERR, 'get_sd', msgstring, source, revision, revdate)
+endif
end function get_sd
@@ -523,12 +534,18 @@
subroutine set_inflate(inflate_handle, inflate)
-! Sets the single inflation value in the type
+! The single real value inflate contains the obs_space inflation value
+! when obs_space inflation is in use and this sets it.
type(adaptive_inflate_type), intent(inout) :: inflate_handle
real(r8), intent(in) :: inflate
-inflate_handle%inflate = inflate
+if(do_obs_inflate(inflate_handle)) then
+ inflate_handle%inflate = inflate
+else
+ write(msgstring, *) 'This routine can only be used with obs_space inflation'
+ call error_handler(E_ERR, 'set_inflate', msgstring, source, revision, revdate)
+endif
end subroutine set_inflate
@@ -536,12 +553,18 @@
subroutine set_sd(inflate_handle, sd)
-! Sets the single sd value in the type
+! The single real value inflate_sd contains the obs_space inflate_sd value
+! when obs_space inflation is in use and this sets it.
type(adaptive_inflate_type), intent(inout) :: inflate_handle
real(r8), intent(in) :: sd
-inflate_handle%sd = sd
+if(do_obs_inflate(inflate_handle)) then
+ inflate_handle%sd = sd
+else
+ write(msgstring, *) 'This routine can only be used with obs_space inflation'
+ call error_handler(E_ERR, 'set_sd', msgstring, source, revision, revdate)
+endif
end subroutine set_sd
@@ -563,9 +586,8 @@
! If unit is -1, it hasn't been opened yet, do it.
if(inflate_handle%obs_diag_unit == -1) then
! Open the file
- inflate_handle%obs_diag_unit = get_unit()
- open(unit = inflate_handle%obs_diag_unit, file = inflate_handle%diag_file_name, &
- action = 'write', form = 'formatted')
+ inflate_handle%obs_diag_unit = open_file(inflate_handle%diag_file_name, &
+ form='formatted', action='write')
endif
! Get the time in days and seconds
@@ -573,6 +595,10 @@
! Write out the time followed by the values
write(inflate_handle%obs_diag_unit, *) days, seconds, inflate_handle%inflate, &
inflate_handle%sd
+
+ ! This code intentionally does not close the file handle because
+ ! this routine will be called multiple times. It is closed in the
+ ! adaptive_inflate_end routine.
endif
end subroutine output_inflate_diagnostics
@@ -616,17 +642,19 @@
! on the final prior/posterior increment pairs to avoid large regression
! error as per stochastic filter algorithms. This might help to avoid
! problems with generating gravity waves in the Bgrid model, for instance.
+
+ ! The following code does not do the sort.
- ! Figure out required sd for random noise being added
- ! Don't allow covariance deflation in this version
- if(inflate > 1.0_r8) then
- rand_sd = sqrt(inflate*var - var)
- ! Add random sample from this noise into the ensemble
- do i = 1, ens_size
- ens(i) = random_gaussian(inflate_handle%ran_seq, ens(i), rand_sd)
- end do
- ! Adjust the mean back to the original value
- ens = ens - (sum(ens) / ens_size - mean)
+ ! Figure out required sd for random noise being added
+ ! Don't allow covariance deflation in this version
+ if(inflate > 1.0_r8) then
+ rand_sd = sqrt(inflate*var - var)
+ ! Add random sample from this noise into the ensemble
+ do i = 1, ens_size
+ ens(i) = random_gaussian(inflate_handle%ran_seq, ens(i), rand_sd)
+ end do
+ ! Adjust the mean back to the original value
+ ens = ens - (sum(ens) / ens_size - mean)
endif
endif
@@ -651,8 +679,8 @@
real(r8) :: new_inflate, new_inflate_sd
-! If the inflate_sd is negative, just keep everything the same
-if(inflate_sd < 0.0_r8) return
+! If the inflate_sd not positive, keep everything the same
+if(inflate_sd <= 0.0_r8) return
! A lower bound on the updated inflation sd and an upper bound
! on the inflation itself are provided in the inflate_handle.
@@ -706,6 +734,10 @@
!!!if(gamma > 0.99_r8) then
if(gamma > 1.01_r8) then
+! FIXME: rectify the comment in the following line (about gamma
+! and 1.0) and the test above, which clearly used to be true for
+! 1.0 but is not anymore.
+
! The solution of the cubic below only works if gamma is 1.0
! Can analytically find the maximum of the product: d/dlambda is a
! cubic polynomial in lambda**2; solve using cubic formula for real root
Modified: DART/branches/development/adaptive_inflate/adaptive_inflate_mod.html
===================================================================
--- DART/branches/development/adaptive_inflate/adaptive_inflate_mod.html 2013-05-30 16:45:12 UTC (rev 6198)
+++ DART/branches/development/adaptive_inflate/adaptive_inflate_mod.html 2013-05-30 17:57:31 UTC (rev 6199)
@@ -51,11 +51,20 @@
<a href="http://www.image.ucar.edu/DAReS">IMAGe/DAReS</a> web page document
the algorithms in detail. The <em class=file>DART/tutorial/section12</em>
chapter has more information.
-<BR>
-<BR>
+</P> <P>
Details on controlling the inflation options are contained in the
documentation for the filter. The filter_nml controls what inflation
options are used.
+</P> <P>
+Inflation flavor 3 (spatially-constant state space) reads and writes
+a restart file that is the full size of the state vector, however it
+takes the first value in the array and replicates that throughout the
+array. This allows one to switch between flavors 2 and 3. Going from
+inflation flavor 3 to 2 the initial value for all items in the state
+vector will be a constant value and will then start to adapt.
+Going from inflation flavor 2 to 3 whatever value is in the array at
+index 1 will be replicated and used for the entire rest of the state
+vector items.
</P>
<!--==================================================================-->
More information about the Dart-dev
mailing list