[Dart-dev] [7103] DART/trunk/obs_sequence: Extended to handle mandatory vertical levels as defined in the obsdef_mask .nc
nancy at ucar.edu
nancy at ucar.edu
Mon Aug 4 14:43:54 MDT 2014
Revision: 7103
Author: thoar
Date: 2014-08-04 14:43:53 -0600 (Mon, 04 Aug 2014)
Log Message:
-----------
Extended to handle mandatory vertical levels as defined in the obsdef_mask.nc
The file input mechanism now uses the utilities_mod:set_filename_list()
routine to unify input filename implementation.
There are some 'verbose' timestamps to see if writing each observation
is too slow for production runs.
Modified Paths:
--------------
DART/trunk/obs_sequence/obs_seq_verify.f90
DART/trunk/obs_sequence/obs_seq_verify.html
DART/trunk/obs_sequence/obs_seq_verify.nml
-------------- next part --------------
Modified: DART/trunk/obs_sequence/obs_seq_verify.f90
===================================================================
--- DART/trunk/obs_sequence/obs_seq_verify.f90 2014-08-04 20:34:37 UTC (rev 7102)
+++ DART/trunk/obs_sequence/obs_seq_verify.f90 2014-08-04 20:43:53 UTC (rev 7103)
@@ -10,19 +10,19 @@
!
! This program creates a netCDF file suitable for forecast evaluation.
!
-! (copy, station, level, ensemble, analysisT, fcstlead)
-! | | | | | |
-! | | | | | +- forecast length : 0,3,6,9,12,...
-! | | | | |
-! | | | | +---------- analysis time/date
-! | | | |
-! | | | +--------------------- ensemble member index
-! | | |
-! | | +----------------------------- vertical level index
-! | |
-! | +------------------------------------- (horizontal) station index
-! |
-! +--------------------------------------------- obs value, prior, obs_err.
+! (analysisT, station, level, copy, ensemble, forecast_lead)
+! | | | | | |
+! | | | | | +- forecast length : 0,3,6,9,12,...
+! | | | | |
+! | | | | +----------- ensemble member index
+! | | | |
+! | | | +------------------ obs value, prior, obs_err
+! | | |
+! | | +------------------------- vertical level index
+! | |
+! | +--------------------------------- (horizontal) station index
+! |
+! +------------------------------------------- analysis time/date
!
! I think the logic of the program should be as follows:
! 1) The list of 'stations' is read from a netCDF file -- the product of obs_seq_coverage.f90
@@ -34,17 +34,11 @@
! 5) On to the next obs_seq.fcst file ... (step 3)
! 6) wrap up ...
!
-! Soyoung's wish list:
-! Now I'm done running filter for 24-hr forecast with 3-hrly observations as an
+! Soyoung's (original) wish list - which has been subsequently modified:
+! "Now I'm done running filter for 24-hr forecast with 3-hrly observations as an
! evaluation mode only, and am ready to hand obs_seq.final over to you for the final
! conversion process.
!
-! be1005en.ucar.edu:/ptmp/syha/wrfruns/ENS_FCST/Verify>
-! -rw-r--r-- 1 syha ncar 1006513852 Nov 19 12:24 Prior_Diag.nc
-! -rw-r--r-- 1 syha ncar 1006513856 Nov 19 12:24 Posterior_Diag.nc
-! -rw-r--r-- 1 syha ncar 111484368 Nov 19 12:24 prior_inflate_restart
-! -rw-r--r-- 1 syha ncar 126848308 Nov 19 12:24 obs_seq.final
-!
! Ideally, in obs_seq_fcst.nc (if I can name it on my own), I would like to have a
! data structure of (copy, station, level, ensemble, date, time) for each variable
! and each obs type, where copy is (observation value, prior observation value
@@ -55,40 +49,46 @@
! same length of forecast lead times. This dimension is the most critical one in the
! statistical verification since this dimension number is the actual sample size to
! tell if we have enough samples to make the verification result statistically
-! significant or not.
+! significant or not."
!
!-----------------------------------------------------------------------
use types_mod, only : r4, r8, digits12, MISSING_R8, 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_obs_values, &
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, get_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, VERTISUNDEF, &
+ vert_is_undef, vert_is_surface, vert_is_pressure
+
use time_manager_mod, only : time_type, set_date, set_time, get_time, print_time, &
set_time_missing, print_date, set_calendar_type, &
operator(>), operator(<), operator(==), &
operator(<=), operator(-), operator(+), operator(/=)
-use utilities_mod, only : get_unit, close_file, register_module, &
+
+use utilities_mod, only : get_unit, close_file, register_module, timestamp, &
file_exist, error_handler, E_ERR, E_WARN, E_MSG, &
initialize_utilities, finalize_utilities, nmlfileunit, &
find_namelist_in_file, check_namelist_read, nc_check, &
- next_file, get_next_filename, find_textfile_dims, &
+ next_file, set_filename_list, find_textfile_dims, &
file_to_text, do_nml_file, do_nml_term
-use sort_mod, only : sort
use typeSizes
use netcdf
@@ -107,29 +107,62 @@
! written to the appropriate slots in the netCDF file.
!---------------------------------------------------------------------
-type station
+type voxel_type
+ integer :: good
integer :: obs_type
- type(location_type) :: location
type(time_type) :: first_time
type(time_type) :: last_time
- integer :: ntimes
- type(time_type), pointer :: times(:)
- integer :: analysisT
- integer, pointer, dimension( :) :: orgqc
- integer, pointer, dimension( :) :: dartqc
- real(r8), pointer, dimension( :) :: observation ! (fcstlead)
- real(r8), pointer, dimension( :) :: obserror ! (fcstlead)
- real(r8), pointer, dimension(:,:) :: forecast ! (enssize, fcstlead) priors
-end type station
+ type(time_type), pointer :: times(:) ! actual obs time closest to nominal
+ type(location_type) :: location
+ real(r8), dimension(3) :: LonLatLev ! easy access to location components
+ integer :: station_id ! index of unique horizontal location
+ integer :: levelindex ! index of mandatory level
+end type voxel_type
-! (copy, station, level, ensemble, analysisT, fcstlead)
+type station_type
+ integer :: obs_type
+ type(location_type) :: location
+end type station_type
-logical, allocatable, dimension(:) :: DesiredStations
-type(station), allocatable, dimension(:) :: stations
-integer :: ensemble_size ! the # of ensemble members in the obs_seq file
-integer :: num_stations ! This is the current number of unique locations
+type nc6Dvar
+ character(len=NF90_MAX_NAME) :: AnalysisDimName = 'analysisT'
+ integer :: AnalysisDimID
+ integer :: AnalysisDimLen
+ character(len=NF90_MAX_NAME) :: CopyDimName = 'copy'
+ integer :: CopyDimID
+ integer :: CopyDimLen
+ character(len=NF90_MAX_NAME) :: StationsDimName = 'station'
+ integer :: StationsDimID
+ integer :: StationsDimLen
+ character(len=NF90_MAX_NAME) :: LevelsDimName = 'level'
+ integer :: LevelsDimID
+ integer :: LevelsDimLen
+ character(len=NF90_MAX_NAME) :: EnsembleDimName = 'ensemble'
+ integer :: EnsembleDimID
+ integer :: EnsembleDimLen
+ character(len=NF90_MAX_NAME) :: ForecastDimName = 'forecast_lead'
+ integer :: ForecastDimID
+ integer :: ForecastDimLen
+end type nc6Dvar
+
+logical, allocatable, dimension(:) :: DesiredStations
+type(station_type), allocatable, dimension(:) :: station
+type(voxel_type), allocatable, dimension(:) :: voxel
+type(nc6Dvar) :: ncmeta
+
+integer :: ensemble_size ! the # of ensemble members in the obs_seq files
+integer :: num_voxels ! The number of unique locations
+integer :: num_stations ! The number of unique horizontal locations
+integer :: num_levels ! number of mandatory levels
integer :: num_verif_times ! number of all possible verification times.
+integer :: num_forecasts ! number of forecasts
+integer, parameter :: NCOPIES = 3 ! obs value, prior, obs error variance
+character(len=metadatalength), dimension(NCOPIES) :: copy_metadata = &
+ (/ 'observation ', &
+ 'forecast ', &
+ 'observation error variance' /)
+
!---------------------------------------------------------------------
! variables associated with the observation
!---------------------------------------------------------------------
@@ -139,9 +172,6 @@
type(obs_def_type) :: obs_def
type(location_type) :: obs_loc
-character(len = 129) :: obs_seq_in_file_name
-character(len = 129), allocatable, dimension(:) :: obs_seq_filenames
-
integer :: flavor
integer :: num_copies, num_qc, num_obs, max_num_obs, obs_seq_file_id
integer :: obtype_integer
@@ -155,23 +185,29 @@
logical :: pre_I_format
logical :: last_ob_flag
+integer, parameter :: MAX_NUM_INPUT_FILES = 500
+integer :: num_input_files
+character(len=256) :: obs_seq_in_file_name
+
!-----------------------------------------------------------------------
! Namelist with (some scalar) default values
!-----------------------------------------------------------------------
-character(len = 129) :: obs_sequence_name = 'obs_seq.final'
-character(len = 129) :: obs_sequence_list = ''
-character(len = 129) :: station_template = 'obsdef_mask.nc'
-character(len = 129) :: netcdf_out = 'obs_seq_fcst.nc'
-character(len = 129) :: calendar = 'gregorian'
-character(len = obstypelength) :: obtype_string
+character(len=256) :: obs_sequences(MAX_NUM_INPUT_FILES) = ''
+character(len=256) :: obs_sequence_list = ''
+character(len=256) :: input_template = 'obsdef_mask.nc'
+character(len=256) :: netcdf_out = 'forecast.nc'
+character(len=129) :: calendar = 'gregorian'
+character(len=obstypelength) :: obtype_string
-logical :: verbose = .false.
-logical :: debug = .false. ! undocumented ... on purpose
+integer :: print_every = 10000
+logical :: verbose = .true.
+logical :: debug = .false.
-namelist /obs_seq_verify_nml/ obs_sequence_name, obs_sequence_list, &
- station_template, netcdf_out, &
- calendar, obtype_string, verbose, debug
+namelist /obs_seq_verify_nml/ obs_sequences, obs_sequence_list, &
+ input_template, netcdf_out, &
+ obtype_string, calendar, &
+ print_every, verbose, debug
!-----------------------------------------------------------------------
! Quantities of interest
@@ -180,30 +216,34 @@
integer :: qc_index ! copy index of the original qc value
integer :: dart_qc_index ! copy index of the DART qc value
integer :: obs_copy_index ! copy index of the observation
-integer :: station_id
-integer :: fcst_T_index ! verification time/forecast lead time index
+integer :: voxel_id
+integer :: fcst_lead_index ! forecast lead time index
integer :: AnalysisIndex ! index of the forecast experiment
+type(time_type) :: analysisT ! valid time of analysis at start of forecast
character(len=metadatalength), dimension(:), allocatable :: module_obs_copy_names
-integer, dimension(:), allocatable :: copy_indices
+integer, dimension(:), allocatable :: prior_copy_indices
integer, dimension(:), allocatable :: forecast_leads
-real(digits12), dimension(:,:), allocatable :: ExperimentTimes
+real(r8), dimension(:), allocatable :: mandatory_level
type(time_type), dimension(:,:), allocatable :: VerifyTimes
+real(r8) :: obs_error_variance
!-----------------------------------------------------------------------
! General purpose variables
!-----------------------------------------------------------------------
-integer :: ifile, nread, ngood
+integer :: ifile, iobs, ngood
integer :: i, io, ncunit
+integer, dimension(6) :: mystart, mycount
type(time_type) :: obs_time
-type(time_type) :: analysisT ! valid time of analysis at start of forecast
-character(len = 129) :: ncName, string1, string2
+character(len=256) :: ncName
+character(len=512) :: 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)
+real(r8), parameter :: OnePa = 1.0_r8
!=======================================================================
! Get the party started
@@ -229,6 +269,15 @@
if (do_nml_file()) write(nmlfileunit, nml=obs_seq_verify_nml)
if (do_nml_term()) write( * , nml=obs_seq_verify_nml)
+num_input_files = set_filename_list(obs_sequences, obs_sequence_list,'obs_seq_verify')
+
+if (debug) then
+ write(*,*)'There are ',num_input_files,' input observation sequence files.'
+ do ifile = 1,num_input_files
+ write(*,*)'file ',ifile,' is ['//trim(obs_sequences(ifile))//']'
+ enddo
+endif
+
string1 = adjustl(obtype_string)
string2 = adjustl(netcdf_out)
obtype_string = trim(string1)
@@ -237,26 +286,21 @@
if (obtype_integer < 1) then
write(string1,*)'obtype_string ',trim(obtype_string),' is unknown. change input.nml'
- call error_handler(E_ERR,'obs_seq_verify',string1,source,revision,revdate)
+ call error_handler(E_ERR,'obs_seq_verify:',string1,source,revision,revdate)
endif
call set_calendar_type(calendar)
-! Check the user input for sanity
-if ((obs_sequence_name /= '') .and. (obs_sequence_list /= '')) then
- write(string1,*)'specify "obs_sequence_name" or "obs_sequence_list"'
- write(string2,*)'set other to an empty string ... i.e. ""'
- call error_handler(E_ERR, 'obs_seq_verify', string1, &
- source, revision, revdate, text2=string2)
-endif
-
-! Determine the number of stations from the input netcdf file,
+! Determine the number of voxels from the input netcdf file.
! Initialize the output netcdf file. Must know the ensemble size.
+! Determine which voxels share the same horizontal location.
+! Each voxel already has a vertical level index associated with it.
-ensemble_size = find_ensemble_size() ! from the obs_seq copy metadata
-num_stations = fill_stations( station_template, ensemble_size )
+call find_ensemble_size() ! from the obs_seq copy metadata of the first file
+num_voxels = fill_voxels( input_template )
+num_stations = relate_voxels_to_stations()
-allocate(copy_indices(ensemble_size), module_obs_copy_names(ensemble_size))
+allocate(prior_copy_indices(ensemble_size), module_obs_copy_names(ensemble_size))
ncName = trim(obtype_string)//'_'//trim(netcdf_out)
@@ -266,74 +310,63 @@
! Prepare the variables
!----------------------------------------------------------------------
-allocate(obs_seq_filenames(1000))
-obs_seq_filenames = 'null'
-
-ObsFileLoop : do ifile=1, size(obs_seq_filenames)
+ObsFileLoop : do ifile=1,num_input_files
!-----------------------------------------------------------------------
+! By construction -
+! each file contains all the observations for an entire forecast.
+! The file NAME specifies the analysis time from which the forecast is started.
- ! Because of the ability to 'cycle' the ObsFileLoop, we need to
- ! destroy and deallocate at the top of the loop.
+ ! Because of the ability to 'cycle' the ObsFileLoop, we need to
+ ! destroy and deallocate at the top of the loop.
- call destroy_obs(obs1)
- call destroy_obs(obs2)
- call destroy_obs_sequence(seq)
+ if (ifile /= 1) then
+ call destroy_obs(obs1)
+ call destroy_obs(obs2)
+ call destroy_obs_sequence(seq)
+ endif
if (allocated( copy_values)) deallocate( copy_values)
if (allocated( qc_values)) deallocate( qc_values)
if (allocated( qc_copy_names)) deallocate( qc_copy_names)
if (allocated(obs_copy_names)) deallocate(obs_copy_names)
- ! Determine the next input filename ...
+ ! Determine the next input filename and check for existence.
- if (obs_sequence_list == '') then
- obs_seq_in_file_name = trim(next_file(obs_sequence_name,ifile))
- else
- obs_seq_in_file_name = trim(get_next_filename(obs_sequence_list,ifile))
- if (obs_seq_in_file_name == '') exit ObsFileLoop
- endif
+ obs_seq_in_file_name = obs_sequences(ifile)
- if ( file_exist(trim(obs_seq_in_file_name)) ) then
+ if ( file_exist(obs_seq_in_file_name) ) then
write(*,*) ! whitespace
write(*,*) ! whitespace
write(string1,*)'opening ', trim(obs_seq_in_file_name)
- call error_handler(E_MSG,'obs_seq_verify:',string1,source,revision,revdate)
+ call error_handler(E_MSG,'obs_seq_verify:',string1)
else
- write(string1,*)trim(obs_seq_in_file_name),&
- ' does not exist. Finishing up.'
- call error_handler(E_MSG,'obs_seq_verify:',string1,source,revision,revdate)
+ write(string1,*)'['//trim(obs_seq_in_file_name)//'] does not exist.'
+ call error_handler(E_ERR,'obs_seq_verify:',string1,source,revision,revdate)
exit ObsFileLoop
endif
- ! Read in information about observation sequence so we can allocate
- ! observations. We need info about how many copies, qc values, etc.
-
- obs_seq_filenames(ifile) = trim(obs_seq_in_file_name)
-
! Determine the analysis time from the file name
call find_analysis_time(obs_seq_in_file_name, AnalysisIndex, analysisT)
- if (verbose) then
- write(string1,*)'analysis ',AnalysisIndex,' date is '
- call print_date(analysisT,string1)
- endif
+ ! Read in information about observation sequence so we can allocate
+ ! observations. We need info about how many copies, qc values, etc.
call read_obs_seq_header(obs_seq_in_file_name, &
num_copies, num_qc, num_obs, max_num_obs, &
obs_seq_file_id, obs_seq_read_format, pre_I_format, &
close_the_file = .true.)
+ if ((num_qc <= 0) .or. (num_copies <= 0)) then
+ write(string1,*)'need at least 1 qc and 1 observation copy'
+ call error_handler(E_ERR,'obs_seq_verify:',string1,source,revision,revdate)
+ endif
+
! Initialize some (individual) observation variables
call init_obs(obs1, num_copies, num_qc) ! First obs in sequence
call init_obs(obs2, num_copies, num_qc)
- if ((num_qc <= 0) .or. (num_copies <=0)) then
- write(string1,*)'need at least 1 qc and 1 observation copy'
- call error_handler(E_ERR,'obs_seq_verify',string1,source,revision,revdate)
- endif
-
allocate( copy_values(num_copies), qc_values(num_qc))
allocate( obs_copy_names(num_copies), qc_copy_names(num_qc))
@@ -347,14 +380,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
! * Find the N copy indices we want
@@ -366,69 +391,76 @@
call read_obs_seq(obs_seq_in_file_name, 0, 0, 0, seq)
do i=1, num_copies
- string1 = trim(get_copy_meta_data(seq,i))
- obs_copy_names(i) = adjustl(string1)
+ obs_copy_names(i) = get_copy_meta_data(seq,i)
enddo
do i=1, num_qc
- string1 = trim(get_qc_meta_data(seq,i))
- qc_copy_names(i) = adjustl(string1)
+ qc_copy_names(i) = get_qc_meta_data(seq,i)
enddo
! Determine which qc copy is original qc and dart qc
- call find_our_copies(seq, obs_copy_index, copy_indices, qc_index, dart_qc_index )
+ call find_our_copies(seq, obs_copy_index, prior_copy_indices)
! The first trip through sets the module_obs_copy_names so we can
! be sure we are stuffing compatible objects into the same slots
if ( ifile == 1 ) then
- module_obs_copy_names = obs_copy_names(copy_indices(1:ensemble_size))
+ module_obs_copy_names = obs_copy_names(prior_copy_indices(1:ensemble_size))
else
! Check to make sure the ensemble members are in the expected copies
do i = 1,ensemble_size
- if ( obs_copy_names(copy_indices(i)) /= module_obs_copy_names(i) ) then
+ if ( obs_copy_names(prior_copy_indices(i)) /= module_obs_copy_names(i) ) then
write(string1,'(''module has '',A)') trim(module_obs_copy_names(i))
write(string2,'(A,'' has '',A)') trim(obs_seq_in_file_name), &
- trim(obs_copy_names(copy_indices(i)))
- call error_handler(E_ERR,'obs_seq_verify', &
- 'mismatch in observation copies',source,revision,revdate,text2=string1,text3=string2)
+ trim(obs_copy_names(prior_copy_indices(i)))
+ call error_handler(E_ERR,'obs_seq_verify:', &
+ 'mismatch in observation copies', source, &
+ revision, revdate, text2=string1,text3=string2)
endif
enddo
endif
+ ngood = 0
+
!--------------------------------------------------------------------
- ! Read the first observation in the sequence
+ ObservationLoop : do iobs = 1,num_obs
!--------------------------------------------------------------------
- ngood = 0
-
- if ( .not. get_first_obs(seq, obs1) ) &
- call error_handler(E_ERR,'obs_seq_verify', &
+ if (iobs == 1) then
+ if ( .not. get_first_obs(seq, obs1) ) &
+ call error_handler(E_ERR,'obs_seq_verify:', &
'No first observation in sequence.', &
source,revision,revdate)
+ else
+ call get_next_obs(seq, obs1, obs2, last_ob_flag)
+ obs1 = obs2
+ endif
- !--------------------------------------------------------------------
- ObservationLoop : do nread = 1,num_obs
- !--------------------------------------------------------------------
+ if ( verbose .and. (mod(iobs,print_every) == 0) ) &
+ write(*,*)'Processing obs ',iobs,' of ',num_obs
- if ( verbose .and. (mod(nread,10000) == 0) ) &
- write(*,*)'Processing obs ',nread,' of ',num_obs
-
call get_obs_def(obs1, obs_def)
flavor = get_obs_kind( obs_def)
obs_time = get_obs_def_time( obs_def)
obs_loc = get_obs_def_location(obs_def)
+ call get_qc(obs1, qc_values)
- ! if (debug .and. any(nread == (/1, 295, 1908, 6265/) )) then
- if (debug .and. (nread == 1)) then
+ ! if (debug .and. any(iobs == (/1, 295, 1908, 6265/) )) then
+ if (debug .and. (iobs == 1)) then
+ write(*,*)'looking for observation types of ',obtype_integer, &
+ get_obs_kind_name(obtype_integer)
+
+ write(*,*)'First observation "type" happens to be ',flavor, &
+ get_obs_kind_name(flavor)
call print_time(obs_time,'First observation time')
call print_date(obs_time,'First observation date')
- call get_qc( obs1, qc_values)
+ call write_location(6,obs_loc,fform='ascii')
+ call write_location(6,obs_loc,fform='ascii',charstring=string1)
+ write(*,*)trim(string1)
+
+ write(*,*)' qc_index and value is ', qc_index, qc_values( qc_index)
+ write(*,*)'dart_qc_index and value is ',dart_qc_index, qc_values(dart_qc_index)
+
call get_obs_values(obs1, copy_values)
- write(*,*)'looking for observation kinds of ',obtype_integer, &
- get_obs_kind_name(obtype_integer)
-
- write(*,*)'First observation "kind" is ',flavor, &
- get_obs_kind_name(flavor)
write(*,*)'observation values are:'
do i = 1,size(copy_values)
write(*,*)i,copy_values(i)
@@ -438,15 +470,6 @@
write(*,*)i,qc_values(i)
enddo
- call write_location(6,obs_loc,fform='ascii')
- call write_location(6,obs_loc,fform='ascii',charstring=string1)
- write(*,*)trim(string1)
-
- if (qc_index > 0) &
- write(*,*)' qc_index and value is ', qc_index, qc_values( qc_index)
- if (dart_qc_index > 0) &
- write(*,*)'dart_qc_index and value is ',dart_qc_index, qc_values(dart_qc_index)
-
endif
! 0) is it the right type of observation
@@ -454,58 +477,47 @@
! 2) Is it one of the times of interest?
! 3) stuff it into the appropriate station structure
- if ( obtype_integer /= flavor ) goto 100
+ ! FIXME ... What do we do with DART QCs of 7
+ if ( flavor /= obtype_integer) cycle ObservationLoop
+ if (qc_values(dart_qc_index) >= 4) cycle ObservationLoop
- station_id = find_station_location(flavor, obs_loc, stations)
- if ( station_id < 1 ) goto 100
+ call match_obs_to_voxel(iobs, flavor, obs_loc, obs_time, voxel_id, fcst_lead_index)
+ if ( fcst_lead_index < 1 ) cycle ObservationLoop
- if ( .not. time_wanted( obs_time, station_id, fcst_T_index)) goto 100
-
- if (debug) write(*,*)'obs ',nread,' is station ',station_id,' at fcst index ',fcst_T_index
-
- call get_qc( obs1, qc_values)
call get_obs_values(obs1, copy_values)
+ obs_error_variance = get_obs_def_error_variance(obs_def)
- if ( qc_index > 0) stations(station_id)%orgqc( fcst_T_index) = qc_values( qc_index)
- if (dart_qc_index > 0) stations(station_id)%dartqc(fcst_T_index) = qc_values(dart_qc_index)
-
- stations(station_id)%obserror( fcst_T_index) = get_obs_def_error_variance(obs_def)
- stations(station_id)%observation(fcst_T_index) = copy_values(obs_copy_index)
- stations(station_id)%forecast(:, fcst_T_index) = copy_values(copy_indices)
-
! if (debug) then ! GIGANTIC OUTPUT - BE CAREFUL
- if (debug .and. (station_id == 1)) then ! LESS GIGANTIC OUTPUT
- do i = 1,size(copy_indices)
- write(*,'(''ob'',3(1x,i6),4(1x,f14.7))') &
- nread, fcst_T_index, copy_indices(i), &
- stations(station_id)%obserror( fcst_T_index),&
- stations(station_id)%observation(fcst_T_index),&
- stations(station_id)%forecast(i, fcst_T_index),&
- copy_values(copy_indices(i))
+ if (debug .and. (voxel(voxel_id)%station_id == 1)) then ! LESS GIGANTIC OUTPUT
+ do i = 1,size(prior_copy_indices)
+ write(*,'('' DEBUG ob '',i10,3(1x,i6),3(1x,f17.7))') &
+ iobs, i, fcst_lead_index, prior_copy_indices(i), &
+ obs_error_variance, &
+ copy_values(obs_copy_index), &
+ copy_values(prior_copy_indices(i))
enddo
write(*,*) ! a little whitespace
endif
- 100 continue
+ ! Determine where to stuff it into the output netCDF file.
+ ! YES ... one observation at a time.
+ ! FIXME allocate a potentially LARGE array instead of writing one-at-a-time
+ ! FIXME ... one for each Analysis Time
- call get_next_obs(seq, obs1, obs2, last_ob_flag)
- if (.not. last_ob_flag) obs1 = obs2
+ call determine_hyperslab_indices(AnalysisIndex, fcst_lead_index, &
+ voxel_id, mystart, mycount)
+ call WriteNetCDF(ncunit, trim(ncName), copy_values, obs_copy_index, &
+ prior_copy_indices, obs_error_variance, qc_values(qc_index), &
+ qc_values(dart_qc_index), mystart, mycount)
+
!--------------------------------------------------------------------
enddo ObservationLoop
!--------------------------------------------------------------------
- ! So by now, all the observations have been stuffed into the 'station'
- ! array - not that we know if the station array is complete ...
- ! Take what we have (i.e. for one analysis time) and push it into
- ! the output forecast netcdf file.
- call WriteNetCDF(ncunit, trim(ncName), ifile)
-
- call reset_stations()
-
enddo ObsFileLoop
-call CloseNetCDF(ncunit, trim(ncName))
+call CloseNetCDF(ncunit, ncName)
!-----------------------------------------------------------------------
! Really, really, done.
@@ -514,440 +526,460 @@
call destroy_obs(obs1)
call destroy_obs(obs2)
call destroy_obs_sequence(seq)
-call destroy_stations(stations)
+! call destroy_stations(station)
+call destroy_voxels()
if (allocated(qc_values)) deallocate(qc_values)
if (allocated(qc_copy_names)) deallocate(qc_copy_names)
if (allocated(copy_values)) deallocate(copy_values)
if (allocated(obs_copy_names)) deallocate(obs_copy_names)
-if (allocated(obs_seq_filenames)) deallocate(obs_seq_filenames)
if (allocated(DesiredStations)) deallocate(DesiredStations)
-call error_handler(E_MSG,'obs_seq_verify','Finished successfully.',source,revision,revdate)
+call error_handler(E_MSG,'obs_seq_verify:','Finished successfully.')
call finalize_utilities()
+
!======================================================================
CONTAINS
!======================================================================
-subroutine reset_stations()
-
-integer :: i
-
-do i = 1,size(stations)
- stations(i)%orgqc = MISSING_I
- stations(i)%dartqc = MISSING_I
- stations(i)%obserror = MISSING_R8
- stations(i)%observation = MISSING_R8
- stations(i)%forecast = MISSING_R8
-enddo
-
-end subroutine reset_stations
-
-
-
-function find_station_location(ObsType, ObsLocation, stationlist) result(station_id)
+function fill_voxels( filename ) result(nvoxels)
!----------------------------------------------------------------------------
-! 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(in) :: stationlist
-integer :: station_id
-
-integer :: i
-real(r8), dimension(3) :: obslocarray, stnlocarray
-real(r8) :: londiff, latdiff
-
-station_id = 0;
-
-FindLoop : do i = 1,num_stations
-
- obslocarray = get_location(ObsLocation) ! returns degrees
- stnlocarray = get_location(stationlist(i)%location)
-
- londiff = abs(obslocarray(1) - stnlocarray(1))
- latdiff = abs(obslocarray(2) - stnlocarray(2))
-
- if ( (londiff <= HALF_METER) .and. &
- (latdiff <= HALF_METER) .and. &
- (ObsType == stationlist(i)%obs_type) ) then
-
- station_id = i
- exit FindLoop
- endif
-
-enddo FindLoop
-
-end function find_station_location
-
-
-!============================================================================
-
-
-function fill_stations( filename, nensemble ) result(nstations)
-!----------------------------------------------------------------------------
! Check to make sure the file exists.
-! Read the 'station' variable (binary value for desired or not).
+! Read the 'voxel' variable (binary value for desired or not).
! count the number of non-zero entries ...
+!
+! Since the 'voxels' coming in are assumed vertically independent, but the output
+! structure 'station' has a vertical component ... we must determine which input
+! 'voxels' are actually part of a station.
!-----------------------------------------------------------------------
character(len=*), intent(in) :: filename
-integer, intent(in) :: nensemble
-integer :: nstations
+integer :: nvoxels
integer :: DimID, VarID
-integer :: nmax, nanalyses, nforecasts, locNdim, strlen
-integer :: ncid, i, j, istation, ndims, mylen
-integer :: my_type
+integer :: nanalyses, locNdim, strlen
+integer :: ncid, i, j, ivoxel, ndims, mylen
integer :: seconds, days
type(location_type) :: myloc
integer, dimension(nf90_max_var_dims) :: dimIDs
-! These temporarily mimic the contents of the netCDF file
+! These mimic the contents of the netCDF file
character(len=obstypelength), dimension(:), allocatable :: obs_types
-integer, dimension(:), allocatable :: allstations
-real(r8), dimension(:,:), allocatable :: locations
+integer, dimension(:), allocatable :: voxel_flag
integer, dimension(:), allocatable :: which_vert
-integer, dimension(:), allocatable :: ntimearray
+integer, dimension(:), allocatable :: voxel_level_index
real(digits12), dimension(:), allocatable :: first_time
real(digits12), dimension(:), allocatable :: last_time
real(digits12), dimension(:,:), allocatable :: ReportTimes
+real(digits12), dimension(:,:), allocatable :: ExperimentTimes
+real(r8), dimension(:,:), allocatable :: locations
+if (verbose) then
+ call timestamp(string1='fill_voxels: reading the voxel information at ',pos='brief')
+endif
+
!-----------------------------------------------------------------------
! Open the netCDF file
call nc_check(nf90_open(trim(filename), NF90_NOWRITE, &
- ncid = ncid), 'fill_stations', 'open '//trim(filename))
+ ncid = ncid), 'fill_voxels', 'open '//trim(filename))
!-----------------------------------------------------------------------
! Determine the maximum number of stations from the station dimension,
! the number of times, and the number of dimensions in the location
! used for integer array indicating yes/no - do we want the station
-call nc_check(nf90_inq_dimid(ncid, 'station', dimid=DimID), &
- 'fill_stations', 'inquire station dimid '//trim(filename))
-call nc_check(nf90_inquire_dimension(ncid, DimID, len=nmax), &
- 'fill_stations', 'inquire station nmax '//trim(filename))
+call nc_check(nf90_inq_dimid(ncid, 'voxel', dimid=DimID), &
+ 'fill_voxels', 'inquire voxel dimid '//trim(filename))
+call nc_check(nf90_inquire_dimension(ncid, DimID, len=nvoxels), &
+ 'fill_voxels', 'inquire voxel nvoxels '//trim(filename))
! used for array of all possible verification times needed
call nc_check(nf90_inq_dimid(ncid, 'time', dimid=DimID), &
- 'fill_stations', 'inquire time dimid '//trim(filename))
+ 'fill_voxels', 'inquire time dimid '//trim(filename))
call nc_check(nf90_inquire_dimension(ncid, DimID, len=num_verif_times), &
- 'fill_stations', 'inquire time len '//trim(filename))
+ 'fill_voxels', 'inquire time len '//trim(filename))
! used for array of the analysis times / 'zero-length' forecast times
call nc_check(nf90_inq_dimid(ncid, 'analysisT', dimid=DimID), &
- 'fill_stations', 'inquire analysisT dimid '//trim(filename))
+ 'fill_voxels', 'inquire analysisT dimid '//trim(filename))
call nc_check(nf90_inquire_dimension(ncid, DimID, len=nanalyses), &
- 'fill_stations', 'inquire analysisT len '//trim(filename))
+ 'fill_voxels', 'inquire analysisT len '//trim(filename))
! used for array of the verification times
call nc_check(nf90_inq_dimid(ncid, 'forecast_lead', dimid=DimID), &
- 'fill_stations', 'inquire forecast_lead dimid '//trim(filename))
-call nc_check(nf90_inquire_dimension(ncid, DimID, len=nforecasts), &
- 'fill_stations', 'inquire forecast_lead len '//trim(filename))
+ 'fill_voxels', 'inquire forecast_lead dimid '//trim(filename))
+call nc_check(nf90_inquire_dimension(ncid, DimID, len=num_forecasts), &
+ 'fill_voxels', 'inquire forecast_lead len '//trim(filename))
call nc_check(nf90_inq_dimid(ncid, 'stringlength', dimid=DimID), &
- 'fill_stations', 'inquire stringlength dimid '//trim(filename))
+ 'fill_voxels', 'inquire stringlength dimid '//trim(filename))
call nc_check(nf90_inquire_dimension(ncid, DimID, len=strlen), &
- 'fill_stations', 'inquire stringlength len '//trim(filename))
+ 'fill_voxels', 'inquire stringlength len '//trim(filename))
call nc_check(nf90_inq_dimid(ncid, 'location', dimid=DimID), &
- 'fill_stations', 'inquire location dimid '//trim(filename))
+ 'fill_voxels', 'inquire location dimid '//trim(filename))
call nc_check(nf90_inquire_dimension(ncid, DimID, len=locNdim), &
- 'fill_stations', 'inquire location len '//trim(filename))
+ 'fill_voxels', 'inquire location len '//trim(filename))
+call nc_check(nf90_inq_dimid(ncid, 'nlevels', dimid=DimID), &
+ 'fill_voxels', 'inquire nlevels dimid '//trim(filename))
+call nc_check(nf90_inquire_dimension(ncid, DimID, len=num_levels), &
+ 'fill_voxels', 'inquire nlevels len '//trim(filename))
+
!-----------------------------------------------------------------------
-! Read the "stations" 1D variable. Each non-zero entry indicates a
-! 'station' we wanted to keep when we ran obs_seq_coverage.f90.
+! Read the "voxel" 1D variable. Each non-zero entry indicates a
+! voxel we wanted to keep when we ran obs_seq_coverage.f90.
-call nc_check(nf90_inq_varid(ncid, 'station', varid=VarID), &
- 'fill_stations', 'inq_varid:station '//trim(filename))
+call nc_check(nf90_inq_varid(ncid, 'voxel', varid=VarID), &
+ 'fill_voxels', 'inq_varid:voxel '//trim(filename))
call nc_check(nf90_inquire_variable(ncid, VarID, ndims=ndims, dimids=dimIDs), &
- 'fill_stations', 'inquire station ndims '//trim(filename))
+ 'fill_voxels', 'inquire voxel ndims '//trim(filename))
if (ndims /= 1) then
- write(string1,*)'station is expected to be a 1D array, it is ',ndims
- call error_handler(E_MSG,'fill_stations:',string1,source,revision,revdate)
+ write(string1,*)'voxel is expected to be a 1D array, it is ',ndims
+ call error_handler(E_MSG,'fill_voxels:',string1)
endif
call nc_check(nf90_inquire_dimension(ncid, dimIDs(1), len=mylen), &
- 'fill_stations', 'station inquire dimid(1) '//trim(filename))
+ 'fill_voxels', 'voxel inquire dimid(1) '//trim(filename))
-if (mylen /= nmax) then
- write(string1,*)'station length expected to be ',nmax,' it is ',mylen
- call error_handler(E_MSG,'fill_stations:',string1,source,revision,revdate)
+if (mylen /= nvoxels) then
+ write(string1,*)'voxel length expected to be ',nvoxels,' it is ',mylen
+ call error_handler(E_MSG,'fill_voxels:',string1)
endif
-allocate(allstations(nmax), obs_types(nmax), which_vert(nmax), &
- ntimearray(nmax), first_time(nmax), last_time(nmax))
-allocate(locations(locNdim,nmax))
-allocate(forecast_leads(nforecasts)) ! how far into the forecast (seconds)
-allocate(ReportTimes(num_verif_times,nmax)) ! observation times that are closest
-allocate(ExperimentTimes(nforecasts,nanalyses)) ! verification times of interest
-allocate(VerifyTimes(nanalyses,nforecasts)) ! verification times of interest
+allocate( voxel_flag(nvoxels), obs_types(nvoxels), which_vert(nvoxels), &
+ first_time(nvoxels), last_time(nvoxels), voxel_level_index(nvoxels))
+allocate(locations(locNdim,nvoxels))
+allocate(ReportTimes(num_verif_times,nvoxels)) ! observation times that are closest
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list