[Dart-dev] [3953] DART/trunk/perfect_model_obs/perfect_model_obs.f90: Allow ' output_interval' to be set to 0 without a core dump.
nancy at ucar.edu
nancy at ucar.edu
Tue Jun 30 13:51:29 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090630/93c2d8b8/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/perfect_model_obs/perfect_model_obs.f90
===================================================================
--- DART/trunk/perfect_model_obs/perfect_model_obs.f90 2009-06-30 19:46:50 UTC (rev 3952)
+++ DART/trunk/perfect_model_obs/perfect_model_obs.f90 2009-06-30 19:51:29 UTC (rev 3953)
@@ -25,7 +25,7 @@
write_obs_seq, get_num_obs, init_obs, assignment(=), &
static_init_obs_sequence, get_num_qc, read_obs_seq_header, &
set_qc_meta_data, get_expected_obs, delete_seq_head, &
- delete_seq_tail, set_obs_def, destroy_obs
+ delete_seq_tail, set_obs_def, destroy_obs, destroy_obs_sequence
use obs_def_mod, only : obs_def_type, get_obs_def_error_variance, &
set_obs_def_error_variance, get_obs_def_time
@@ -77,6 +77,8 @@
integer :: obs_window_seconds = -1
integer :: tasks_per_model_advance = 1
integer :: output_interval = 1
+integer :: print_every_nth_obs = 0
+
character(len = 129) :: restart_in_file_name = 'perfect_ics', &
restart_out_file_name = 'perfect_restart', &
obs_seq_in_file_name = 'obs_seq.in', &
@@ -91,8 +93,10 @@
obs_seq_in_file_name, obs_seq_out_file_name, &
adv_ens_command, tasks_per_model_advance, &
obs_window_days, obs_window_seconds, silence, &
- trace_execution, output_timestamps
+ trace_execution, output_timestamps, &
+ print_every_nth_obs
+
!------------------------------------------------------------------------------
! Doing this allows independent scoping for subroutines in main program file
@@ -117,7 +121,7 @@
integer :: j, iunit, time_step_number, obs_seq_file_id
integer :: cnum_copies, cnum_qc, cnum_obs, cnum_max
integer :: additional_qc, additional_copies
-integer :: ierr, io, istatus, num_obs_in_set
+integer :: ierr, io, istatus, num_obs_in_set, nth_obs
integer :: model_size, key_bounds(2), num_qc, last_key_used
real(r8) :: true_obs(1), obs_value(1), qc(1), errvar
@@ -140,17 +144,25 @@
if (do_nml_file()) write(nmlfileunit, nml=perfect_model_obs_nml)
if (do_nml_term()) write( * , nml=perfect_model_obs_nml)
+! set the level of output
+call set_trace(trace_execution, output_timestamps, silence)
+
+call trace_message('Perfect_model start')
+call timestamp_message('Perfect_model start')
+
! Don't let this run with more than one task; just a waste of resource
if(task_count() > 1) then
write(msgstring, *) 'Only use one mpi process here: ', task_count(), ' were requested'
call error_handler(E_ERR, 'perfect_main', msgstring, &
source, revision, revdate)
endif
-call task_sync()
-! set the level of output
-call set_trace(trace_execution, output_timestamps, silence)
+! Default to printing nothing
+nth_obs = -1
+call trace_message('Before setting up space for observations')
+call timestamp_message('Before observation setup')
+
! Find out how many data copies are in the obs_sequence
call read_obs_seq_header(obs_seq_in_file_name, cnum_copies, cnum_qc, cnum_obs, cnum_max, &
obs_seq_file_id, obs_seq_read_format, pre_I_format, close_the_file = .true.)
@@ -185,17 +197,24 @@
call set_copy_meta_data(seq, 1, copy_meta_data(1))
call set_copy_meta_data(seq, 2, copy_meta_data(2))
+call timestamp_message('After observation setup')
+call trace_message('After setting up space for observations')
+
! Initialize the model now that obs_sequence is all set up
model_size = get_model_size()
write(msgstring,*)'Model size = ',model_size
call error_handler(E_MSG,'perfect_main',msgstring)
! Set up the ensemble storage and read in the restart file
+call trace_message('Before reading in ensemble restart file')
call perfect_read_restart(ens_handle, model_size)
+call trace_message('After reading in ensemble restart file')
+call trace_message('Before initializing output diagnostic file')
state_meta(1) = 'true state'
! Set up output of truth for state
StateUnit = init_diag_output('True_State', 'true state from control', 1, state_meta)
+call trace_message('After initializing output diagnostic file')
! Initialize a repeatable random sequence for perturbations
call init_random_seq(random_seq)
@@ -208,6 +227,8 @@
write(msgstring, *) 'number of qc values is ',num_qc
call error_handler(E_MSG,'perfect_main',msgstring)
+call trace_message('Before trimming obs seq if start/stop time specified')
+
! Need to find first obs with appropriate time, delete all earlier ones
if(first_obs_seconds > 0 .or. first_obs_days > 0) then
first_obs_time = set_time(first_obs_seconds, first_obs_days)
@@ -230,6 +251,8 @@
endif
endif
+call trace_message('After trimming obs seq if start/stop time specified')
+
! Time step number is used to do periodic diagnostic output
time_step_number = -1
window_time = set_time(0,0)
@@ -239,56 +262,93 @@
time_step_number = time_step_number + 1
write(msgstring , '(A,I5)') 'Main evaluation loop, starting iteration', time_step_number
- if (.not.silence) call error_handler(E_MSG,'', ' ')
- if (.not.silence) call error_handler(E_MSG,'perfect_model_obs:', msgstring)
+ call trace_message(' ', ' ', -1)
+ call trace_message(msgstring, 'perfect_model_obs: ', -1)
+ call trace_message('Before move_ahead checks time of data and next obs')
+
! Get the model to a good time to use a next set of observations
call move_ahead(ens_handle, 1, seq, last_key_used, window_time, &
key_bounds, num_obs_in_set, curr_ens_time, next_ens_time)
if(key_bounds(1) < 0) then
- call error_handler(E_MSG,'perfect_model_obs:', 'No more obs to evaluate, exiting main loop')
+ call trace_message('No more obs to evaluate, exiting main loop', 'perfect_model_obs:', -1)
exit AdvanceTime
endif
+ call trace_message('After move_ahead checks time of data and next obs')
+
if (curr_ens_time /= next_ens_time) then
- if (.not.silence) call error_handler(E_MSG,'perfect_model_obs:', 'Ready to run model to advance data time')
+ call trace_message('Ready to run model to advance data time', 'perfect_model_obs:', -1)
+ call print_ens_time(ens_handle, 'Ensemble data time before advance')
+ call trace_message('Before advance_state called to run model')
+ call timestamp_message('Before model advance', sync=.true.)
+
call advance_state(ens_handle, 1, next_ens_time, async, &
adv_ens_command, tasks_per_model_advance)
+
+ call timestamp_message('After model advance', sync=.true.)
+ call trace_message('After advance_state called to run model')
+ call print_ens_time(ens_handle, 'Ensemble data time after advance')
else
- if (.not.silence) call error_handler(E_MSG,'perfect_model_obs:', 'Model does not need to run; data already at required time')
+ call trace_message('Model does not need to run; data already at required time', 'perfect_model_obs:', -1)
endif
+ call trace_message('Before setup for next group of observations')
+ write(msgstring, '(A,I7)') 'Number of observations to be evaluated', &
+ num_obs_in_set
+ call trace_message(msgstring)
+ call print_obs_time(seq, key_bounds(1), 'Time of first observation in window')
+ call print_obs_time(seq, key_bounds(2), 'Time of last observation in window')
+
! Allocate storage for observation keys for this part of sequence
allocate(keys(num_obs_in_set))
! Get all the keys associated with this set of observations
call get_time_range_keys(seq, key_bounds, num_obs_in_set, keys)
+ call trace_message('After setup for next group of observations')
+
! Output the true state to the netcdf file
- if(time_step_number / output_interval * output_interval == time_step_number) &
+ if((output_interval > 0) .and. &
+ (time_step_number / output_interval * output_interval == time_step_number)) then
+
+ call trace_message('Before updating truth diagnostics file')
call aoutput_diagnostics(StateUnit, ens_handle%time(1), ens_handle%vars(:, 1), 1)
+ call trace_message('After updating truth diagnostics file')
+ endif
write(msgstring, '(A,I8,A)') 'Ready to evaluate up to', size(keys), ' observations'
- if (.not.silence) call error_handler(E_MSG,'perfect_model_obs:', msgstring)
+ call trace_message(msgstring, 'perfect_model_obs:', -1)
- ! How many observations in this set
- write(msgstring, *) 'num_obs_in_set is ', num_obs_in_set
- call error_handler(E_DBG,'perfect_main',msgstring,source,revision,revdate)
-
! Compute the forward observation operator for each observation in set
do j = 1, num_obs_in_set
+
+ ! Some compilers do not like mod by 0, so test first.
+ if (print_every_nth_obs > 0) nth_obs = mod(j, print_every_nth_obs)
+
+ ! If requested, print out a message every Nth observation
+ ! to indicate progress is being made and to allow estimates
+ ! of how long the assim will take.
+ if (nth_obs == 0) then
+ write(msgstring, '(A,1x,I8,1x,A,I8)') 'Processing observation ', j, &
+ ' of ', num_obs_in_set
+ call trace_message(msgstring, 'perfect_model_obs:', -1)
+ ! or if you want timestamps:
+ ! call timestamp(msgstring, pos="debug")
+ endif
+
! Compute the observations from the state
call get_expected_obs(seq, keys(j:j), ens_handle%vars(:, 1), &
true_obs(1:1), istatus, assimilate_this_ob, evaluate_this_ob)
- ! If observation is not being evaluated or assimilated, skip it
- ! Ends up setting a 1000 qc field so observation is not used again.
! Get the observational error covariance (diagonal at present)
! Generate the synthetic observations by adding in error samples
call get_obs_from_key(seq, keys(j), obs)
call get_obs_def(obs, obs_def)
+ ! If observation is not being evaluated or assimilated, skip it
+ ! Ends up setting a 1000 qc field so observation is not used again.
if(istatus == 0 .and. (assimilate_this_ob .or. evaluate_this_ob)) then
! DEBUG: try this out to see if it's useful. if the incoming error
! variance is negative, add no noise to the values, but do switch the
@@ -329,24 +389,41 @@
end do AdvanceTime
-call error_handler(E_MSG,'perfect_model_obs:', 'End of main filter evaluation loop, starting cleanup')
+call trace_message('End of main evaluation loop, starting cleanup', 'perfect_model_obs:', -1)
! properly dispose of the diagnostics files
+call trace_message('Before finalizing diagnostics file')
ierr = finalize_diag_output(StateUnit)
+call trace_message('After finalizing diagnostics file')
! Write out the sequence
+call trace_message('Before writing output sequence file')
call write_obs_seq(seq, obs_seq_out_file_name)
+call trace_message('After writing output sequence file')
! Output a restart file if requested
-if(output_restart) &
+if(output_restart) then
+ call trace_message('Before writing state restart file')
call write_ensemble_restart(ens_handle, restart_out_file_name, 1, 1, &
force_single_file = .true.)
+ call trace_message('After writing state restart file')
+endif
+call trace_message('Before ensemble and obs memory cleanup')
+
! Release storage for ensemble
call end_ensemble_manager(ens_handle)
-call error_handler(E_MSG,'perfect_main','FINISHED',source,revision,revdate)
+! Free up the observation kind and obs sequence
+call destroy_obs(obs)
+call destroy_obs_sequence(seq)
+call trace_message('After ensemble and obs memory cleanup')
+call trace_message('Perfect_model done')
+call timestamp_message('Perfect_model done')
+
+!call error_handler(E_MSG,'perfect_main','FINISHED',source,revision,revdate)
+
! closes the log file.
call timestamp(string1=source,string2=revision,string3=revdate,pos='end')
@@ -476,7 +553,7 @@
if (sync) call task_sync()
endif
-call timestamp(trim(msg), pos='debug')
+call timestamp(trim(msg), pos='brief')
end subroutine timestamp_message
@@ -495,8 +572,8 @@
if (get_my_num_copies(ens_handle) < 1) return
call get_ensemble_time(ens_handle, 1, mtime)
-call print_time(mtime, msg, logfileunit)
-call print_time(mtime, msg)
+call print_time(mtime, ' p_m_o trace: '//msg, logfileunit)
+call print_time(mtime, ' p_m_o trace: '//msg)
end subroutine print_ens_time
@@ -519,8 +596,8 @@
call get_obs_from_key(seq, key, obs)
call get_obs_def(obs, obs_def)
mtime = get_obs_def_time(obs_def)
-call print_time(mtime, msg, logfileunit)
-call print_time(mtime, msg)
+call print_time(mtime, ' p_m_o trace: '//msg, logfileunit)
+call print_time(mtime, ' p_m_o trace: '//msg)
call destroy_obs(obs)
end subroutine print_obs_time
More information about the Dart-dev
mailing list