[Dart-dev] [7007] DART/trunk: Now handles mandatory pressure levels as well as ' surface' or 'undefined'.
nancy at ucar.edu
nancy at ucar.edu
Tue Jun 3 10:12:40 MDT 2014
Revision: 7007
Author: thoar
Date: 2014-06-03 10:12:40 -0600 (Tue, 03 Jun 2014)
Log Message:
-----------
Now handles mandatory pressure levels as well as 'surface' or 'undefined'.
Modified Paths:
--------------
DART/trunk/models/mpas_atm/work/input.nml
DART/trunk/models/mpas_atm/work/path_names_obs_seq_coverage
DART/trunk/obs_sequence/obs_seq_coverage.f90
-------------- next part --------------
Modified: DART/trunk/models/mpas_atm/work/input.nml
===================================================================
--- DART/trunk/models/mpas_atm/work/input.nml 2014-06-02 15:41:43 UTC (rev 7006)
+++ DART/trunk/models/mpas_atm/work/input.nml 2014-06-03 16:12:40 UTC (rev 7007)
@@ -297,7 +297,7 @@
&obs_seq_to_netcdf_nml
- obs_sequence_name = 'obs_seq.raobT1.final'
+ obs_sequence_name = 'obs_seq.final'
obs_sequence_list = ''
append_to_netcdf = .false.
lonlim1 = 0.0
@@ -390,7 +390,7 @@
&obs_seq_coverage_nml
obs_sequence_list = 'obs_coverage_list.txt'
obs_sequence_name = ''
- obs_of_interest = 'RADIOSONDE_TEMPERATURE'
+ obs_of_interest = 'RADIOSONDE_TEMPERATURE', 'RADIOSONDE_U_WIND_COMPONENT', 'RADIOSONDE_V_WIND_COMPONENT'
textfile_out = 'obsdef_mask.txt'
netcdf_out = 'obsdef_mask.nc'
calendar = 'Gregorian'
@@ -405,7 +405,7 @@
latlim1 = -90.0
latlim2 = 90.0
verbose = .true.
- debug = .true.
+ debug = .false.
/
# selections_file is a list of obs_defs output
Modified: DART/trunk/models/mpas_atm/work/path_names_obs_seq_coverage
===================================================================
--- DART/trunk/models/mpas_atm/work/path_names_obs_seq_coverage 2014-06-02 15:41:43 UTC (rev 7006)
+++ DART/trunk/models/mpas_atm/work/path_names_obs_seq_coverage 2014-06-03 16:12:40 UTC (rev 7007)
@@ -12,6 +12,5 @@
obs_sequence/obs_seq_coverage.f90
obs_sequence/obs_sequence_mod.f90
random_seq/random_seq_mod.f90
-sort/sort_mod.f90
time_manager/time_manager_mod.f90
utilities/utilities_mod.f90
Modified: DART/trunk/obs_sequence/obs_seq_coverage.f90
===================================================================
--- DART/trunk/obs_sequence/obs_seq_coverage.f90 2014-06-02 15:41:43 UTC (rev 7006)
+++ DART/trunk/obs_sequence/obs_seq_coverage.f90 2014-06-03 16:12:40 UTC (rev 7007)
@@ -13,52 +13,44 @@
!
! The observation sequence file only contains lat/lon/level/which_vert,
! so this is all we have to work with.
-!
-! From: Soyoung Ha
-! Date: November 5, 2010 12:00:22 PM MDT
-! To: Tim Hoar, Nancy Collins
-! Subject: obs_seq.out on bluefire
-!
-! Hi Tim and Nancy,
-!
-! I just found that my current DART codes were compiled with r8, not r4.
-! But the MADIS obs netcdf files have lat/lon/elevation information in "float"
-! and only time information in "double".
-! So, although I see the "double" accuracy in my obs_seq.out,
-! I guess those were artificially generated by DART.
-! My obs_seq.out files which were preprocessed using wrf_dart_obs_preprocess
-! every 3 hrs are all available in be:/ptmp/syha/SFCDA/OBS_SEQ/Jun08/Preproc_3hrly/
-! for a month-long period of Jun 2008.
-! Just for your information, I put an ascii header file named "20080601_0000.head.txt"
-! in the same place, which is what I've got from doing "ncdump -h" for the original
-! MADIS metar netcdf file valid at the initial cycle.
-
!-----------------------------------------------------------------------
-use types_mod, only : r4, r8, digits12, MISSING_R8, MISSING_R4, &
+use types_mod, only : r4, r8, digits12, MISSING_R8, MISSING_R4, MISSING_I, &
metadatalength, obstypelength
+
use obs_sequence_mod, only : read_obs_seq, obs_type, obs_sequence_type, get_first_obs, &
get_obs_def, get_copy_meta_data, &
get_next_obs, init_obs, init_obs_sequence, &
assignment(=), get_num_copies, get_num_qc, get_qc, &
static_init_obs_sequence, destroy_obs_sequence, destroy_obs, &
read_obs_seq_header, get_qc_meta_data
+
use obs_def_mod, only : obs_def_type, get_obs_def_time, get_obs_kind, write_obs_def, &
get_obs_def_location, set_obs_def_time, &
set_obs_def_location, set_obs_def_kind, set_obs_def_error_variance
+
use obs_kind_mod, only : max_obs_kinds, get_obs_kind_name, get_obs_kind_index, &
write_obs_kind
+
use location_mod, only : location_type, get_location, set_location_missing, &
write_location, operator(/=), operator(==), &
set_location, is_location_in_region, query_location, &
nc_write_location_atts, nc_get_location_varids, &
- nc_write_location, LocationDims
+ nc_write_location, LocationDims, &
+ vert_is_undef, VERTISUNDEF, &
+ vert_is_surface, VERTISSURFACE, &
+ vert_is_level, VERTISLEVEL, &
+ vert_is_pressure, VERTISPRESSURE, &
+ vert_is_height, VERTISHEIGHT, &
+ vert_is_scale_height, VERTISSCALEHEIGHT
+
use time_manager_mod, only : time_type, set_date, set_time, get_time, &
set_calendar_type, get_calendar_string, &
print_time, print_date, &
operator(+), operator(-), operator(<), operator(>), &
operator(==), operator(/=), operator(<=), operator(/), &
operator(>=), operator(*)
+
use utilities_mod, only : get_unit, close_file, register_module, &
file_exist, error_handler, E_ERR, E_WARN, E_MSG, &
initialize_utilities, nmlfileunit, finalize_utilities, &
@@ -80,6 +72,7 @@
!---------------------------------------------------------------------
!---------------------------------------------------------------------
+! FIXME ... "voxel" is a much better term than "station".
type station
integer :: obs_type
type(location_type) :: location
@@ -98,7 +91,9 @@
integer :: num_out_stat ! total number of desired stations found
integer :: num_out_total ! total number of desired locations * times found
-integer, parameter :: MAX_OBS_INPUT_TYPES = 500 ! lazy, just going big
+integer, parameter :: MAX_OBS_NAMES_IN_NAMELIST = 500 ! lazy, just going big
+integer, parameter :: NUM_MANDATORY_LEVELS = 14 ! number of mandatory pressure levels
+real(r8), parameter :: OnePa = 1.0_r8 ! 1 Pa ... tolerance for vertical comparisons
!---------------------------------------------------------------------
! variables associated with the observation
@@ -125,7 +120,7 @@
character(len = 129) :: obs_sequence_list = 'obs_coverage_list.txt'
character(len = 129) :: obs_sequence_name = ''
-character(len = obstypelength) :: obs_of_interest(MAX_OBS_INPUT_TYPES) = ''
+character(len = obstypelength) :: obs_of_interest(MAX_OBS_NAMES_IN_NAMELIST) = ''
character(len = 129) :: textfile_out = 'obsdef_mask.txt'
character(len = 129) :: netcdf_out = 'obsdef_mask.nc'
character(len = 129) :: calendar = 'Gregorian'
@@ -143,6 +138,7 @@
logical :: verbose = .false.
logical :: debug = .false. ! undocumented ... on purpose
+
namelist /obs_seq_coverage_nml/ obs_sequence_list, obs_sequence_name, &
obs_of_interest, textfile_out, netcdf_out, calendar, &
first_analysis, last_analysis, forecast_length_days, &
@@ -151,6 +147,21 @@
verbose, debug
!-----------------------------------------------------------------------
+! General purpose variables
+!-----------------------------------------------------------------------
+
+integer :: ifile, iobs, ngood
+integer :: i, j, io, ncunit
+
+type(time_type) :: obs_time, no_time, last_possible_time
+
+character(len = 129) :: ncName, string1, string2, string3
+
+! ~# of degrees for 1/2 meter at Earth equator
+! 360 deg-earth/(40000 km-earth * 1000m-km)
+real(r8), parameter :: HALF_METER = 180.0_r8 / (40000.0_r8 * 1000.0_r8)
+
+!-----------------------------------------------------------------------
! Quantities of interest
!-----------------------------------------------------------------------
@@ -162,7 +173,7 @@
character(len=metadatalength), allocatable, dimension(:) :: module_qc_copy_names
character(len=metadatalength), allocatable, dimension(:) :: qc_copy_names
-integer, dimension(max_obs_kinds) :: obs_kinds_inds = 0
+integer, dimension(max_obs_kinds) :: obs_type_inds = 0
real(r8), allocatable, dimension(:) :: qc_values
type(time_type), allocatable, dimension(:) :: all_verif_times
type(time_type), allocatable, dimension(:,:) :: verification_times
@@ -174,21 +185,8 @@
integer :: num_verification_times ! number of verification times - total
integer :: nT_minimum ! will settle for this many verif times - total
-!-----------------------------------------------------------------------
-! General purpose variables
-!-----------------------------------------------------------------------
+real(r8), dimension(NUM_MANDATORY_LEVELS) :: mandatory_levels = MISSING_R8 ! pressure level (hPa)
-integer :: ifile, nread, ngood
-integer :: i, j, io, ncunit
-
-type(time_type) :: obs_time, no_time, last_possible_time
-
-character(len = 129) :: ncName, string1, string2, string3
-
-! ~# of degrees for 1/2 meter at Earth equator
-! 360 deg-earth/(40000 km-earth * 1000m-km)
-real(r8), parameter :: HALF_METER = 180.0_r8 / (40000.0_r8 * 1000.0_r8)
-
!=======================================================================
! Get the party started
!=======================================================================
@@ -218,7 +216,7 @@
if (temporal_coverage_percent < 100.0_r8) then
write(string1,*)'namelist: temporal_coverage_percent (',temporal_coverage_percent,&
') must be == 100.0 for now.)'
-! call error_handler(E_ERR, 'obs_seq_coverage', string1, source, revision, revdate)
+ call error_handler(E_MSG, 'obs_seq_coverage', string1, source, revision, revdate)
endif
if ((obs_sequence_name /= '') .and. (obs_sequence_list /= '')) then
@@ -230,16 +228,17 @@
call set_calendar_type(calendar)
call get_calendar_string(calendar)
+call SetPressureLevels()
-minl = set_location( (/ lonlim1, latlim1, 0.0_r8, 1.0_r8 /)) ! vertical unimportant
-maxl = set_location( (/ lonlim2, latlim2, 0.0_r8, 1.0_r8 /)) ! vertical unimportant
+minl = set_location(lonlim1, latlim1, 0.0_r8, VERTISUNDEF)
+maxl = set_location(lonlim2, latlim2, 0.0_r8, VERTISUNDEF)
! Determine if the desired observation types exist
-TypeLoop : do i = 1,MAX_OBS_INPUT_TYPES
+TypeLoop : do i = 1,MAX_OBS_NAMES_IN_NAMELIST
if ( (len(obs_of_interest(i)) == 0) .or. &
- (obs_of_interest(i) == "") ) exit TypeLoop
+ (obs_of_interest(i) == "") ) exit TypeLoop
string2 = adjustl(obs_of_interest(i))
@@ -252,10 +251,18 @@
if (verbose) write(*,*)trim(string2),' is type ',flavor_of_interest
- obs_kinds_inds(flavor_of_interest) = 1 ! indicate that we want this one
+ obs_type_inds(flavor_of_interest) = 1 ! indicate that we want this one
enddo TypeLoop
+if (all(obs_type_inds < 1)) then
+ write(string1,*)' No desired observation types of interest specified.'
+ write(string2,*)' Check the value of "obs_of_interest".'
+ call error_handler(E_ERR, 'obs_seq_coverage', string1, &
+ source, revision, revdate, text2=string2)
+endif
+
+
! Set the verification time array (global storage)
call set_required_times(first_analysis, last_analysis, &
@@ -371,14 +378,6 @@
write(*,*)
endif
- if (num_obs <= 1) then
- last_ob_flag = .TRUE.
- else
- last_ob_flag = .FALSE.
- endif
-
- if ( debug ) write(*,*)'num obs/last_ob_flag are ',num_obs, last_ob_flag
-
!--------------------------------------------------------------------
! Read the entire observation sequence - allocates 'seq' internally
!--------------------------------------------------------------------
@@ -444,24 +443,25 @@
endif
- !--------------------------------------------------------------------
- ! Read the first observation in the sequence
- !--------------------------------------------------------------------
-
ngood = 0
- if ( .not. get_first_obs(seq, obs1) ) &
- call error_handler(E_ERR,'obs_seq_coverage', &
- 'No first observation in sequence.', &
- source,revision,revdate)
-
!--------------------------------------------------------------------
- ObservationLoop : do nread = 1,num_obs
+ ObservationLoop : do iobs = 1,num_obs
!--------------------------------------------------------------------
- if ( verbose .and. (mod(nread,10000) == 0) ) &
- write(*,*)'Processing obs ',nread,' of ',num_obs
+ if (iobs == 1) then
+ if ( .not. get_first_obs(seq, obs1) ) &
+ call error_handler(E_ERR,'obs_seq_coverage', &
+ 'No first observation in sequence.', &
+ source,revision,revdate)
+ else
+ call get_next_obs(seq, obs1, obs2, last_ob_flag)
+ obs1 = obs2
+ endif
+ if ( verbose .and. (mod(iobs,10000) == 0) ) &
+ write(*,*)'Processing obs ',iobs,' of ',num_obs
+
call get_obs_def(obs1, obs_def)
flavor = get_obs_kind( obs_def)
obs_time = get_obs_def_time( obs_def)
@@ -469,7 +469,7 @@
call get_qc( obs1, qc_values)
- if (verbose .and. (nread == 1)) then
+ if (verbose .and. (iobs == 1)) then
call print_time(obs_time,'First observation time')
call print_date(obs_time,'First observation date')
write(*,*) ! whitespace
@@ -478,37 +478,35 @@
if (obs_time > last_possible_time) exit ObsFileLoop
!-----------------------------------------------------------------
+ ! * reject if dart_qc exists and QC is undesirable
+ ! * reject if not a type we want [tracked in obs_type_inds(:)]
! * reject if not in desired region
- ! * reject if not a type we want [tracked in obs_kinds_inds(:)]
- ! * reject if dart_qc exists and is a 4 ...
+ ! * reject if outside vertical region of interest
!-----------------------------------------------------------------
- if ( .not. is_location_in_region(obs_loc,minl,maxl) ) goto 100
- if ( obs_kinds_inds(flavor) <= 0) goto 100
if ( dart_qc_index > 0 ) then
- if (qc_values(dart_qc_index) == 4) goto 100
+ if (qc_values(dart_qc_index) >= 4) cycle ObservationLoop
endif
- ngood = ngood + 1
+ if ( obs_type_inds(flavor) <= 0 ) cycle ObservationLoop
+ if ( .not. is_location_in_region(obs_loc,minl,maxl) ) cycle ObservationLoop
+
+ if ( .not. vertically_desired(obs_loc) ) cycle ObservationLoop
+
+ ngood = ngood + 1
+
! determine if obs is a new location or time at an existing loc
- station_id = find_station_location(flavor, obs_loc, stations)
+ station_id = find_station_location(flavor, obs_loc)
if ( station_id < 1 ) then
- station_id = add_new_station(flavor, obs_loc, stations)
+ station_id = add_new_station(flavor, obs_loc)
endif
- if ( time_is_wanted( obs_time, station_id, stations, timeindex) ) &
- call update_time( obs_time, station_id, stations, timeindex)
+ if ( time_is_wanted( obs_time, station_id, timeindex) ) &
+ call update_time( obs_time, station_id, timeindex)
- 100 continue
-
- call get_next_obs(seq, obs1, obs2, last_ob_flag)
- if (.not. last_ob_flag) obs1 = obs2
-
- if (debug) write(*,*)'obs(',nread,') last_ob_flag ',last_ob_flag
-
!--------------------------------------------------------------------
enddo ObservationLoop
!--------------------------------------------------------------------
@@ -556,8 +554,7 @@
if (num_out_stat < 1) then
write(string1,*)'No location had at least ',nT_minimum,' reporting times.'
- call error_handler(E_ERR, 'obs_seq_coverage', string1, source, revision, &
- revdate, text2=string2)
+ call error_handler(E_ERR, 'obs_seq_coverage', string1, source, revision, revdate)
endif
! Output the file of desired observation locations and times.
@@ -590,36 +587,49 @@
!======================================================================
-function find_station_location(ObsType, ObsLocation, stationlist) result(station_id)
+function find_station_location(ObsType, ObsLocation) result(station_id)
! Simply try to find a matching lat/lon for an observation type
! The lons/lats get yanked around "a lot" - being converted from ASCII radians
! to r8 degrees to r8 radians to r8 degrees and then checked for "equality".
! So - we're actually just checking to see if the lat/lon is within something
! like 500 cm either direction. Seems like a reasonable definition of 'match'.
!
-integer, intent(in) :: ObsType
-type(location_type), intent(in) :: ObsLocation
-type(station), dimension(:), intent(inout) :: stationlist
-integer :: station_id
+! By this point the ObsLocation has one of the mandatory vertical levels.
+integer, intent(in) :: ObsType
+type(location_type), intent(in) :: ObsLocation
+integer :: station_id
+
integer :: i
real(r8), dimension(3) :: obslocarray, stnlocarray
-real(r8) :: londiff, latdiff
+real(r8) :: londiff, latdiff, hgtdiff
-station_id = 0;
+station_id = 0
+if (num_stations == 0) return
+
+obslocarray = get_location(ObsLocation)
+
FindLoop : do i = 1,num_stations
- obslocarray = get_location(ObsLocation)
- stnlocarray = get_location(stationlist(i)%location)
+ if (ObsType /= stations(i)%obs_type) cycle FindLoop
+ stnlocarray = get_location(stations(i)%location)
+
londiff = abs(obslocarray(1) - stnlocarray(1))
latdiff = abs(obslocarray(2) - stnlocarray(2))
+ if (vert_is_pressure(ObsLocation) ) then
+ hgtdiff = abs(obslocarray(3) - stnlocarray(3))
+ else
+ ! if the vertical is undefined or surface, then
+ ! the vertical will always match.
+ hgtdiff = 0.0_r8
+ endif
+
if ( (londiff <= HALF_METER) .and. &
(latdiff <= HALF_METER) .and. &
- (ObsType == stationlist(i)%obs_type) ) then
-
+ (hgtdiff < OnePa) ) then ! within 1 Pa is close enough
station_id = i
exit FindLoop
endif
@@ -632,21 +642,19 @@
!============================================================================
-function add_new_station(ObsType, ObsLocation, stationlist) result(station_id)
+function add_new_station(ObsType, ObsLocation) result(station_id)
! Ugh ... if a new location is found, add it. If the stationlist does not have
! enough space, must copy the info to a temporary list, deallocate/reallocate
! copy the info back, and deallocate the temporary list. Ugh.
-integer, intent(in) :: ObsType
-type(location_type), intent(in) :: ObsLocation
-type(station), allocatable, dimension(:), intent(inout) :: stationlist
-integer :: station_id
+integer, intent(in) :: ObsType
+type(location_type), intent(in) :: ObsLocation
+integer :: station_id
+type(station), allocatable, dimension(:) :: templist
integer :: i
-type(station), allocatable, dimension(:) :: templist
-
if ( num_stations >= max_stations ) then ! need to make room
if (verbose) write(*,*)'Doubling number of possible stations from ', &
@@ -666,43 +674,43 @@
! Copy the information to the temporary space.
DupLoop1 : do i = 1,num_stations
- templist(i)%obs_type = stationlist(i)%obs_type
- templist(i)%location = stationlist(i)%location
- templist(i)%first_time = stationlist(i)%first_time
- templist(i)%last_time = stationlist(i)%last_time
- templist(i)%ntimes = stationlist(i)%ntimes
+ templist(i)%obs_type = stations(i)%obs_type
+ templist(i)%location = stations(i)%location
+ templist(i)%first_time = stations(i)%first_time
+ templist(i)%last_time = stations(i)%last_time
+ templist(i)%ntimes = stations(i)%ntimes
! Make sure the time array is the right size, then copy.
if (associated(templist(i)%times)) then
deallocate( templist(i)%times )
nullify( templist(i)%times )
endif
- allocate( templist(i)%times( size(stationlist(i)%times) ) )
- templist(i)%times = stationlist(i)%times
+ allocate( templist(i)%times( size(stations(i)%times) ) )
+ templist(i)%times = stations(i)%times
enddo DupLoop1
! Deallocate the stations, double the array length,
! allocate the new space.
- call destroy_stations(stationlist)
+ call destroy_stations(stations)
max_stations = 2 * max_stations
- call initialize_stations(max_stations, stationlist)
+ call initialize_stations(max_stations, stations)
! Copy the information BACK to the new space.
DupLoop2 : do i = 1,num_stations
- stationlist(i)%obs_type = templist(i)%obs_type
- stationlist(i)%location = templist(i)%location
- stationlist(i)%first_time = templist(i)%first_time
- stationlist(i)%last_time = templist(i)%last_time
- stationlist(i)%ntimes = templist(i)%ntimes
+ stations(i)%obs_type = templist(i)%obs_type
+ stations(i)%location = templist(i)%location
+ stations(i)%first_time = templist(i)%first_time
+ stations(i)%last_time = templist(i)%last_time
+ stations(i)%ntimes = templist(i)%ntimes
- if (associated(stationlist(i)%times)) then
- deallocate( stationlist(i)%times )
- nullify( stationlist(i)%times )
+ if (associated(stations(i)%times)) then
+ deallocate( stations(i)%times )
+ nullify( stations(i)%times )
endif
- allocate( stationlist(i)%times( size(templist(i)%times) ) )
- stationlist(i)%times = templist(i)%times
+ allocate( stations(i)%times( size(templist(i)%times) ) )
+ stations(i)%times = templist(i)%times
enddo DupLoop2
@@ -711,15 +719,26 @@
endif
-! Add the new station information
+! Add the new station information.
+! Create station with nominal (mandatory) vertical value.
+! The vertically_desired() routine replaces the ob vertical value with
+! the closest mandatory value, so we're good.
num_stations = num_stations + 1
station_id = num_stations
-if (debug) write(*,*)'Adding station ',station_id,' for type ',ObsType
+stations(station_id)%obs_type = ObsType
+stations(station_id)%location = ObsLocation
-stationlist(station_id)%obs_type = ObsType
-stationlist(station_id)%location = ObsLocation
+if (debug) then
+ call write_location(0,ObsLocation,'ascii',string1)
+ call write_location(0,stations(station_id)%location,'ascii',string2)
+ write(*,*)
+ write(*,*)'Added station ',station_id,' for type ',ObsType
+ write(*,*)'observation location', trim(string1)
+ write(*,*)'voxel location', trim(string2)
+ write(*,*)
+endif
end function add_new_station
@@ -727,17 +746,16 @@
!============================================================================
-function time_is_wanted(ObsTime, stationid, stationlist, timeindex)
+function time_is_wanted(ObsTime, stationid, timeindex)
! The station has a list of the observation times closest to the
! verification times. Determine if the observation time is closer to
! the verification time than what we already have.
-type(time_type), intent(in) :: ObsTime
-integer, intent(in) :: stationid
-type(station), dimension(:), intent(in) :: stationlist
-integer, intent(out) :: timeindex
-logical :: time_is_wanted
+type(time_type), intent(in) :: ObsTime
+integer, intent(in) :: stationid
+integer, intent(out) :: timeindex
+logical :: time_is_wanted
type(time_type) :: stndelta, obdelta
integer :: i
@@ -756,12 +774,12 @@
if (obdelta >= half_stride) cycle TimeLoop
! we must be close now ...
- stndelta = stationlist(stationid)%times(i) - all_verif_times(i)
+ stndelta = stations(stationid)%times(i) - all_verif_times(i)
! Check to see if the observation is closer to the verification time
! than the one we have.
if (obdelta < stndelta) then
- if (debug) call print_time(stationlist(stationid)%times(i),'replacing ')
+ if (debug) call print_time(stations(stationid)%times(i),'replacing ')
if (debug) call print_time(ObsTime,'with this observation time')
timeindex = i
time_is_wanted = .TRUE.
@@ -776,38 +794,37 @@
!============================================================================
-subroutine update_time(ObsTime, stationid, stationlist, timeindex)
+subroutine update_time(ObsTime, stationid, timeindex)
! The station has a list of the observation times closest to the
! verification times.
! Add a new time to the station registry.
-type(time_type), intent(in) :: ObsTime
-integer, intent(in) :: stationid
-type(station), dimension(:), intent(inout) :: stationlist
-integer, intent(in) :: timeindex
+type(time_type), intent(in) :: ObsTime
+integer, intent(in) :: stationid
+integer, intent(in) :: timeindex
! Update stuff that seems like a good idea,
! but I don't really know if I'll use it ...
-if ( stationlist(stationid)%ntimes == 0 ) then
- stationlist(stationid)%first_time = ObsTime
- stationlist(stationid)%last_time = ObsTime
+if ( stations(stationid)%ntimes == 0 ) then
+ stations(stationid)%first_time = ObsTime
+ stations(stationid)%last_time = ObsTime
endif
-if ( stationlist(stationid)%first_time > ObsTime ) &
- stationlist(stationid)%first_time = ObsTime
-if ( stationlist(stationid)%last_time < ObsTime ) &
- stationlist(stationid)%last_time = ObsTime
+if ( stations(stationid)%first_time > ObsTime ) &
+ stations(stationid)%first_time = ObsTime
+if ( stations(stationid)%last_time < ObsTime ) &
+ stations(stationid)%last_time = ObsTime
if (debug) write(*,*)'Stuffing time into station ',stationid,' at timestep ', timeindex
! as long as ntimes /= 0 we are OK.
! When the stations get written to the netCDF file, count the
! number of non-zero times in the times array for a real count.
-stationlist(stationid)%ntimes = stationlist(stationid)%ntimes + 1
+stations(stationid)%ntimes = stations(stationid)%ntimes + 1
! Stuff the time in the appropriate slot ... finally.
-stationlist(stationid)%times(timeindex) = ObsTime
+stations(stationid)%times(timeindex) = ObsTime
end subroutine update_time
@@ -823,6 +840,7 @@
integer :: LineLenDimID, nlinesDimID, stringDimID
integer :: TimeDimID, StationsDimID, FcstDimID, VerifyDimID
integer :: VarID, FcstVarID, VerifVarID, ExperimentVarID
+integer :: nlevDimID, plevelVarID
character(len=8) :: crdate ! needed by F90 DATE_AND_TIME intrinsic
character(len=10) :: crtime ! needed by F90 DATE_AND_TIME intrinsic
@@ -847,9 +865,6 @@
call nc_check(nf90_create(path = trim(fname), cmode = nf90_clobber, &
ncid = ncid), 'obs_seq_coverage:InitNetCDF', 'create '//trim(fname))
-if (debug) write(string1,*)trim(fname), ' is fortran unit ',ncid
-call error_handler(E_MSG,'InitNetCDF',string1,source,revision,revdate)
-
!----------------------------------------------------------------------------
! Write Global Attributes (mostly namelist input, that sort of thing)
!----------------------------------------------------------------------------
@@ -878,9 +893,9 @@
! Write all desired observation types.
! As a sanity check - do it from our working array.
ntypes = 0
-TYPELOOP : do i = 1,size(obs_kinds_inds)
+TYPELOOP : do i = 1,size(obs_type_inds)
- if (obs_kinds_inds(i) < 1) cycle TYPELOOP
+ if (obs_type_inds(i) < 1) cycle TYPELOOP
ntypes = ntypes + 1
@@ -990,6 +1005,19 @@
call nc_check(nf90_put_att(ncid, ExperimentVarID, 'cols', 'each verification time'), &
'InitNetCDF', 'experiment:put_att cols')
+! the number of levels
+
+call nc_check(nf90_def_dim(ncid=ncid, &
+ name='nlevels', len= NUM_MANDATORY_LEVELS, dimid=nlevDimID), &
+ 'InitNetCDF', 'def_dim:nlevels '//trim(fname))
+call nc_check(nf90_def_var(ncid=ncid, name='plevel', xtype=nf90_real, &
+ dimids = (/ nlevDimID /), varid=plevelVarID), &
+ 'InitNetCDF', 'plevel:def_var')
+call nc_check(nf90_put_att(ncid, plevelVarID, 'long_name', 'mandatory pressure levels'), &
+ 'InitNetCDF', 'plevel:long_name')
+call nc_check(nf90_put_att(ncid, plevelVarID, 'units', 'Pa'), &
+ 'InitNetCDF', 'plevel:units')
+
! write all namelist quantities
call find_textfile_dims('input.nml', nlines, linelen)
@@ -1167,6 +1195,8 @@
call nc_check(nf90_put_var(ncid, ExperimentVarID, experiment_Tr8 ), &
'InitNetCDF', 'put_var:verification_times')
+call nc_check(nf90_put_var(ncid, plevelVarID, mandatory_levels ), &
+ 'InitNetCDF', 'put_var:plevel')
!----------------------------------------------------------------------------
! Finish up ...
!----------------------------------------------------------------------------
@@ -1181,7 +1211,7 @@
!============================================================================
-Subroutine WriteNetCDF(ncid, fname, stations)
+subroutine WriteNetCDF(ncid, fname, stations)
integer, intent(in) :: ncid
character(len=*), intent(in) :: fname
type(station), dimension(:), intent(in) :: stations
@@ -1296,13 +1326,13 @@
call nc_check(nf90_sync( ncid), 'WriteNetCDF', 'sync '//trim(fname))
-end Subroutine WriteNetCDF
+end subroutine WriteNetCDF
!============================================================================
-Subroutine CloseNetCDF(ncid, fname)
+subroutine CloseNetCDF(ncid, fname)
integer, intent(in) :: ncid
character(len=*), intent(in) :: fname
@@ -1311,34 +1341,28 @@
call nc_check(nf90_sync( ncid), 'CloseNetCDF', 'sync '//trim(fname))
call nc_check(nf90_close(ncid), 'CloseNetCDF', 'close '//trim(fname))
-end Subroutine CloseNetCDF
+end subroutine CloseNetCDF
!============================================================================
-subroutine initialize_stations(Nstations, stations)
+subroutine initialize_stations(Nstations, mystations)
integer, intent(in) :: Nstations
-type(station), allocatable, dimension(:), intent(out) :: stations
+type(station), allocatable, dimension(:), intent(out) :: mystations
-integer :: i,N
+integer :: i
-allocate(stations(Nstations))
-N = size(stations)
+allocate(mystations(Nstations))
-if (N /= Nstations) then
- write(string1,*)'wanted ',Nstations,' allocated ',N
- call error_handler(E_ERR, 'initialize_stations', string1, source, revision, revdate)
-endif
-
-do i = 1,N
- stations(i)%obs_type = 0
- stations(i)%location = set_location_missing()
- stations(i)%ntimes = 0
- allocate( stations(i)%times( num_verification_times ) )
- stations(i)%first_time = no_time
- stations(i)%last_time = no_time
- stations(i)%times = no_time
+do i = 1,Nstations
+ mystations(i)%obs_type = 0
+ mystations(i)%location = set_location_missing()
+ mystations(i)%ntimes = 0
+ allocate( mystations(i)%times( num_verification_times ) )
+ mystations(i)%first_time = no_time
+ mystations(i)%last_time = no_time
+ mystations(i)%times = no_time
enddo
end subroutine initialize_stations
@@ -1347,21 +1371,21 @@
!============================================================================
-subroutine destroy_stations(stations)
-type(station), allocatable, dimension(:), intent(inout) :: stations
+subroutine destroy_stations(mystations)
+type(station), allocatable, dimension(:), intent(inout) :: mystations
integer :: i,N
-N = size(stations)
+N = size(mystations)
do i = 1,N
- if (associated(stations(i)%times)) then
- deallocate( stations(i)%times )
- nullify( stations(i)%times )
+ if (associated(mystations(i)%times)) then
+ deallocate( mystations(i)%times )
+ nullify( mystations(i)%times )
endif
enddo
-if (allocated(stations)) deallocate(stations)
+if (allocated(mystations)) deallocate(mystations)
end subroutine destroy_stations
@@ -1390,7 +1414,7 @@
! and times fit the requirements.
write(iunit,*)'num_definitions ',num_out_total
-call write_obs_kind(iunit, fform='formatted', use_list=obs_kinds_inds)
+call write_obs_kind(iunit, fform='formatted', use_list=obs_type_inds)
call set_obs_def_error_variance(obs_def, MISSING_R8)
write(*,*) ! whitespace
@@ -1401,17 +1425,17 @@
write(*,*) ! whitespace
if ( debug ) then
- TYPELOOP : do i = 1,size(obs_kinds_inds)
- if (obs_kinds_inds(i) < 1) cycle TYPELOOP
+ TYPELOOP : do i = 1,size(obs_type_inds)
+ if (obs_type_inds(i) < 1) cycle TYPELOOP
string2 = adjustl(get_obs_kind_name(i))
- write(*,*)'i,obs_kind_inds(i)',i,obs_kinds_inds(i),trim(string2)
+ write(*,*)'i,obs_type_inds(i)',i,obs_type_inds(i),trim(string2)
enddo TYPELOOP
write(*,*)
endif
-Summary : do i = 1,num_stations
+Selections : do i = 1,num_stations
- if ( .not. DesiredStations(i) ) cycle Summary
+ if ( .not. DesiredStations(i) ) cycle Selections
call set_obs_def_kind( obs_def, stations(i)%obs_type)
call set_obs_def_location(obs_def, stations(i)%location)
@@ -1424,13 +1448,14 @@
enddo TimeLoop
if (verbose) then
- call get_time(stations(i)%first_time,sec1,day1)
- call get_time(stations(i)%last_time, secN,dayN)
- write(*,'(''station '',i6,'' has '',i3,'' obs between ['',i7,1x,i5,'' and '',i7,1x,i5,'']'')') &
+ call get_time(stations(i)%first_time,sec1,day1)
+ call get_time(stations(i)%last_time, secN,dayN)
+ write(*,'(''station '',i6,'' has '',i3,'' obs between ['',&
+ &i7,1x,i5,'' and '',i7,1x,i5,'']'')') &
i,stations(i)%ntimes,day1,sec1,dayN,secN
endif
-enddo Summary
+enddo Selections
close(iunit)
@@ -1481,7 +1506,6 @@
type(time_type) :: T1, TN, TimeMax, thistime, nexttime
integer :: i, j, nsteps, seconds, days
-real(digits12) :: steps
verification_stride = set_time(v_int_seconds,0) ! global variable
half_stride = verification_stride / 2
@@ -1501,9 +1525,8 @@
! Check to make sure the forecast length is a multiple
! of the verification interval.
-steps = flen / verification_stride ! time_manager_mod:time_real_divide
-nsteps = nint(steps)
-nexttime = verification_stride * nsteps ! time_manager_mod:time_scalar_divide
+nsteps = flen / verification_stride ! time_manager_mod:time_divide
+nexttime = verification_stride * nsteps ! time_manager_mod:time_scalar_mult
if (nexttime /= flen) then
call print_time(verification_stride,'verification stride')
@@ -1533,9 +1556,8 @@
if (TN /= T1) then
thistime = TN - T1
- steps = thistime / verification_stride ! time_manager_mod:time_real_divide
- nsteps = nint(steps)
- nexttime = verification_stride * nsteps ! time_manager_mod:time_scalar_divide
+ nsteps = thistime / verification_stride
+ nexttime = verification_stride * nsteps
if (nexttime /= thistime) then
call print_time(verification_stride,'verification stride')
@@ -1562,8 +1584,7 @@
TimeMax = TN + flen ! The time of the last verification
thistime = TimeMax - T1
-steps = thistime / verification_stride
-nsteps = nint(steps)
+nsteps = thistime / verification_stride
nexttime = verification_stride * nsteps
if (nexttime /= thistime) then
@@ -1577,9 +1598,6 @@
num_verification_times = nsteps+1 ! SET GLOBAL VALUE
if (verbose) write(*,*)'There are ',num_verification_times,' total times we need observations.'
-if (allocated(all_verif_times)) deallocate(all_verif_times)
-if (allocated( verification_times)) deallocate( verification_times)
-if (allocated( experiment_Tr8)) deallocate( experiment_Tr8)
allocate(all_verif_times(num_verification_times))
allocate(verification_times(num_analyses,num_verify_per_fcst))
allocate(experiment_Tr8(num_verify_per_fcst,num_analyses)) ! TRANSPOSED for netCDF
@@ -1592,8 +1610,8 @@
! Use those desired times to fill the times for each forecast experiment
do i = 1,num_analyses ! SET GLOBAL VALUE
- verification_times(i, 1:num_verify_per_fcst) = &
- all_verif_times(i:(i+num_verify_per_fcst-1))
+ verification_times(i, 1:num_verify_per_fcst) = &
+ all_verif_times(i:(i+num_verify_per_fcst-1))
enddo
! Use those to create the matching array of real(digits12)
@@ -1645,11 +1663,11 @@
enddo QCMetaDataLoop
if ( qcindex < 0 ) then
- write(string1,*)'metadata:Quality Control copyindex not found'
+ write(string1,*)'metadata:source Quality Control copyindex not found'
call error_handler(E_MSG,'find_qc_indices',string1,source,revision,revdate)
endif
if ( dart_qcindex < 0 ) then
- write(string1,*)'metadata:DART quality control copyindex not found'
+ write(string1,*)'metadata:DART Quality Control copyindex not found'
call error_handler(E_MSG,'find_qc_indices',string1,source,revision,revdate)
endif
@@ -1667,6 +1685,98 @@
!============================================================================
+subroutine setPressureLevels()
+
+! We are using the mandatory levels as defined in:
+! http://www.meteor.wisc.edu/~hopkins/aos100/raobdoc.htm
+!
+! There are 14 MANDATORY pressure levels for radiosonde observations:
+! 1000, 925, 850, 700, 500, 400, 300, 250, 200, 150, 100, 70, 50, and 10 (hPa)
+
+mandatory_levels = (/ 1000.0_r8, 925.0_r8, 850.0_r8, 700.0_r8, 500.0_r8, &
+ 400.0_r8, 300.0_r8, 250.0_r8, 200.0_r8, 150.0_r8, &
+ 100.0_r8, 70.0_r8, 50.0_r8, 10.0_r8/)
+
+mandatory_levels = mandatory_levels * 100.0_r8 ! convert hPa to Pa
+
+end subroutine setPressureLevels
+
+
+!======================================================================
+
+
+function vertically_desired(location)
+
+! For obs on pressure levels ... how close to the mandatory levels is close
+! enough? A quick examination of 3000+ radiosonde obs revealed that all the
+! pressures were recorded with a precision of 10 Pa (e.g. 110.0, 120.0, ...)
+!
+! So if the observation is within 1 Pa of a mandatory level, that's good enough.
+
+type(location_type), intent(inout) :: location
+logical :: vertically_desired
+
+integer :: iz
+real(r8), dimension(4) :: obslocarray
+real(r8) :: zdist
+
+if (vert_is_undef( location) .or. &
+ vert_is_surface(location)) then
+
+ vertically_desired = .true.
+
+elseif ( vert_is_pressure(location) ) then
+
+ obslocarray(1:3) = get_location(location)
+ obslocarray( 4 ) = query_location(location,'which_vert')
+
+ iz = FindClosestPressureLevel(obslocarray(3))
+
+ zdist = abs(obslocarray(3) - mandatory_levels(iz))
+
+ if ( zdist < OnePa ) then
+
+ obslocarray(3) = mandatory_levels(iz)
+ location = set_location( obslocarray )
+ vertically_desired = .true.
+
+ else ! not close enough
+
+ vertically_desired = .false.
+
+ endif
+
+else
+
+ ! do not support height, scale_height, or level
+ vertically_desired = .false.
+
+endif
+
+end function vertically_desired
+
+
+!======================================================================
+
+
+function FindClosestPressureLevel(obslevel)
+real(r8), intent(in) :: obslevel
+integer :: FindClosestPressureLevel
+
+real(r8), dimension(NUM_MANDATORY_LEVELS) :: distances
+integer, dimension(1) :: iz
+
+distances = abs(obslevel - mandatory_levels)
+iz = minloc(distances)
+
+FindClosestPressureLevel = iz(1)
+
+end function FindClosestPressureLevel
+
+
+!======================================================================
+
+
end program obs_seq_coverage
! <next few lines under version control, do not edit>
More information about the Dart-dev
mailing list