[Dart-dev] [4721] DART/trunk/system_simulation: Minor cleanup of full_error. f90 -- ensemble size is
nancy at ucar.edu
nancy at ucar.edu
Tue Feb 15 15:47:43 MST 2011
Revision: 4721
Author: nancy
Date: 2011-02-15 15:47:43 -0700 (Tue, 15 Feb 2011)
Log Message:
-----------
Minor cleanup of full_error.f90 -- ensemble size is
a namelist option, output goes both to stdout and directly
via a formatted write to an output file in exactly the
right format for input to filter if sampling error
correction option enabled. minimal html doc page
describing namelist and expected filenames, and pointer
to submitted paper as reference. add missing input.nml
in work directory for namelists, and a quickbuild to compile
everything in work directory.
Modified Paths:
--------------
DART/trunk/system_simulation/full_error.f90
Added Paths:
-----------
DART/trunk/system_simulation/system_simulation.html
DART/trunk/system_simulation/work/input.nml
DART/trunk/system_simulation/work/quickbuild.csh
-------------- next part --------------
Modified: DART/trunk/system_simulation/full_error.f90
===================================================================
--- DART/trunk/system_simulation/full_error.f90 2011-02-15 22:44:11 UTC (rev 4720)
+++ DART/trunk/system_simulation/full_error.f90 2011-02-15 22:47:43 UTC (rev 4721)
@@ -10,9 +10,21 @@
! $Revision$
! $Date$
-use types_mod, only : r8
+! Correct covariances for fixed ensemble sizes.
+! See Anderson, J. L., 2011: Localization and Sampling Error Correction
+! in Ensemble Kalman Filter Data Assimilation.
+! Submitted for publication, Jan 2011. Contact author.
+
+! This version of the program reads the ensemble size and base filename
+! for the output from a namelist.
+
+use types_mod, only : r8
+use utilities_mod, only : open_file, close_file, error_handler, &
+ find_namelist_in_file, check_namelist_read, &
+ do_nml_file, do_nml_term, nmlfileunit, &
+ initialize_utilities, finalize_utilities, E_ERR
use random_seq_mod, only : random_seq_type, init_random_seq, random_gaussian, &
- twod_gaussians, random_uniform
+ twod_gaussians, random_uniform
implicit none
@@ -22,33 +34,67 @@
revision = "$Revision$", &
revdate = "$Date$"
-integer, parameter :: num_times = 1, num_points = 100000000
-!!!integer, parameter :: num_times = 1, num_points = 10000000
-integer, parameter :: ens_size = 80
+integer, parameter :: num_times = 1
+integer, parameter :: num_samples = 100000000
+
+
+! ---------------
+! namelist items
+
+integer :: ens_size = 80
+character(len=128) :: output_filename = 'full_error'
+
+namelist /full_error_nml/ ens_size, output_filename
+
+
type (random_seq_type) :: ran_id
-real(r8) :: pairs(2, ens_size), temp(ens_size)
+real(r8), allocatable :: pairs(:,:), temp(:)
+
real(r8) :: zero_2(2) = 0.0, cov(2, 2)
real(r8) :: t_correl, correl_mean, sample_correl, alpha, beta
real(r8) :: s_mean(2), s_var(2), reg_mean, reg_sd, t_sd1, t_sd2, true_correl_mean
real(r8) :: tcorrel_sum(201), reg_sum(201), reg_2_sum(201)
-integer i, j, k, bin_num, bin_count(201)
+integer :: i, j, k, bin_num, bin_count(201)
+integer :: iunit, io
+
+character(len=132) :: outname, errstring
+character(len=16) :: formstring
+
+!
+! start of executable code
+!
+
+call initialize_utilities('full_error')
+
+! Read the namelist entry
+call find_namelist_in_file("input.nml", "full_error_nml", iunit, .false.)
+read(iunit, nml = full_error_nml, iostat = io)
+call check_namelist_read(iunit, io, "full_error_nml", .false.)
+
+! Record the namelist values used for the run
+if (do_nml_file()) write(nmlfileunit, nml=full_error_nml)
+if (do_nml_term()) write( * , nml=full_error_nml)
+
+call close_file(iunit)
+
call init_random_seq(ran_id)
-
write(*, *) 'stats for ensemble size ', ens_size
-bin_count = 0
+allocate(pairs(2, ens_size), temp(ens_size))
+
+bin_count = 0
tcorrel_sum = 0.0_r8
-reg_sum = 0.0_r8
-reg_2_sum = 0.0_r8
+reg_sum = 0.0_r8
+reg_2_sum = 0.0_r8
! Uniformly distributed sequence of true correlations
-do j = 1, num_points
+do j = 1, num_samples
! Assume that true correlation is uniform [-1,1]
- t_correl = -1.0_r8 + 2.0_r8 * ((j - 1) / (num_points - 1.0_r8))
- if(j > num_points) t_correl = 0.0_r8
+ t_correl = -1.0_r8 + 2.0_r8 * ((j - 1) / (num_samples - 1.0_r8))
+ if(j > num_samples) t_correl = 0.0_r8
do k = 1, num_times
@@ -56,8 +102,10 @@
t_sd2 = 1.0_r8
! Generate the covariance matrix for this correlation
- cov(1, 1) = t_sd1**2; cov(2, 2) = t_sd2**2
- cov(1, 2) = t_correl * t_sd1 * t_sd2; cov(2, 1) = cov(1, 2)
+ cov(1, 1) = t_sd1**2
+ cov(2, 2) = t_sd2**2
+ cov(1, 2) = t_correl * t_sd1 * t_sd2
+ cov(2, 1) = cov(1, 2)
! Loop to generate an ensemble size sample from this correl
! Generate a random sample of size ens_size from something with this correlation
@@ -68,7 +116,8 @@
! Compute the sample correlation
call comp_correl(pairs, ens_size, sample_correl)
- ! Bin statistics for each 0.01 sample correlation bin; start with just mean true correlation
+ ! Bin statistics for each 0.01 sample correlation bin;
+ ! start with just mean true correlation
bin_num = (sample_correl + 1.0_r8) / 0.01_r8 + 1
! Keep a sum of the sample_correl and squared to compute standard deviation
tcorrel_sum(bin_num) = tcorrel_sum(bin_num) + t_correl
@@ -89,26 +138,66 @@
end do
end do
+! make the size of the integer output .X, .XX, .XXX, etc
+! depending on the number of decimal digits in the value.
+if (ens_size < 10) then
+ formstring = '(2A,I1)'
+else if (ens_size < 100) then
+ formstring = '(2A,I2)'
+else if (ens_size < 1000) then
+ formstring = '(2A,I3)'
+else if (ens_size < 10000) then
+ formstring = '(2A,I4)'
+else if (ens_size < 100000) then
+ formstring = '(2A,I5)'
+else
+ formstring = '(2A,I8)'
+endif
-do i = 1, 201
- if(bin_count(i) > 0) then
- ! Compute the standard deviation of the true correlations
- true_correl_mean = tcorrel_sum(i) / bin_count(i)
- reg_mean = reg_sum(i) / bin_count(i)
- reg_sd = sqrt((reg_2_sum(i) - bin_count(i) * reg_mean**2) / (bin_count(i) - 1))
+! filename
+write(outname, formstring) trim(output_filename), '.', ens_size
- if(reg_sd <= 0.0_r8) then
- alpha = 1.0_r8
- else
- !!!beta = reg_mean**2 / reg_sd**2
- ! CRAZY IDEA???
- beta = reg_mean**2 / (reg_sd**2 * (1.0_r8 + 1.0_r8 / ens_size))
- alpha = beta / (1.0_r8 + beta)
- endif
+! text file, overwrite existing file if present
+iunit = open_file(outname, 'formatted', 'write')
- write(*, *) 'bin, count, mean ', i, bin_count(i), true_correl_mean, alpha
+! always generate exactly 200 entries in the output file
+do i = 1, 200
+ ! must have at least 2 counts to compute a std dev
+ if(bin_count(i) <= 1) then
+ write(errstring, *) 'Bin ', i, ' has ', bin_count(i), ' counts'
+ call error_handler(E_ERR,'full_error', errstring, &
+ source, revision, revdate, text2="All bins must have at least 2 counts")
endif
+
+ ! Compute the standard deviation of the true correlations
+ true_correl_mean = tcorrel_sum(i) / bin_count(i)
+ reg_mean = reg_sum(i) / bin_count(i)
+ reg_sd = sqrt((reg_2_sum(i) - bin_count(i) * reg_mean**2) / (bin_count(i) - 1))
+
+ if(reg_sd <= 0.0_r8) then
+ alpha = 1.0_r8
+ else
+ !!!beta = reg_mean**2 / reg_sd**2
+ ! Correct for bias in the standard deviation for very small ensembles, too
+ beta = reg_mean**2 / (reg_sd**2 * (1.0_r8 + 1.0_r8 / ens_size))
+ alpha = beta / (1.0_r8 + beta)
+ endif
+
+ write(*, *) 'bin, count, mean ', i, bin_count(i), true_correl_mean, alpha
+ write(iunit, 10) i, bin_count(i), true_correl_mean, alpha
end do
+
+10 format (I4,I9,2G25.16)
+
+call close_file(iunit)
+
+! print out just to stdout how many counts, if any, fell beyond bin 200
+write(*, *) 'bin, count, mean ', 201, bin_count(201), 0, 0
+
+deallocate(pairs, temp)
+
+call finalize_utilities()
+
end program full_error
Added: DART/trunk/system_simulation/system_simulation.html
===================================================================
--- DART/trunk/system_simulation/system_simulation.html (rev 0)
+++ DART/trunk/system_simulation/system_simulation.html 2011-02-15 22:47:43 UTC (rev 4721)
@@ -0,0 +1,219 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN"
+ "http://www.w3.org/TR/html4/strict.dtd">
+<HTML>
+<HEAD>
+<TITLE>system simulation programs</TITLE>
+<link rel="stylesheet" type="text/css" href="../doc/html/doc.css">
+</HEAD>
+<BODY>
+<A NAME="TOP"></A>
+
+<H1>PROGRAM <em class=program>full_error and others</em></H1>
+
+<table border=0 summary="" cellpadding=5>
+<tr>
+ <td valign=middle>
+ <img src="../doc/html/Dartboard9.png" alt="DART project logo" height=70 />
+ </td>
+ <td>
+ <P>Jump to <a href="../index.html">DART Documentation Main Index</a><br />
+ <small><small>version information for this file: <br />
+ <!-- version tag follows, do not edit -->
+ $Id: filter.html 4522 2010-10-13 17:05:42Z thoar $</small></small>
+ </P></td>
+</tr>
+</table>
+
+<A HREF="#Namelist">NAMELIST</A> /
+<A HREF="#Modules">MODULES</A> /
+<A HREF="#FilesUsed">FILES</A> /
+<A HREF="#References">REFERENCES</A> /
+<A HREF="#Errors">ERRORS</A> /
+<A HREF="#FuturePlans">PLANS</A> /
+<A HREF="#Legalese">TERMS OF USE</A>
+
+<H2>Overview</H2>
+
+<P>
+A collection of standalone programs for simulating
+various properties of ensembles. The program of
+most interest here is <em class=file>full_error.f90</em>
+which generates the auxiliary information needed when
+using sampling error correction in the filter.
+</P>
+<P>
+Talk to Jeff Anderson about the other programs in
+this directory.
+</P>
+
+<!--==================================================================-->
+<!--============== DESCRIPTION OF A NAMELIST ========================-->
+<!--==================================================================-->
+
+<A NAME="Namelist"></A>
+<HR>
+<H2>NAMELIST</H2>
+<P>
+<A HREF="#Namelist"> <em class=code>&full_error_nml</em> </A>
+is always read from file <em class=file>input.nml</em>.
+</P>
+<P>We adhere to the F90 standard of starting a namelist with an ampersand
+'&' and terminating with a slash '/' for all our namelist input.
+Character strings that contain a '/' must be
+enclosed in quotes to prevent them from prematurely terminating the namelist.
+The namelist declaration is:
+</P>
+<div class=namelist>
+<pre>
+<em class=call>namelist / full_error_nml / </em>
+ ens_size,
+ output_filename
+</pre>
+</div>
+
+<div class=indent1>
+<!-- Description -->
+
+<TABLE border=0 cellpadding=3 width=100%>
+<TR><TH align=left>Contents </TH>
+ <TH align=left>Type </TH>
+ <TH align=left>Description </TH></TR>
+
+<TR><!--contents--><TD valign=top>ens_size</TD>
+ <!-- type --><TD valign=top>integer</TD>
+ <!--descript--><TD>Number of ensemble members to compute
+ sampling error for. This computation is independent
+ of any of the model details; it only depends on the
+ count of ensemble members.
+ Default: 80</TD></TR>
+
+<TR><!--contents--><TD valign=top>output_filename</TD>
+ <!-- type --><TD valign=top>character(len=128)</TD>
+ <!--descript--><TD>Basename for the output file written
+ by this program. The program will append
+ <em class=code>.N</em> to the name, where N is the
+ ensemble size, creating filenames
+ of the form "base.8", "base.32", etc. The base name
+ should be "full_error" to match the code in the filter
+ program which reads this file.
+ Default: "full_error"</TD></TR>
+
+</TABLE>
+
+</div>
+<br />
+
+
+
+
+
+</P>
+
+<!--==================================================================-->
+
+<A NAME="Modules"></A>
+<HR>
+<H2>MODULES USED</H2>
+<PRE>
+types_mod
+utilities_mod
+random_seq_mod
+</PRE>
+
+<!--==================================================================-->
+<!-- Describe the Files Used by this module. -->
+<!--==================================================================-->
+
+<A NAME="FilesUsed"></A>
+<HR>
+<H2>FILES</H2>
+<P>
+The full_error program creates files which should be called
+<em class=file>full_error.N</em> to match the filenames expected
+by the filter program.
+</P>
+
+<!--==================================================================-->
+<!-- Cite references, if need be. -->
+<!--==================================================================-->
+
+<A NAME="References"></A>
+<HR>
+<H2>REFERENCES</H2>
+<ul>
+<li>Anderson, J. L., 2011:
+Localization and Sampling Error Correction
+in Ensemble Kalman Filter Data Assimilation.
+Submitted for publication, Jan 2011.
+Contact author.</li>
+
+</ul>
+<br />
+
+<!--==================================================================-->
+<!-- Describe all the error conditions and codes. -->
+<!--==================================================================-->
+
+<A NAME="Errors"></A>
+<HR>
+<H2>ERROR CODES and CONDITIONS</H2>
+<div class="errors">
+<TABLE border=1 cellspacing=1 cellpadding=10 width=100%>
+<TR><TH>Routine</TH><TH>Message</TH><TH>Comment</TH></TR>
+
+<TR><!-- routine --><TD VALIGN=top>full_error</TD>
+ <!-- message --><TD VALIGN=top>cannot handle task counts > 99999</TD>
+ <!-- comment --><TD VALIGN=top>Ensemble size must be less than 100,000. </TD></TR>
+
+<TR><!-- routine --><TD VALIGN=top>full_error</TD>
+ <!-- message --><TD VALIGN=top>empty bin</TD>
+ <!-- comment --><TD VALIGN=top>The sample size must be large enough for
+ all bins to have counts</TD></TR>
+
+</TABLE>
+</div>
+
+<H2>KNOWN BUGS</H2>
+<P>
+none
+</P>
+
+<!--==================================================================-->
+<!-- Descibe Future Plans. -->
+<!--==================================================================-->
+
+<A NAME="FuturePlans"></A>
+<HR>
+<H2>FUTURE PLANS</H2>
+<P>
+none at this time.
+</P>
+
+<!--==================================================================-->
+<!-- Legalese & Metadata -->
+<!--==================================================================-->
+
+<A NAME="Legalese"></A>
+<HR>
+<H2>Terms of Use</H2>
+
+<P>
+DART software - Copyright © 2004 - 2010 UCAR.<br>
+This open source software is provided by UCAR, "as is",<br>
+without charge, subject to all terms of use at<br>
+<a href="http://www.image.ucar.edu/DAReS/DART/DART_download">
+http://www.image.ucar.edu/DAReS/DART/DART_download</a>
+</P>
+
+<TABLE border=0 cellpadding=0 width=100% summary="">
+<TR><TD valign=top>Contact: </TD><TD> Jeff Anderson </TD></TR>
+<TR><TD valign=top>Revision: </TD><TD> $Revision: 4222 $ </TD></TR>
+<TR><TD valign=top>Source: </TD><TD> $URL: https://proxy.subversion.ucar.edu/DAReS/DART/trunk/filter/wakeup_filter.html $ </TD></TR>
+<TR><TD valign=top>Change Date: </TD><TD> $Date: 2010-01-20 11:27:41 -0700 (Wed, 20 Jan 2010) $ </TD></TR>
+<TR><TD valign=top>Change history: </TD><TD> try "svn log" or "svn diff" </TD></TR>
+</TABLE>
+
+<!--==================================================================-->
+
+</BODY>
+</HTML>
Added: DART/trunk/system_simulation/work/input.nml
===================================================================
--- DART/trunk/system_simulation/work/input.nml (rev 0)
+++ DART/trunk/system_simulation/work/input.nml 2011-02-15 22:47:43 UTC (rev 4721)
@@ -0,0 +1,7 @@
+&utilities_nml
+/
+
+&full_error_nml
+ ens_size = 10,
+ output_filename = "final_full",
+/
Added: DART/trunk/system_simulation/work/quickbuild.csh
===================================================================
--- DART/trunk/system_simulation/work/quickbuild.csh (rev 0)
+++ DART/trunk/system_simulation/work/quickbuild.csh 2011-02-15 22:47:43 UTC (rev 4721)
@@ -0,0 +1,44 @@
+#!/bin/csh
+#
+# DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+# provided by UCAR, "as is", without charge, subject to all terms of use at
+# http://www.image.ucar.edu/DAReS/DART/DART_download
+#
+# $Id: quickbuild.csh 4504 2010-09-24 21:05:09Z nancy $
+#
+# This script compiles all executables in this directory.
+
+set MODEL = "system_simulation"
+
+@ n = 0
+
+#----------------------------------------------------------------------
+# Build all the single-threaded targets
+#----------------------------------------------------------------------
+
+foreach TARGET ( mkmf_* )
+
+ set PROG = `echo $TARGET | sed -e 's#mkmf_##'`
+
+ @ n = $n + 1
+ echo
+ echo "---------------------------------------------------"
+ echo "${MODEL} build number ${n} is ${PROG}"
+ \rm -f ${PROG}
+ csh $TARGET || exit $n
+ make || exit $n
+
+end
+
+# clean up. comment this out if you want to keep the .o and .mod files around
+\rm -f *.o *.mod input.nml.*_default
+
+echo "Success: All DART programs compiled."
+
+exit 0
+
+# <next few lines under version control, do not edit>
+# $URL: https://proxy.subversion.ucar.edu/DAReS/DART/trunk/observations/gps/work/quickbuild.csh $
+# $Revision: 4504 $
+# $Date: 2010-09-24 15:05:09 -0600 (Fri, 24 Sep 2010) $
+
Property changes on: DART/trunk/system_simulation/work/quickbuild.csh
___________________________________________________________________
Added: svn:executable
+ *
More information about the Dart-dev
mailing list