[Dart-dev] [4742] DART/trunk: Allow localization distances to be specified by-type.
nancy at ucar.edu
nancy at ucar.edu
Tue Feb 22 13:50:20 MST 2011
Revision: 4742
Author: nancy
Date: 2011-02-22 13:50:20 -0700 (Tue, 22 Feb 2011)
Log Message:
-----------
Allow localization distances to be specified by-type.
In the wrf model_mod, new option to localize in units of Scale Height,
and compute vertical height as geometric and not geopotential meters.
Details of the WRF model_mod changes include:
- heights in meters are now computed as geometric height and not geopotential
height. to check backwards compatibility the original geopotential heights
can be computed by setting 'use_geopotential_height' to true.
- add scale_height option and conversions between other
vertical units. thanks to ryan torn for this code.
- rename the 'vert_interpolate()' subroutine to 'vert_convert()'
since it does no interpolation but does do conversion.
- rewrite the vert_convert routine to use nested select
statements instead of if/then/else's. add comments, and
return sooner in the routine if no conversion is needed.
- change 'vert_coord' in the derived type to 'localization_coord'
since that is more descriptive of what that variable is used for.
- add a real perturbation routine but leave it commented out by default.
see the comments in the code if you want to use it.
if you have no observations using height in meters (e.g. the vertical coord
is in pressure) and if you localize in pressure, this version should give
identical results to the previous version.
the changes to assim_tools and the location files are all for the
'different localization distances by-type' option. see the html file
for details on how to use it.
Modified Paths:
--------------
DART/trunk/assim_tools/assim_tools_mod.f90
DART/trunk/assim_tools/assim_tools_mod.html
DART/trunk/assim_tools/assim_tools_mod.nml
DART/trunk/location/annulus/location_mod.f90
DART/trunk/location/column/location_mod.f90
DART/trunk/location/oned/location_mod.f90
DART/trunk/location/oned/location_mod.html
DART/trunk/location/threed_sphere/location_mod.f90
DART/trunk/location/threed_sphere/location_mod.html
DART/trunk/location/threed_sphere/location_mod.nml
DART/trunk/location/threed_sphere/test/path_names_location_test
DART/trunk/location/threed_sphere/test/test.in
DART/trunk/location/twod/location_mod.f90
DART/trunk/location/twod_sphere/location_mod.f90
DART/trunk/models/wrf/model_mod.f90
DART/trunk/models/wrf/work/input.nml
-------------- next part --------------
Modified: DART/trunk/assim_tools/assim_tools_mod.f90
===================================================================
--- DART/trunk/assim_tools/assim_tools_mod.f90 2011-02-22 20:25:29 UTC (rev 4741)
+++ DART/trunk/assim_tools/assim_tools_mod.f90 2011-02-22 20:50:20 UTC (rev 4742)
@@ -12,7 +12,7 @@
! A variety of operations required by assimilation.
-use types_mod, only : r8, digits12, PI
+use types_mod, only : r8, digits12, PI, missing_r8
use utilities_mod, only : file_exist, get_unit, check_namelist_read, do_output, &
find_namelist_in_file, register_module, error_handler, &
E_ERR, E_MSG, nmlfileunit, do_nml_file, do_nml_term, &
@@ -28,6 +28,9 @@
use obs_def_mod, only : obs_def_type, get_obs_def_location, get_obs_def_time, &
get_obs_def_error_variance, get_obs_kind
+use obs_kind_mod, only : get_num_obs_kinds, get_obs_kind_index
+
+
use cov_cutoff_mod, only : comp_cov_factor
use reg_factor_mod, only : comp_reg_factor
@@ -63,14 +66,17 @@
! Indicates if module initialization subroutine has been called yet
logical :: module_initialized = .false.
+
integer :: print_timestamps = 0
integer :: print_trace_details = 0
-
! True if random sequence needs to be initialized
logical :: first_inc_ran_call = .true.
type (random_seq_type) :: inc_ran_seq
+integer :: num_types = 0
+real(r8), allocatable :: cutoff_list(:)
+logical :: has_special_cutoffs
character(len = 129) :: errstring
! Need to read in table for off-line based sampling correction and store it
@@ -93,27 +99,43 @@
! 4 = particle filter
! 5 = random draw from posterior
! 6 = deterministic draw from posterior with fixed kurtosis
+! 8 = Rank Histogram Filter (see Anderson 2011)
+!
+! special_localization_obs_types -> Special treatment for the specified observation types
+! special_localization_cutoffs -> Different cutoff value for each specified obs type
+!
integer :: filter_kind = 1
real(r8) :: cutoff = 0.2_r8
logical :: sort_obs_inc = .false.
logical :: spread_restoration = .false.
logical :: sampling_error_correction = .false.
integer :: adaptive_localization_threshold = -1
+real(r8) :: adaptive_cutoff_floor = 0.0_r8
integer :: print_every_nth_obs = 0
-logical :: output_localization_diagnostics = .false.
+
+! since this is in the namelist, it has to have a fixed size.
+integer, parameter :: MAX_ITEMS = 300
+character(len = 129) :: special_localization_obs_types(MAX_ITEMS)
+real(r8) :: special_localization_cutoffs(MAX_ITEMS)
+
+logical :: output_localization_diagnostics = .false.
character(len = 129) :: localization_diagnostics_file = "localization_diagnostics"
+
! Following only relevant for filter_kind = 8
logical :: rectangular_quadrature = .true.
logical :: gaussian_likelihood_tails = .false.
+
! Not in the namelist; this var disables the experimental
! linear and spherical case code in the adaptive localization
! sections. to try out the alternatives, set this to .false.
logical :: only_area_adapt = .true.
namelist / assim_tools_nml / filter_kind, cutoff, sort_obs_inc, &
- spread_restoration, sampling_error_correction, adaptive_localization_threshold, &
+ spread_restoration, sampling_error_correction, &
+ adaptive_localization_threshold, adaptive_cutoff_floor, &
print_every_nth_obs, rectangular_quadrature, gaussian_likelihood_tails, &
- output_localization_diagnostics, localization_diagnostics_file
+ output_localization_diagnostics, localization_diagnostics_file, &
+ special_localization_obs_types, special_localization_cutoffs
!============================================================================
@@ -123,10 +145,17 @@
subroutine assim_tools_init()
-integer :: iunit, io
+integer :: iunit, io, i, j
+integer :: num_special_cutoff, type_index
call register_module(source, revision, revdate)
+! give these guys initial values at runtime *before* we read
+! in the namelist. this is to help detect how many items are
+! actually given in the namelist.
+special_localization_obs_types(:) = 'null'
+special_localization_cutoffs(:) = missing_r8
+
! Read the namelist entry
call find_namelist_in_file("input.nml", "assim_tools_nml", iunit)
read(iunit, nml = assim_tools_nml, iostat = io)
@@ -138,22 +167,81 @@
! FOR NOW, can only do spread restoration with filter option 1 (need to extend this)
if(spread_restoration .and. .not. filter_kind == 1) then
- write(errstring, *) 'cant combine spread_restoration and filter_kind ', filter_kind
+ write(errstring, *) 'cannot combine spread_restoration and filter_kind ', filter_kind
call error_handler(E_ERR,'assim_tools_init:', errstring, source, revision, revdate)
endif
+! allocate a list in all cases - even the ones where there is only
+! a single cutoff value. note that in spite of the name these
+! are specific types (e.g. RADIOSONDE_TEMPERATURE, AIRCRAFT_TEMPERATURE)
+! because that's what get_close() is passed.
+num_types = get_num_obs_kinds()
+allocate(cutoff_list(num_types))
+cutoff_list(:) = cutoff
+has_special_cutoffs = .false.
+
+! Go through special-treatment observation kinds, if any.
+num_special_cutoff = 0
+j = 0
+do i = 1, MAX_ITEMS
+ if(special_localization_obs_types(i) == 'null') exit
+ if(special_localization_cutoffs(i) == MISSING_R8) then
+ write(errstring, *) 'cutoff value', i, ' is uninitialized.'
+ call error_handler(E_ERR,'assim_tools_init:', &
+ 'special cutoff namelist for types and distances do not match', &
+ source, revision, revdate, &
+ text2='kind = '//trim(special_localization_obs_types(i)), &
+ text3=trim(errstring))
+ endif
+ j = j + 1
+enddo
+num_special_cutoff = j
+
+if (num_special_cutoff > 0) has_special_cutoffs = .true.
+
+do i = 1, num_special_cutoff
+ type_index = get_obs_kind_index(special_localization_obs_types(i))
+ if (type_index < 0) then
+ write(errstring, *) 'unrecognized TYPE_ in the special localization namelist:'
+ call error_handler(E_ERR,'assim_tools_init:', errstring, source, revision, revdate, &
+ text2=trim(special_localization_obs_types(i)))
+ endif
+ cutoff_list(type_index) = special_localization_cutoffs(i)
+end do
+
if (do_output()) then
write(errstring, '(A,F18.6)') 'The cutoff namelist value is ', cutoff
call error_handler(E_MSG,'assim_tools_init:', errstring)
- write(errstring, '(A,F18.6)') 'cutoff = localization half-width parameter, so'
+ write(errstring, '(A)') 'cutoff is the localization half-width parameter,'
call error_handler(E_MSG,'assim_tools_init:', errstring)
- write(errstring, '(A,F18.6)') 'the effective localization radius is ', cutoff*2.0_r8
+ write(errstring, '(A,F18.6)') 'so the effective localization radius is ', cutoff*2.0_r8
call error_handler(E_MSG,'assim_tools_init:', errstring)
+ if (has_special_cutoffs) then
+ call error_handler(E_MSG, '', '')
+ call error_handler(E_MSG,'assim_tools_init:','Observations with special localization treatment:')
+ call error_handler(E_MSG,'assim_tools_init:','(type name, specified cutoff distance, effective localization radius)')
+
+ do i = 1, num_special_cutoff
+ write(errstring, '(A32,F18.6,F18.6)') special_localization_obs_types(i), &
+ special_localization_cutoffs(i), special_localization_cutoffs(i)*2.0_r8
+ call error_handler(E_MSG,'assim_tools_init:', errstring)
+ end do
+ call error_handler(E_MSG,'assim_tools_init:','all other observation types will use the default cutoff distance')
+ call error_handler(E_MSG, '', '')
+ endif
+
+
if(adaptive_localization_threshold > 0) then
write(errstring, '(A,I10,A)') 'Using adaptive localization, threshold ', &
adaptive_localization_threshold, ' obs'
call error_handler(E_MSG,'assim_tools_init:', errstring)
+ if(adaptive_cutoff_floor > 0) then
+ write(errstring, '(A,I10,A)') 'Minimum cutoff will not go below ', &
+ adaptive_cutoff_floor
+ call error_handler(E_MSG,'assim_tools_init:', 'Using adaptive localization cutoff floor.', &
+ text2=errstring)
+ endif
endif
if(output_localization_diagnostics) then
@@ -199,7 +287,7 @@
real(r8) :: last_close_obs_dist(obs_ens_handle%my_num_vars)
real(r8) :: last_close_state_dist(ens_handle%my_num_vars)
-integer :: my_num_obs, i, j, owner, owners_index, my_num_state
+integer :: my_num_obs, i, j, ispc, owner, owners_index, my_num_state
integer :: my_obs(obs_ens_handle%my_num_vars), my_state(ens_handle%my_num_vars)
integer :: this_obs_key, obs_mean_index, obs_var_index
integer :: grp_beg(num_groups), grp_end(num_groups), grp_size, grp_bot, grp_top, group
@@ -231,7 +319,6 @@
logical :: local_obs_inflate
logical :: get_close_buffering
-
! Initialize assim_tools_module if needed
if(.not. module_initialized) then
call assim_tools_init
@@ -321,14 +408,30 @@
end do
endif
+! the get_close_maxdist_init() call is quick - it just allocates space.
+! it's the get_close_obs_init() call that takes time. correct me if
+! i'm confused but the obs list is different each time, so it can't
+! be skipped. i guess the state space obs_init is redundant and we
+! could just do it once the first time - but knowing when to
+! deallocate the space is a question.
+
! NOTE THESE COULD ONLY BE DONE ONCE PER RUN!!! FIGURE THIS OUT.
! The computations in the two get_close_maxdist_init are redundant
+
! Initialize the method for getting state variables close to a given ob on my process
-call get_close_maxdist_init(gc_state, 2.0_r8*cutoff)
+if (has_special_cutoffs) then
+ call get_close_maxdist_init(gc_state, 2.0_r8*cutoff, 2.0_r8*cutoff_list)
+else
+ call get_close_maxdist_init(gc_state, 2.0_r8*cutoff)
+endif
call get_close_obs_init(gc_state, my_num_state, my_state_loc)
! Initialize the method for getting obs close to a given ob on my process
-call get_close_maxdist_init(gc_obs, 2.0_r8*cutoff)
+if (has_special_cutoffs) then
+ call get_close_maxdist_init(gc_obs, 2.0_r8*cutoff, 2.0_r8*cutoff_list)
+else
+ call get_close_maxdist_init(gc_obs, 2.0_r8*cutoff)
+endif
call get_close_obs_init(gc_obs, my_num_obs, my_obs_loc)
if (get_close_buffering) then
@@ -502,20 +605,23 @@
endif
endif
+ ! set the cutoff default
+ cutoff_rev = cutoff_list(base_obs_kind)
+
! For adaptive localization, need number of other obs close to the chosen observation
- cutoff_rev = cutoff
if(adaptive_localization_threshold > 0) then
! this does a cross-task sum, so all tasks must make this call.
total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_kind, &
- close_obs_dist, cutoff*2.0_r8)
+ close_obs_dist, cutoff_rev*2.0_r8)
! Want expected number of close observations to be reduced to some threshold;
! accomplish this by cutting the size of the cutoff distance.
if(total_num_close_obs > adaptive_localization_threshold) then
- cutoff_rev = revised_distance(cutoff*2.0_r8, adaptive_localization_threshold, &
- total_num_close_obs, base_obs_loc) / 2.0_r8
+ cutoff_rev = revised_distance(cutoff_rev*2.0_r8, adaptive_localization_threshold, &
+ total_num_close_obs, base_obs_loc, &
+ adaptive_cutoff_floor*2.0_r8) / 2.0_r8
if ( output_localization_diagnostics ) then
@@ -545,7 +651,8 @@
call write_location(-1, base_obs_loc, charstring=base_loc_text)
write(localization_unit,'(i8,1x,i5,1x,i8,1x,A,2(f14.5,1x,i10))') i, secs, days, &
- trim(base_loc_text), cutoff, total_num_close_obs, cutoff_rev, rev_num_close_obs
+ trim(base_loc_text), cutoff_list(base_obs_kind), total_num_close_obs, &
+ cutoff_rev, rev_num_close_obs
endif
endif
@@ -557,7 +664,7 @@
! this does a cross-task sum, so all tasks must make this call.
total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_kind, &
- close_obs_dist, cutoff*2.0_r8)
+ close_obs_dist, cutoff_rev*2.0_r8)
if (my_task_id() == 0) then
call get_obs_def(observation, obs_def)
@@ -566,7 +673,7 @@
call write_location(-1, base_obs_loc, charstring=base_loc_text)
write(localization_unit,'(i8,1x,i5,1x,i8,1x,A,f14.5,1x,i10)') i, secs, days, &
- trim(base_loc_text), cutoff, total_num_close_obs
+ trim(base_loc_text), cutoff_rev, total_num_close_obs
endif
endif
@@ -1414,8 +1521,8 @@
write(correction_file_name, 41) 'final_full.', ens_size
else
write(errstring,*)'Trying to use ',ens_size,' model states -- too many.'
- call error_handler(E_MSG,'get_correction_from_file',errstring,source,revision,revdate)
- call error_handler(E_ERR,'get_correction_from_file','Use less than 10000 ensemble.',source,revision,revdate)
+ call error_handler(E_ERR,'get_correction_from_file','Use less than 10000 ens members.',&
+ source,revision,revdate, text2=errstring)
11 format(a11, i1)
21 format(a11, i2)
@@ -2178,10 +2285,11 @@
!--------------------------------------------------------------------
-function revised_distance(orig_dist, newcount, oldcount, base)
+function revised_distance(orig_dist, newcount, oldcount, base, cutfloor)
real(r8), intent(in) :: orig_dist
integer, intent(in) :: newcount, oldcount
type(location_type), intent(in) :: base
+ real(r8), intent(in) :: cutfloor
real(r8) :: revised_distance
@@ -2195,6 +2303,10 @@
if (only_area_adapt) then
revised_distance = orig_dist * sqrt(real(newcount, r8) / oldcount)
+
+ ! allow user to set a minimum cutoff, so even if there are very dense
+ ! observations the cutoff distance won't go below this floor.
+ if (revised_distance < cutfloor) revised_distance = cutfloor
return
endif
@@ -2241,6 +2353,10 @@
source, revision, revdate)
endif
+! allow user to set a minimum cutoff, so even if there are very dense
+! observations the cutoff distance won't go below this floor.
+if (revised_distance < cutfloor) revised_distance = cutfloor
+
end function revised_distance
!--------------------------------------------------------------------
Modified: DART/trunk/assim_tools/assim_tools_mod.html
===================================================================
--- DART/trunk/assim_tools/assim_tools_mod.html 2011-02-22 20:25:29 UTC (rev 4741)
+++ DART/trunk/assim_tools/assim_tools_mod.html 2011-02-22 20:50:20 UTC (rev 4742)
@@ -4,6 +4,7 @@
<HEAD>
<TITLE>module assim_tools_mod</TITLE>
<link rel="stylesheet" type="text/css" href="../doc/html/doc.css" />
+<link href="../doc/html/dart.ico" rel="shortcut icon" />
</HEAD>
<BODY>
<A NAME="TOP"></A>
@@ -58,7 +59,7 @@
<li> 5 = Random draw from posterior (talk to Jeff before using)
<li> 6 = Deterministic draw from posterior with fixed kurtosis (ditto)
<li> 7 = Boxcar kernel filter
-<li> 8 = Rank histogram filter (see Anderson 2011)
+<li> 8 = Rank histogram filter (see Anderson 2010)
</ul>
Most users use type 1, the EAKF.
</P>
@@ -221,7 +222,7 @@
the following namelist items allow this.
<br /> <br />
Optional list of observation types (e.g. "RADAR_REFLECTIVITY",
-"AIRS_RADIANCE") which will use a different cutoff distance.
+"AIRS_TEMPERATURE") which will use a different cutoff distance.
Any observation types not listed here will use the standard cutoff
distance (set by the 'cutoff' namelist value). This is only
implemented for the threed_sphere location module (the one
@@ -577,24 +578,47 @@
<A NAME="References"></A>
<HR>
<H2>REFERENCES</H2>
-<OL>
-<LI>Anderson, Jeffrey L. "A Local Least Squares Framework for Ensemble Filtering"
-April 2003, MWR 131, pp 634-642,
-doi: 10.1175/1520-0493(2003)131<0634:ALLSFF>2.0.CO;2</LI>
-<LI> Anderson, J., Collins, N., 2007,
-"Scalable Implementations of Ensemble Filter Algorithms for Data Assimilation"
-Journal of Atmospheric and Oceanic Technology, 24, 1452-1463,
-doi: 10.1175/JTECH2049.1</LI>
-<LI>Anderson, Jeffrey L. "A Non-Gaussian Ensemble Filter Update for Data Assimilation"
-November 2010, MWR 139, pp 4186-4198,
-DOI: 10.1175/2010MWR3253.1</LI>
-<LI>Anderson, J. L. (2011),
-"Localization and Sampling Error Correction
-in Ensemble Kalman Filter Data Assimilation"
-Submitted for publication, Jan 2011.
+<UL>
+<li>Anderson, J. L., 2001:
+An Ensemble Adjustment Kalman Filter for Data Assimilation.
+<span style="font-style: italic;">Mon. Wea. Rev.</span>,
+<span style="font-weight: bold;">129</span>, 2884-2903.<br />
+<a href="http://dx.doi.org/10.1175/1520-0493(2001)129<2884:AEAKFF>2.0.CO;2"
+target="_blank" >
+doi: 10.1175/1520-0493(2001)129<2884:AEAKFF>2.0.CO;2</a> </li>
+<br />
+<li>Anderson, J. L., 2003:
+A Local Least Squares Framework for Ensemble Filtering.
+<span style="font-style: italic;">Mon. Wea. Rev.</span>,
+<span style="font-weight: bold;">131</span>, 634-642.<br />
+<a href="http://dx.doi.org/10.1175/1520-0493(2003)131<0634:ALLSFF>2.0.CO;2"
+target="_blank" >
+doi: 10.1175/1520-0493(2003)131<0634:ALLSFF>2.0.CO;2</a></li>
+<br />
+<li>Anderson, J., Collins, N., 2007:
+Scalable Implementations of Ensemble Filter Algorithms for Data Assimilation.
+<span style="font-style: italic;">Journal of Atmospheric and Oceanic Technology</span>,
+<span style="font-weight: bold;">24</span>, 1452-1463.<br />
+<a href="http://dx.doi.org/10.1175/JTECH2049.1"
+target="_blank" >doi: 10.1175/JTECH2049.1</a></li>
+<br />
+<li>Anderson, J. L., 2010:
+A Non-Gaussian Ensemble Filter Update for Data Assimilation.
+<span style="font-style: italic;">Mon. Wea. Rev.</span>,
+<span style="font-weight: bold;">139</span>, 4186-4198.<br />
+<a href="http://dx.doi.org/10.1175/2010MWR3253.1"
+target="_blank" >
+doi: 10.1175/2010MWR3253.1</a></lI>
+<br />
+<LI>Anderson, J. L., 2011:,
+Localization and Sampling Error Correction
+in Ensemble Kalman Filter Data Assimilation.
+<span style="font-style: italic;">Submitted for publication</span>, Jan 2011.
Contact author.</LI>
-</OL>
+</UL>
+<br />
+
<!--==================================================================-->
<!-- Describe all the error conditions and codes. -->
<!--==================================================================-->
@@ -612,8 +636,8 @@
</TR>
<TR><!-- routine --><TD VALIGN=top>obs_increment</TD>
- <!-- message --><TD VALIGN=top>Illegal value of filter_kind in assim_tools namelist [1-4 OK]</TD>
- <!-- comment --><TD VALIGN=top>Only 1-6 currently supported.</TD>
+ <!-- message --><TD VALIGN=top>Illegal value of filter_kind in assim_tools namelist [1-8 OK]</TD>
+ <!-- comment --><TD VALIGN=top>Only 1-4,8 currently supported.</TD>
</TR>
<TR><!-- routine --><TD VALIGN=top>obs_increment_eakf <BR>
Modified: DART/trunk/assim_tools/assim_tools_mod.nml
===================================================================
--- DART/trunk/assim_tools/assim_tools_mod.nml 2011-02-22 20:25:29 UTC (rev 4741)
+++ DART/trunk/assim_tools/assim_tools_mod.nml 2011-02-22 20:50:20 UTC (rev 4742)
@@ -5,9 +5,15 @@
spread_restoration = .false.,
sampling_error_correction = .false.,
adaptive_localization_threshold = -1,
+ adaptive_cutoff_floor = -1.0,
output_localization_diagnostics = .false.,
localization_diagnostics_file = "localization_diagnostics",
print_every_nth_obs = 0,
rectangular_quadrature = .true.,
- gaussian_likelihood_tails = .false. /
+ gaussian_likelihood_tails = .false.,
+/
+# specify these in the same order, the same number of items
+# special_localization_obs_types = "",
+# special_localization_cutoffs = -1,
+
Modified: DART/trunk/location/annulus/location_mod.f90
===================================================================
--- DART/trunk/location/annulus/location_mod.f90 2011-02-22 20:25:29 UTC (rev 4741)
+++ DART/trunk/location/annulus/location_mod.f90 2011-02-22 20:50:20 UTC (rev 4742)
@@ -716,10 +716,11 @@
!----------------------------------------------------------------------------
-subroutine get_close_maxdist_init(gc, maxdist)
+subroutine get_close_maxdist_init(gc, maxdist, maxdist_list)
type(get_close_type), intent(inout) :: gc
real(r8), intent(in) :: maxdist
+real(r8), intent(in), optional :: maxdist_list(:)
! Set the maximum distance in the structure
gc%maxdist = maxdist
Modified: DART/trunk/location/column/location_mod.f90
===================================================================
--- DART/trunk/location/column/location_mod.f90 2011-02-22 20:25:29 UTC (rev 4741)
+++ DART/trunk/location/column/location_mod.f90 2011-02-22 20:50:20 UTC (rev 4742)
@@ -574,10 +574,11 @@
!----------------------------------------------------------------------------
-subroutine get_close_maxdist_init(gc, maxdist)
+subroutine get_close_maxdist_init(gc, maxdist, maxdist_list)
type(get_close_type), intent(inout) :: gc
real(r8), intent(in) :: maxdist
+real(r8), intent(in), optional :: maxdist_list(:)
! Set the maximum distance in the structure
gc%maxdist = maxdist
Modified: DART/trunk/location/oned/location_mod.f90
===================================================================
--- DART/trunk/location/oned/location_mod.f90 2011-02-22 20:25:29 UTC (rev 4741)
+++ DART/trunk/location/oned/location_mod.f90 2011-02-22 20:50:20 UTC (rev 4742)
@@ -513,11 +513,13 @@
!----------------------------------------------------------------------------
-subroutine get_close_maxdist_init(gc, maxdist)
+subroutine get_close_maxdist_init(gc, maxdist, maxdist_list)
type(get_close_type), intent(inout) :: gc
real(r8), intent(in) :: maxdist
+real(r8), intent(in), optional :: maxdist_list(:)
+
! Set the maximum distance in the structure
gc%maxdist = maxdist
Modified: DART/trunk/location/oned/location_mod.html
===================================================================
--- DART/trunk/location/oned/location_mod.html 2011-02-22 20:25:29 UTC (rev 4741)
+++ DART/trunk/location/oned/location_mod.html 2011-02-22 20:50:20 UTC (rev 4742)
@@ -4,6 +4,7 @@
<HEAD>
<TITLE>module location_mod</TITLE>
<link rel="stylesheet" type="text/css" href="../../doc/html/doc.css" />
+<link href="../../doc/html/dart.ico" rel="shortcut icon" />
</HEAD>
<BODY>
<A NAME="TOP"></A>
@@ -469,10 +470,12 @@
<A NAME="get_close_maxdist_init"></A>
<br>
<div class=routine>
-<em class=call> call get_close_maxdist_init(gc,maxdist) </em>
+<em class=call> call get_close_maxdist_init(gc,maxdist
+<em class=optionalcode>, [maxdist_list]</em>) </em>
<pre>
type(get_close_type), intent(inout) :: <em class=code>gc</em>
real(r8), intent(in) :: <em class=code>maxdist</em>
+real(r8), intent(in), optional :: <em class=optionalcode>maxdist_list(:)</em>
</pre>
</div>
@@ -493,6 +496,8 @@
<TD>Data for efficiently finding close locations.</TD></TR>
<TR><TD valign=top><em class=code>maxdist</em></TD>
<TD>Anything closer than this distance is a close location.</TD></TR>
+<TR><TD valign=top><em class=optionalcode>maxdist_list</em></TD>
+ <TD>Ignored for this location type.</TD></TR>
</TABLE>
</div>
Modified: DART/trunk/location/threed_sphere/location_mod.f90
===================================================================
--- DART/trunk/location/threed_sphere/location_mod.f90 2011-02-22 20:25:29 UTC (rev 4741)
+++ DART/trunk/location/threed_sphere/location_mod.f90 2011-02-22 20:50:20 UTC (rev 4742)
@@ -30,6 +30,7 @@
check_namelist_read, do_output, do_nml_file, &
do_nml_term, is_longitude_between
use random_seq_mod, only : random_seq_type, init_random_seq, random_uniform
+use obs_kind_mod, only : get_num_obs_kinds, get_obs_kind_name
use mpi_utilities_mod, only : my_task_id, task_count
implicit none
@@ -43,9 +44,9 @@
operator(==), operator(/=), get_dist, get_close_obs_destroy, &
nc_write_location_atts, nc_get_location_varids, nc_write_location, &
vert_is_height, vert_is_pressure, vert_is_undef, vert_is_level, &
- vert_is_surface, has_vertical_localization, &
+ vert_is_surface, vert_is_scale_height, has_vertical_localization, &
VERTISUNDEF, VERTISSURFACE, VERTISLEVEL, VERTISPRESSURE, &
- VERTISHEIGHT, print_get_close_type, horiz_dist_only
+ VERTISHEIGHT, VERTISSCALEHEIGHT, print_get_close_type, horiz_dist_only
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
@@ -56,11 +57,12 @@
! The possible values for the location_type%which_vert component.
! These are intended to be PRIVATE to this module. Do not make public.
-integer, parameter :: VERTISUNDEF = -2 ! has no vertical location (undefined)
-integer, parameter :: VERTISSURFACE = -1 ! surface value
-integer, parameter :: VERTISLEVEL = 1 ! by level
-integer, parameter :: VERTISPRESSURE = 2 ! by pressure
-integer, parameter :: VERTISHEIGHT = 3 ! by height
+integer, parameter :: VERTISUNDEF = -2 ! has no vertical location (undefined)
+integer, parameter :: VERTISSURFACE = -1 ! surface value
+integer, parameter :: VERTISLEVEL = 1 ! by level
+integer, parameter :: VERTISPRESSURE = 2 ! by pressure
+integer, parameter :: VERTISHEIGHT = 3 ! by height
+integer, parameter :: VERTISSCALEHEIGHT = 4 ! by scale height
type location_type
private
@@ -71,21 +73,24 @@
! Type to facilitate efficient computation of observations close to a given location
type get_close_type
private
- integer :: num
- real(r8) :: maxdist
- integer, pointer :: lon_offset(:, :) ! (nlat, nlat);
- integer, pointer :: obs_box(:) ! (nobs); List of obs indices in boxes
- integer, pointer :: count(:, :) ! (nlon, nlat); # of obs in each box
- integer, pointer :: start(:, :) ! (nlon, nlat); Start of list of obs in this box
- real(r8) :: bot_lat, top_lat ! Bottom and top latitudes of latitude boxes
- real(r8) :: bot_lon, top_lon ! Bottom and top longitudes of longitude boxes
- real(r8) :: lon_width, lat_width ! Width of boxes in lon and lat
- logical :: lon_cyclic ! Do boxes wraparound in longitude?
+ integer :: num
+ real(r8) :: maxdist
+ integer, pointer :: lon_offset(:, :) ! (nlat, nlat);
+ integer, pointer :: obs_box(:) ! (nobs); List of obs indices in boxes
+ integer, pointer :: count(:, :) ! (nlon, nlat); # of obs in each box
+ integer, pointer :: start(:, :) ! (nlon, nlat); Start of list of obs in this box
+ real(r8) :: bot_lat, top_lat ! Bottom and top latitudes of latitude boxes
+ real(r8) :: bot_lon, top_lon ! Bottom and top longitudes of longitude boxes
+ real(r8) :: lon_width, lat_width ! Width of boxes in lon and lat
+ logical :: lon_cyclic ! Do boxes wraparound in longitude?
+ integer :: num_types ! if > 0, cutoffs per type
+ real(r8), pointer :: special_maxdist(:)
+ real(r8) :: original_maxdist
end type get_close_type
type(random_seq_type) :: ran_seq
-logical :: ran_seq_init = .false.
-logical, save :: module_initialized = .false.
+logical :: ran_seq_init = .false.
+logical, save :: module_initialized = .false.
integer, parameter :: LocationDims = 3
character(len = 129), parameter :: LocationName = "loc3Dsphere"
@@ -95,7 +100,8 @@
character(len = 129) :: errstring
! Global storage for vertical distance normalization factors
-real(r8) :: vert_normalization(3)
+real(r8) :: vert_normalization(4)
+real(r8), allocatable :: special_vert_norm(:,:) ! if doing per-type
! Global storage for fast approximate sin and cosine lookups
! PAR For efficiency for small cases might want to fill tables as needed
@@ -109,40 +115,46 @@
!-----------------------------------------------------------------
! Namelist with default values
-! horiz_dist_only == .true. -> Only the great circle horizontal distance is
-! computed in get_dist.
-! horiz_dist_only == .false. -> Square root of sum of squared horizontal and
-! normalized vertical dist computed in get_dist
-! vert_normalization_pressure Number pascals that give a distance equivalent
-! to one radian in horizontal
-! vert_normalization_height -> Number meters that give a distance equivalent
-! to one radian in horizontal
-! vert_normalization_level -> Number levels that give a distance equivalent
-! to one radian in horizontal
-! approximate_distance -> Use a faster table lookup for the trig math.
-! Works well for global models and large areas,
-! and improves performance. For smaller regions
-! might introduce banding, so leave .false.
-! nlon -> Number longitude boxes for get_close_obs
-! nlon MUST BE ODD
-! nlat -> Number latitude boxes for get_close_obs
+! horiz_dist_only == .true. -> Only the great circle horizontal distance is
+! computed in get_dist.
+! horiz_dist_only == .false. -> Square root of sum of squared horizontal and
+! normalized vertical dist computed in get_dist
+! vert_normalization_pressure -> Number pascals that give a distance equivalent
+! to one radian in horizontal
+! vert_normalization_height -> Number meters that give a distance equivalent
+! to one radian in horizontal
+! vert_normalization_level -> Number levels that give a distance equivalent
+! to one radian in horizontal
+! vert_normalization_scale_height -> Number scale heights that give a distance
+! equivalent to one radian in horizontal
+! approximate_distance -> Use a faster table lookup for the trig math.
+! Works well for global models and large areas,
+! and improves performance. For smaller regions
+! might introduce banding, so leave .false.
+! nlon -> Number longitude boxes for get_close_obs
+! nlon MUST BE ODD
+! nlat -> Number latitude boxes for get_close_obs
-logical :: horiz_dist_only = .true.
-real(r8) :: vert_normalization_pressure = 100000.0_r8
-real(r8) :: vert_normalization_height = 10000.0_r8
-real(r8) :: vert_normalization_level = 20.0_r8
-logical :: approximate_distance = .false.
-integer :: nlon = 71
-integer :: nlat = 36
-logical :: output_box_info = .false.
+logical :: horiz_dist_only = .true.
+real(r8) :: vert_normalization_pressure = 100000.0_r8
+real(r8) :: vert_normalization_height = 10000.0_r8
+real(r8) :: vert_normalization_level = 20.0_r8
+real(r8) :: vert_normalization_scale_height = 5.0_r8
+logical :: maintain_original_vert = .false.
+logical :: approximate_distance = .false.
+integer :: nlon = 71
+integer :: nlat = 36
+logical :: output_box_info = .false.
! should be fixed in the code - but leave this here for now.
! search for this variable below in the code for more info
-logical :: num_tasks_insensitive = .false.
-integer :: print_box_level = 0
+logical :: num_tasks_insensitive = .false.
+integer :: print_box_level = 0
-namelist /location_nml/ horiz_dist_only, vert_normalization_pressure, &
- vert_normalization_height, vert_normalization_level, approximate_distance, &
- nlon, nlat, output_box_info, print_box_level, num_tasks_insensitive
+namelist /location_nml/ horiz_dist_only, vert_normalization_pressure, &
+ vert_normalization_height, vert_normalization_level, &
+ vert_normalization_scale_height, approximate_distance, nlon, nlat, &
+ output_box_info, print_box_level, num_tasks_insensitive, &
+ maintain_original_vert
!-----------------------------------------------------------------
@@ -190,6 +202,7 @@
vert_normalization(1) = vert_normalization_level
vert_normalization(2) = vert_normalization_pressure
vert_normalization(3) = vert_normalization_height
+vert_normalization(4) = vert_normalization_scale_height
if (horiz_dist_only) then
call error_handler(E_MSG,'location_mod:', &
@@ -199,12 +212,14 @@
call error_handler(E_MSG,'location_mod:', &
'Including vertical separation when computing distances:', &
source,revision,revdate)
- write(str1,'(A,f17.5)') ' # pascals ~ 1 horiz radian: ', vert_normalization_pressure
+ write(str1,'(A,f17.5)') ' # pascals ~ 1 horiz radian: ', vert_normalization_pressure
call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
- write(str1,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', vert_normalization_height
+ write(str1,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', vert_normalization_height
call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
- write(str1,'(A,f17.5)') ' # model levels ~ 1 horiz radian: ', vert_normalization_level
+ write(str1,'(A,f17.5)') ' # model levels ~ 1 horiz radian: ', vert_normalization_level
call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
+ write(str1,'(A,f17.5)') ' # scale heights ~ 1 horiz radian: ', vert_normalization_scale_height
+ call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
endif
! Set up a lookup table for cos and sin for approximate but fast distances
@@ -240,12 +255,17 @@
! If set to false, this will always do full 3d distance if possible. If set to
! true it will never do the full 3d distance. At present asking to do a vertical
! distance computation for incompatible vertical location types results
-! in a fatal error.
+! in a fatal error unless one of the vertical types is UNDEFINED.
-! The kinds are available if a more sophisticated distance computation is required.
-! This generic code does not use the kinds, but several model_mods intercept this
-! call and compute their own distances based on kind, and so the interface needs
-! to match how the higher level call will be made.
+! In spite of the names, the 3rd and 4th argument are actually specific types
+! (e.g. RADIOSONDE_TEMPERATURE, AIRCRAFT_TEMPERATURE). The types are part of
+! the interface in case user-code wants to do a more sophisticated distance
+! calculation based on the base or target types. In the usual case this
+! code still doesn't use the types, but there's an undocumented feature that
+! allows you to maintain the original vertical normalization even when
+! changing the cutoff distance in the horizontal. For that to work we
+! do need to know the type, and we use the type of loc1 to control it.
+!
type(location_type), intent(in) :: loc1, loc2
integer, optional, intent(in) :: kind1, kind2
@@ -292,8 +312,8 @@
endif
endif
-! Now compute a vertical distance if requested,
-! Highest priority is optional no_vert argument, test it first
+! Now compute a vertical distance if requested. Highest priority is
+! the optional no_vert argument, so test it first.
if(present(no_vert)) then
comp_h_only = no_vert
! Namelist horizontal only has second highest priority
@@ -321,8 +341,18 @@
! Compute the difference and divide by the appropriate normalization factor
! Normalization factor computes relative distance in vertical compared to one radian
- vert_dist = abs(loc1%vloc - loc2%vloc) / vert_normalization(loc1%which_vert)
+ ! This is new - if per-type localization distances given, use the kind of loc1
+ ! to determine the vertical mapping distance. it defaults to the 4 standard ones,
+ ! but can be specified separately if desired.
+ ! note that per-type vertical conversion factors are a global here, and were set
+ ! by the last call to gc_init that specified per/type cutoffs.
+ if (allocated(special_vert_norm)) then
+ vert_dist = abs(loc1%vloc - loc2%vloc) / special_vert_norm(loc1%which_vert, kind1)
+ else
+ vert_dist = abs(loc1%vloc - loc2%vloc) / vert_normalization(loc1%which_vert)
+ endif
+
! Spherical distance shape is computed here, other flavors can be computed
get_dist = sqrt(get_dist**2 + vert_dist**2)
endif
@@ -456,8 +486,8 @@
set_location_single%lon = lon * DEG2RAD
set_location_single%lat = lat * DEG2RAD
-if(which_vert < VERTISUNDEF .or. which_vert == 0 .or. which_vert > VERTISHEIGHT) then
- write(errstring,*)'which_vert (',which_vert,') must be one of -2, -1, 1, 2, or 3'
+if(which_vert < VERTISUNDEF .or. which_vert == 0 .or. which_vert > VERTISSCALEHEIGHT) then
+ write(errstring,*)'which_vert (',which_vert,') must be one of -2, -1, 1, 2, 3, or 4'
call error_handler(E_ERR,'set_location', errstring, source, revision, revdate)
endif
@@ -631,6 +661,8 @@
write(charstring, '(A,1X,F12.7,A)') trim(string1), loc%vloc / 100.0_r8, ' hPa'
case (VERTISHEIGHT)
write(charstring, '(A,1X,F12.7,A)') trim(string1), loc%vloc / 1000.0_r8, ' km'
+ case (VERTISSCALEHEIGHT)
+ write(charstring, '(A,1X,F12.7,A)') trim(string1), loc%vloc, ' scale ht'
case default
write(errstring, *) 'unrecognized key for vertical type: ', loc%which_vert
call error_handler(E_ERR, 'write_location', errstring, source, revision, revdate)
@@ -662,10 +694,10 @@
endif
! Now read the location data value
read(locfile, *)read_location%lon, read_location%lat, &
- read_location%vloc, read_location%which_vert
+ read_location%vloc, read_location%which_vert
else
read(locfile)read_location%lon, read_location%lat, &
- read_location%vloc, read_location%which_vert
+ read_location%vloc, read_location%which_vert
endif
end function read_location
@@ -701,6 +733,7 @@
write(*, *)VERTISLEVEL,' --> model level'
write(*, *)VERTISPRESSURE,' --> pressure'
write(*, *)VERTISHEIGHT,' --> height'
+write(*, *)VERTISSCALEHEIGHT,' --> scale height'
100 read(*, *) location%which_vert
if(location%which_vert == VERTISLEVEL ) then
@@ -720,6 +753,9 @@
!location%vloc = 100.0 * location%vloc ! only applies to pressure
write(*, *) 'Vertical coordinate height'
read(*, *) location%vloc
+else if(location%which_vert == VERTISSCALEHEIGHT ) then
+ write(*, *) 'Vertical coordinate scale height (-ln(p/ps))'
+ read(*, *) location%vloc
else if(location%which_vert == VERTISUNDEF ) then
! a valid floating point value, but should be unused.
location%vloc = MISSING_R8
@@ -851,6 +887,8 @@
'nc_write_location_atts', 'which_vert:VERTISPRESSURE')
call nc_check(nf90_put_att(ncFileID, VarID, 'VERTISHEIGHT', VERTISHEIGHT), &
'nc_write_location_atts', 'which_vert:VERTISHEIGHT')
+call nc_check(nf90_put_att(ncFileID, VarID, 'VERTISSCALEHEIGHT', VERTISSCALEHEIGHT), &
+ 'nc_write_location_atts', 'which_vert:VERTISSCALEHEIGHT')
! If we made it to here without error-ing out ... we're good.
@@ -1072,25 +1110,93 @@
type(get_close_type), intent(inout) :: gc
deallocate(gc%obs_box, gc%lon_offset, gc%count, gc%start)
+if (gc%num_types > 0) deallocate(gc%special_maxdist)
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list