[Dart-dev] [6332] DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.f90 : minor updates:
nancy at ucar.edu
nancy at ucar.edu
Tue Jul 30 14:09:00 MDT 2013
Revision: 6332
Author: nancy
Date: 2013-07-30 14:09:00 -0600 (Tue, 30 Jul 2013)
Log Message:
-----------
minor updates:
1) when superob'ing aircraft and satwind obs ignore those which
are within the superob radius of either pole. (the code which
computes the average lat/lon doesn't handle wrapping over the poles.)
print a count if any obs are ignored so user is aware.
2) fix handling of obs close to the prime meridian when superob'ing
aircraft and satwind obs. now even at high latitudes the average
longitude will be computed consistently.
3) use incoming obs count or computed obs count for size of temporary
obs sequences.
Modified Paths:
--------------
DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.f90
-------------- next part --------------
Modified: DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.f90
===================================================================
--- DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.f90 2013-07-30 17:12:57 UTC (rev 6331)
+++ DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.f90 2013-07-30 20:09:00 UTC (rev 6332)
@@ -170,14 +170,14 @@
end if
! create obs sequences for different obs types
-call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_rawin)
-call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_sfc)
-call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_acars)
-call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_satwnd)
-call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_prof)
-call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_gpsro)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_rawin)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_sfc)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_acars)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_satwnd)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_prof)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_gpsro)
call create_new_obs_seq(num_copies, num_qc, 100, seq_tc)
-call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_other)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_other)
! read input obs_seq file, divide into platforms
call read_and_parse_input_seq(file_name_input, real_nx, real_ny, obs_boundary, &
@@ -1359,14 +1359,14 @@
real(r8), intent(in) :: hdist, vdist
integer :: num_copies, num_qc, nloc, k, locdex, obs_kind, n, &
- num_obs
-logical :: last_obs
+ num_obs, poleward_obs
+logical :: last_obs, close_to_greenwich, pole_check
real(r8) :: nuwnd, latu, lonu, preu, uwnd, erru, qcu, nvwnd, latv, &
lonv, prev, vwnd, errv, qcv, ntmpk, latt, lont, pret, &
tmpk, errt, qct, nqvap, latq, lonq, preq, qvap, errq, &
dwpt, errd, qcd, ndwpt, latd, lond, pred, relh, errr, &
- qcr, nrelh, latr, lonr, prer, qcq, obs_dist, &
- xyz_loc(3), obs_val(1), qc_val(1)
+ qcr, nrelh, latr, lonr, prer, qcq, obs_dist, &
+ xyz_loc(3), obs_val(1), qc_val(1), lon_degree_limit
type(location_type) :: obs_loc
type(obs_def_type) :: obs_def
type(obs_type) :: obs, prev_obs
@@ -1394,7 +1394,7 @@
call init_obs(obs, num_copies, num_qc)
call init_obs(prev_obs, num_copies, num_qc)
-last_obs = .false. ; nloc = 0
+last_obs = .false. ; nloc = 0 ; poleward_obs = 0
if ( .not. get_first_obs(seq, obs) ) last_obs = .true.
! loop over all observations in sequence, add to ACARS observation type
@@ -1420,10 +1420,19 @@
if ( locdex < 1 ) then ! create new observation location type
+ ! test if we are within hdist of either pole, and punt for now on those
+ ! obs because we can't accurately average points that wrap the poles.
+ ! (count up obs here and print later)
+ if (pole_check(xyz_loc(1), xyz_loc(2), hdist)) then
+ poleward_obs = poleward_obs + 1
+ goto 200
+ endif
+
nloc = nloc + 1
locdex = nloc
+
+
airobs(locdex)%lon = xyz_loc(1)
- if ( airobs(locdex)%lon < 5.0_r8 ) airobs(locdex)%lon = airobs(locdex)%lon + 360.0_r8
airobs(locdex)%lat = xyz_loc(2)
airobs(locdex)%pressure = xyz_loc(3)
airobs(locdex)%obs_loc = obs_loc
@@ -1477,11 +1486,18 @@
end if
+200 continue ! come here to skip this obs
+
prev_obs = obs
call get_next_obs(seq, prev_obs, obs, last_obs)
end do
+if (poleward_obs > 0) then
+ write(6, *) 'WARNING: skipped ', poleward_obs, ' of ', poleward_obs+nloc, ' aircraft obs because'
+ write(6, *) 'they were within ', hdist, ' KM of the poles (the superobs distance).'
+endif
+
call destroy_obs_sequence(seq)
call create_new_obs_seq(num_copies, num_qc, num_obs, seq)
call init_obs(obs, num_copies, num_qc)
@@ -1504,6 +1520,11 @@
if ( airobs(k)%lat /= missing_r8 ) then ! create initial superob
+ call superob_location_check(airobs(k)%lon, airobs(k)%lat, hdist, &
+ close_to_greenwich, lon_degree_limit)
+ if (close_to_greenwich) call wrap_lon(airobs(k)%lon, airobs(k)%lon - lon_degree_limit, &
+ airobs(k)%lon + lon_degree_limit)
+
if ( airobs(k)%uwnd /= missing_r8 ) then
nuwnd = nuwnd + 1.0_r8
latu = latu + airobs(k)%lat
@@ -1572,6 +1593,9 @@
obs_dist = get_dist(airobs(k)%obs_loc, airobs(n)%obs_loc, 2, 2, .true.) * earth_radius
if ( obs_dist <= hdist .and. abs(airobs(k)%pressure-airobs(n)%pressure) <= vdist ) then
+ if (close_to_greenwich) call wrap_lon(airobs(n)%lon, airobs(k)%lon - lon_degree_limit, &
+ airobs(k)%lon + lon_degree_limit)
+
if ( airobs(n)%uwnd /= missing_r8 ) then
nuwnd = nuwnd + 1.0_r8
latu = latu + airobs(n)%lat
@@ -1655,7 +1679,7 @@
end if
- if ( nvwnd > 0.0_r8 ) then ! write mieridional wind superob
+ if ( nvwnd > 0.0_r8 ) then ! write meridional wind superob
latv = latv / nvwnd
lonv = lonv / nvwnd
@@ -1771,10 +1795,11 @@
real(r8), intent(in) :: hdist, vdist
integer :: num_copies, num_qc, nloc, k, locdex, obs_kind, n, &
- num_obs
-logical :: last_obs
+ num_obs, poleward_obs
+logical :: last_obs, close_to_greenwich, pole_check
real(r8) :: nwnd, lat, lon, pres, uwnd, erru, qcu, vwnd, &
- errv, qcv, obs_dist, xyz_loc(3), obs_val(1), qc_val(1)
+ errv, qcv, obs_dist, xyz_loc(3), obs_val(1), qc_val(1), &
+ lon_degree_limit
type(location_type) :: obs_loc
type(obs_def_type) :: obs_def
@@ -1801,7 +1826,7 @@
call init_obs(obs, num_copies, num_qc)
call init_obs(prev_obs, num_copies, num_qc)
-last_obs = .false. ; nloc = 0
+last_obs = .false. ; nloc = 0 ; poleward_obs = 0
if ( .not. get_first_obs(seq, obs) ) last_obs = .true.
! loop over satellite winds, create list
@@ -1828,10 +1853,19 @@
if ( locdex < 1 ) then ! create new observation type
+ ! test if we are within hdist of either pole, and punt for now on those
+ ! obs because we can't accurately average points that wrap the poles.
+ ! (hdist is radius, in KM, of region of interest.)
+ if (pole_check(xyz_loc(1), xyz_loc(2), hdist)) then
+ ! count up obs here and print later
+ poleward_obs = poleward_obs + 1
+ goto 200
+ endif
+
nloc = nloc + 1
locdex = nloc
+
satobs(locdex)%lon = xyz_loc(1)
- if ( satobs(locdex)%lon < 5.0_r8 ) satobs(locdex)%lon = satobs(locdex)%lon + 360.0_r8
satobs(locdex)%lat = xyz_loc(2)
satobs(locdex)%pressure = xyz_loc(3)
satobs(locdex)%obs_loc = obs_loc
@@ -1856,11 +1890,18 @@
end if
+200 continue ! come here to skip this obs
+
prev_obs = obs
call get_next_obs(seq, prev_obs, obs, last_obs)
end do
+if (poleward_obs > 0) then
+ write(6, *) 'WARNING: skipped ', poleward_obs, ' of ', poleward_obs+nloc, ' satwind obs because'
+ write(6, *) 'they were within ', hdist, ' KM of the poles (the superobs distance).'
+endif
+
! create new sequence
call destroy_obs_sequence(seq)
call create_new_obs_seq(num_copies, num_qc, num_obs, seq)
@@ -1870,6 +1911,11 @@
if ( satobs(k)%uwnd /= missing_r8 .and. satobs(k)%vwnd /= missing_r8 ) then
+ call superob_location_check(satobs(k)%lon, satobs(k)%lat, hdist, &
+ close_to_greenwich, lon_degree_limit)
+ if (close_to_greenwich) call wrap_lon(satobs(k)%lon, satobs(k)%lon - lon_degree_limit, &
+ satobs(k)%lon + lon_degree_limit)
+
nwnd = 1.0_r8
lat = satobs(k)%lat
lon = satobs(k)%lon
@@ -1889,6 +1935,9 @@
obs_dist = get_dist(satobs(k)%obs_loc, satobs(n)%obs_loc, 2, 2, .true.) * earth_radius
if ( obs_dist <= hdist .and. abs(satobs(k)%pressure-satobs(n)%pressure) <= vdist ) then
+ if (close_to_greenwich) call wrap_lon(satobs(n)%lon, satobs(k)%lon - lon_degree_limit, &
+ satobs(k)%lon + lon_degree_limit)
+
nwnd = nwnd + 1.0_r8
lat = lat + satobs(n)%lat
lon = lon + satobs(n)%lon
@@ -1974,6 +2023,175 @@
return
end function surface_obs_check
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! pole_check - determine if we are within km_dist of either pole.
+! function returns true if so, false if not.
+!
+! lon - longitude in degrees
+! lat - latitude in degrees
+! km_dist - horizontal superob radius in kilometers
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+function pole_check(lon, lat, km_dist)
+
+use types_mod, only : r8, earth_radius
+use location_mod, only : location_type, get_dist, set_location, VERTISUNDEF
+
+implicit none
+
+real(r8), intent(in) :: lon, lat, km_dist
+logical :: pole_check
+
+type(location_type) :: thisloc, pole
+
+
+! create a point at this lon/lat, and at the nearest pole
+thisloc = set_location(lon, lat, 0.0_r8, VERTISUNDEF)
+if (lat >= 0) then
+ pole = set_location(0.0_r8, 90.0_r8, 0.0_r8, VERTISUNDEF)
+else
+ pole = set_location(0.0_r8, -90.0_r8, 0.0_r8, VERTISUNDEF)
+endif
+
+! are we within km_distance of that pole?
+if ( get_dist(thisloc, pole, 1, 1, .true.) * earth_radius <= km_dist ) then
+ pole_check = .true.
+else
+ pole_check = .false.
+endif
+
+return
+end function pole_check
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! superob_location_check - determine if a point is close to the greenwich
+! longitude based on the given latitude and distance.
+! if this location is close enough to longitude=0
+! we have to take action in order not to average points
+! at 0 and 360 and end up with 180 by mistake. as long
+! as we treat all points consistently when doing an average
+! the exact determination of 'close enough' isn't critical.
+! (the minimum and maximum extents in longitude for a given
+! point and radius is left as an exercise for the reader.
+! hint, they are not along the same latitude circle.)
+!
+! lon - longitude in degrees (input)
+! lat - latitude in degrees (input)
+! km_dist - horizontal superob radius in kilometers (input)
+! near_greenwich - returns true if the given lon/lat is potentially within
+! km_dist of longitude 0 (output)
+! lon_degree_limit - number of degrees along a latitude circle that the
+! km_dist equates to, plus a tolerance (output)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine superob_location_check(lon, lat, km_dist, near_greenwich, lon_degree_limit)
+
+use types_mod, only : r8, PI, earth_radius, RAD2DEG, DEG2RAD
+use location_mod, only : location_type, get_dist, set_location, VERTISUNDEF
+
+implicit none
+
+real(r8), intent(in) :: lon, lat, km_dist
+logical, intent(out) :: near_greenwich
+real(r8), intent(out) :: lon_degree_limit
+
+real(r8) :: lat_radius
+real(r8), parameter :: fudge_factor = 1.2_r8 ! add a flat 20%
+
+
+! the problem with trying to superob in a circle that's specified as a
+! radius of kilometers is that it isn't parallel with longitude lines as
+! you approach the poles. also when averaging lon values near the prime
+! meridian some values are < 360 and some are > 0 but are close in space.
+! simple suggestion for all but the highest latitudes:
+! if dist between lat/lon and lat/0 is < hdist, add 360 to any values >= 0.
+! if the final averaged point >= 360 subtract 360 before setting the location.
+! this still isn't good enough as you get closer to the poles; there
+! lat/lon averages are a huge pain. hdist could be shorter across the
+! pole and therefore a lon that is 180 degrees away could still be
+! 'close enough' to average in the same superob box. probably the
+! best suggestion for this case is to convert lat/lon into 3d cartesian
+! coords, average the x/y/z's separately, and then convert back.
+! (no code like this has been implemented here.)
+!
+! obs_dist_in_km = earth_radius_in_km * dist_in_radians
+! (which is what get_dist(loc1, loc2, 0, 0, .true.) returns)
+!
+! dist_in_radians = obs_dist_in_km / earth_radius_in_km
+!
+
+! figure out how far in degrees km_dist is at this latitude
+! if we traveled along the latitude line. and since on a sphere
+! the actual closest point to a longitude line isn't found by following
+! a latitude circle, add a percentage. as long as all points are
+! treated consistently it doesn't matter if the degree value is exact.
+
+lat_radius = earth_radius * cos(lat*DEG2RAD)
+lon_degree_limit = ((km_dist / lat_radius) * RAD2DEG) * fudge_factor
+
+! are we within 'lon_degree_limit' of the greenwich line?
+if (lon <= lon_degree_limit .or. (360.0_r8 - lon) <= lon_degree_limit) then
+ near_greenwich = .true.
+else
+ near_greenwich = .false.
+endif
+
+return
+end subroutine superob_location_check
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! wrap_lon - update the incoming longitude possibly + 360 degrees if
+! the given limits define a region that crosses long=0.
+! all values should be in units of degrees. 'lon' value
+! should be between westlon and eastlon.
+!
+! lon - longitude to update, returns either unchanged or + 360
+! westlon - westernmost longitude of region in degrees
+! eastlon - easternmost longitude of region in degrees
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine wrap_lon(lon, westlon, eastlon)
+
+use types_mod, only : r8
+
+! uniform way to treat longitude ranges, in degrees, on a globe.
+! adds 360 to the incoming lon if region crosses longitude 0 and
+! given point is east of lon=0.
+
+real(r8), intent(inout) :: lon
+real(r8), intent(in) :: westlon, eastlon
+
+real(r8) :: westl, eastl
+real(r8), parameter :: circumf = 360.0_r8
+
+! ensure the region boundaries and target point are between 0 and 360.
+! the modulo() function handles negative values ok; mod() does not.
+westl = modulo(westlon, circumf)
+eastl = modulo(eastlon, circumf)
+lon = modulo(lon, circumf)
+
+! if the 'region' is the entire globe you can return now.
+if (westl == eastl) return
+
+! here's where the magic happens:
+! normally the western boundary longitude (westl) has a smaller magnitude than
+! the eastern one (eastl). but westl will be larger than eastl if the region
+! of interest crosses the prime meridian. e.g. westl=100, eastl=120 doesn't
+! cross it, westl=340, eastl=10 does. for regions crossing lon=0, a target lon
+! west of lon=0 should not be changed; a target lon east of lon=0 needs +360 degrees.
+! e.g. lon=350 stays unchanged; lon=5 becomes lon=365.
+
+if (westl > eastl .and. lon <= eastl) lon = lon + circumf
+
+return
+end subroutine wrap_lon
+
+
+
! <next few lines under version control, do not edit>
! $URL$
! $Id$
More information about the Dart-dev
mailing list