[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