[Dart-dev] [5671] DART/branches/development/location: alternate location type donated by alex reinecke from the NRL.
nancy at ucar.edu
nancy at ucar.edu
Mon Apr 9 10:43:32 MDT 2012
Revision: 5671
Author: nancy
Date: 2012-04-09 10:43:31 -0600 (Mon, 09 Apr 2012)
Log Message:
-----------
alternate location type donated by alex reinecke from the NRL.
has not been extensively tested.
Added Paths:
-----------
DART/branches/development/location/channel/
DART/branches/development/location/channel/location_mod.f90
DART/branches/development/location/channel/location_mod.html
DART/branches/development/location/channel/location_mod.nml
DART/branches/development/location/channel/test/
DART/branches/development/location/channel/test/input.nml
DART/branches/development/location/channel/test/location_test_file
DART/branches/development/location/channel/test/mkmf_location_test
DART/branches/development/location/channel/test/mkmf_location_test3
DART/branches/development/location/channel/test/path_names_location_test
DART/branches/development/location/channel/test/path_names_location_test3
DART/branches/development/location/channel/test/test.in
-------------- next part --------------
Added: DART/branches/development/location/channel/location_mod.f90
===================================================================
--- DART/branches/development/location/channel/location_mod.f90 (rev 0)
+++ DART/branches/development/location/channel/location_mod.f90 2012-04-09 16:43:31 UTC (rev 5671)
@@ -0,0 +1,2218 @@
+! DART software - Copyright 2004 - 2011 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
+
+module location_mod
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+! Implements location interfaces for a 3d channel in X,Y,Z where X is periodic,
+! Y has walls (limited domain), and Z is infinite
+
+use types_mod, only : r8, MISSING_R8, MISSING_I, PI, RAD2DEG, DEG2RAD
+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 mpi_utilities_mod, only : my_task_id, task_count
+
+implicit none
+private
+
+public :: location_type, get_location, set_location, &
+ set_location_missing, is_location_in_region, &
+ write_location, read_location, interactive_location, query_location, &
+ LocationDims, LocationName, LocationLName, get_close_obs, &
+ get_close_maxdist_init, get_close_obs_init, get_close_type, &
+ 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, vert_is_scale_height, has_vertical_localization, &
+ print_get_close_type, find_nearest
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+type location_type
+ private
+ real(r8) :: x, y, z
+end type location_type
+
+! This version supports both regularly spaced boxes, and octree division
+! of the space. for octrees, divide each dim in half until N numbers of filled
+! boxes, or octree reaches some depth? give some threshold where you don't
+! divide a box with less than N points in it?
+
+! contrast with kD-trees (divide along dimensions, not points), and there are
+! two types of octrees - PR (point region) where the regions split at an
+! explicit point, vs MX tree where the split is defined to be at the center
+! of the region.
+
+! if the underlying geometry is spherical, there will be many many empty boxes
+! if we uniformly divide up space, and worse, existing locations will be
+! clustered in a few boxes.
+
+
+! fortran doesn't let you make arrays of pointers, but you can make a
+! derived type containing a pointer, and then make arrays of that derived type.
+! i'm sure if i think about this hard enough i'll figure out why this is so,
+! but for now i'll just believe the great google which tells me it's this way.
+type octree_ptr
+ private
+ type(octree_type), pointer :: p
+end type octree_ptr
+
+type octree_type
+ private
+ integer :: count ! count in this cube, -1 for non-terminal cube
+ integer, pointer :: index(:) ! list of indices in this cube, count long
+ type(octree_ptr), allocatable :: children(:,:,:) ! subcubes
+ type(octree_type), pointer :: parent ! who made you
+ type(location_type) :: llb ! xyz of lower left bottom
+ type(location_type) :: split ! xyz of split point
+ type(location_type) :: urt ! xyz of upper right top
+end type octree_type
+
+type box_type
+ private
+ integer, pointer :: obs_box(:) ! (nobs); List of obs indices in boxes
+ integer, pointer :: count(:, :, :) ! (nx, ny, nz); # of obs in each box
+ integer, pointer :: start(:, :, :) ! (nx, ny, nz); Start of list of obs in this box
+ real(r8) :: bot_x, top_x ! extents in x, y, z
+ real(r8) :: bot_y, top_y
+ real(r8) :: bot_z, top_z
+ real(r8) :: x_width, y_width, z_width ! widths of boxes in x,y,z
+ real(r8) :: nboxes_x, nboxes_y, nboxes_z ! based on maxdist how far to search
+end type box_type
+
+! Type to facilitate efficient computation of observations close to a given location
+type get_close_type
+ private
+ integer :: num
+ real(r8) :: maxdist
+ type(box_type) :: box
+ type(octree_type) :: root
+end type get_close_type
+
+type(random_seq_type) :: ran_seq
+logical :: ran_seq_init = .false.
+logical, save :: module_initialized = .false.
+
+integer, parameter :: LocationDims = 3
+character(len = 129), parameter :: LocationName = "loc3Dchan"
+character(len = 129), parameter :: LocationLName = &
+ "threed channel locations: x, y, z"
+
+character(len = 512) :: errstring
+
+real(r8) :: radius ! used only for converting points on a sphere into x,y,z and back
+
+! If maxdist stays the same, don't need to do box distance calculations
+integer :: last_maxdist = -1.0
+
+!-----------------------------------------------------------------
+! Namelist with default values
+
+! FIXME: give these better defaults? or put in namelist
+! right now they are recomputed based on the setting of nboxes;
+! each is about the cube root of nboxes, but if the bounds are
+! very different distances in each dim, we may want better control
+! over the ratios of each.
+integer :: nx = 10
+integer :: ny = 10
+integer :: nz = 10
+
+logical :: output_box_info = .false.
+integer :: print_box_level = 0
+! tuning options
+integer :: nboxes = 1000 ! suggestion for max number of nodes
+integer :: maxdepth = 4 ! suggestion for max tree depth
+integer :: filled = 10 ! threshold at which you quit splitting
+logical :: use_octree = .false. ! if false, use regular boxes, true = octree
+
+! Option for verification using exhaustive search
+logical :: compare_to_correct = .true. ! normally false
+
+namelist /location_nml/ &
+ filled, nboxes, maxdepth, use_octree, &
+ compare_to_correct, output_box_info, print_box_level
+
+!-----------------------------------------------------------------
+
+interface operator(==); module procedure loc_eq; end interface
+interface operator(/=); module procedure loc_ne; end interface
+
+interface set_location
+ module procedure set_location_single
+ module procedure set_location_array
+end interface set_location
+
+contains
+
+!----------------------------------------------------------------------------
+
+subroutine initialize_module
+
+! things which need doing exactly once.
+
+integer :: iunit, io, i
+character(len=129) :: str1
+
+if (module_initialized) return
+
+call register_module(source, revision, revdate)
+module_initialized = .true.
+
+! Read the namelist entry
+call find_namelist_in_file("input.nml", "location_nml", iunit)
+read(iunit, nml = location_nml, iostat = io)
+call check_namelist_read(iunit, io, "location_nml")
+
+! Write the namelist values to the log file
+
+if(do_nml_file()) write(nmlfileunit, nml=location_nml)
+if(do_nml_term()) write( * , nml=location_nml)
+
+if (filled < 1) then
+ write(errstring,*)'filled sets limit for number of points per box. must be >= 1'
+ call error_handler(E_ERR, 'set_location', errstring, source, revision, revdate)
+endif
+
+end subroutine initialize_module
+
+!----------------------------------------------------------------------------
+
+function get_dist(loc1, loc2, kind1, kind2)
+
+! returns the distance between 2 locations
+
+! 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
+real(r8) :: get_dist
+
+real(r8) :: x_dif, y_dif, z_dif
+
+if ( .not. module_initialized ) call initialize_module
+
+x_dif = loc1%x - loc2%x
+y_dif = loc1%y - loc2%y
+z_dif = loc1%z - loc2%z
+
+get_dist = sqrt(x_dif * x_dif + y_dif * y_dif + z_dif * z_dif)
+
+end function get_dist
+
+!---------------------------------------------------------------------------
+
+function loc_eq(loc1,loc2)
+
+! Interface operator used to compare two locations.
+! Returns true only if all components are 'the same' to within machine
+! precision.
+
+type(location_type), intent(in) :: loc1, loc2
+logical :: loc_eq
+
+if ( .not. module_initialized ) call initialize_module
+
+loc_eq = .false.
+
+if ( abs(loc1%x - loc2%x ) > epsilon(loc1%x ) ) return
+if ( abs(loc1%y - loc2%y ) > epsilon(loc1%y ) ) return
+if ( abs(loc1%z - loc2%z ) > epsilon(loc1%z ) ) return
+
+loc_eq = .true.
+
+end function loc_eq
+
+!---------------------------------------------------------------------------
+
+function loc_ne(loc1,loc2)
+
+! Interface operator used to compare two locations.
+! Returns true if locations are not identical to machine precision.
+
+type(location_type), intent(in) :: loc1, loc2
+logical :: loc_ne
+
+if ( .not. module_initialized ) call initialize_module
+
+loc_ne = (.not. loc_eq(loc1,loc2))
+
+end function loc_ne
+
+!---------------------------------------------------------------------------
+
+function get_location(loc)
+
+! Given a location type return the x,y,z coordinates
+
+type(location_type), intent(in) :: loc
+real(r8), dimension(3) :: get_location
+
+if ( .not. module_initialized ) call initialize_module
+
+get_location(1) = loc%x
+get_location(2) = loc%y
+get_location(3) = loc%z
+
+end function get_location
+
+!---------------------------------------------------------------------------
+
+function set_location_single(x, y, z)
+
+! Puts the x, y, z into a location datatype.
+
+real(r8), intent(in) :: x, y, z
+type (location_type) :: set_location_single
+
+if ( .not. module_initialized ) call initialize_module
+
+set_location_single%x = x
+set_location_single%y = y
+set_location_single%z = z
+
+end function set_location_single
+
+!----------------------------------------------------------------------------
+
+function set_location_array(list)
+
+! location semi-independent interface routine
+! given 3 float numbers, call the underlying set_location routine
+
+real(r8), intent(in) :: list(:)
+type (location_type) :: set_location_array
+
+if ( .not. module_initialized ) call initialize_module
+
+if (size(list) < 3) then
+ write(errstring,*)'requires 3 input values'
+ call error_handler(E_ERR, 'set_location', errstring, source, revision, revdate)
+endif
+
+set_location_array = set_location_single(list(1), list(2), list(3))
+
+end function set_location_array
+
+!----------------------------------------------------------------------------
+
+function set_location_missing()
+
+! Initialize a location type to indicate the contents are unset.
+
+type (location_type) :: set_location_missing
+
+if ( .not. module_initialized ) call initialize_module
+
+set_location_missing%x = MISSING_R8
+set_location_missing%y = MISSING_R8
+set_location_missing%z = MISSING_R8
+
+end function set_location_missing
+
+!---------------------------------------------------------------------------
+
+function query_location(loc, attr)
+
+! Returns the value of the attribute
+
+type(location_type), intent(in) :: loc
+character(len=*), optional, intent(in) :: attr
+real(r8) :: query_location
+
+if ( .not. module_initialized ) call initialize_module
+
+! before you change any of the code in this subroutine,
+! check out the extensive comments in the threed_sphere
+! version of this routine. then proceed with caution.
+
+if (.not. present(attr)) then
+ query_location = loc%x
+ return
+endif
+
+select case(attr)
+ case ('x','X')
+ query_location = loc%x
+ case ('y','Y')
+ query_location = loc%y
+ case ('z','Z')
+ query_location = loc%z
+ case default
+ call error_handler(E_ERR, 'query_location:', &
+ 'Only "X","Y","Z" are legal attributes to request from location', &
+ source, revision, revdate)
+end select
+
+end function query_location
+
+!----------------------------------------------------------------------------
+
+subroutine write_location(locfile, loc, fform, charstring)
+
+! Writes a location to a file.
+! most recent change: adding the optional charstring option. if present,
+! locfile is ignored, and a pretty-print formatting is done into charstring.
+
+integer, intent(in) :: locfile
+type(location_type), intent(in) :: loc
+character(len = *), intent(in), optional :: fform
+character(len = *), intent(out), optional :: charstring
+
+integer :: charlength
+logical :: writebuf
+character(len=129) :: string1
+
+10 format(1X,3(G25.16,1X))
+
+if ( .not. module_initialized ) call initialize_module
+
+! writing to a file (normal use) or to a character buffer?
+writebuf = present(charstring)
+
+! output file; test for ascii or binary, write what's asked, and return
+if (.not. writebuf) then
+ if (ascii_file_format(fform)) then
+ write(locfile, '(''loc3Dchan'')' )
+ write(locfile, 10) loc%x, loc%y, loc%z
+ else
+ write(locfile) loc%x, loc%y, loc%z
+ endif
+ return
+endif
+
+! you only get here if you're writing to a buffer and not
+! to a file, and you can't have binary format set.
+if (.not. ascii_file_format(fform)) then
+ call error_handler(E_ERR, 'write_location', &
+ 'Cannot use string buffer with binary format', &
+ source, revision, revdate)
+endif
+
+! format the location to be human-friendly
+
+! the output can be no longer than this
+charlength = 70
+
+if (len(charstring) < charlength) then
+ write(errstring, *) 'charstring buffer must be at least ', charlength, ' chars long'
+ call error_handler(E_ERR, 'write_location', errstring, source, revision, revdate)
+endif
+
+! format into the outout string
+write(charstring, '(A,3(G20.8,1X))') 'X/Y/Z: ', loc%x, loc%y, loc%z
+
+end subroutine write_location
+
+!----------------------------------------------------------------------------
+
+function read_location(locfile, fform)
+
+! Reads a location from a file that was written by write_location.
+! See write_location for additional discussion.
+
+integer, intent(in) :: locfile
+character(len = *), intent(in), optional :: fform
+type(location_type) :: read_location
+
+character(len=8) :: header
+
+if ( .not. module_initialized ) call initialize_module
+
+if (ascii_file_format(fform)) then
+ read(locfile, '(A8)' ) header
+ if(header /= 'loc3Dchan') then
+ write(errstring,*)'Expected location header "loc3Dchan" in input file, got ', header
+ call error_handler(E_ERR, 'read_location', errstring, source, revision, revdate)
+ endif
+ ! Now read the location data value
+ read(locfile, *)read_location%x, read_location%y, read_location%z
+else
+ read(locfile)read_location%x, read_location%y, read_location%z
+endif
+
+end function read_location
+
+!--------------------------------------------------------------------------
+
+subroutine interactive_location(location, set_to_default)
+
+! Allows for interactive input of a location. Also gives option of selecting
+! a uniformly distributed random location.
+
+type(location_type), intent(out) :: location
+logical, intent(in), optional :: set_to_default
+
+real(r8) :: v(3), minv, maxv
+character(len=1) :: l(3)
+integer :: i, r
+
+if ( .not. module_initialized ) call initialize_module
+
+! If set_to_default is true, then just zero out and return
+if(present(set_to_default)) then
+ if(set_to_default) then
+ location%x = 0.0
+ location%y = 0.0
+ location%z = 0.0
+ return
+ endif
+endif
+
+l(1) = 'X'
+l(2) = 'Y'
+l(3) = 'Z'
+
+! prompt for an explicit location or a random one.
+! if random, generate all 3 x/y/z values randomly.
+! if you want to make some combination of x/y/z random
+! and specify others, you would have to move this read
+! into the loop.
+
+r = 1
+do while (r > 0)
+
+ write(*, *) 'Input 0 to specify a value for the location, or'
+ write(*, *) '-1 for a uniformly distributed random location'
+ read(*, *) r
+
+ if (r > 0) write(*, *) 'Please input 0 or -1 for selection'
+enddo
+
+do i = 1, 3
+ if (r == 0) then
+ write(*, *) 'Input value for ', l(i)
+ read (*,*) v(i)
+
+ else
+ ! Need to make sure random sequence is initialized
+
+ if(.not. ran_seq_init) then
+ call init_random_seq(ran_seq)
+ ran_seq_init = .TRUE.
+ endif
+
+ write(*, *) 'Input minimum ', l(i), ' value '
+ read(*, *) minv
+
+ write(*, *) 'Input maximum ', l(i), ' value '
+ read(*, *) maxv
+
+ v(i) = random_uniform(ran_seq) * (maxv-minv) + minv
+
+ write(*, *) 'random location is ', v(i)
+
+ endif
+
+enddo
+
+location%x = v(1)
+location%y = v(2)
+location%z = v(3)
+
+end subroutine interactive_location
+
+!----------------------------------------------------------------------------
+
+function nc_write_location_atts( ncFileID, fname, ObsNumDimID ) result (ierr)
+
+! Writes the "location module" -specific attributes to a netCDF file.
+
+use typeSizes
+use netcdf
+
+integer, intent(in) :: ncFileID ! handle to the netcdf file
+character(len=*), intent(in) :: fname ! file name (for printing purposes)
+integer, intent(in) :: ObsNumDimID ! handle to the dimension that grows
+integer :: ierr
+
+integer :: LocDimID
+integer :: VarID
+
+if ( .not. module_initialized ) call initialize_module
+
+ierr = -1 ! assume things will fail ...
+
+! define the rank/dimension of the location information
+call nc_check(nf90_def_dim(ncid=ncFileID, name='location', len=LocationDims, &
+ dimid = LocDimID), 'nc_write_location_atts', 'def_dim:location '//trim(fname))
+
+! Define the observation location variable and attributes
+
+call nc_check(nf90_def_var(ncid=ncFileID, name='location', xtype=nf90_double, &
+ dimids=(/ LocDimID, ObsNumDimID /), varid=VarID), &
+ 'nc_write_location_atts', 'location:def_var')
+
+call nc_check(nf90_put_att(ncFileID, VarID, 'description', &
+ 'location coordinates'), 'nc_write_location_atts', 'location:description')
+call nc_check(nf90_put_att(ncFileID, VarID, 'location_type', &
+ trim(LocationName)), 'nc_write_location_atts', 'location:location_type')
+call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', &
+ trim(LocationLName)), 'nc_write_location_atts', 'location:long_name')
+call nc_check(nf90_put_att(ncFileID, VarID, 'storage_order', &
+ 'X Y Z'), 'nc_write_location_atts', 'location:storage_order')
+call nc_check(nf90_put_att(ncFileID, VarID, 'units', &
+ 'X Y Z'), 'nc_write_location_atts', 'location:units')
+
+ierr = 0
+
+end function nc_write_location_atts
+
+!----------------------------------------------------------------------------
+
+subroutine nc_get_location_varids( ncFileID, fname, LocationVarID, WhichVertVarID )
+
+! Return the LocationVarID and WhichVertVarID variables from a given netCDF file.
+!
+! ncFileId the netcdf file descriptor
+! fname the name of the netcdf file (for error messages only)
+! LocationVarID the integer ID of the 'location' variable in the netCDF file
+! WhichVertVarID the integer ID of the 'which_vert' variable in the netCDF file
+
+use typeSizes
+use netcdf
+
+integer, intent(in) :: ncFileID ! handle to the netcdf file
+character(len=*), intent(in) :: fname ! file name (for printing purposes)
+integer, intent(out) :: LocationVarID, WhichVertVarID
+
+if ( .not. module_initialized ) call initialize_module
+
+call nc_check(nf90_inq_varid(ncFileID, 'location', varid=LocationVarID), &
+ 'nc_get_location_varids', 'inq_varid:location '//trim(fname))
+
+WhichVertVarID = -1
+
+end subroutine nc_get_location_varids
+
+!----------------------------------------------------------------------------
+
+subroutine nc_write_location(ncFileID, LocationVarID, loc, obsindex, WhichVertVarID)
+
+! Writes a SINGLE location to the specified netCDF variable and file.
+! The LocationVarID and WhichVertVarID must be the values returned from
+! the nc_get_location_varids call.
+
+use typeSizes
+use netcdf
+
+integer, intent(in) :: ncFileID, LocationVarID
+type(location_type), intent(in) :: loc
+integer, intent(in) :: obsindex
+integer, intent(in) :: WhichVertVarID
+
+real(r8), dimension(LocationDims) :: locations
+integer, dimension(1) :: intval
+
+if ( .not. module_initialized ) call initialize_module
+
+locations = get_location( loc )
+
+call nc_check(nf90_put_var(ncFileID, LocationVarId, locations, &
+ start=(/ 1, obsindex /), count=(/ LocationDims, 1 /) ), &
+ 'nc_write_location', 'put_var:location')
+
+end subroutine nc_write_location
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_obs_init(gc, num, obs)
+
+! Initializes part of get_close accelerator that depends on the particular obs
+
+type(get_close_type), intent(inout) :: gc
+integer, intent(in) :: num
+type(location_type), intent(in) :: obs(num)
+
+if (use_octree) then
+ call get_close_init_otree(gc, num, obs)
+else
+ call get_close_init_boxes(gc, num, obs)
+endif
+
+end subroutine get_close_obs_init
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_init_boxes(gc, num, obs)
+
+! Initializes part of get_close accelerator that depends on the particular obs
+
+type(get_close_type), intent(inout) :: gc
+integer, intent(in) :: num
+type(location_type), intent(in) :: obs(num)
+
+integer :: i, j, k, cum_start, l
+integer :: x_box(num), y_box(num), z_box(num)
+integer :: tstart(nx, ny, nz)
+
+if ( .not. module_initialized ) call initialize_module
+
+! Allocate storage for obs number dependent part
+allocate(gc%box%obs_box(num))
+gc%box%obs_box(:) = -1
+
+! Set the value of num_obs in the structure
+gc%num = num
+
+! If num == 0, no point in going any further.
+if (num == 0) return
+
+! FIXME: compute nx, ny, nz from nboxes? or put in namelist
+nx = nint(real(nboxes, r8)**0.33333) ! roughly cube root
+ny = nint(real(nboxes, r8)**0.33333) ! roughly cube root
+nz = nint(real(nboxes, r8) / real(nx * ny, r8)) ! whatever is left
+
+! Determine where the boxes should be for this set of obs and maxdist
+call find_box_ranges(gc, obs, num)
+
+! Begin by computing the number of observations in each box in x,y,z
+gc%box%count = 0
+do i = 1, num
+
+!print *, i, obs(i)%x, obs(i)%y, obs(i)%z
+ x_box(i) = floor((obs(i)%x - gc%box%bot_x) / gc%box%x_width) + 1
+ if(x_box(i) > nx) x_box(i) = nx
+ if(x_box(i) < 1) x_box(i) = 1
+
+ y_box(i) = floor((obs(i)%y - gc%box%bot_y) / gc%box%y_width) + 1
+ if(y_box(i) > ny) y_box(i) = ny
+ if(y_box(i) < 1) y_box(i) = 1
+
+ z_box(i) = floor((obs(i)%z - gc%box%bot_z) / gc%box%z_width) + 1
+ if(z_box(i) > nz) z_box(i) = nz
+ if(z_box(i) < 1) z_box(i) = 1
+
+ gc%box%count(x_box(i), y_box(i), z_box(i)) = gc%box%count(x_box(i), y_box(i), z_box(i)) + 1
+!print *, 'adding count to box ', x_box(i), y_box(i), z_box(i), &
+! gc%box%count(x_box(i), y_box(i), z_box(i))
+end do
+
+! Figure out where storage for each boxes members should begin
+cum_start = 1
+do i = 1, nx
+ do j = 1, ny
+ do k = 1, nz
+ gc%box%start(i, j, k) = cum_start
+ cum_start = cum_start + gc%box%count(i, j, k)
+ end do
+ end do
+end do
+
+! Now we know how many are in each box, get a list of which are in each box
+tstart = gc%box%start
+do i = 1, num
+ gc%box%obs_box(tstart(x_box(i), y_box(i), z_box(i))) = i
+ tstart(x_box(i), y_box(i), z_box(i)) = tstart(x_box(i), y_box(i), z_box(i)) + 1
+end do
+
+do i = 1, nx
+ do j = 1, ny
+ do k = 1, nz
+if (gc%box%count(i,j,k) > 0) print *, i,j,k, gc%box%count(i,j,k), gc%box%start(i,j,k)
+ do l=1, gc%box%count(i,j,k)
+!print *, l, gc%box%obs_box(l)
+ enddo
+ end do
+ end do
+end do
+
+
+! info on how well the boxes are working. by default print nothing.
+! set print_box_level to higher values to get more and more detail.
+! user info should be level 1; 2 and 3 should be for debug only.
+! special for grid-decomposition debugging; set print level to -8.
+if (output_box_info) then
+ ! if this task normally prints, call the print routine.
+ ! 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)
+ 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 set_output(.false.)
+ endif
+endif
+
+end subroutine get_close_init_boxes
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_init_otree(gc, num, locs)
+
+! Octree version
+
+! Initializes part of get_close accelerator that depends on the particular obs
+
+type(get_close_type), intent(inout), target :: gc
+integer, intent(in) :: num
+type(location_type), intent(in) :: locs(num)
+
+integer :: i, j, k, n
+
+type(octree_type), pointer :: r, c
+real(r8) :: xl, xu, yl, yu, zl, zu
+
+if ( .not. module_initialized ) call initialize_module
+
+r => gc%root
+
+! Set the value of num_obs in the structure
+gc%num = num
+
+! If num == 0, no point in going any further.
+if (num == 0) return
+
+r%count = num
+
+! need to include space outside the limits of the initialization set,
+! so points outside boundary but closer than maxdist will still match.
+r%llb = set_location(minval(locs(:)%x)-gc%maxdist, &
+ minval(locs(:)%y)-gc%maxdist, &
+ minval(locs(:)%z)-gc%maxdist)
+r%urt = set_location(maxval(locs(:)%x)+gc%maxdist, &
+ maxval(locs(:)%y)+gc%maxdist, &
+ maxval(locs(:)%z)+gc%maxdist)
+! for now, split is midpoint in each dim. this does NOT have to be the case.
+r%split = set_location((r%urt%x + r%llb%x) / 2.0_r8, &
+ (r%urt%y + r%llb%y) / 2.0_r8, &
+ (r%urt%z + r%llb%z) / 2.0_r8)
+
+! initially everyone is in the original list
+allocate(r%index(num))
+do i=1,num
+ r%index(i) = i
+enddo
+
+! recursion starts here
+if (r%count > filled) then
+ call split_tree(r)
+ call move_to_children(r, locs)
+endif
+
+! info on how well the boxes are working. by default print nothing.
+! set print_box_level to higher values to get more and more detail.
+! user info should be level 1; 2 and 3 should be for debug only.
+! special for grid-decomposition debugging; set print level to -8.
+if (output_box_info) then
+ ! if this task normally prints, call the print routine.
+ ! 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)
+ 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 set_output(.false.)
+ endif
+endif
+
+end subroutine get_close_init_otree
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_obs_destroy(gc)
+
+type(get_close_type), intent(inout) :: gc
+
+if (use_octree) then
+ call get_close_destroy_otree(gc)
+else
+ call get_close_destroy_boxes(gc)
+endif
+
+end subroutine get_close_obs_destroy
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_destroy_boxes(gc)
+
+type(get_close_type), intent(inout) :: gc
+
+deallocate(gc%box%obs_box, gc%box%count, gc%box%start)
+
+end subroutine get_close_destroy_boxes
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_destroy_otree(gc)
+
+type(get_close_type), intent(inout) :: gc
+
+! FIXME: do a depth-first search and deallocate
+! the index() arrays, then deallocate the octree structs
+! one by one.
+
+! gc%root -> children until unallocated
+
+end subroutine get_close_destroy_otree
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_maxdist_init(gc, maxdist)
+
+type(get_close_type), intent(inout) :: gc
+real(r8), intent(in) :: maxdist
+
+character(len=129) :: str1
+integer :: i
+
+! set the default value.
+gc%maxdist = maxdist
+!print *, 'setting maxdist to ', maxdist
+
+if (.not. use_octree) then
+ ! Allocate the storage for the grid dependent boxes
+ allocate(gc%box%count(nx,ny,nz), gc%box%start(nx,ny,nz))
+ gc%box%count = -1
+ gc%box%start = -1
+endif
+
+end subroutine get_close_maxdist_init
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_obs(gc, base_obs_loc, base_obs_kind, obs, obs_kind, &
+ num_close, close_ind, dist)
+
+! FIXME: these work on any locations. the names of these args should be:
+! gc, base_loc, base_kind, locs, locs_kinds, ...
+
+type(get_close_type), intent(in) :: gc
+type(location_type), intent(in) :: base_obs_loc, obs(:)
+integer, intent(in) :: base_obs_kind, obs_kind(:)
+integer, intent(out) :: num_close, close_ind(:)
+real(r8), optional, intent(out) :: dist(:)
+
+if (use_octree) then
+ call get_close_otree(gc, base_obs_loc, base_obs_kind, obs, obs_kind, &
+ num_close, close_ind, dist)
+else
+ call get_close_boxes(gc, base_obs_loc, base_obs_kind, obs, obs_kind, &
+ num_close, close_ind, dist)
+endif
+
+end subroutine get_close_obs
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_boxes(gc, base_obs_loc, base_obs_kind, obs, obs_kind, &
+ num_close, close_ind, dist)
+
+! FIXME: these work on any locations. the names of these args should be:
+! gc, base_loc, base_kind, locs, locs_kinds, ...
+
+type(get_close_type), intent(in) :: gc
+type(location_type), intent(in) :: base_obs_loc, obs(:)
+integer, intent(in) :: base_obs_kind, obs_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 :: x_box, y_box, z_box, i, j, k, l
+integer :: start_x, end_x, start_y, end_y, start_z, end_z
+integer :: n_in_box, st, t_ind
+real(r8) :: this_dist, this_maxdist
+
+! Variables needed for comparing against correct case.
+! these could be large - make them allocatable
+! and only allocate them if needed.
+integer :: cnum_close
+integer, allocatable :: cclose_ind(:)
+real(r8), allocatable :: cdist(:)
+
+! First, set the intent out arguments to a missing value
+num_close = 0
+close_ind = -99
+if(present(dist)) dist = -1e38_r8 ! big but negative
+this_dist = 1e38_r8 ! something big and positive.
+
+! the list of locations in the obs() argument must be the same
+! as the list of locations passed into get_close_obs_init(), so
+! gc%num and size(obs) better be the same. if the list changes,
+! you have to destroy the old gc and init a new one.
+if (size(obs) /= gc%num) then
+ write(errstring,*)'obs() array must match one passed to get_close_obs_init()'
+ call error_handler(E_ERR, 'get_close_obs', errstring, source, revision, revdate)
+endif
+
+! If num == 0, no point in going any further.
+if (gc%num == 0) return
+
+this_maxdist = gc%maxdist
+
+
+! For validation, it is useful to be able to compare against exact
+! exhaustive search
+if(compare_to_correct) then
+ call exhaustive_collect(gc, base_obs_loc, obs, &
+ cnum_close, cclose_ind, cdist)
+endif
+
+
+! Begin by figuring out which box the base loc is in
+x_box = floor((base_obs_loc%x - gc%box%bot_x) / gc%box%x_width) + 1
+y_box = floor((base_obs_loc%y - gc%box%bot_y) / gc%box%y_width) + 1
+z_box = floor((base_obs_loc%z - gc%box%bot_z) / gc%box%z_width) + 1
+
+! If it is not in any box, then it is more than the maxdist away from everybody
+if(x_box > nx .or. x_box < 1 .or. x_box < 0) return
+if(y_box > ny .or. y_box < 1 .or. y_box < 0) return
+if(z_box > nz .or. z_box < 1 .or. z_box < 0) return
+
+! figure out how many boxes need searching
+! FIXME: if we support a variable maxdist, nboxes_X will need to
+! be computed on the fly here instead of precomputed at init time.
+!print *, 'nboxes x, y, z = ', gc%box%nboxes_x, gc%box%nboxes_y, gc%box%nboxes_z
+!print *, 'base_loc in box ', x_box, y_box, z_box
+
+start_x = x_box - gc%box%nboxes_x
+if (start_x < 1) start_x = 1
+end_x = x_box + gc%box%nboxes_x
+if (end_x > nx) end_x = nx
+
+start_y = y_box - gc%box%nboxes_y
+if (start_y < 1) start_y = 1
+end_y = y_box + gc%box%nboxes_y
+if (end_y > ny) end_y = ny
+
+start_z = z_box - gc%box%nboxes_z
+if (start_z < 1) start_z = 1
+end_z = z_box + gc%box%nboxes_z
+if (end_z > nz) end_z = nz
+
+!print *, 'looping from '
+!print *, 'x: ', start_x, end_x
+!print *, 'y: ', start_y, end_y
+!print *, 'z: ', start_z, end_z
+
+! Next, loop through each box that is close to this box
+do i = start_x, end_x
+ do j = start_y, end_y
+ do k = start_z, end_z
+
+ ! Box to search is i,j,k
+ n_in_box = gc%box%count(i, j, k)
+ st = gc%box%start(i,j,k)
+
+
+ ! Loop to check how close all obs in the box are; add those that are close
+ do l = 1, n_in_box
+
+ t_ind = gc%box%obs_box(st - 1 + l)
+!print *, 'l, t_ind = ', l, t_ind
+
+ ! Only compute distance if dist is present
+ if(present(dist)) then
+ this_dist = get_dist(base_obs_loc, obs(t_ind))
+!print *, 'this_dist = ', this_dist
+ ! If this obs' distance is less than cutoff, add it in list
+ if(this_dist <= this_maxdist) then
+ num_close = num_close + 1
+ close_ind(num_close) = t_ind
+ dist(num_close) = this_dist
+ endif
+ else
+ ! Dist isn't present; add this ob to list without computing distance
+ num_close = num_close + 1
+ close_ind(num_close) = t_ind
+ endif
+
+ end do
+ end do
+ end do
+end do
+
+
+! Verify by comparing to exhaustive search
+if(compare_to_correct) then
+ call exhaustive_report(cnum_close, num_close, cclose_ind, close_ind, cdist, dist)
+endif
+
+
+end subroutine get_close_boxes
+
+!----------------------------------------------------------------------------
+
+subroutine get_close_otree(gc, base_obs_loc, base_obs_kind, obs, obs_kind, &
+ num_close, close_ind, dist)
+
+! FIXME: these work on any locations. the names of these args should be:
+! gc, base_loc, base_kind, locs, locs_kinds, ...
+
+type(get_close_type), intent(in), target :: gc
+type(location_type), intent(in) :: base_obs_loc, obs(:)
+integer, intent(in) :: base_obs_kind, obs_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
+
+type(octree_type), pointer :: r, c
+
+integer :: x_box, y_box, z_box, i, j, k, l, n
+integer :: start_x, end_x, start_y, end_y, start_z, end_z
+integer :: n_in_box, st, t_ind
+real(r8) :: this_dist
+
+! Variables needed for comparing against correct case.
+! these could be large - make them allocatable
+! and only allocate them if needed.
+integer :: cnum_close
+integer, allocatable :: cclose_ind(:)
+real(r8), allocatable :: cdist(:)
+
+! First, set the intent out arguments to a missing value
+num_close = 0
+close_ind = -99
+if(present(dist)) dist = -1e38_r8 ! big but negative
+this_dist = 1e38_r8 ! something big and positive.
+
+! the list of locations in the obs() argument must be the same
+! as the list of locations passed into get_close_obs_init(), so
+! gc%num and size(obs) better be the same. if the list changes,
+! you have to destroy the old gc and init a new one.
+if (size(obs) /= gc%num) then
+ write(errstring,*)'obs() array must match one passed to get_close_obs_init()'
+ call error_handler(E_ERR, 'get_close_obs', errstring, source, revision, revdate)
+endif
+
+! If num == 0, no point in going any further.
+if (gc%num == 0) return
+
+
+! For validation, it is useful to be able to compare against exact
+! exhaustive search
+if(compare_to_correct) then
+ call exhaustive_collect(gc, base_obs_loc, obs, &
+ cnum_close, cclose_ind, cdist)
+endif
+
+
+! revised plan:
+! find min/max extents in each dim, +/- maxdist from base
+! and only descend lower in a branch if the region intersects
+! the range that this branch covers. maxdist doesn't have to
+! be a constant in this case - it can be specified per search.
+
+call collect_nearby(gc%root, base_obs_loc, gc%maxdist, obs, base_obs_kind, obs_kind, &
+ num_close, close_ind, dist)
+
+
+! Verify by comparing to exhaustive search
+if(compare_to_correct) then
+ call exhaustive_report(cnum_close, num_close, cclose_ind, close_ind, cdist, dist)
+endif
+
+
+end subroutine get_close_otree
+
+!--------------------------------------------------------------------------
+
+subroutine find_box_ranges(gc, locs, num)
+
+! Finds boundaries for x,y,z boxes.
+! FIXME: ways boxes could be divided:
+! - evenly along each axis
+! - octree-like, divide each axis so roughly half the points are
+! on each side of the dividing plane.
+! - about 100 other schemes
+
+type(get_close_type), intent(inout) :: gc
+integer, intent(in) :: num
+type(location_type), intent(in) :: locs(num)
+
+logical :: old_out
+
+
+! FIXME: this space could be very sparse
+
+gc%box%bot_x = minval(locs(:)%x)
+gc%box%bot_y = minval(locs(:)%y)
+gc%box%bot_z = minval(locs(:)%z)
+
+gc%box%top_x = maxval(locs(:)%x)
+gc%box%top_y = maxval(locs(:)%y)
+gc%box%top_z = maxval(locs(:)%z)
+
+gc%box%x_width = (gc%box%top_x - gc%box%bot_x) / nx
+gc%box%y_width = (gc%box%top_y - gc%box%bot_y) / ny
+gc%box%z_width = (gc%box%top_z - gc%box%bot_z) / nz
+
+! FIXME: compute a sphere of radius maxdist and see how
+! many boxes in x, y, z that would include.
+gc%box%nboxes_x = aint((gc%maxdist + (gc%box%x_width-1)) / gc%box%x_width)
+gc%box%nboxes_y = aint((gc%maxdist + (gc%box%y_width-1)) / gc%box%y_width)
+gc%box%nboxes_z = aint((gc%maxdist + (gc%box%z_width-1)) / gc%box%z_width)
+
+
+!if(compare_to_correct) then
+! old_out = do_output()
+! call set_output(.true.)
+! write(errstring, *) 'x bot, top, width, nboxes ', gc%box%bot_x, gc%box%top_x, gc%box%x_width, gc%box%nboxes_x
+! call error_handler(E_MSG, 'find_box_ranges', errstring)
+! write(errstring, *) 'y bot, top, width, nboxes ', gc%box%bot_y, gc%box%top_y, gc%box%y_width, gc%box%nboxes_y
+! call error_handler(E_MSG, 'find_box_ranges', errstring)
+! write(errstring, *) 'z bot, top, width, nboxes ', gc%box%bot_z, gc%box%top_z, gc%box%z_width, gc%box%nboxes_z
+! call error_handler(E_MSG, 'find_box_ranges', errstring)
+! call set_output(old_out)
+!endif
+
+end subroutine find_box_ranges
+
+!----------------------------------------------------------------------------
+
+recursive subroutine split_tree(r)
+ type(octree_type), pointer :: r
+
+integer :: i, j, k
+real(r8) :: xl, xu, yl, yu, zl, zu
+type(octree_type), pointer :: c
+
+
+allocate(r%children(2,2,2))
+
+do i=1,2
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list