[Dart-dev] [7754] DART/trunk/location/threed_sphere/location_mod.f90: optimized the code when specifying per-type cutoff values;
nancy at ucar.edu
nancy at ucar.edu
Tue Mar 24 16:49:44 MDT 2015
Revision: 7754
Author: nancy
Date: 2015-03-24 16:49:43 -0600 (Tue, 24 Mar 2015)
Log Message:
-----------
optimized the code when specifying per-type cutoff values;
creates a separate gc structure for each different distance
so each search is as specific as possible. compares bitwise
with previous versions of the locations mod.
Modified Paths:
--------------
DART/trunk/location/threed_sphere/location_mod.f90
-------------- next part --------------
Modified: DART/trunk/location/threed_sphere/location_mod.f90
===================================================================
--- DART/trunk/location/threed_sphere/location_mod.f90 2015-03-24 20:53:16 UTC (rev 7753)
+++ DART/trunk/location/threed_sphere/location_mod.f90 2015-03-24 22:49:43 UTC (rev 7754)
@@ -19,14 +19,14 @@
! Note that for now, lev = -1 represents a surface quantity independent
! of vertical discretization as required for Bgrid surface pressure.
-use types_mod, only : r8, MISSING_R8, MISSING_I, PI, RAD2DEG, DEG2RAD
+use types_mod, only : r8, MISSING_R8, MISSING_I, PI, RAD2DEG, DEG2RAD, OBSTYPELENGTH
use utilities_mod, only : register_module, error_handler, E_ERR, ascii_file_format, &
nc_check, E_MSG, open_file, close_file, set_output, &
logfileunit, nmlfileunit, find_namelist_in_file, &
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 obs_kind_mod, only : get_num_obs_kinds, get_obs_kind_name, get_obs_kind_index
use mpi_utilities_mod, only : my_task_id, task_count
implicit none
@@ -66,24 +66,50 @@
integer :: which_vert ! determines if by level, height, pressure, ...
end type location_type
-! Type to facilitate efficient computation of observations close to a given location
+! Derived type to facilitate efficient computation of locations close to a given observation.
+
+type get_close_type_by_type
+ private
+ integer :: num
+ real(r8) :: maxdist ! furthest separation between "close" locations
+ integer, allocatable :: lon_offset(:, :) ! (nlat, nlat), lon box indices searched (nlat 2x IS correct)
+ integer, allocatable :: loc_box(:) ! (nloc), List of loc indices in boxes
+ integer, allocatable :: count(:, :) ! (nlon, nlat), # of loc in each box
+ integer, allocatable :: start(:, :) ! (nlon, nlat), Start of list of loc 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?
+end type get_close_type_by_type
+
+! Support more than a single cutoff distance. nt is the count of
+! distinct cutoffs, which are selected per specific observation type.
+! The map associates the incoming location type with the
+! corresponding gtt index. There are only as many close_types
+! as there are distinct cutoff distances; if more than one specific
+! type has the same cutoff distances they share the type.
+
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_types ! if > 0, cutoffs per type
- real(r8), pointer :: special_maxdist(:)
- real(r8) :: original_maxdist
+ integer :: nt ! number of distinct cutoffs
+ integer, allocatable :: type_to_cutoff_map(:) ! mapping of types to index
+ type(get_close_type_by_type), allocatable :: gtt(:) ! array of close_types by type
end type get_close_type
+! Horizontal localization/cutoff values are passed in by the caller.
+! The Vertical normalization values are globals; are specified by namelist
+! here, and apply to all specific types unless a 'special list' is also specified
+! that overrides the default values.
+
+logical :: has_special_vertical_norms = .false.
+integer :: num_special_vert_norms = 0
+
+! Some calls include information about the type or kind of the location.
+! If the location refers to an observation it is possible to have both
+! a specific type and a generic kind. If the location refers to a
+! state vector item, it only has a generic kind. Some routines have
+! a specific type and a generic kind as arguments; look carefully at
+! the argument names before using.
+
type(random_seq_type) :: ran_seq
logical :: ran_seq_init = .false.
logical, save :: module_initialized = .false.
@@ -93,16 +119,31 @@
character(len = 129), parameter :: LocationLName = &
"threed sphere locations: lon, lat, vertical"
-character(len = 512) :: msgstring
+character(len = 512) :: msgstring, msgstring1, msgstring2
! Global storage for vertical distance normalization factors
-real(r8) :: vert_normalization(4)
-real(r8), allocatable :: special_vert_norm(:,:) ! if doing per-type
+! The 4 below is for the 4 vertical units (pressure, level, height,
+! scale height). undefined and surface don't need vert distances.
+! NOTE: Code that uses VERT_TYPE_COUNT depends on pressure, level,
+! height, and scale height having actual values between 1 and 4, or
+! this code will break.
+integer, parameter :: VERT_TYPE_COUNT = 4
+real(r8) :: vert_normalization(VERT_TYPE_COUNT)
+real(r8), allocatable :: per_type_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
-real(r8) :: my_sin(-630:630), my_cos(-630:630), my_acos(-1000:1000)
+! Also these could be larger for more accurate computations, if needed.
+! 630 is 2 * PI rounded up, times 100.
+real(r8), parameter :: SINCOS_DELTA = 100.0_r8
+integer, parameter :: SINCOS_LIMIT = 630
+real(r8), parameter :: ACOS_DELTA = 1000.0_r8
+integer, parameter :: ACOS_LIMIT = 1000
+real(r8) :: my_sin(-SINCOS_LIMIT:SINCOS_LIMIT)
+real(r8) :: my_cos(-SINCOS_LIMIT:SINCOS_LIMIT)
+real(r8) :: my_acos(-ACOS_LIMIT:ACOS_LIMIT)
+
! Option for verification using exhaustive search
logical :: COMPARE_TO_CORRECT = .false. ! normally false
@@ -139,16 +180,24 @@
integer :: nlat = 36
logical :: output_box_info = .false.
integer :: print_box_level = 0
-! obsolete now - code fixed. leave in for backwards compatibility
-! but does nothing now.
-logical :: num_tasks_insensitive = .false.
+integer, parameter :: MAX_ITEMS = 500
+character(len=OBSTYPELENGTH) :: special_vert_normalization_obs_types(MAX_ITEMS)
+real(r8) :: special_vert_normalization_pressures(MAX_ITEMS)
+real(r8) :: special_vert_normalization_heights(MAX_ITEMS)
+real(r8) :: special_vert_normalization_levels(MAX_ITEMS)
+real(r8) :: special_vert_normalization_scale_heights(MAX_ITEMS)
+
+
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
+ output_box_info, print_box_level, maintain_original_vert, &
+ special_vert_normalization_obs_types, special_vert_normalization_pressures, &
+ special_vert_normalization_heights, special_vert_normalization_levels, &
+ special_vert_normalization_scale_heights
+
!-----------------------------------------------------------------
interface operator(==); module procedure loc_eq; end interface
@@ -167,14 +216,22 @@
! things which need doing exactly once.
-integer :: iunit, io, i
-character(len=129) :: str1
+integer :: iunit, io, i, v, k, typecount, type_index
+
if (module_initialized) return
call register_module(source, revision, revdate)
module_initialized = .true.
+! give these initial values before reading them from the namelist.
+special_vert_normalization_obs_types(:) = 'null'
+special_vert_normalization_pressures(:) = missing_r8
+special_vert_normalization_heights(:) = missing_r8
+special_vert_normalization_levels(:) = missing_r8
+special_vert_normalization_scale_heights(:) = missing_r8
+
+
! Read the namelist entry
call find_namelist_in_file("input.nml", "location_nml", iunit)
read(iunit, nml = location_nml, iostat = io)
@@ -185,24 +242,99 @@
if(do_nml_file()) write(nmlfileunit, nml=location_nml)
if(do_nml_term()) write( * , nml=location_nml)
-! deprecated namelist item
-if (num_tasks_insensitive) then
- call error_handler(E_MSG, 'location_mod:', &
- 'WARNING: namelist item "num_tasks_insensitive" is deprecated and will be removed in the future')
+call error_handler(E_MSG, 'location_mod:', 'using code with optimized cutoffs', &
+ source, revision, revdate)
+
+! Stop supporting this now that you can specify per-type vertical normalization factors
+if (maintain_original_vert) then
+ write(msgstring, '(A)') 'namelist item "maintain_original_vert" is deprecated'
+ write(msgstring1, '(A)') 'use the per-type vertical normalization settings now'
+ write(msgstring2, '(A)') 'remove it from the namelist, or set it to .false. to continue'
+ call error_handler(E_ERR, 'location_mod', msgstring, &
+ source, revision, revdate, text2=msgstring1, text3=msgstring2)
+
endif
-! Make sure that the number of longitudes, nlon, for get_close_obs is odd
+! Simple error checking for nlon, nlat to be sure we can use them for allocations.
+if (nlon < 1 .or. nlat < 1) then
+ write(msgstring, '(A,I5,A,I5)') 'from the namelist, nlon is ', nlon, '; nlat is ', nlat
+ write(msgstring2, '(A)') 'both must be >= 1, and nlon must be odd'
+ call error_handler(E_ERR, 'location_mod', msgstring, &
+ source, revision, revdate, text2=msgstring2)
+endif
+
+! One additional restriction; number of longitudes, nlon, for get_close_obs must be odd
if((nlon / 2) * 2 == nlon) then
- call error_handler(E_ERR, 'initialize_module', 'nlon must be odd', &
- source, revision, revdate)
+ write(msgstring, '(A,I5,A)') 'nlon is ', nlon, '. Must be odd'
+ call error_handler(E_ERR, 'location_mod', msgstring, source, revision, revdate)
endif
! Copy the normalization factors in the vertical into an array
-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
+vert_normalization(VERTISLEVEL) = vert_normalization_level
+vert_normalization(VERTISPRESSURE) = vert_normalization_pressure
+vert_normalization(VERTISHEIGHT) = vert_normalization_height
+vert_normalization(VERTISSCALEHEIGHT) = vert_normalization_scale_height
+! If the user has set a namelist to make different vertical normalization factors
+! when computing the localization distances, allocate and set that up here.
+! It overrides the defaults above.
+
+if (special_vert_normalization_obs_types(1) /= 'null' .or. &
+ special_vert_normalization_pressures(1) /= missing_r8 .or. &
+ special_vert_normalization_heights(1) /= missing_r8 .or. &
+ special_vert_normalization_levels(1) /= missing_r8 .or. &
+ special_vert_normalization_scale_heights(1) /= missing_r8) then
+
+ ! FIXME: add code to check for mismatched length lists. are we going to force
+ ! users to specify all 4 values for any obs type that is not using the defaults?
+
+ typecount = get_num_obs_kinds() ! ignore function name, this is specific type count
+ allocate(per_type_vert_norm(VERT_TYPE_COUNT, typecount))
+
+ ! Set the defaults for all specific types not listed in the special list
+ per_type_vert_norm(VERTISLEVEL, :) = vert_normalization_level
+ per_type_vert_norm(VERTISPRESSURE, :) = vert_normalization_pressure
+ per_type_vert_norm(VERTISHEIGHT, :) = vert_normalization_height
+ per_type_vert_norm(VERTISSCALEHEIGHT, :) = vert_normalization_scale_height
+
+ ! Go through special-treatment observation kinds, if any.
+ num_special_vert_norms = 0
+ k = 0
+ do i = 1, MAX_ITEMS
+ if(special_vert_normalization_obs_types(i) == 'null') exit
+ k = k + 1
+ enddo
+ num_special_vert_norms = k
+
+ if (num_special_vert_norms > 0) has_special_vertical_norms = .true.
+
+ do i = 1, num_special_vert_norms
+ type_index = get_obs_kind_index(special_vert_normalization_obs_types(i))
+ if (type_index < 0) then
+ write(msgstring, *) 'unrecognized TYPE_ in the special vertical normalization namelist:'
+ call error_handler(E_ERR,'location_mod:', msgstring, source, revision, revdate, &
+ text2=trim(special_vert_normalization_obs_types(i)))
+ endif
+
+ per_type_vert_norm(VERTISLEVEL, i) = special_vert_normalization_levels(i)
+ per_type_vert_norm(VERTISPRESSURE, i) = special_vert_normalization_pressures(i)
+ per_type_vert_norm(VERTISHEIGHT, i) = special_vert_normalization_heights(i)
+ per_type_vert_norm(VERTISSCALEHEIGHT, i) = special_vert_normalization_scale_heights(i)
+
+ enddo
+
+ if (any(per_type_vert_norm == missing_r8)) then
+ write(msgstring, *) 'one or more special vertical normalization values is uninitialized.'
+ call error_handler(E_ERR,'location_mod:', &
+ 'special vert normalization value namelist requires all 4 values per type', &
+ source, revision, revdate, &
+ text2=trim(msgstring))
+ endif
+
+endif
+
+! if the namelist control says we are only basing the localization on
+! distances in the horizontal, log that in the dart log file.
if (horiz_dist_only) then
call error_handler(E_MSG,'location_mod:', &
'Ignoring vertical when computing distances; horizontal only', &
@@ -211,26 +343,49 @@
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
- call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
- 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
- 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)
+ write(msgstring,'(A,f17.5)') ' # pascals ~ 1 horiz radian: ', vert_normalization_pressure
+ call error_handler(E_MSG,'location_mod:',msgstring,source,revision,revdate)
+ write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', vert_normalization_height
+ call error_handler(E_MSG,'location_mod:',msgstring,source,revision,revdate)
+ write(msgstring,'(A,f17.5)') ' # model levels ~ 1 horiz radian: ', vert_normalization_level
+ call error_handler(E_MSG,'location_mod:',msgstring,source,revision,revdate)
+ write(msgstring,'(A,f17.5)') ' # scale heights ~ 1 horiz radian: ', vert_normalization_scale_height
+ call error_handler(E_MSG,'location_mod:',msgstring,source,revision,revdate)
+
+ do i = 1, num_special_vert_norms
+ if ((per_type_vert_norm(VERTISLEVEL, i) /= vert_normalization_level) .or. &
+ (per_type_vert_norm(VERTISPRESSURE, i) /= vert_normalization_pressure) .or. &
+ (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) .or. &
+ (per_type_vert_norm(VERTISSCALEHEIGHT, i) /= vert_normalization_scale_height)) then
+
+ write(msgstring,'(2A)') 'Altering vertical normalization for type ', trim(special_vert_normalization_obs_types(i))
+ call error_handler(E_MSG,'location_mod:',msgstring,source,revision,revdate)
+ write(msgstring,'(A,f17.5)') ' # pascals ~ 1 horiz radian: ', &
+ per_type_vert_norm(VERTISPRESSURE, i)
+ call error_handler(E_MSG,'location_mod:',msgstring,source,revision,revdate)
+ write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', &
+ per_type_vert_norm(VERTISHEIGHT, i)
+ call error_handler(E_MSG,'location_mod:',msgstring,source,revision,revdate)
+ write(msgstring,'(A,f17.5)') ' # model levels ~ 1 horiz radian: ', &
+ per_type_vert_norm(VERTISLEVEL, i)
+ call error_handler(E_MSG,'location_mod:',msgstring,source,revision,revdate)
+ write(msgstring,'(A,f17.5)') ' # scale heights ~ 1 horiz radian: ', &
+ per_type_vert_norm(VERTISSCALEHEIGHT, i)
+ call error_handler(E_MSG,'location_mod:',msgstring,source,revision,revdate)
+ endif
+ enddo
endif
! Set up a lookup table for cos and sin for approximate but fast distances
! Don't worry about rounding errors as long as one gives more weight
! Really only need tables half this size, too (sin from -pi/2 to pi/2, cos only +)
if(approximate_distance) then
- do i = -630, 630
- my_cos(i) = cos(i / 100.0_r8)
- my_sin(i) = sin(i / 100.0_r8)
+ do i = -SINCOS_LIMIT, SINCOS_LIMIT
+ my_cos(i) = cos(i / SINCOS_DELTA)
+ my_sin(i) = sin(i / SINCOS_DELTA)
end do
- do i = -1000, 1000
- my_acos(i) = acos(i / 1000.0_r8)
+ do i = -ACOS_LIMIT, ACOS_LIMIT
+ my_acos(i) = acos(i / ACOS_DELTA)
end do
call error_handler(E_MSG,'location_mod:', &
'Using table-lookup approximation for distance computations', &
@@ -286,11 +441,11 @@
if(approximate_distance) then
! Option 1: Use table lookup; faster but less accurate
- lat1_ind = int(loc1%lat*100.0_r8)
- lat2_ind = int(loc2%lat*100.0_r8)
- lon_ind = int(lon_dif*100.0_r8)
- temp = int(1000.0_r8 * (my_sin(lat2_ind) * my_sin(lat1_ind) + &
- my_cos(lat2_ind) * my_cos(lat1_ind) * my_cos(lon_ind)))
+ lat1_ind = int(loc1%lat*SINCOS_DELTA)
+ lat2_ind = int(loc2%lat*SINCOS_DELTA)
+ lon_ind = int(lon_dif *SINCOS_DELTA)
+ temp = int(ACOS_DELTA * (my_sin(lat2_ind) * my_sin(lat1_ind) + &
+ my_cos(lat2_ind) * my_cos(lat1_ind) * my_cos(lon_ind)))
get_dist = my_acos(temp)
else
! Option 2: Use pre-defined trig functions: accurate but slow
@@ -350,9 +505,13 @@
! 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, type1)
+ ! by namelist values. they apply to all calls to get_dist() based on the obs type.
+ if (allocated(per_type_vert_norm)) then
+ if (.not. present(type1)) then
+ write(msgstring, *) 'obs type required in get_dist() if doing per-type vertical normalization'
+ call error_handler(E_MSG, 'get_dist', msgstring, source, revision, revdate)
+ endif
+ vert_dist = abs(loc1%vloc - loc2%vloc) / per_type_vert_norm(loc1%which_vert, type1)
else
vert_dist = abs(loc1%vloc - loc2%vloc) / vert_normalization(loc1%which_vert)
endif
@@ -428,42 +587,6 @@
!---------------------------------------------------------------------------
-!FIXME: this routine does not exist in any other locations module, and i cannot
-! find any code which uses it. I propose we obsolete this code. get_location()
-! returns the contents of a locations type in a location-independent routine.
-function get_location_lon(loc)
-
-! Given a location type, return the longitude. Values stored in radians but
-! returned in degrees.
-
-type(location_type), intent(in) :: loc
-real(r8) :: get_location_lon
-
-if ( .not. module_initialized ) call initialize_module
-
-get_location_lon = loc%lon * RAD2DEG
-
-end function get_location_lon
-
-!---------------------------------------------------------------------------
-
-!FIXME: ditto for comment above. this should be deprecated.
-function get_location_lat(loc)
-
-! Given a location type, return the latitude. Values stored in radians but
-! returned in degrees.
-
-type(location_type), intent(in) :: loc
-real(r8) :: get_location_lat
-
-if ( .not. module_initialized ) call initialize_module
-
-get_location_lat = loc%lat * RAD2DEG
-
-end function get_location_lat
-
-!----------------------------------------------------------------------------
-
function set_location_single(lon, lat, vert_loc, which_vert)
! Puts the given longitude, latitude, and vertical location
@@ -964,126 +1087,135 @@
!----------------------------------------------------------------------------
-subroutine get_close_obs_init(gc, num, obs)
+!----------------------------------------------------------------------------
+
+subroutine get_close_obs_init(gc, num, locs)
-! Initializes part of get_close accelerator that depends on the particular obs
+! Initializes part of get_close accelerator that depends on the particular loc
type(get_close_type), intent(inout) :: gc
integer, intent(in) :: num
-type(location_type), intent(in) :: obs(num)
+type(location_type), intent(in) :: locs(num)
integer :: blat_ind, tlat_ind
integer :: bot_tlat_ind, top_tlat_ind
real(r8) :: base_lat(2), target_lat(2), del_lon, cos_del_lon
real(r8) :: max_del_lon
-integer :: i, j, cum_start, lon_box(num), lat_box(num), tstart(nlon, nlat), tj, bj
+integer :: i, j, n, cum_start, lon_box(num), lat_box(num), tstart(nlon, nlat), tj, bj
-if ( .not. module_initialized ) call initialize_module
+if ( .not. module_initialized ) call initialize_module()
-! Allocate storage for obs number dependent part
-allocate(gc%obs_box(num))
-gc%obs_box(:) = -1
+! store the location counts in all derived types
+do i=1, gc%nt
+ gc%gtt(i)%num = num
+enddo
-! Set the value of num_obs in the structure
-gc%num = num
-
-! If num == 0, no point in going any further.
+! If there are no locs to operate on, no point in going any further.
if (num == 0) return
-! Determine where the boxes should be for this set of obs and maxdist
-call find_box_ranges(gc, obs, num)
+do i=1, gc%nt
+ ! Allocate storage for locs number dependent part
+ allocate(gc%gtt(i)%loc_box(num))
+ gc%gtt(i)%loc_box(:) = -1
+enddo
-! Figure out which boxes are close to a box on a given latitude circle
-! MIGHT AVOID DOING THIS WITH A COPY ROUTINE: HAVE SAME BOXES OFTEN
-do blat_ind = 1, nlat
- ! Search from east side of base block
- ! Start searching out, have to look for closest point in box being checked
- ! Only have to search latitude boxes that are within maximum distance
- bot_tlat_ind = blat_ind - floor(gc%maxdist / gc%lat_width ) - 1
- if(bot_tlat_ind < 1) bot_tlat_ind = 1
- top_tlat_ind = blat_ind + floor(gc%maxdist / gc%lat_width) + 1
- if(top_tlat_ind > nlat) top_tlat_ind = nlat
- do tlat_ind = bot_tlat_ind, top_tlat_ind
- ! Spherical geometry can be tricky
- ! We want to find the MINIMUM distance between two lat/lon boxes
- ! This distance is known to be corners of the boxes. It is also known
- ! to be between the corners such that the longitude difference between
- ! the corners is a minimum. HOWEVER, determining whether it is between
- ! the closest latitudes or not is a non-trivial computation. Hence,
- ! since this isn't done much, we just check all four possible combinations
- ! of latitude and pick the one that gives the closest distance.
- do j = 1, 2
- base_lat(j) = gc%bot_lat + (blat_ind - 2 + j) * gc%lat_width
- target_lat(j) = gc%bot_lat + (tlat_ind - 2 + j) * gc%lat_width
- end do
+do n=1, gc%nt
+ ! Determine where the boxes should be for this set of locs and maxdist
+ call find_box_ranges(gc%gtt(n), num, locs)
- ! If the max distance > PI, then everybody is close.
- ! Do a test for something slightly less than PI to avoid round-off error.
- ! Set max_del_lon to something much larger than PI since it doesn't matter.
- if(gc%maxdist > PI - 0.0001_r8) then
- max_del_lon = 2.0 * PI
- else
- ! Find the maximum longitude offset for the different possible latitudes edges
- max_del_lon = 0.0_r8
- do tj = 1, 2
- do bj = 1, 2
- ! Compute the lon offset directly by inverting distance
- cos_del_lon = (cos(gc%maxdist) - sin(base_lat(bj)) * sin(target_lat(tj))) / &
- (cos(base_lat(bj)) * cos(target_lat(tj)))
- if(cos_del_lon < -1.0_r8) then
- del_lon = PI
- else if(cos_del_lon > 1.0_r8) then
- del_lon = 0.0_r8
- else
- del_lon = acos(cos_del_lon)
- endif
- if(del_lon > max_del_lon) max_del_lon = del_lon
-
+ ! Figure out which boxes are close to a box on a given latitude circle
+ ! MIGHT AVOID DOING THIS WITH A COPY ROUTINE: HAVE SAME BOXES OFTEN
+ do blat_ind = 1, nlat
+ ! Search from east side of base block
+ ! Start searching out, have to look for closest point in box being checked
+ ! Only have to search latitude boxes that are within maximum distance
+ bot_tlat_ind = blat_ind - floor(gc%gtt(n)%maxdist / gc%gtt(n)%lat_width) - 1
+ if(bot_tlat_ind < 1) bot_tlat_ind = 1
+ top_tlat_ind = blat_ind + floor(gc%gtt(n)%maxdist / gc%gtt(n)%lat_width) + 1
+ if(top_tlat_ind > nlat) top_tlat_ind = nlat
+ do tlat_ind = bot_tlat_ind, top_tlat_ind
+ ! Spherical geometry can be tricky
+ ! We want to find the MINIMUM distance between two lat/lon boxes
+ ! This distance is known to be corners of the boxes. It is also known
+ ! to be between the corners such that the longitude difference between
+ ! the corners is a minimum. HOWEVER, determining whether it is between
+ ! the closest latitudes or not is a non-trivial computation. Hence,
+ ! since this isn't done much, we just check all four possible combinations
+ ! of latitude and pick the one that gives the closest distance.
+ do j = 1, 2
+ base_lat(j) = gc%gtt(n)%bot_lat + (blat_ind - 2 + j) * gc%gtt(n)%lat_width
+ target_lat(j) = gc%gtt(n)%bot_lat + (tlat_ind - 2 + j) * gc%gtt(n)%lat_width
+ end do
+
+ ! If the max distance > PI, then everybody is close.
+ ! Do a test for something slightly less than PI to avoid round-off error.
+ ! Set max_del_lon to something much larger than PI since it doesn't matter.
+ if(gc%gtt(n)%maxdist > PI - 0.0001_r8) then
+ max_del_lon = 2.0 * PI
+ else
+ ! Find the maximum longitude offset for the different possible latitudes edges
+ max_del_lon = 0.0_r8
+ do tj = 1, 2
+ do bj = 1, 2
+ ! Compute the lon offset directly by inverting distance
+ cos_del_lon = (cos(gc%gtt(n)%maxdist) - sin(base_lat(bj)) * sin(target_lat(tj))) / &
+ (cos(base_lat(bj)) * cos(target_lat(tj)))
+ if(cos_del_lon < -1.0_r8) then
+ del_lon = PI
+ else if(cos_del_lon > 1.0_r8) then
+ del_lon = 0.0_r8
+ else
+ del_lon = acos(cos_del_lon)
+ endif
+ if(del_lon > max_del_lon) max_del_lon = del_lon
+
+ end do
end do
- end do
- endif
-
- ! Compute the number of boxes to search in longitude for maximum del_lon
- gc%lon_offset(blat_ind, tlat_ind) = floor(max_del_lon / gc%lon_width) + 1
- ! Watch for roundoff leading to a search of more offsets than exist
- if(gc%lon_offset(blat_ind, tlat_ind) > nlon / 2) &
- gc%lon_offset(blat_ind, tlat_ind) = nlon / 2
-
+ endif
+
+ ! Compute the number of boxes to search in longitude for maximum del_lon
+ gc%gtt(n)%lon_offset(blat_ind, tlat_ind) = floor(max_del_lon / gc%gtt(n)%lon_width) + 1
+ ! Watch for roundoff leading to a search of more offsets than exist
+ if(gc%gtt(n)%lon_offset(blat_ind, tlat_ind) > nlon / 2) &
+ gc%gtt(n)%lon_offset(blat_ind, tlat_ind) = nlon / 2
+
+ end do
end do
-end do
-! Begin by computing the number of observations in each box in lat/lon
-gc%count = 0
-do i = 1, num
- lon_box(i) = get_lon_box(gc, obs(i)%lon)
- if(lon_box(i) < 0 .or. lon_box(i) > nlon) then
- write(msgstring, *) 'Contact Dart Developers: this error should not happen'
- call error_handler(E_MSG, 'get_close_obs_init', msgstring, source, revision, revdate)
- write(msgstring, *) 'obs outside grid boxes, index value:', lon_box(i)
- call error_handler(E_ERR, 'get_close_obs_init', msgstring, source, revision, revdate)
- endif
-
- lat_box(i) = floor((obs(i)%lat - gc%bot_lat) / gc%lat_width) + 1
- if(lat_box(i) > nlat) lat_box(i) = nlat
- if(lat_box(i) < 1) lat_box(i) = 1
-
- gc%count(lon_box(i), lat_box(i)) = gc%count(lon_box(i), lat_box(i)) + 1
-end do
-
-! Figure out where storage for each boxes members should begin
-cum_start = 1
-do i = 1, nlon
- do j = 1, nlat
- gc%start(i, j) = cum_start
- cum_start = cum_start + gc%count(i, j)
+ ! Begin by computing the number of locations in each box in lat/lon
+ gc%gtt(n)%count = 0
+ do i = 1, num
+ lon_box(i) = get_lon_box(gc%gtt(n), locs(i)%lon)
+ if(lon_box(i) < 0 .or. lon_box(i) > nlon) then
+ write(msgstring, *) 'Contact Dart Developers: this error should not happen'
+ call error_handler(E_MSG, 'get_close_obs_init', msgstring, source, revision, revdate)
+ write(msgstring, *) 'loc outside grid boxes, index value:', lon_box(i)
+ call error_handler(E_ERR, 'get_close_obs_init', msgstring, source, revision, revdate)
+ endif
+
+ lat_box(i) = floor((locs(i)%lat - gc%gtt(n)%bot_lat) / gc%gtt(n)%lat_width) + 1
+ if(lat_box(i) > nlat) lat_box(i) = nlat
+ if(lat_box(i) < 1) lat_box(i) = 1
+
+ gc%gtt(n)%count(lon_box(i), lat_box(i)) = gc%gtt(n)%count(lon_box(i), lat_box(i)) + 1
end do
-end do
+
+ ! Figure out where storage for each boxes members should begin
+ cum_start = 1
+ do i = 1, nlon
+ do j = 1, nlat
+ gc%gtt(n)%start(i, j) = cum_start
+ cum_start = cum_start + gc%gtt(n)%count(i, j)
+ end do
+ end do
+
+ ! Now we know how many are in each box, get a list of which are in each box
+ tstart = gc%gtt(n)%start
+ do i = 1, num
+ gc%gtt(n)%loc_box(tstart(lon_box(i), lat_box(i))) = i
+ tstart(lon_box(i), lat_box(i)) = tstart(lon_box(i), lat_box(i)) + 1
+ end do
-! Now we know how many are in each box, get a list of which are in each box
-tstart = gc%start
-do i = 1, num
- gc%obs_box(tstart(lon_box(i), lat_box(i))) = i
- tstart(lon_box(i), lat_box(i)) = tstart(lon_box(i), lat_box(i)) + 1
end do
! info on how well the boxes are working. by default print nothing.
@@ -1095,12 +1227,12 @@
! if print level > 2, set all tasks to print and call print.
! then reset the status to off again.
if (do_output()) then
- call print_get_close_type(gc, print_box_level)
+ call print_get_close_type(gc, 1, print_box_level)
else if (print_box_level >= 2 .or. print_box_level < 0) then
! print status was false, but turn on temporarily
! to output box info from all tasks.
call set_output(.true.)
- call print_get_close_type(gc, print_box_level)
+ call print_get_close_type(gc, 1, print_box_level)
call set_output(.false.)
endif
endif
@@ -1113,12 +1245,17 @@
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)
-! since this is a global, keep it around. it is always allocated
-! the same size and can be reused for any gc.
-!if (allocated(special_vert_norm)) deallocate(special_vert_norm)
+integer :: i
+! order matters here. do all arrays inside the gtt types
+! before releasing the array of types
+do i=1, gc%nt
+ if (allocated(gc%gtt(i)%loc_box)) deallocate(gc%gtt(i)%loc_box)
+ deallocate(gc%gtt(i)%lon_offset, gc%gtt(i)%count, gc%gtt(i)%start)
+enddo
+deallocate(gc%type_to_cutoff_map)
+deallocate(gc%gtt)
+
end subroutine get_close_obs_destroy
!----------------------------------------------------------------------------
@@ -1129,83 +1266,61 @@
real(r8), intent(in) :: maxdist
real(r8), intent(in), optional :: maxdist_list(:)
-character(len=129) :: str1
-integer :: i
-! a bit of a hack - it only prints the first time this routine is
-! called. in filter it's called twice with the same args and this
-! code is counting on that.
-logical, save :: firsttime = .true.
+integer :: i, typecount, distcount, mapnum
+real(r8), allocatable :: distlist(:)
-! Allocate the storage for the grid dependent boxes
-allocate(gc%lon_offset(nlat, nlat), gc%count(nlon, nlat), gc%start(nlon, nlat))
-gc%lon_offset = -1
-gc%count = -1
-gc%start = -1
+! Support per-loc-type localization more efficiently.
+typecount = get_num_obs_kinds() ! ignore function name, this is specific type count
+allocate(gc%type_to_cutoff_map(typecount))
-! set the default value. if there is a list and any distance in the list
-! is larger, use that instead. (the boxes need to be calculated based on
-! the largest possible distance).
-gc%maxdist = maxdist
-gc%original_maxdist = maxdist
-gc%num_types = 0
-
if (present(maxdist_list)) then
- gc%num_types = get_num_obs_kinds()
- if (size(maxdist_list) .ne. gc%num_types) then
- write(msgstring,'(A,I8,A,I8)')'maxdist_list len must equal number of kinds, ', &
- size(maxdist_list), ' /= ', gc%num_types
+ if (size(maxdist_list) .ne. typecount) then
+ write(msgstring,'(A,I8,A,I8)')'maxdist_list len must equal number of specific types, ', &
+ size(maxdist_list), ' /= ', typecount
call error_handler(E_ERR, 'get_close_maxdist_init', msgstring, source, revision, revdate)
endif
- allocate(gc%special_maxdist(gc%num_types))
- if (.not.allocated(special_vert_norm)) &
- allocate(special_vert_norm(4, gc%num_types)) ! fill from namelist here, or assim_tools?
-
- gc%special_maxdist(:) = maxdist_list(:)
- gc%maxdist = maxval(gc%special_maxdist)
+
+ allocate(distlist(typecount))
+ call distinct_values(maxdist_list, distcount, distlist, gc%type_to_cutoff_map)
+ gc%nt = distcount
+ if (gc%nt <= 0) then
+ write(msgstring,'(A)')'error getting count of distinct cutoff dists; should not happen'
+ call error_handler(E_ERR, 'get_close_maxdist_init', msgstring, source, revision, revdate)
+ endif
+else
+ gc%nt = 1
+ gc%type_to_cutoff_map(:) = 1
+endif
- ! by default, the vertical changes along with the horizontal to keep
- ! the aspect ratio of the elipsoid constant.
- special_vert_norm(1, :) = vert_normalization_level
- special_vert_norm(2, :) = vert_normalization_pressure
- special_vert_norm(3, :) = vert_normalization_height
- special_vert_norm(4, :) = vert_normalization_scale_height
+! make a gtt type array for each different cutoff distance
+! (just 1 type is the most common case.)
+allocate(gc%gtt(gc%nt))
- ! keep the original vertical distance constant even as you change
- ! the horizontal localization distance.
- if (maintain_original_vert) then
- special_vert_norm(1, :) = vert_normalization_level * &
- (gc%original_maxdist / maxdist_list(:))
- special_vert_norm(2, :) = vert_normalization_pressure * &
- (gc%original_maxdist / maxdist_list(:))
- special_vert_norm(3, :) = vert_normalization_height * &
- (gc%original_maxdist / maxdist_list(:))
- special_vert_norm(4, :) = vert_normalization_scale_height * &
- (gc%original_maxdist / maxdist_list(:))
- if (firsttime) then
- do i = 1, gc%num_types
- if (abs(gc%original_maxdist - maxdist_list(i)) < 1.0e-14_r8) cycle
-
- write(str1,'(2A)') 'Altering vertical normalization for type ', get_obs_kind_name(i)
- call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
- write(str1,'(A,f17.5)') ' # pascals ~ 1 horiz radian: ', special_vert_norm(2, i)
- call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
- write(str1,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', special_vert_norm(3, i)
- call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
- write(str1,'(A,f17.5)') ' # model levels ~ 1 horiz radian: ', special_vert_norm(1, i)
- call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
- write(str1,'(A,f17.5)') ' # scale heights ~ 1 horiz radian: ', special_vert_norm(4, i)
- call error_handler(E_MSG,'location_mod:',str1,source,revision,revdate)
- enddo
- firsttime = .false.
- endif
- endif
+if (present(maxdist_list)) then
+ do i=1, gc%nt
+ gc%gtt(i)%maxdist = distlist(i)
+ enddo
+else
+ ! no per-type settings, everyone uses same distance
+ gc%gtt(1)%maxdist = maxdist
endif
+if (present(maxdist_list)) deallocate(distlist) ! temp storage
+
+! Allocate the storage for the grid dependent boxes
+do i=1, gc%nt
+ allocate(gc%gtt(i)%count(nlon, nlat), gc%gtt(i)%start(nlon, nlat))
+ allocate(gc%gtt(i)%lon_offset(nlat, nlat))
+ gc%gtt(i)%lon_offset = -1
+ gc%gtt(i)%count = -1
+ gc%gtt(i)%start = -1
+enddo
+
end subroutine get_close_maxdist_init
!----------------------------------------------------------------------------
-subroutine get_close_obs(gc, base_obs_loc, base_obs_type, obs, obs_kind, &
+subroutine get_close_obs(gc, base_obs_loc, base_obs_type, locs, loc_kind, &
num_close, close_ind, dist)
! The specific type of the base observation, plus the generic kinds list
@@ -1213,57 +1328,78 @@
! distance computation is needed.
type(get_close_type), intent(in) :: gc
-type(location_type), intent(in) :: base_obs_loc, obs(:)
-integer, intent(in) :: base_obs_type, obs_kind(:)
+type(location_type), intent(in) :: base_obs_loc, locs(:)
+integer, intent(in) :: base_obs_type, loc_kind(:)
integer, intent(out) :: num_close, close_ind(:)
real(r8), optional, intent(out) :: dist(:)
! If dist is NOT present, just find everybody in a box, put them in the list,
! but don't compute any distances
-integer :: lon_box, lat_box, i, j, k, n_lon, lon_ind, n_in_box, st, t_ind
+integer :: lon_box, lat_box, i, j, k, n_lon, lon_ind, n_in_box, st, t_ind, bt
real(r8) :: this_dist, this_maxdist
! Variables needed for comparing against correct case
-integer :: cnum_close, cclose_ind(size(obs))
-real(r8) :: cdist(size(obs))
+integer :: cnum_close, cclose_ind(size(locs))
+real(r8) :: cdist(size(locs))
+
! First, set the intent out arguments to a missing value
num_close = 0
close_ind = -99
if(present(dist)) dist = -99.0_r8
this_dist = 999999.0_r8 ! something big.
-! the list of locations in the obs() argument must be the same
+! map from type index to gtt index
+if (base_obs_type < 1 .or. base_obs_type > size(gc%type_to_cutoff_map)) then
+ write(msgstring,'(A,I8)')'base_obs_type out of range, is ', base_obs_type
+ write(msgstring1,'(A,2I8)')'must be between ', 1, size(gc%type_to_cutoff_map)
+ call write_location (0, base_obs_loc, charstring=msgstring2)
+ call error_handler(E_ERR, 'get_close_obs', msgstring, source, revision, revdate, &
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list