[Dart-dev] [3809] DART/trunk: Fixed errors in the obs_sequence_tool, the ncep observation
nancy at ucar.edu
nancy at ucar.edu
Mon Apr 13 10:21:35 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090413/0d67d824/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/diagnostics/threed_sphere/obs_diag.f90
===================================================================
--- DART/trunk/diagnostics/threed_sphere/obs_diag.f90 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/diagnostics/threed_sphere/obs_diag.f90 2009-04-13 16:21:33 UTC (rev 3809)
@@ -30,8 +30,9 @@
get_obs_def_location, get_obs_kind, get_obs_name
use obs_kind_mod, only : max_obs_kinds, get_obs_kind_var_type, get_obs_kind_name, &
KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT
-use location_mod, only : location_type, get_location, set_location_missing, &
- write_location, operator(/=), &
+use location_mod, only : location_type, get_location, set_location_missing, &
+ write_location, operator(/=), is_location_in_region, &
+ set_location, &
vert_is_undef, VERTISUNDEF, &
vert_is_surface, VERTISSURFACE, &
vert_is_level, VERTISLEVEL, &
@@ -165,6 +166,7 @@
real(r8), dimension(MaxRegions) :: lonlim1= MISSING_R8, lonlim2= MISSING_R8
real(r8), dimension(MaxRegions) :: latlim1= MISSING_R8, latlim2= MISSING_R8
character(len = stringlength), dimension(MaxRegions) :: reg_names = 'null'
+type(location_type), dimension(MaxRegions) :: min_loc, max_loc
real(r8):: rat_cri = 5000.0_r8 ! QC ratio
real(r8):: input_qc_threshold = 4.0_r8 ! maximum NCEP QC factor
@@ -331,6 +333,9 @@
call SetScaleFactors(scale_factor, logfileunit) ! for plotting purposes
+call SetRegionLimits(Nregions, lonlim1, lonlim2, latlim1, latlim2, &
+ min_loc, max_loc)
+
if (verbose) then
write(*,*)'pressure levels = ',plevel( 1:Nplevels)
write(*,*)'pressure interfaces = ',plevel_edges(1:Nplevels+1)
@@ -926,8 +931,7 @@
Areas : do iregion =1, Nregions
- keeper = InRegion( obslon, obslat, lonlim1(iregion), lonlim2(iregion), &
- latlim1(iregion), latlim2(iregion))
+ keeper = is_location_in_region( obs_loc, min_loc(iregion), max_loc(iregion) )
if ( .not. keeper ) cycle Areas
!-----------------------------------------------------------
@@ -2119,29 +2123,29 @@
- Function InRegion( lon, lat, lon1, lon2, lat1, lat2 )
- ! InRegion handles regions that 'wrap' in longitude
- ! if the easternmost longitude of the region is > 360.0
- ! For example; lon1 == 320, lon = 10, lon2 == 380 --> .true.
- !
- logical :: InRegion
- real(r8), intent(in) :: lon, lat, lon1, lon2, lat1, lat2
- real(r8) :: mylon
+ Subroutine SetRegionLimits(Nregions, lonlim1, lonlim2, latlim1, latlim2, &
+ min_loc, max_loc)
- InRegion = .false.
+ ! Set the min and max location_types for each region
+
+ integer, intent(in) :: Nregions
+ real(r8), dimension(*), intent(in) :: lonlim1, lonlim2, latlim1, latlim2
+ type(location_type), dimension(*), intent(out) :: min_loc, max_loc
- if( (lon >= lon1) .and. (lon <= lon2) .and. &
- (lat >= lat1) .and. (lat <= lat2) ) InRegion = .true.
+ integer :: i
+ real(r8) :: lon
- if( lon2 > 360.0_r8 ) then ! checking the wraparound case
- mylon = lon + 360.0_r8
- if( (mylon >= lon1) .and. (mylon <= lon2) .and. &
- ( lat >= lat1) .and. ( lat <= lat2) ) then
- InRegion = .true.
- endif
- endif
+ ! only set the horizontal limits; vertical will be ignored.
+ ! modulo() return is always positive, even if input is negative.
+ ! this is not true for mod().
+ do i = 1, Nregions
+ lon = modulo(lonlim1(i), 360.0_r8)
+ min_loc(i) = set_location(lon, latlim1(i), 0.0_r8, VERTISUNDEF)
+ lon = modulo(lonlim2(i), 360.0_r8)
+ max_loc(i) = set_location(lon, latlim2(i), 0.0_r8, VERTISUNDEF)
+ enddo
- end Function InRegion
+ end Subroutine SetRegionLimits
Modified: DART/trunk/location/threed_sphere/location_mod.f90
===================================================================
--- DART/trunk/location/threed_sphere/location_mod.f90 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/location/threed_sphere/location_mod.f90 2009-04-13 16:21:33 UTC (rev 3809)
@@ -27,7 +27,7 @@
use types_mod, only : r8, PI, RAD2DEG, DEG2RAD, MISSING_R8, MISSING_I
use utilities_mod, only : register_module, error_handler, E_ERR, E_MSG, &
logfileunit, nmlfileunit, find_namelist_in_file, &
- check_namelist_read, do_output
+ check_namelist_read, do_output, is_longitude_between
use random_seq_mod, only : random_seq_type, init_random_seq, random_uniform
implicit none
@@ -1678,20 +1678,24 @@
function is_location_in_region(loc, minl, maxl)
!----------------------------------------------------------------------------
!
-! Returns true if the given location is between the other two.
+! 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)
implicit none
logical :: is_location_in_region
type(location_type), intent(in) :: loc, minl, maxl
-
character(len=129) :: errstring
if ( .not. module_initialized ) call initialize_module
-!if ((minl%which_vert /= maxl%which_vert) .or. &
-! (minl%which_vert /= loc%which_vert)) then
+! 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(errstring,*)'which_vert (',loc%which_vert,') must be same in all args'
! call error_handler(E_ERR, 'is_location_in_region', errstring, source, revision, revdate)
!endif
@@ -1700,10 +1704,17 @@
! set to success only at the bottom after all tests have passed.
is_location_in_region = .false.
-if ((loc%lon < minl%lon) .or. (loc%lon > maxl%lon)) return
+! latitude: we do not allow wrap of rectangular regions over the poles.
if ((loc%lat < minl%lat) .or. (loc%lat > maxl%lat)) return
-!if ((loc%vloc < minl%vloc) .or. (loc%vloc > maxl%vloc)) 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
Modified: DART/trunk/ncep_obs/create_real_obs.f90
===================================================================
--- DART/trunk/ncep_obs/create_real_obs.f90 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/ncep_obs/create_real_obs.f90 2009-04-13 16:21:33 UTC (rev 3809)
@@ -97,10 +97,6 @@
if (do_output()) write(logfileunit, nml=ncepobs_nml)
if (do_output()) write( * , nml=ncepobs_nml)
-lon1 = min(max(lon1,0.0_r8),360.0_r8)
-lon2 = min(max(lon2,0.0_r8),360.0_r8)
-if ( lon1 > lon2 ) lon2 = lon2 + 360.0_r8
-
! Loop through the days interested.
do ii = 1, tot_days
Modified: DART/trunk/ncep_obs/real_obs_mod.f90
===================================================================
--- DART/trunk/ncep_obs/real_obs_mod.f90 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/ncep_obs/real_obs_mod.f90 2009-04-13 16:21:33 UTC (rev 3809)
@@ -21,8 +21,9 @@
get_date, set_time, GREGORIAN
use utilities_mod, only : get_unit, open_file, close_file, file_exist, &
register_module, error_handler, &
- E_ERR, E_MSG, timestamp
-use location_mod, only : location_type, set_location, VERTISPRESSURE, VERTISSURFACE
+ E_ERR, E_MSG, timestamp, is_longitude_between
+use location_mod, only : location_type, set_location, &
+ VERTISPRESSURE, VERTISSURFACE
use obs_sequence_mod, only : init_obs_sequence, init_obs, insert_obs_in_seq, &
set_obs_values, set_qc, obs_sequence_type, obs_type, &
copy_obs, set_copy_meta_data, set_qc_meta_data, set_obs_def
@@ -111,7 +112,7 @@
integer :: obs_unit
integer :: obs_prof, obs_kind, obs_kind_gen, which_vert, iqc, obstype, pc
real (r8) :: obs_err, lon, lat, lev, zob, time, pre_time, rcount, zob2
-real (r8) :: vloc, obs_value, aqc, var2, lon2c, lonc
+real (r8) :: vloc, obs_value, aqc, var2
real (r8) :: bin_beg, bin_end
@@ -208,11 +209,9 @@
cycle obsloop
endif
- lonc = lon
- if (lon2 > 360.0_r8 .and. lon < 180.0_r8) lonc = lon + 360.0_r8
-
! reject observations outside the bounding box
- if(lat < lat1 .or. lat > lat2 .or. lonc < lon1 .or. lonc > lon2) then
+ if(lat < lat1 .or. lat > lat2 .or. &
+ .not. is_longitude_between(lon, lon1, lon2)) then
iskip = iskip + 1
cycle obsloop
endif
@@ -398,30 +397,24 @@
var2, aqc, obs_kind, which_vert, seconds, days)
- if(obs_num == 1) then ! for the first observation
+ if(obs_num == 1) then ! the first observation, no prev_obs yet
- call insert_obs_in_seq(real_obs_sequence, obs)
- call copy_obs(prev_obs, obs)
- pre_time = time
+ call insert_obs_in_seq(real_obs_sequence, obs)
- else ! not the first observation
+ else if(time >= pre_time) then ! same time or later than previous obs
- if(time == pre_time) then ! same time as previous observation
+ call insert_obs_in_seq(real_obs_sequence, obs, prev_obs)
- call insert_obs_in_seq(real_obs_sequence, obs, prev_obs)
- call copy_obs(prev_obs, obs)
- pre_time = time
+ else ! earlier time, must start search from start of list
- else ! not the same time
+ call insert_obs_in_seq(real_obs_sequence, obs)
- call insert_obs_in_seq(real_obs_sequence, obs)
- call copy_obs(prev_obs, obs)
- pre_time = time
+ endif
- endif
+ ! save for next iteration
+ call copy_obs(prev_obs, obs)
+ pre_time = time
- endif
-
end do obsloop
200 continue
Modified: DART/trunk/obs_sequence/obs_sequence_mod.f90
===================================================================
--- DART/trunk/obs_sequence/obs_sequence_mod.f90 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/obs_sequence/obs_sequence_mod.f90 2009-04-13 16:21:33 UTC (rev 3809)
@@ -22,8 +22,8 @@
! FOR THESE TYPES THAT HAVE COPY SUBROUTINES.
use types_mod, only : r8, DEG2RAD, MISSING_R8
-use location_mod, only : location_type, interactive_location
- !%! location_inside
+use location_mod, only : location_type, interactive_location, &
+ is_location_in_region
use obs_def_mod, only : obs_def_type, get_obs_def_time, read_obs_def, &
write_obs_def, destroy_obs_def, copy_obs_def, &
interactive_obs_def, get_obs_def_location, &
@@ -34,8 +34,8 @@
use time_manager_mod, only : time_type, operator(>), operator(<), &
operator(>=), operator(/=), set_time, &
operator(-), operator(+), operator(==)
-use utilities_mod, only : get_unit, close_file, &
- register_module, error_handler, &
+use utilities_mod, only : get_unit, close_file, &
+ register_module, error_handler, &
find_namelist_in_file, check_namelist_read, &
E_ERR, E_WARN, E_MSG, nmlfileunit, do_output
@@ -1915,12 +1915,6 @@
type(location_type) :: location
logical :: is_this_last, inside, first_obs
-! FIXME:
-write(msg_string, *) 'Function not implemented yet'
-call error_handler(E_ERR, 'select_obs_by_location', msg_string, &
- source, revision, revdate)
-return
-! FIXME
! Initialize an observation type with appropriate size
call init_obs(obs, get_num_copies(seq), get_num_qc(seq))
@@ -1946,10 +1940,7 @@
location = get_obs_def_location(obs_def)
! each diff locations mod has a different one of these
- !%! FIXME:
- !%! inside = location_inside(location, min_loc, max_loc)
- inside = .false.
- !%! FIXME
+ inside = is_location_in_region(location, min_loc, max_loc)
! same code as delete/keep by obstype; do any code fixes both places
if (.not. inside) then
Modified: DART/trunk/obs_sequence/obs_sequence_tool.f90
===================================================================
--- DART/trunk/obs_sequence/obs_sequence_tool.f90 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/obs_sequence/obs_sequence_tool.f90 2009-04-13 16:21:33 UTC (rev 3809)
@@ -17,8 +17,8 @@
use utilities_mod, only : timestamp, register_module, initialize_utilities, &
find_namelist_in_file, check_namelist_read, &
error_handler, E_ERR, E_MSG, nmlfileunit
-use location_mod, only : location_type, get_location, &
- LocationName !! , vert_is_height !! set_location2
+use location_mod, only : location_type, get_location, set_location2, &
+ LocationName !! , vert_is_height
use obs_def_mod, only : obs_def_type, get_obs_def_time, get_obs_kind, &
get_obs_def_location
use obs_kind_mod, only : max_obs_kinds, get_obs_kind_name, get_obs_kind_index
@@ -37,9 +37,8 @@
get_num_key_range, delete_obs_by_typelist, &
get_obs_key, &
delete_obs_from_seq, get_next_obs_from_key, &
- delete_obs_by_qc, delete_obs_by_copy
-
- !%! select_obs_by_location ! not in repository yet
+ delete_obs_by_qc, delete_obs_by_copy, &
+ select_obs_by_location
implicit none
! <next few lines under version control, do not edit>
@@ -76,7 +75,7 @@
character(len = 32) :: obs_types(max_obs_input_types)
logical :: restrict_by_obs_type
integer :: num_obs_input_types
-logical :: restrict_by_location, restrict_by_latlon
+logical :: restrict_by_location
logical :: restrict_by_qc, restrict_by_copy, restrict_by_height
@@ -119,23 +118,15 @@
min_gps_height
!----------------------------------------------------------------
-! Start of the routine.
-! This routine basically opens the second observation sequence file
-! and 'inserts' itself into the first observation sequence file.
-! Each observation in the second file is independently inserted
-! or appended into the first file.
-! Inserting takes time, appending is much faster when possible.
+! Start of the program:
+!
+! Process each input observation sequence file in turn, optionally
+! selecting observations to insert into the output sequence file.
!----------------------------------------------------------------
call obs_seq_modules_used()
-! if you are not using a gregorian cal, set this to false
-! if anyone cares, we can add calendar integer to list
-if (gregorian_cal) then
- call set_calendar_type(GREGORIAN)
-endif
-
-! Initialize input obs_seq filenames and obs_types
+! Initialize input obs_seq filenames and obs_types before reading namelist.
do i = 1, max_num_input_files
filename_seq(i) = ""
enddo
@@ -156,7 +147,17 @@
! Record the namelist values used for the run ...
write(nmlfileunit, nml=obs_sequence_tool_nml)
+!if (do_nml_file()) write(nmlfileunit, nml=obs_sequence_tool_nml)
+!if (do_nml_term()) write( * , nml=obs_sequence_tool_nml)
+! if you are not using a gregorian cal, set this to false in the namelist.
+! if users need it, we could add a calendar type integer to the namelist,
+! if users want to specify a particular calendar which is not gregorian.
+! (earlier versions of this file had the test before the namelist read - duh.)
+if (gregorian_cal) then
+ call set_calendar_type(GREGORIAN)
+endif
+
! See if the user is restricting the obs types to be processed, and set up
! the values if so.
num_obs_input_types = 0
@@ -171,19 +172,28 @@
endif
! See if the user is restricting the obs locations to be processed, and set up
-! the values if so.
-!!if ((minval(min_box).ne.missing_r8) .or. (maxval(min_box).ne.missing_r8) .or. &
-!! (minval(max_box).ne.missing_r8) .or. (maxval(max_box).ne.missing_r8)) then
-!! restrict_by_location = .true.
-!! restrict_by_latlon = .false.
-!! min_loc = set_location2(min_box)
-!! max_loc = set_location2(max_box)
-!!else
+! the values if so. Note that if any of the values are not missing_r8, all of
+! them must have good (and consistent) values. i don't have code in here yet
+! to idiotproof this, but i should add it.
+if ((minval(min_box).ne.missing_r8) .or. (maxval(min_box).ne.missing_r8) .or. &
+ (minval(max_box).ne.missing_r8) .or. (maxval(max_box).ne.missing_r8)) then
+ restrict_by_location = .true.
+ min_loc = set_location2(min_box)
+ max_loc = set_location2(max_box)
+else
+ restrict_by_location = .false.
+endif
+
+! 3d sphere box check - locations module dependent, but an important one.
if ((min_lat /= -90.0_r8) .or. (max_lat /= 90.0_r8) .or. &
(min_lon /= 0.0_r8) .or. (max_lon /= 360.0_r8)) then
- ! 3d sphere box check - locations module dependent, but an important one.
+ ! do not allow namelist to set BOTH min/max box and by lat/lon.
+ if (restrict_by_location) then
+ call error_handler(E_ERR,'obs_sequence_tool', &
+ 'use either lat/lon box or min/max box but not both', &
+ source,revision,revdate)
+ endif
restrict_by_location = .true.
- restrict_by_latlon = .true.
if (trim(LocationName) /= 'loc3Dsphere') then
call error_handler(E_ERR,'obs_sequence_tool', &
'can only use lat/lon box with 3d sphere locations', &
@@ -223,18 +233,37 @@
source,revision,revdate)
endif
- ! handle wrap in lon; at the end of this block the min is [0,360) and
- ! max is [0, 720) and greater than min. modulo() must be used here and not
- ! mod() -- the result of modulo() is positive even if the input is negative.
- min_lon = modulo(min_lon,360.0_r8)
- max_lon = modulo(max_lon,360.0_r8)
- if (min_lon > max_lon) max_lon = max_lon + 360.0_r8
-
+ ! it is risky to allow == operations on floating point numbers.
+ ! if someone wants to select only obs exactly on a particular lon,
+ ! force them to do an interval [min-epsilon, max-epsilon].
+ ! move this test BEFORE the modulo so that if afterwards they are
+ ! the same (e.g. 0 and 360) that's ok -- it will mean all longitudes
+ ! will be accepted.
if (min_lon == max_lon) then
call error_handler(E_ERR,'obs_sequence_tool', &
- 'min_lon cannot equal max_lon', &
+ 'min_lon cannot exactly equal max_lon', &
source,revision,revdate)
endif
+
+ ! do not test for min < max, because lon can wrap around 0. ensure that
+ ! the values after here are [0, 360) and the compare routine knows about wrap.
+ ! must use modulo() and not just mod() to handle negative values correctly;
+ ! the result of modulo() is positive even if the input is negative.
+ min_lon = modulo(min_lon,360.0_r8)
+ max_lon = modulo(max_lon,360.0_r8)
+
+ ! create real location_type objects, knowing we are running with the
+ ! 3d sphere locations module.
+ min_box(1) = min_lon
+ min_box(2) = min_lat
+ min_box(3) = 0.0_r8
+ min_box(4) = -2.0_r8 !! FIXME: VERTISUNDEF, but not all loc mods have it
+ min_loc = set_location2(min_box)
+ max_box(1) = max_lon
+ max_box(2) = max_lat
+ max_box(3) = 0.0_r8
+ max_box(4) = -2.0_r8 !! FIXME: VERTISUNDEF, but not all loc mods have it
+ max_loc = set_location2(max_box)
else
restrict_by_location = .false.
endif
@@ -683,18 +712,13 @@
endif
endif
- ! optionally select obs in a lat/lon box
- ! the more generic is a box in any number of dimensions which matches
- ! the locations module which is determined at link time. but for now...
+ ! optionally select only obs inside a bounding box
+ ! most common expected use is a lat/lon box with the 3d sphere locations
+ ! mod, but this should work with any locations mod if the namelist uses
+ ! the right number of corners appropriate for the dimensions of the location
+ ! when setting min_box and max_box.
if (restrict_by_location) then
- if (restrict_by_latlon) then
- call select_obs_by_latlon(min_lon, max_lon, min_lat, max_lat, &
- seq, all_gone)
- else
- !%! call select_obs_by_location(min_loc, max_loc, seq, all_gone)
- msgstring = 'Select by general bounding box not implemented yet'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring)
- endif
+ call select_obs_by_location(min_loc, max_loc, seq, all_gone)
if(all_gone) then
if (print_msg) then
msgstring = 'Skipping: no obs in ' // trim(seqfilename) // &
@@ -1023,116 +1047,7 @@
end subroutine print_metadata
-!---------------------------------------------------------------------
-!---------------------------------------------------------------------
-! WARNING:
-! this routine should be in obs_sequence_mod.f90 but is has one line which
-! is dependent on the locations module being the 3d sphere; if it is
-! something else, this code will not work. i have a first cut at a loc-indep
-! routine but not all the kinks are worked out yet so put this code here for
-! now. to use merge with any other locations module, comment out
-! the calls to get_location() and do not select by bounding box in the nml.
-subroutine select_obs_by_latlon(min_lon, max_lon, min_lat, max_lat, &
- seq, all_gone)
-
-! Delete all observations in the sequence which are outside the bounding box.
-! If there are no obs left afterwards return that the sequence is all_gone.
-
-real(r8), intent(in) :: min_lon, max_lon, min_lat, max_lat
-type(obs_sequence_type), intent(inout) :: seq
-logical, intent(out) :: all_gone
-
-type(obs_def_type) :: obs_def
-type(obs_type) :: obs, prev_obs
-integer :: i, key
-type(location_type) :: location
-logical :: out_of_range, is_this_last, inside, first_obs
-real(r8) :: ll(3), lat, lon
-
-
-! Initialize an observation type with appropriate size
-call init_obs(obs, get_num_copies(seq), get_num_qc(seq))
-call init_obs(prev_obs, get_num_copies(seq), get_num_qc(seq))
-
-! Iterate entire sequence, deleting obs which are not in the box.
-! First, make sure there are obs to delete, and initialize first obs.
-if(.not. get_first_obs(seq, obs)) then
- all_gone = .true.
- call destroy_obs(obs)
- call destroy_obs(prev_obs)
- return
-endif
-
-first_obs = .true.
-prev_obs = obs
-
-! This is going to be O(n), n=num obs in seq
-is_this_last = .false.
-allobs : do while (.not. is_this_last)
-
- call get_obs_def(obs, obs_def)
- location = get_obs_def_location(obs_def)
-
- ! each diff locations mod has a different one of these
- ! FIXME: this next line is what makes this locations dependent.
- ! this also is ignoring the vertical for now.
- ll = get_location(location)
- lon = ll(1)
- lat = ll(2)
-
- ! if wrap in longitude, lon now (0,720]
- !if (lon < max_lon-360.0_r8) lon = lon + 360.0
- if (max_lon >= 360.0) lon = lon + 360.0
-
- ! box test.
- if ((lon < min_lon) .or. (lon > max_lon) .or. &
- (lat < min_lat) .or. (lat > max_lat)) then
- inside = .false.
- else
- inside = .true.
- endif
-
- ! same code as delete/keep by obstype; do any code fixes both places
- if (.not. inside) then
- if (first_obs) then
- call delete_obs_from_seq(seq, obs)
- if(.not. get_first_obs(seq, obs)) exit allobs
- else
-!print *, 'going to del obs key ', obs%key
-!print *, 'prev key is ', prev_obs%key
- call delete_obs_from_seq(seq, obs)
- ! cannot simply use prev_obs; cached copy out of sync with seq one
- key = get_obs_key(prev_obs)
- call get_next_obs_from_key(seq, key, obs, is_this_last)
-!print *, 'next obs now is key ', obs%key
- endif
- else
-!print *, 'no del, keep this obs key ', obs%key
- first_obs = .false.
- prev_obs = obs
-!print *, 'prev obs now is key ', prev_obs%key
-!print *, 'obs was key ', obs%key
- call get_next_obs(seq, prev_obs, obs, is_this_last)
-!print *, 'obs now is key ', obs%key
- endif
-
-end do allobs
-
-! Figure out if there are no more obs left in the sequence.
-if(.not. get_first_obs(seq, obs)) then
- all_gone = .true.
-else
- all_gone = .false.
-endif
-
-! Done. delete temp storage and return.
-call destroy_obs(obs)
-call destroy_obs(prev_obs)
-
-end subroutine select_obs_by_latlon
-
-
!---------------------------------------------------------------------
subroutine select_gps_by_height(min_height, seq, all_gone)
Modified: DART/trunk/observations/quikscat/convert_L2b.f90
===================================================================
--- DART/trunk/observations/quikscat/convert_L2b.f90 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/observations/quikscat/convert_L2b.f90 2009-04-13 16:21:33 UTC (rev 3809)
@@ -83,10 +83,6 @@
write(nmlfileunit, nml=convert_L2b_nml)
write( * , nml=convert_L2b_nml)
-lon1 = min(max(lon1,0.0_r8),360.0_r8)
-lon2 = min(max(lon2,0.0_r8),360.0_r8)
-if ( lon1 > lon2 ) lon2 = lon2 + 360.0_r8
-
call create_output_filename(l2b_file, output_name)
datafile = trim( datadir)//'/'//trim(l2b_file)
dartfile = trim(outputdir)//'/'//trim(output_name)
Modified: DART/trunk/observations/quikscat/quikscat_JPL_mod.f90
===================================================================
--- DART/trunk/observations/quikscat/quikscat_JPL_mod.f90 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/observations/quikscat/quikscat_JPL_mod.f90 2009-04-13 16:21:33 UTC (rev 3809)
@@ -25,7 +25,7 @@
use utilities_mod, only : get_unit, open_file, close_file, file_exist, &
register_module, error_handler, &
- E_ERR, E_MSG, timestamp
+ E_ERR, E_MSG, timestamp, is_longitude_between
use location_mod, only : location_type, set_location, VERTISHEIGHT
@@ -164,7 +164,7 @@
integer :: which_vert, uobstype, vobstype
real(r4) :: speed, dir
-real(r8) :: lonc, lon, lat, vloc
+real(r8) :: lon, lat, vloc
real(r8) :: u_obs, v_obs, u_var, v_var
real(r8) :: aqc
real(r8) :: sintheta, costheta, dirvar, speedvar
@@ -240,12 +240,9 @@
endif
! reject observations outside the bounding box (allowing wrapping)
- lonc = lon
- if (lon2 > 360.0_r8 .and. lon < 180.0_r8) lonc = lon + 360.0_r8
+ if(( lat < lat1) .or. ( lat > lat2 ) .or. &
+ (.not. is_longitude_between(lon, lon1, lon2))) cycle wvcloop
- if(( lat < lat1) .or. ( lat > lat2) .or. &
- (lonc < lon1) .or. (lonc > lon2)) cycle wvcloop
-
! QuikSCAT uses the oceanographic/flow convention ...
! 0.0 is TOWARD the north - in direct contradiction to
! atmospheric convention. Convert by adding 180 modulo 360
Modified: DART/trunk/utilities/utilities_mod.f90
===================================================================
--- DART/trunk/utilities/utilities_mod.f90 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/utilities/utilities_mod.f90 2009-04-13 16:21:33 UTC (rev 3809)
@@ -13,7 +13,7 @@
!-----------------------------------------------------------------------
!
-! A collection of simple useful programs.
+! A collection of simple useful routines:
!
! file_exist Function that returns if a given
! file name exists
@@ -68,6 +68,17 @@
! a single text variable ... to record in the
! netcdf output files.
!
+! is_longitude_between checks whether a given longitude is between
+! the two given limits, starting at the first and
+! going EAST until reaching the second. the end
+! points are included. if min=max, all points are
+! considered inside. there is no rejection of the
+! input values based on range; they are all converted
+! to [0-360) by calling modulo() before starting.
+! default is degrees, but there is an optional
+! argument to select radians instead.
+!
+!
! nsc start 31jan07
! idea - add some unit number routine here?
! you can extract the filename associated with a fortran unit number
@@ -79,7 +90,7 @@
!logical :: is_named
!integer :: rc
!
-!inquire(ncFileID, named=is_named, name=filename, iostat=rc)
+!inquire(unitnum, named=is_named, name=filename, iostat=rc)
!print *, 'is_named =', is_named, 'name = ', trim(filename)
!if ((rc /= 0) .or. (.not. is_named)) filename = 'unknown file'
!
@@ -87,7 +98,7 @@
!
!-----------------------------------------------------------------------
-use types_mod, only : r8
+use types_mod, only : r8, PI
use netcdf
implicit none
@@ -103,10 +114,12 @@
integer, parameter :: E_DBG = -1, E_MSG = 0, E_WARN = 1, E_ERR = 2
integer, parameter :: DEBUG = -1, MESSAGE = 0, WARNING = 1, FATAL = 2
+real(r8), parameter :: TWOPI = PI * 2.0_r8
+
public :: file_exist, get_unit, open_file, close_file, timestamp, &
register_module, error_handler, to_upper, &
nc_check, logfileunit, nmlfileunit, &
- find_textfile_dims, file_to_text, &
+ find_textfile_dims, file_to_text, is_longitude_between, &
initialize_utilities, finalize_utilities, dump_unit_attributes, &
find_namelist_in_file, check_namelist_read, &
set_tasknum, set_output, do_output, &
@@ -1256,6 +1269,48 @@
end subroutine file_to_text
+!#######################################################################
+
+function is_longitude_between (lon, minlon, maxlon, doradians)
+
+! uniform way to treat longitude ranges, in degrees, on a globe.
+! returns true if lon is between min and max, starting at min
+! and going EAST until reaching max. wraps across 0 longitude.
+! if min == max, all points are inside. includes edges.
+! if optional arg doradians is true, do computation in radians
+! between 0 and 2*PI instead of 360.
+
+real(r8), intent(in) :: lon, minlon, maxlon
+logical, intent(in), optional :: doradians
+logical :: is_longitude_between
+
+real(r8) :: minl, maxl, lon2, circumf
+
+circumf = 360.0_r8
+if (present(doradians)) then
+ if (doradians) circumf = TWOPI
+endif
+
+minl = modulo(minlon, circumf)
+maxl = modulo(maxlon, circumf)
+
+if (minl == maxl) then
+ is_longitude_between = .true. ! entire globe
+ return
+endif
+
+lon2 = modulo(lon, circumf)
+
+if (minl > maxl) then
+ maxl = maxl + circumf
+ if (lon2 < minl) lon2 = lon2 + circumf
+endif
+
+is_longitude_between = ((lon2 >= minl) .and. (lon2 <= maxl))
+
+
+end function is_longitude_between
+
!=======================================================================
! End of utilities_mod
!=======================================================================
Modified: DART/trunk/utilities/utilities_mod.html
===================================================================
--- DART/trunk/utilities/utilities_mod.html 2009-04-10 21:41:29 UTC (rev 3808)
+++ DART/trunk/utilities/utilities_mod.html 2009-04-13 16:21:33 UTC (rev 3809)
@@ -104,6 +104,9 @@
<TR><TD> </TD><TD><A HREF="#dump_unit_attributes">dump_unit_attributes</A></TD></TR>
<TR><TD> </TD><TD><A HREF="#find_namelist_in_file">find_namelist_in_file</A></TD></TR>
<TR><TD> </TD><TD><A HREF="#check_namelist_read">check_namelist_read</A></TD></TR>
+<TR><TD> </TD><TD><A HREF="#find_textfile_dims">find_textfile_dims</A></TD></TR>
+<TR><TD> </TD><TD><A HREF="#file_to_text">file_to_text</A></TD></TR>
+<TR><TD> </TD><TD><A HREF="#is_longitude_between">is_longitude_between</A></TD></TR>
<TR><TD> </TD><TD><A HREF="#set_tasknum">set_tasknum</A></TD></TR>
<TR><TD> </TD><TD><A HREF="#set_output">set_output</A></TD></TR>
<TR><TD> </TD><TD><A HREF="#do_output">do_output</A></TD></TR>
@@ -401,6 +404,95 @@
<BR>
<!--============= DESCRIPTION OF A SUBROUTINE =======================-->
+ <A NAME="find_textfile_dims"></A>
+ <P></P><HR><P></P>
+ <div class=routine>
+ <em class=call> call find_textfile_dims (fname, nlines, linelen)</em>
+ <pre>
+ character(len=*), intent (IN) :: <em class=code>fname</em>
+ integer, intent (OUT) :: <em class=code>nlines</em>
+ integer, intent (OUT) :: <em class=code>linelen</em>
+ </pre></div>
+ <H3 class=indent1>Description</H3>
+ <P>
+ Determines the number of lines and maximum line length of an
+ ASCII text file.
+ </P>
+ <TABLE width=100% border=0 summary="" celpadding=3>
+ <TR><TD valign=top><em class=code>fname</em></TD>
+ <TD>input, character string file name</TD></TR>
+ <TR><TD valign=top><em class=code>nlines</em></TD>
+ <TD>output, number of lines in the file</TD></TR>
+ <TR><TD valign=top><em class=code>linelen</em></TD>
+ <TD>output, length of longest line in the file</TD></TR>
+ </TABLE>
+ <BR>
+
+<!--============= DESCRIPTION OF A SUBROUTINE =======================-->
+ <A NAME="file_to_text"></A>
+ <P></P><HR><P></P>
+ <div class=routine>
+ <em class=call> call file_to_text (fname, textblock)</em>
+ <pre>
+ character(len=*), intent (IN) :: <em class=code>fname</em>
+ character(len=*), dimension(:), intent (OUT) :: <em class=code>textblock</em>
+ </pre></div>
+ <H3 class=indent1>Description</H3>
+ <P>
+ Opens the given filename and reads ASCII text lines into a
+ character array.
+ </P>
+ <TABLE width=100% border=0 summary="" celpadding=3>
+ <TR><TD valign=top><em class=code>fname</em></TD>
+ <TD>input, character string file name</TD></TR>
+ <TR><TD valign=top><em class=code>textblock</em></TD>
+ <TD>output, character array of text in the file</TD></TR>
+ </TABLE>
+ <BR>
+
+<!--============= DESCRIPTION OF A FUNCTION ========================-->
+ <A NAME="is_longitude_between"></A>
+ <P></P><HR><P></P>
+ <div class=routine>
+ <em class=call> var = is_longitude_between(lon, minlon, maxlon
+ <em class=optionalcode>[, doradians]</em>) </em>
+ <pre>
+ real(r8), intent(in) :: <em class=code>lon</em>
+ real(r8), intent(in) :: <em class=code>minlon</em>
+ real(r8), intent(in) :: <em class=code>maxlon</em>
+ logical, intent(in), optional :: <em class=optionalcode>doradians</em>
+ logical :: <em class=code>is_longitude_between</em>
+ </pre></div>
+ <H3 class=indent1>Description</H3>
+ <P>
+ Uniform way to test longitude ranges, in degrees, on a globe.
+ Returns true if lon is between min and max, starting at min
+ and going EAST until reaching max. Wraps across 0 longitude.
+ If min equals max, all points are inside. Includes endpoints.
+ If optional arg doradians is true, do computation in radians
+ between 0 and 2*PI instead of default 360. There is no rejection
+ of input values based on range; they are all converted to a known
+ range by calling modulo() first.
+ </P>
+ <TABLE width=100% border=0 summary="" celpadding=3>
+ <TR><TD valign=top><em class=code>var </em></TD>
+ <TD>True if lon is between min and max.</TD></TR>
+ <TR><TD valign=top><em class=code>lon </em></TD>
+ <TD>Location to test.</TD></TR>
+ <TR><TD valign=top><em class=code>minlon </em></TD>
+ <TD>Minimum longitude. Region will start here and go east.</TD></TR>
+ <TR><TD valign=top><em class=code>maxlon </em></TD>
+ <TD>Maximum longitude. Region will end here.</TD></TR>
+ <TR><TD valign=top><em class=optionalcode>doradians </em></TD>
+ <TD>Optional argument. Default computations are in degrees. If
+ this argument is specified and is .true., do the computation
+ in radians, and wrap across the globe at 2 * PI. All inputs
+ must then be specified in radians.</TD></TR>
+ </TABLE>
+ <BR>
+
+
+<!--============= DESCRIPTION OF A SUBROUTINE =======================-->
<A NAME="to_upper"></A>
<P></P><HR><P></P>
<div class=routine>
More information about the Dart-dev
mailing list