[Dart-dev] [7756] DART/trunk/location/threed_sphere/location_mod.f90: last commit of this batch; simply rearrange the order
nancy at ucar.edu
nancy at ucar.edu
Tue Mar 24 17:36:35 MDT 2015
Revision: 7756
Author: nancy
Date: 2015-03-24 17:36:35 -0600 (Tue, 24 Mar 2015)
Log Message:
-----------
last commit of this batch; simply rearrange the order
of the subroutines. all routines that deal only with the
locations structure are first; all that deal with get_close
are last. the get_close routines now are in the same order
as their usage.
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 23:28:29 UTC (rev 7755)
+++ DART/trunk/location/threed_sphere/location_mod.f90 2015-03-24 23:36:35 UTC (rev 7756)
@@ -147,6 +147,17 @@
real(r8) :: my_cos(-SINCOS_LIMIT:SINCOS_LIMIT)
real(r8) :: my_acos(-ACOS_LIMIT:ACOS_LIMIT)
+! Tolerance for the top latitude boundary test. All locations which
+! are located on a box boundary are added to the bin on the larger
+! side of the boundary (e.g. for north-south latitudes, it rounds
+! towards the north). But when a value falls exactly on the edge
+! of the last box, technically it is inside the region but would be
+! rounded up and outside the region unless handled specially.
+! this tolerance below is used to determine if a value is within
+! the range of the last box boundary and if so, the location is
+! included in the last box. In particular, for global grids this
+! preserves locations which are at exactly 90.0 degrees latitude.
+real(r8), parameter :: EDGE_TOLERANCE = 100.0_r8 * epsilon(0.0_r8)
! Option for verification using exhaustive search
logical :: COMPARE_TO_CORRECT = .false. ! normally false
@@ -924,7 +935,6 @@
maxlon = maxlon * DEG2RAD
! Longitude is random from minlon to maxlon
-! location%lon = random_uniform(ran_seq) * 2.0_r8 * PI
location%lon = random_uniform(ran_seq) * (maxlon-minlon) + minlon
write(*, *) 'Input minimum latitude (-90.0 to 90.0)'
@@ -936,7 +946,6 @@
maxlat = sin(maxlat * DEG2RAD)
! Latitude must be area weighted
-! location%lat = asin(random_uniform(ran_seq) * 2.0_r8 - 1.0_r8)
location%lat = asin(random_uniform(ran_seq) * (maxlat-minlat) + minlat)
write(*, *) 'random location is ', location%lon / DEG2RAD, &
@@ -1091,11 +1100,238 @@
end subroutine nc_write_location
!----------------------------------------------------------------------------
+
+function is_location_in_region(loc, minl, maxl)
+
+! Returns true if the given location is inside the rectangular
+! region defined by minl as the lower left, maxl the upper right.
+! test is inclusive; values on the edges are considered inside.
+! Periodic in longitude (box can cross the 2PI -> 0 line)
+
+logical :: is_location_in_region
+type(location_type), intent(in) :: loc, minl, maxl
+
+if ( .not. module_initialized ) call initialize_module()
+
+! maybe could use VERTISUNDEF in the minl and maxl args to indicate
+! we want to test only in horizontal? and if not, vtypes must match?
+!if ( (minl%which_vert /= maxl%which_vert) .or. &
+! ((minl%which_vert /= loc%which_vert).and.(minl%which_vert /= VERTISUNDEF))) then
+! write(msgstring,*)'which_vert (',loc%which_vert,') must be same in all args'
+! call error_handler(E_ERR, 'is_location_in_region', msgstring, source, revision, revdate)
+!endif
+
+! assume failure and return as soon as we are confirmed right.
+! set to success only at the bottom after all tests have passed.
+is_location_in_region = .false.
+
+! latitude: we do not allow wrap of rectangular regions over the poles.
+if ((loc%lat < minl%lat) .or. (loc%lat > maxl%lat)) return
+
+! use common routine in utilities module to do all the wrapping
+if (.not. is_longitude_between(loc%lon, minl%lon, maxl%lon, doradians=.TRUE.)) return
+
+! once we decide what to do about diff vert units, this is the test.
+!if ((minl%which_vert .ne. VERTISUNDEF) .and.
+! (loc%vloc < minl%vloc) .or. (loc%vloc > maxl%vloc)) return
+
+is_location_in_region = .true.
+
+end function is_location_in_region
+
+!---------------------------------------------------------------------------
+
+function vert_is_undef(loc)
+
+! Given a location, return true if vertical coordinate is undefined, else false
+
+logical :: vert_is_undef
+type(location_type), intent(in) :: loc
+
+if ( .not. module_initialized ) call initialize_module()
+
+if(loc%which_vert == VERTISUNDEF) then
+ vert_is_undef = .true.
+else
+ vert_is_undef = .false.
+endif
+
+end function vert_is_undef
+
+!---------------------------------------------------------------------------
+
+function vert_is_surface(loc)
+
+! Given a location, return true if vertical coordinate is surface, else false
+
+logical :: vert_is_surface
+type(location_type), intent(in) :: loc
+
+if ( .not. module_initialized ) call initialize_module()
+
+if(loc%which_vert == VERTISSURFACE) then
+ vert_is_surface = .true.
+else
+ vert_is_surface = .false.
+endif
+
+end function vert_is_surface
+
+!---------------------------------------------------------------------------
+
+function vert_is_pressure(loc)
+
+! Given a location, return true if vertical coordinate is pressure, else false
+
+logical :: vert_is_pressure
+type(location_type), intent(in) :: loc
+
+if ( .not. module_initialized ) call initialize_module()
+
+if(loc%which_vert == VERTISPRESSURE) then
+ vert_is_pressure = .true.
+else
+ vert_is_pressure = .false.
+endif
+
+end function vert_is_pressure
+
+!---------------------------------------------------------------------------
+
+function vert_is_height(loc)
+
+! Given a location, return true if vertical coordinate is height, else false
+
+logical :: vert_is_height
+type(location_type), intent(in) :: loc
+
+if ( .not. module_initialized ) call initialize_module()
+
+if(loc%which_vert == VERTISHEIGHT ) then
+ vert_is_height = .true.
+else
+ vert_is_height = .false.
+endif
+
+end function vert_is_height
+
+!---------------------------------------------------------------------------
+
+function vert_is_level(loc)
+
+! Given a location, return true if vertical coordinate is level, else false
+
+logical :: vert_is_level
+type(location_type), intent(in) :: loc
+
+if ( .not. module_initialized ) call initialize_module()
+
+if(loc%which_vert == VERTISLEVEL) then
+ vert_is_level = .true.
+else
+ vert_is_level = .false.
+endif
+
+end function vert_is_level
+
+!---------------------------------------------------------------------------
+
+function vert_is_scale_height(loc)
+
+! Given a location, return true if vertical coordinate is scale height, else false
+
+logical :: vert_is_scale_height
+type(location_type), intent(in) :: loc
+
+if ( .not. module_initialized ) call initialize_module()
+
+if(loc%which_vert == VERTISSCALEHEIGHT ) then
+ vert_is_scale_height = .true.
+else
+ vert_is_scale_height = .false.
+endif
+
+end function vert_is_scale_height
+
+!---------------------------------------------------------------------------
+
+function has_vertical_localization()
+
+! Return the (opposite) namelist setting for horiz_dist_only.
+
+logical :: has_vertical_localization
+
+if ( .not. module_initialized ) call initialize_module()
+
+has_vertical_localization = .not. horiz_dist_only
+
+end function has_vertical_localization
+
+
!----------------------------------------------------------------------------
+!----------------------------------------------------------------------------
! get close routines
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
+subroutine get_close_maxdist_init(gc, maxdist, maxdist_list)
+
+type(get_close_type), intent(inout) :: gc
+real(r8), intent(in) :: maxdist
+real(r8), intent(in), optional :: maxdist_list(:)
+
+integer :: i, typecount, distcount, mapnum
+real(r8), allocatable :: distlist(:)
+
+! 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))
+
+if (present(maxdist_list)) then
+ 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(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
+
+! make a gtt type array for each different cutoff distance
+! (just 1 type is the most common case.)
+allocate(gc%gtt(gc%nt))
+
+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_init(gc, num, locs)
@@ -1269,66 +1505,6 @@
!----------------------------------------------------------------------------
-subroutine get_close_maxdist_init(gc, maxdist, maxdist_list)
-
-type(get_close_type), intent(inout) :: gc
-real(r8), intent(in) :: maxdist
-real(r8), intent(in), optional :: maxdist_list(:)
-
-integer :: i, typecount, distcount, mapnum
-real(r8), allocatable :: distlist(:)
-
-! 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))
-
-if (present(maxdist_list)) then
- 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(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
-
-! make a gtt type array for each different cutoff distance
-! (just 1 type is the most common case.)
-allocate(gc%gtt(gc%nt))
-
-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, locs, loc_kind, &
num_close, close_ind, dist)
@@ -2292,175 +2468,6 @@
end subroutine print_get_close_type
!----------------------------------------------------------------------------
-
-function is_location_in_region(loc, minl, maxl)
-
-! Returns true if the given location is inside the rectangular
-! region defined by minl as the lower left, maxl the upper right.
-! test is inclusive; values on the edges are considered inside.
-! Periodic in longitude (box can cross the 2PI -> 0 line)
-
-logical :: is_location_in_region
-type(location_type), intent(in) :: loc, minl, maxl
-
-if ( .not. module_initialized ) call initialize_module()
-
-! maybe could use VERTISUNDEF in the minl and maxl args to indicate
-! we want to test only in horizontal? and if not, vtypes must match?
-!if ( (minl%which_vert /= maxl%which_vert) .or. &
-! ((minl%which_vert /= loc%which_vert).and.(minl%which_vert /= VERTISUNDEF))) then
-! write(msgstring,*)'which_vert (',loc%which_vert,') must be same in all args'
-! call error_handler(E_ERR, 'is_location_in_region', msgstring, source, revision, revdate)
-!endif
-
-! assume failure and return as soon as we are confirmed right.
-! set to success only at the bottom after all tests have passed.
-is_location_in_region = .false.
-
-! latitude: we do not allow wrap of rectangular regions over the poles.
-if ((loc%lat < minl%lat) .or. (loc%lat > maxl%lat)) return
-
-! use common routine in utilities module to do all the wrapping
-if (.not. is_longitude_between(loc%lon, minl%lon, maxl%lon, doradians=.TRUE.)) return
-
-! once we decide what to do about diff vert units, this is the test.
-!if ((minl%which_vert .ne. VERTISUNDEF) .and.
-! (loc%vloc < minl%vloc) .or. (loc%vloc > maxl%vloc)) return
-
-is_location_in_region = .true.
-
-end function is_location_in_region
-
-!---------------------------------------------------------------------------
-
-function vert_is_undef(loc)
-
-! Given a location, return true if vertical coordinate is undefined, else false
-
-logical :: vert_is_undef
-type(location_type), intent(in) :: loc
-
-if ( .not. module_initialized ) call initialize_module()
-
-if(loc%which_vert == VERTISUNDEF) then
- vert_is_undef = .true.
-else
- vert_is_undef = .false.
-endif
-
-end function vert_is_undef
-
-!---------------------------------------------------------------------------
-
-function vert_is_surface(loc)
-
-! Given a location, return true if vertical coordinate is surface, else false
-
-logical :: vert_is_surface
-type(location_type), intent(in) :: loc
-
-if ( .not. module_initialized ) call initialize_module()
-
-if(loc%which_vert == VERTISSURFACE) then
- vert_is_surface = .true.
-else
- vert_is_surface = .false.
-endif
-
-end function vert_is_surface
-
-!---------------------------------------------------------------------------
-
-function vert_is_pressure(loc)
-
-! Given a location, return true if vertical coordinate is pressure, else false
-
-logical :: vert_is_pressure
-type(location_type), intent(in) :: loc
-
-if ( .not. module_initialized ) call initialize_module()
-
-if(loc%which_vert == VERTISPRESSURE) then
- vert_is_pressure = .true.
-else
- vert_is_pressure = .false.
-endif
-
-end function vert_is_pressure
-
-!---------------------------------------------------------------------------
-
-function vert_is_height(loc)
-
-! Given a location, return true if vertical coordinate is height, else false
-
-logical :: vert_is_height
-type(location_type), intent(in) :: loc
-
-if ( .not. module_initialized ) call initialize_module()
-
-if(loc%which_vert == VERTISHEIGHT ) then
- vert_is_height = .true.
-else
- vert_is_height = .false.
-endif
-
-end function vert_is_height
-
-!---------------------------------------------------------------------------
-
-function vert_is_level(loc)
-
-! Given a location, return true if vertical coordinate is level, else false
-
-logical :: vert_is_level
-type(location_type), intent(in) :: loc
-
-if ( .not. module_initialized ) call initialize_module()
-
-if(loc%which_vert == VERTISLEVEL) then
- vert_is_level = .true.
-else
- vert_is_level = .false.
-endif
-
-end function vert_is_level
-
-!---------------------------------------------------------------------------
-
-function vert_is_scale_height(loc)
-
-! Given a location, return true if vertical coordinate is scale height, else false
-
-logical :: vert_is_scale_height
-type(location_type), intent(in) :: loc
-
-if ( .not. module_initialized ) call initialize_module()
-
-if(loc%which_vert == VERTISSCALEHEIGHT ) then
- vert_is_scale_height = .true.
-else
- vert_is_scale_height = .false.
-endif
-
-end function vert_is_scale_height
-
-!---------------------------------------------------------------------------
-
-function has_vertical_localization()
-
-! Return the (opposite) namelist setting for horiz_dist_only.
-
-logical :: has_vertical_localization
-
-if ( .not. module_initialized ) call initialize_module()
-
-has_vertical_localization = .not. horiz_dist_only
-
-end function has_vertical_localization
-
-
-!----------------------------------------------------------------------------
! end of location/threed_sphere/location_mod.f90
!----------------------------------------------------------------------------
More information about the Dart-dev
mailing list