[Dart-dev] [7092] DART/trunk: The obs_seq_coverage routine now handles mandatory pressre levels as well
nancy at ucar.edu
nancy at ucar.edu
Tue Jul 29 16:35:31 MDT 2014
Revision: 7092
Author: thoar
Date: 2014-07-29 16:35:31 -0600 (Tue, 29 Jul 2014)
Log Message:
-----------
The obs_seq_coverage routine now handles mandatory pressre levels as well
as surface and undefined vertical level quantities.
It writes out 'voxels' where the vertical thickness of the voxel is 1Pa for
all pressure-level observations. There are 14 mandatory pressure levels.
Only those 14 are possible.
Modified Paths:
--------------
DART/trunk/models/mpas_atm/work/input.nml
DART/trunk/obs_sequence/obs_seq_coverage.f90
DART/trunk/obs_sequence/obs_seq_coverage.html
DART/trunk/obs_sequence/obs_seq_coverage.nml
-------------- next part --------------
Modified: DART/trunk/models/mpas_atm/work/input.nml
===================================================================
--- DART/trunk/models/mpas_atm/work/input.nml 2014-07-28 20:29:08 UTC (rev 7091)
+++ DART/trunk/models/mpas_atm/work/input.nml 2014-07-29 22:35:31 UTC (rev 7092)
@@ -130,8 +130,8 @@
/
&obs_kind_nml
- assimilate_these_obs_types = 'RADIOSONDE_TEMPERATURE'
- 'RADIOSONDE_U_WIND_COMPONENT'
+ assimilate_these_obs_types = 'RADIOSONDE_TEMPERATURE',
+ 'RADIOSONDE_U_WIND_COMPONENT',
'RADIOSONDE_V_WIND_COMPONENT'
/
@@ -158,35 +158,35 @@
highest_obs_pressure_mb = -1.0
/
-! 'theta', 'KIND_POTENTIAL_TEMPERATURE'
-! 'surface_pressure', 'KIND_SURFACE_PRESSURE'
-! 'uReconstructZonal', 'KIND_U_WIND_COMPONENT'
-! 'uReconstructMeridional','KIND_V_WIND_COMPONENT'
-! 'u', 'KIND_EDGE_NORMAL_SPEED'
-! 'w', 'KIND_VERTICAL_VELOCITY'
-! 'rho', 'KIND_DENSITY'
-! 'qv', 'KIND_VAPOR_MIXING_RATIO'
-! 'qc', 'KIND_CLOUDWATER_MIXING_RATIO'
-! 'qr', 'KIND_RAINWATER_MIXING_RATIO'
-! 'qi', 'KIND_ICE_MIXING_RATIO'
-! 'qs', 'KIND_SNOW_MIXING_RATIO'
-! 'qg', 'KIND_GRAUPEL_MIXING_RATIO'
-! 'rho', 'KIND_DENSITY'
-! 'salinity', 'KIND_SALINITY'
+! 'theta', 'KIND_POTENTIAL_TEMPERATURE',
+! 'surface_pressure', 'KIND_SURFACE_PRESSURE',
+! 'uReconstructZonal', 'KIND_U_WIND_COMPONENT',
+! 'uReconstructMeridional','KIND_V_WIND_COMPONENT',
+! 'u', 'KIND_EDGE_NORMAL_SPEED',
+! 'w', 'KIND_VERTICAL_VELOCITY',
+! 'rho', 'KIND_DENSITY',
+! 'qv', 'KIND_VAPOR_MIXING_RATIO',
+! 'qc', 'KIND_CLOUDWATER_MIXING_RATIO',
+! 'qr', 'KIND_RAINWATER_MIXING_RATIO',
+! 'qi', 'KIND_ICE_MIXING_RATIO',
+! 'qs', 'KIND_SNOW_MIXING_RATIO',
+! 'qg', 'KIND_GRAUPEL_MIXING_RATIO',
+! 'rho', 'KIND_DENSITY',
+! 'salinity', 'KIND_SALINITY',
! 'temperature', 'KIND_TEMPERATURE'
&mpas_vars_nml
mpas_state_variables =
- 'uReconstructZonal', 'KIND_U_WIND_COMPONENT'
- 'uReconstructMeridional','KIND_V_WIND_COMPONENT'
- 'theta', 'KIND_POTENTIAL_TEMPERATURE'
- 'surface_pressure', 'KIND_SURFACE_PRESSURE'
- 'w', 'KIND_VERTICAL_VELOCITY'
- 'rho', 'KIND_DENSITY'
+ 'uReconstructZonal', 'KIND_U_WIND_COMPONENT',
+ 'uReconstructMeridional','KIND_V_WIND_COMPONENT',
+ 'theta', 'KIND_POTENTIAL_TEMPERATURE',
+ 'surface_pressure', 'KIND_SURFACE_PRESSURE',
+ 'w', 'KIND_VERTICAL_VELOCITY',
+ 'rho', 'KIND_DENSITY',
'qv', 'KIND_VAPOR_MIXING_RATIO'
- mpas_state_bounds = 'qv','0.0','NULL','CLAMP'
- 'qc','0.0','NULL','CLAMP'
+ mpas_state_bounds = 'qv','0.0','NULL','CLAMP',
+ 'qc','0.0','NULL','CLAMP',
'qr','0.0','NULL','CLAMP'
/
@@ -215,11 +215,12 @@
output_obs_def_mod_file = '../../../obs_def/obs_def_mod.f90'
input_obs_kind_mod_file = '../../../obs_kind/DEFAULT_obs_kind_mod.F90'
output_obs_kind_mod_file = '../../../obs_kind/obs_kind_mod.f90'
- input_files = '../../../obs_def/obs_def_reanalysis_bufr_mod.f90'
- '../../../obs_def/obs_def_altimeter_mod.f90'
- '../../../obs_def/obs_def_gts_mod.f90'
- '../../../obs_def/obs_def_metar_mod.f90'
- '../../../obs_def/obs_def_gps_mod.f90'
+ input_files = '../../../obs_def/obs_def_reanalysis_bufr_mod.f90',
+ '../../../obs_def/obs_def_altimeter_mod.f90',
+ '../../../obs_def/obs_def_gts_mod.f90',
+ '../../../obs_def/obs_def_metar_mod.f90',
+ '../../../obs_def/obs_def_gps_mod.f90',
+ '../../../obs_def/obs_def_dew_point_mod.f90'
/
&obs_sequence_tool_nml
@@ -394,12 +395,12 @@
textfile_out = 'obsdef_mask.txt'
netcdf_out = 'obsdef_mask.nc'
calendar = 'Gregorian'
- first_analysis = 2008, 8, 1, 18, 0, 0
- last_analysis = 2008, 8, 1, 18, 0, 0
- forecast_length_days = 0
- forecast_length_seconds = 21600
- verification_interval_seconds = 21600
- temporal_coverage_percent = 50.0
+ first_analysis = 2008, 6, 1, 00, 0, 0
+ last_analysis = 2008, 6, 1, 12, 0, 0
+ forecast_length_days = 3
+ forecast_length_seconds = 0
+ verification_interval_seconds = 43200
+ temporal_coverage_percent = 90.0
lonlim1 = 0.0
lonlim2 = 360.0
latlim1 = -90.0
Modified: DART/trunk/obs_sequence/obs_seq_coverage.f90
===================================================================
--- DART/trunk/obs_sequence/obs_seq_coverage.f90 2014-07-28 20:29:08 UTC (rev 7091)
+++ DART/trunk/obs_sequence/obs_seq_coverage.f90 2014-07-29 22:35:31 UTC (rev 7092)
@@ -72,20 +72,22 @@
!---------------------------------------------------------------------
!---------------------------------------------------------------------
-type voxel
+type voxel_type
integer :: obs_type
type(location_type) :: location
+ integer :: levelindx
type(time_type) :: first_time
type(time_type) :: last_time
integer :: ntimes
type(time_type), pointer :: times(:)
-end type voxel
+end type voxel_type
-logical, allocatable, dimension(:) :: Desiredvoxels
-type(voxel), allocatable, dimension(:) :: voxels
-integer :: num_voxels ! This is the current number of unique locations
-integer :: max_voxels ! This is the largest possible number of uniq locs
-integer :: voxel_id ! the index (into voxels) of an existing location
+logical, allocatable, dimension(:) :: Desiredvoxels
+type(voxel_type), allocatable, dimension(:) :: voxels
+
+integer :: num_voxels ! This is the current number of unique locations
+integer :: max_voxels ! This is the largest possible number of uniq locs
+integer :: voxel_id ! the index (into voxels) of an existing location
integer :: timeindex ! the index (into the time array of a voxel)
integer :: num_out_stat ! total number of desired voxels found
integer :: num_out_total ! total number of desired locations * times found
@@ -183,6 +185,7 @@
integer :: num_verify_per_fcst
integer :: num_verification_times ! number of verification times - total
integer :: nT_minimum ! will settle for this many verif times - total
+integer :: ilev ! index of mandatory level
real(r8), dimension(NUM_MANDATORY_LEVELS) :: mandatory_levels = MISSING_R8 ! pressure level (hPa)
@@ -261,33 +264,27 @@
source, revision, revdate, text2=string2)
endif
-
! Set the verification time array (global storage)
call set_required_times(first_analysis, last_analysis, &
forecast_length_days, forecast_length_seconds, &
verification_interval_seconds, temporal_coverage_percent)
-if (verbose) then
+write(*,*) ! whitespace
+write(*,*)'The analysis times (the start of the forecasts) are:'
+do i=1,size(verification_times,1)
+ write(string1,*)'analysis # ',i,' at '
+ call print_date(verification_times(i,1),trim(string1))
+enddo
+write(*,*) ! whitespace
+write(*,*)'At least',nT_minimum,' observations are required from the following:'
+do i=1,num_verification_times
+ write(string1,*)'verification # ',i,' at '
+ call print_date(all_verif_times(i),trim(string1))
+ call print_time(all_verif_times(i),trim(string1))
+enddo
+write(*,*) ! whitespace
- write(*,*) ! whitespace
- write(*,*)'The analysis times (the start of the forecasts) are:'
- do i=1,size(verification_times,1)
- write(string1,*)'analysis # ',i,' at '
- call print_date(verification_times(i,1),trim(string1))
- enddo
-
- write(*,*) ! whitespace
- write(*,*)'At least',nT_minimum,' observations times are required during:'
- do i=1,num_verification_times
- write(string1,*)'verification # ',i,' at '
- call print_date(all_verif_times(i),trim(string1))
- enddo
-
-
- write(*,*) ! whitespace
-endif
-
last_possible_time = all_verif_times(num_verification_times) + half_stride
! Allocate a hunk of voxels. If we fill this up, we will
@@ -491,7 +488,7 @@
if ( .not. is_location_in_region(obs_loc,minl,maxl) ) cycle ObservationLoop
- if ( .not. vertically_desired(obs_loc) ) cycle ObservationLoop
+ if ( .not. vertically_desired(obs_loc, ilev) ) cycle ObservationLoop
ngood = ngood + 1
@@ -500,9 +497,14 @@
voxel_id = find_voxel_location(flavor, obs_loc)
if ( voxel_id < 1 ) then
- voxel_id = add_new_voxel(flavor, obs_loc)
+ voxel_id = add_new_voxel(flavor, obs_loc, ilev)
endif
+ if (debug) then
+ call write_location(0,obs_loc,'ascii',string1)
+ write(*,*)'observation',iobs,'is voxel',voxel_id,trim(string1)
+ endif
+
if ( time_is_wanted( obs_time, voxel_id, timeindex) ) &
call update_time( obs_time, voxel_id, timeindex)
@@ -575,7 +577,7 @@
if (allocated(module_obs_copy_names)) deallocate(module_obs_copy_names)
if (allocated(module_qc_copy_names )) deallocate(module_qc_copy_names )
if (allocated(obs_seq_filenames)) deallocate(obs_seq_filenames)
-if (allocated(Desiredvoxels)) deallocate(Desiredvoxels)
+if (allocated(Desiredvoxels)) deallocate(Desiredvoxels)
call error_handler(E_MSG,'obs_seq_coverage','Finished successfully.',source,revision,revdate)
call finalize_utilities()
@@ -641,7 +643,7 @@
!============================================================================
-function add_new_voxel(ObsType, ObsLocation) result(voxel_id)
+function add_new_voxel(ObsType, ObsLocation, ilevel) result(voxel_id)
! Ugh ... if a new location is found, add it. If the voxellist does not have
! enough space, must copy the info to a temporary list, deallocate/reallocate
@@ -649,9 +651,10 @@
integer, intent(in) :: ObsType
type(location_type), intent(in) :: ObsLocation
+integer, intent(in) :: ilevel
integer :: voxel_id
-type(voxel), allocatable, dimension(:) :: templist
+type(voxel_type), allocatable, dimension(:) :: templist
integer :: i
if ( num_voxels >= max_voxels ) then ! need to make room
@@ -675,6 +678,7 @@
templist(i)%obs_type = voxels(i)%obs_type
templist(i)%location = voxels(i)%location
+ templist(i)%levelindx = voxels(i)%levelindx
templist(i)%first_time = voxels(i)%first_time
templist(i)%last_time = voxels(i)%last_time
templist(i)%ntimes = voxels(i)%ntimes
@@ -700,6 +704,7 @@
voxels(i)%obs_type = templist(i)%obs_type
voxels(i)%location = templist(i)%location
+ voxels(i)%levelindx = templist(i)%levelindx
voxels(i)%first_time = templist(i)%first_time
voxels(i)%last_time = templist(i)%last_time
voxels(i)%ntimes = templist(i)%ntimes
@@ -726,8 +731,9 @@
num_voxels = num_voxels + 1
voxel_id = num_voxels
-voxels(voxel_id)%obs_type = ObsType
-voxels(voxel_id)%location = ObsLocation
+voxels(voxel_id)%obs_type = ObsType
+voxels(voxel_id)%location = ObsLocation
+voxels(voxel_id)%levelindx = ilevel
if (debug) then
call write_location(0,ObsLocation,'ascii',string1)
@@ -736,6 +742,7 @@
write(*,*)'Added voxel ',voxel_id,' for type ',ObsType
write(*,*)'observation location', trim(string1)
write(*,*)'voxel location', trim(string2)
+ write(*,*)'voxel level', voxels(voxel_id)%levelindx
write(*,*)
endif
@@ -1009,13 +1016,13 @@
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, &
+call nc_check(nf90_def_var(ncid=ncid, name='mandatory_level', xtype=nf90_real, &
dimids = (/ nlevDimID /), varid=plevelVarID), &
- 'InitNetCDF', 'plevel:def_var')
+ 'InitNetCDF', 'mandatory_level:def_var')
call nc_check(nf90_put_att(ncid, plevelVarID, 'long_name', 'mandatory pressure levels'), &
- 'InitNetCDF', 'plevel:long_name')
+ 'InitNetCDF', 'mandatory_level:long_name')
call nc_check(nf90_put_att(ncid, plevelVarID, 'units', 'Pa'), &
- 'InitNetCDF', 'plevel:units')
+ 'InitNetCDF', 'mandatory_level:units')
! write all namelist quantities
@@ -1063,6 +1070,18 @@
call error_handler(E_ERR,'InitNetCDF',string1,source,revision,revdate)
endif
+! Define the mandatory level corresponding to each voxel
+
+call nc_check(nf90_def_var(ncid=ncid, name='voxel_level_index', xtype=nf90_int, &
+ dimids=(/ voxelsDimID /), varid=VarID), &
+ 'InitNetCDF', 'voxel_level_index:def_var')
+call nc_check(nf90_put_att(ncid, VarID, 'long_name', &
+ 'index of the mandatory level of this voxel'), &
+ 'InitNetCDF', 'voxel_level_index:put_att long_name')
+call nc_check(nf90_put_att(ncid, VarID, 'valid_range', &
+ (/ 1, NUM_MANDATORY_LEVELS /) ), &
+ 'InitNetCDF', 'voxel_level_index:put_att valid_range')
+
! Define the number of observation times
call nc_check(nf90_def_var(ncid=ncid, name='ntimes', xtype=nf90_int, &
@@ -1196,6 +1215,7 @@
call nc_check(nf90_put_var(ncid, plevelVarID, mandatory_levels ), &
'InitNetCDF', 'put_var:plevel')
+
!----------------------------------------------------------------------------
! Finish up ...
!----------------------------------------------------------------------------
@@ -1213,14 +1233,14 @@
subroutine WriteNetCDF(ncid, fname, voxels)
integer, intent(in) :: ncid
character(len=*), intent(in) :: fname
-type(voxel), dimension(:), intent(in) :: voxels
+type(voxel_type), dimension(:), intent(in) :: voxels
integer :: DimID, ntimes, voxelindex, days, secs, i
integer, dimension(1) :: istart, icount
integer :: voxelVarID, TimeVarID, NTimesVarID, &
T1VarID, TNVarID, ObsTypeVarID, &
- LocationVarID, WhichVertVarID
+ LocationVarID, WhichVertVarID, IlevVarID
real(digits12), allocatable, dimension(:) :: mytimes
integer, dimension(num_voxels) :: gooduns ! Cray compiler likes this better
@@ -1260,6 +1280,9 @@
call nc_check(nf90_inq_varid(ncid, 'last_time', varid=TNVarID), &
'WriteNetCDF', 'inq_varid:last_time '//trim(fname))
+call nc_check(nf90_inq_varid(ncid, 'voxel_level_index', varid=IlevVarID), &
+ 'WriteNetCDF', 'inq_varid:voxel_level_index '//trim(fname))
+
call nc_get_location_varids(ncid, fname, LocationVarID, WhichVertVarID)
allocate(mytimes(ntimes))
@@ -1296,6 +1319,9 @@
call nc_check(nf90_put_var(ncid, NTimesVarId, (/ voxels(voxelindex)%ntimes /), &
start=istart, count=icount), 'WriteNetCDF', 'put_var:ntimes')
+ call nc_check(nf90_put_var(ncid, IlevVarId, (/ voxels(voxelindex)%levelindx /), &
+ start=istart, count=icount), 'WriteNetCDF', 'put_var:voxel_level_index')
+
!----------------------------------------------------------------------------
! time : fill, write
!----------------------------------------------------------------------------
@@ -1348,7 +1374,7 @@
subroutine initialize_voxels(Nvoxels, myvoxels)
integer, intent(in) :: Nvoxels
-type(voxel), allocatable, dimension(:), intent(out) :: myvoxels
+type(voxel_type), allocatable, dimension(:), intent(out) :: myvoxels
integer :: i
@@ -1371,7 +1397,7 @@
subroutine destroy_voxels(myvoxels)
-type(voxel), allocatable, dimension(:), intent(inout) :: myvoxels
+type(voxel_type), allocatable, dimension(:), intent(inout) :: myvoxels
integer :: i,N
@@ -1594,9 +1620,18 @@
call error_handler(E_ERR,'set_required_times',string1,source,revision,revdate)
endif
-num_verification_times = nsteps+1 ! SET GLOBAL VALUE
-if (verbose) write(*,*)'There are ',num_verification_times,' total times we need observations.'
+num_verification_times = nsteps + 1 ! SET GLOBAL VALUE
+nT_minimum = nint(num_verification_times * coverage_pcnt/100.0_r8) ! SET GLOBAL VALUE
+if (verbose) then
+ write(string1,*)'Need ',num_verification_times, &
+ ' verification times for 100.00% coverage.'
+ write(string2,'(''at least '',i5,'' verification times for'',&
+ & f7.2,''% coverage.'')') nT_minimum, coverage_pcnt
+
+ call error_handler(E_MSG,'set_required_times',string1,text2=string2)
+endif
+
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
@@ -1622,18 +1657,7 @@
enddo
enddo
-nT_minimum = nint(num_verification_times * coverage_pcnt/100.0_r8) ! SET GLOBAL VALUE
-if (debug) then
- write(*,*) ! whitespace
- write(*,*)'At least',nT_minimum,' observations times are required at:'
- do i=1,num_verification_times
- write(string1,*)'verification # ',i,' at '
- call print_date(all_verif_times(i),trim(string1))
- enddo
- write(*,*) ! whitespace
-endif
-
end subroutine set_required_times
@@ -1704,7 +1728,7 @@
!======================================================================
-function vertically_desired(location)
+function vertically_desired(location,ilevel)
! 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
@@ -1713,16 +1737,20 @@
! So if the observation is within 1 Pa of a mandatory level, that's good enough.
type(location_type), intent(inout) :: location
+integer, intent(out) :: ilevel
logical :: vertically_desired
integer :: iz
real(r8), dimension(4) :: obslocarray
real(r8) :: zdist
+ilevel = MISSING_I
+
if (vert_is_undef( location) .or. &
vert_is_surface(location)) then
vertically_desired = .true.
+ ilevel = 1
elseif ( vert_is_pressure(location) ) then
@@ -1738,6 +1766,7 @@
obslocarray(3) = mandatory_levels(iz)
location = set_location( obslocarray )
vertically_desired = .true.
+ ilevel = iz
else ! not close enough
Modified: DART/trunk/obs_sequence/obs_seq_coverage.html
===================================================================
--- DART/trunk/obs_sequence/obs_seq_coverage.html 2014-07-28 20:29:08 UTC (rev 7091)
+++ DART/trunk/obs_sequence/obs_seq_coverage.html 2014-07-29 22:35:31 UTC (rev 7092)
@@ -41,7 +41,7 @@
sequence files to determine which observation locations report
frequently enough to be useful for a verification study. The big picture
is to be able to pare down a large set of observations into a compact
-observation sequence file to run through <em class=program>filter</em>
+observation sequence file to run through <a href="../filter/filter.html">filter</a>
with all of the intended observation types flagged as
<em class=option>evaluate_only</em>. DART's forward operators then
get applied and all the forecasts are preserved in a
@@ -89,6 +89,7 @@
through time boils down to the simple statement that if the observation is within
about 500cm of a previous observation, they are treated as co-located observations.
For some very high resolution applications, this may be insufficient, but there it is.
+For observations at pressure levels, see the <a href="#Usage">Word about vertical levels</a>.
<br />
<br />
The only complicated part of determining the verification network is the
@@ -185,24 +186,24 @@
<div class=namelist>
<pre>
&obs_seq_coverage_nml
- obs_sequence_list = 'obs_coverage_list.txt',
- obs_sequence_name = '',
- obs_of_interest = '',
- textfile_out = 'obsdef_mask.txt',
- netcdf_out = 'obsdef_mask.nc',
- calendar = 'Gregorian',
- first_analysis = 2003, 1, 1, 0, 0, 0,
- last_analysis = 2003, 1, 2, 0, 0, 0,
- forecast_length_days = 1,
- forecast_length_seconds = 0,
- verification_interval_seconds = 21600,
- temporal_coverage_percent = 100.0,
- lonlim1 = missing_r8,
- lonlim2 = missing_r8,
- latlim1 = missing_r8,
- latlim2 = missing_r8,
- verbose = .false.,
- debug = .false.,
+ obs_sequence_list = 'obs_coverage_list.txt'
+ obs_sequence_name = ''
+ obs_of_interest = ''
+ textfile_out = 'obsdef_mask.txt'
+ netcdf_out = 'obsdef_mask.nc'
+ calendar = 'Gregorian'
+ first_analysis = 2003, 1, 1, 0, 0, 0
+ last_analysis = 2003, 1, 2, 0, 0, 0
+ forecast_length_days = 1
+ forecast_length_seconds = 0
+ verification_interval_seconds = 21600
+ temporal_coverage_percent = 100.0
+ lonlim1 = -888888.0
+ lonlim2 = -888888.0
+ latlim1 = -888888.0
+ latlim2 = -888888.0
+ verbose = .false.
+ debug = .false.
/
</pre>
</div>
@@ -210,7 +211,7 @@
<br />
<br />
-<P>Note that 'missing_r8' is not a valid number.
+<P>Note that -888888.0 is not a useful number.
To use the defaults delete these lines from the namelist,
or set them to 0.0, 360.0 and -90.0, 90.0.
</P><P>
@@ -320,18 +321,20 @@
<TD> integer </TD>
<TD>The number of seconds between each verification.
<UL style="list-style: none;">
-<LI> 1 h == 3600s</LI>
-<LI> 2 h == 7120s</LI>
-<LI> 3 h == 10800s</LI>
-<LI> 6 h == 21600s</LI>
+<LI> 1 h == 3600s</LI>
+<LI> 2 h == 7120s</LI>
+<LI> 3 h == 10800s</LI>
+<LI> 6 h == 21600s</LI>
+<LI> 12 h == 43200s</LI>
</UL>
</TD></TR>
<TR><TD> temporal_coverage_percent </TD>
<TD> real </TD>
- <TD>In the future, it may be possible to specify that you do not need an
-observation at <strong>every</strong> time. Presently, this is
-<strong>required</strong> to be 100%
+ <TD> While it is possible to specify that you do not need an
+observation at <strong>every</strong> time, it makes the most sense.
+This is not actually <strong>required</strong> to be 100% but 100%
+results in the most robust comparison.
</TD></TR>
<TR><TD> lonlim1 </TD>
@@ -379,24 +382,24 @@
<P>For example:
<pre>
&obs_seq_coverage_nml
- obs_sequence_list = 'obs_coverage_list.txt',
- obs_sequence_name = '',
+ obs_sequence_list = 'obs_coverage_list.txt'
+ obs_sequence_name = ''
obs_of_interest = 'METAR_U_10_METER_WIND',
- 'METAR_V_10_METER_WIND',
- textfile_out = 'obsdef_mask.txt',
- netcdf_out = 'obsdef_mask.nc',
- calendar = 'Gregorian',
- first_analysis = 2003, 1, 1, 0, 0, 0,
- last_analysis = 2003, 1, 2, 0, 0, 0,
- forecast_length_days = 1,
- forecast_length_seconds = 0,
- verification_interval_seconds = 21600,
- temporal_coverage_percent = 100.0,
- lonlim1 = 0.0,
- lonlim2 = 360.0,
- latlim1 = -90.0,
- latlim2 = 90.0,
- verbose = .false.,
+ 'METAR_V_10_METER_WIND'
+ textfile_out = 'obsdef_mask.txt'
+ netcdf_out = 'obsdef_mask.nc'
+ calendar = 'Gregorian'
+ first_analysis = 2003, 1, 1, 0, 0, 0
+ last_analysis = 2003, 1, 2, 0, 0, 0
+ forecast_length_days = 1
+ forecast_length_seconds = 0
+ verification_interval_seconds = 21600
+ temporal_coverage_percent = 100.0
+ lonlim1 = 0.0
+ lonlim2 = 360.0
+ latlim1 = -90.0
+ latlim2 = 90.0
+ verbose = .false.
/
</pre>
@@ -464,9 +467,18 @@
<br />
<br />
There is no requirement on the reporting time/frequence of the
-candidate stations. Once the verification times have been defined,
+candidate voxels. Once the verification times have been defined,
the observation <strong>closest in time</strong> to the verification
-time is selected, the others are ignored.
+time is selected, the others are ignored. Only observations within half
+the verification interval are eligible to be considered "close".
+<br />
+<br />
+<strong>A word about vertical levels.</strong>
+If the desired observation type has UNDEFINED or SURFACE for the vertical coordinate system,
+there is no concern about trying to match the vertical. If the desired observation types use
+PRESSURE; the following 14 levels are used as the standard levels: 1000, 925, 850, 700, 500,
+400, 300, 250, 200, 150, 100, 70, 50, 10 (all hPa). <strong>No other vertical coordinate system
+is supported.</strong>
</P>
<a NAME="example48x6"></a>
@@ -479,7 +491,7 @@
file for a single forecast. All the required input observation sequence
filenames will be contained in a file referenced by the
<em class=option>obs_sequence_list</em> variable. We'll also restrict
-the observations to a specific rectangular (in Lat/Lon) region.
+the observations to a specific rectangular (in Lat/Lon) region at a particular level.
It is convenient to turn on the verbose option the first time
to get a feel for the logic. Here are the namelist settings if you want to
verify the METAR_U_10_METER_WIND and METAR_V_10_METER_WIND observations over the
@@ -488,24 +500,24 @@
<div class=routine>
<pre>
&obs_seq_coverage_nml
- obs_sequence_name = '',
- obs_sequence_list = <em class=input>'obs_file_list.txt'</em>,
+ obs_sequence_name = ''
+ obs_sequence_list = <em class=input>'obs_file_list.txt'</em>
obs_of_interest = <em class=input>'METAR_U_10_METER_WIND'</em>,
- <em class=input>'METAR_V_10_METER_WIND'</em>,
- textfile_out = 'obsdef_mask.txt',
- netcdf_out = 'obsdef_mask.nc',
- calendar = 'Gregorian',
- first_analysis = <em class=input> 2008, 6, 8, 18, 0, 0 </em>,
- last_analysis = <em class=input> 2008, 6, 8, 18, 0, 0 </em>,
- forecast_length_days = <em class=input>2</em>,
- forecast_length_seconds = <em class=input>0</em>,
- verification_interval_seconds = <em class=input>21600</em>,
- temporal_coverage_percent = 100.0,
- lonlim1 = <em class=input> 0.0</em>,
- lonlim2 = <em class=input> 360.0</em>,
- latlim1 = <em class=input> -90.0</em>,
- latlim2 = <em class=input> 90.0</em>,
- verbose = <em class=input>.true.</em>,
+ <em class=input>'METAR_V_10_METER_WIND'</em>
+ textfile_out = 'obsdef_mask.txt'
+ netcdf_out = 'obsdef_mask.nc'
+ calendar = 'Gregorian'
+ first_analysis = <em class=input> 2008, 6, 8, 18, 0, 0 </em>
+ last_analysis = <em class=input> 2008, 6, 8, 18, 0, 0 </em>
+ forecast_length_days = <em class=input>2</em>
+ forecast_length_seconds = <em class=input>0</em>
+ verification_interval_seconds = <em class=input>21600</em>
+ temporal_coverage_percent = 100.0
+ lonlim1 = <em class=input> 0.0</em>
+ lonlim2 = <em class=input> 360.0</em>
+ latlim1 = <em class=input> -90.0</em>
+ latlim2 = <em class=input> 90.0</em>
+ verbose = <em class=input>.true.</em>
/
</pre>
</div>
@@ -600,7 +612,7 @@
Processing obs 80000 of 84691
obs_seq_coverage doneDONEdoneDONE does not exist. Finishing up.
- There were 442 stations matching the input criterion.
+ There were 442 voxels matching the input criterion.
...
</pre>
</div>
@@ -615,7 +627,7 @@
contains the metadata for all candidate locations as well as a lot of
information about the desired verification times. It is possible to explore
<em class=file>obsdef_mask.nc</em> to review the selection criteria to
-include observations/"stations" that do not perfectly match the original selection
+include observations/"voxels" that do not perfectly match the original selection
criteria.
<br />
<br />
@@ -627,18 +639,19 @@
<pre>
netcdf obsdef_mask {
dimensions:
- stations = UNLIMITED ; // (512 currently)
+ voxel = UNLIMITED ; // (512 currently)
time = 9 ;
analysisT = 1 ;
forecast_lead = 9 ;
+ nlevels = 14 ;
linelen = 129 ;
nlines = 446 ;
stringlength = 32 ;
location = 3 ;
variables:
- int stations(stations) ;
- stations:long_name = "desired station flag" ;
- stations:description = "1 == good station" ;
+ int voxel(voxel) ;
+ voxel:long_name = "desired voxel flag" ;
+ voxel:description = "1 == good voxel" ;
double time(time) ;
time:long_name = "verification time" ;
time:units = "days since 1601-1-1" ;
@@ -656,17 +669,20 @@
verification_times:calendar = "GREGORIAN" ;
verification_times:rows = "each forecast" ;
verification_times:cols = "each verification time" ;
+ float mandatory_level(nlevels) ;
+ mandatory_level:long_name = "mandatory pressure levels" ;
+ mandatory_level:units = "Pa" ;
char namelist(nlines, linelen) ;
namelist:long_name = "input.nml contents" ;
- char obs_type(stations, stringlength) ;
- obs_type:long_name = "observation type string at this station" ;
- double location(stations, location) ;
+ char obs_type(voxel, stringlength) ;
+ obs_type:long_name = "observation type string at this voxel" ;
+ double location(voxel, location) ;
location:description = "location coordinates" ;
location:location_type = "loc3Dsphere" ;
location:long_name = "threed sphere locations: lon, lat, vertical" ;
location:storage_order = "Lon Lat Vertical" ;
location:units = "degrees degrees which_vert" ;
- int which_vert(stations) ;
+ int which_vert(voxel) ;
which_vert:long_name = "vertical coordinate system code" ;
which_vert:VERTISUNDEF = -2 ;
which_vert:VERTISSURFACE = -1 ;
@@ -674,17 +690,17 @@
which_vert:VERTISPRESSURE = 2 ;
which_vert:VERTISHEIGHT = 3 ;
which_vert:VERTISSCALEHEIGHT = 4 ;
- int ntimes(stations) ;
- ntimes:long_name = "number of observation times at this station" ;
- double first_time(stations) ;
- first_time:long_name = "first valid observation time at this station" ;
+ int ntimes(voxel) ;
+ ntimes:long_name = "number of observation times at this voxel" ;
+ double first_time(voxel) ;
+ first_time:long_name = "first valid observation time at this voxel" ;
first_time:units = "days since 1601-1-1" ;
first_time:calendar = "GREGORIAN" ;
- double last_time(stations) ;
- last_time:long_name = "last valid observation time at this station" ;
+ double last_time(voxel) ;
+ last_time:long_name = "last valid observation time at this voxel" ;
last_time:units = "days since 1601-1-1" ;
last_time:calendar = "GREGORIAN" ;
- double ReportTime(stations, time) ;
+ double ReportTime(voxel, time) ;
ReportTime:long_name = "time of observation" ;
ReportTime:units = "days since 1601-1-1" ;
ReportTime:calendar = "GREGORIAN" ;
@@ -713,22 +729,22 @@
</pre>
</div>
<P>
-The first thing to note is that there are more stations (512) than reported
-during the run-time output (442). Typically, there will be many more stations
+The first thing to note is that there are more voxels (512) than reported
+during the run-time output (442). Typically, there will be many more voxels
in the netCDF file than will meet the selection criteria - but this is just
-an example. Some of the stations in the netCDF file do not meet the selection
+an example. Some of the voxels in the netCDF file do not meet the selection
criteria - meaning they do not have observations at all 9 required times.
Furthermore, there are 512 locations for ALL of the desired observation types.
In keeping with the DART philosophy of scalar observations, each observation
-type gets a separate station. There are <strong>not</strong> 512
+type gets a separate voxel. There are <strong>not</strong> 512
METAR_U_10_METER_WIND observations and 512 METAR_V_10_METER_WIND observations.
There are N METAR_U_10_METER_WIND observations and M METAR_V_10_METER_WIND
observations where N+M = 512. And only 442 of them have observations at all the
times required for the verification. Dump the <em class=option>obs_type</em> variable
-to see what station has what observation type.
+to see what voxel has what observation type.
<br />
<br />
-The <em class=option>stations</em> variable is fundamentally a flag that
+The <em class=option>voxel</em> variable is fundamentally a flag that
indicates if the station has all of the desired verification times. Combine
that information with the <em class=option>obs_type</em> and
<em class=option>location</em> to determine where your verifications of
@@ -837,7 +853,7 @@
<div class="top">[<a href="#">top</a>]</div><hr />
<H2>FUTURE PLANS</H2>
<P>
-Relax the restriction requiring 100.0% temporal coverage.<br />
+<em class=removed>Relax the restriction requiring 100.0% temporal coverage.</em><br />
Sensibly require that we only require against observations that DART can compute.
i.e. the prior forward operator must complete successfully - almost all cases
of the operator failing are extrapolation issues; the observation is outside
Modified: DART/trunk/obs_sequence/obs_seq_coverage.nml
===================================================================
--- DART/trunk/obs_sequence/obs_seq_coverage.nml 2014-07-28 20:29:08 UTC (rev 7091)
+++ DART/trunk/obs_sequence/obs_seq_coverage.nml 2014-07-29 22:35:31 UTC (rev 7092)
@@ -11,23 +11,23 @@
# 21600 secs == 6 hours
&obs_seq_coverage_nml
- obs_sequence_list = 'obs_coverage_list.txt',
- obs_sequence_name = '',
- obs_of_interest = '',
- textfile_out = 'obsdef_mask.txt',
- netcdf_out = 'obsdef_mask.nc',
- calendar = 'Gregorian',
- first_analysis = 2003, 1, 1, 0, 0, 0,
- last_analysis = 2003, 1, 2, 0, 0, 0,
- forecast_length_days = 1,
- forecast_length_seconds = 0,
- verification_interval_seconds = 21600,
- temporal_coverage_percent = 100.0,
- lonlim1 = missing_r8,
- lonlim2 = missing_r8,
- latlim1 = missing_r8,
- latlim2 = missing_r8,
- verbose = .false.,
- debug = .false.,
+ obs_sequence_list = 'obs_coverage_list.txt'
+ obs_sequence_name = ''
+ obs_of_interest = ''
+ textfile_out = 'obsdef_mask.txt'
+ netcdf_out = 'obsdef_mask.nc'
+ calendar = 'Gregorian'
+ first_analysis = 2003, 1, 1, 0, 0, 0
+ last_analysis = 2003, 1, 2, 0, 0, 0
+ forecast_length_days = 1
+ forecast_length_seconds = 0
+ verification_interval_seconds = 21600
+ temporal_coverage_percent = 100.0
+ lonlim1 = -888888.0
+ lonlim2 = -888888.0
+ latlim1 = -888888.0
+ latlim2 = -888888.0
+ verbose = .false.
+ debug = .false.
/
More information about the Dart-dev
mailing list