[Dart-dev] [6222] DART/branches/development: The obs_diag program is refactored.
nancy at ucar.edu
nancy at ucar.edu
Sat Jun 1 12:08:16 MDT 2013
Revision: 6222
Author: thoar
Date: 2013-06-01 12:08:15 -0600 (Sat, 01 Jun 2013)
Log Message:
-----------
The obs_diag program is refactored. It also provides the ability
for users to input level edges -or- the traditional midpoints.
This makes it a lot easier to work with heights and depths.
The make_fake program is a prototype to create a synthetic obs_seq.final
with known characteristics so it is possible to TEST obs_diag with
an obs_sequence file that has all the observation types and coordinates
needed.
Modified Paths:
--------------
DART/branches/development/diagnostics/threed_sphere/obs_diag.f90
DART/branches/development/diagnostics/threed_sphere/obs_diag.nml
DART/branches/development/models/bgrid_solo/work/input.nml
Added Paths:
-----------
DART/branches/development/observations/utilities/make_fake_obs.f90
DART/branches/development/observations/utilities/threed_sphere/input.make_fake.nml
DART/branches/development/observations/utilities/threed_sphere/mkmf_make_fake_obs
DART/branches/development/observations/utilities/threed_sphere/path_names_make_fake_obs
-------------- next part --------------
Modified: DART/branches/development/diagnostics/threed_sphere/obs_diag.f90
===================================================================
--- DART/branches/development/diagnostics/threed_sphere/obs_diag.f90 2013-06-01 17:38:29 UTC (rev 6221)
+++ DART/branches/development/diagnostics/threed_sphere/obs_diag.f90 2013-06-01 18:08:15 UTC (rev 6222)
@@ -17,19 +17,20 @@
! All 'possible' obs_kinds are treated separately.
!-----------------------------------------------------------------------
-! In Atmospheric Science, 'spread' has units of standard deviations ...
-!
+! In Atmospheric Science, 'spread' has units of standard deviations ...
+!
! I should rename some of the variables I use as variances to reflect this.
-! 'priorspred' should really be 'priorvar' since you have to accumulate variances
+! 'priorspred' should really be 'priorvar' since you have to accumulate variances
! the math is correct as it is, but the variable names don't make it easy ...
-use types_mod, only : r4, r8, digits12, MISSING_R8, MISSING_R4, metadatalength
+use types_mod, only : r4, r8, digits12, MISSING_R8, MISSING_R4, MISSING_I, &
+ metadatalength
use obs_sequence_mod, only : read_obs_seq, obs_type, obs_sequence_type, get_first_obs, &
get_obs_from_key, get_obs_def, get_copy_meta_data, &
get_obs_time_range, get_time_range_keys, get_num_obs, &
get_next_obs, get_num_times, get_obs_values, init_obs, &
assignment(=), get_num_copies, static_init_obs_sequence, &
- get_qc, destroy_obs_sequence, read_obs_seq_header, &
+ get_qc, destroy_obs_sequence, read_obs_seq_header, &
get_last_obs, destroy_obs, get_num_qc, get_qc_meta_data
use obs_def_mod, only : obs_def_type, get_obs_def_error_variance, get_obs_def_time, &
get_obs_def_location, get_obs_kind
@@ -37,12 +38,13 @@
KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT
use location_mod, only : location_type, get_location, set_location_missing, &
write_location, operator(/=), is_location_in_region, &
- set_location, query_location, &
- vert_is_undef, VERTISUNDEF, &
- vert_is_surface, VERTISSURFACE, &
- vert_is_level, VERTISLEVEL, &
- vert_is_pressure, VERTISPRESSURE, &
- vert_is_height, VERTISHEIGHT
+ set_location, query_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, print_time, &
set_calendar_type, print_date, GREGORIAN, &
operator(*), operator(+), operator(-), &
@@ -52,7 +54,7 @@
file_exist, error_handler, E_ERR, E_WARN, E_MSG, &
initialize_utilities, logfileunit, nmlfileunit, &
find_namelist_in_file, check_namelist_read, &
- nc_check, do_nml_file, do_nml_term, timestamp, &
+ nc_check, do_nml_file, do_nml_term, finalize_utilities, &
next_file, get_next_filename, find_textfile_dims, &
file_to_text
use sort_mod, only : sort
@@ -70,6 +72,7 @@
revdate = '$Date$'
!---------------------------------------------------------------------
+! Some basic parameters
!---------------------------------------------------------------------
integer, parameter :: MaxLevels = 50
@@ -89,11 +92,11 @@
character(len = 129) :: obs_seq_in_file_name
character(len = 129), allocatable, dimension(:) :: obs_seq_filenames
-character(len = stringlength), dimension(MaxTrusted) :: trusted_obsname = 'null'
+character(len = stringlength), dimension(MaxTrusted) :: trusted_list = 'null'
! Storage with fixed size for observation space diagnostics
real(r8), dimension(1) :: prior_mean, posterior_mean, prior_spread, posterior_spread
-real(r8) :: pr_mean, po_mean ! same as above, without useless dimension
+real(r8) :: pr_mean, po_mean ! same as above, without useless dimension
real(r8) :: pr_sprd, po_sprd ! same as above, without useless dimension
! We are treating winds as a vector pair, but we are handling the
@@ -113,7 +116,7 @@
integer :: obs_index, prior_mean_index, posterior_mean_index
integer :: prior_spread_index, posterior_spread_index
-integer :: flavor, wflavor ! THIS IS THE (global) 'KIND' in the obs_def_mod list.
+integer :: flavor, wflavor ! THIS IS THE (global) 'KIND' in the obs_def_mod list.
integer :: num_copies, num_qc, num_obs, max_num_obs, obs_seq_file_id
integer :: num_obs_types
@@ -132,7 +135,7 @@
integer, allocatable, dimension(:) :: keys
integer, allocatable, dimension(:) :: ens_copy_index
-logical :: out_of_range, is_there_one, keeper
+logical :: out_of_range, keeper
!---------------------------------------------------------------------
! variables associated with quality control
@@ -141,7 +144,7 @@
! Most frequently represents the value NCEP assigned to their
! observations.
!
-! dart_qc_index
+! dart_qc_index
! 0 observation assimilated
! 1 observation evaluated only
! --- everything above this means the prior and posterior are OK
@@ -181,13 +184,16 @@
integer, dimension(6) :: time_to_skip = (/ 0, 0, 1, 0, 0, 0 /)
integer :: max_num_bins = 1000 ! maximum number of temporal bins to consider
-real(r8), dimension(MaxLevels) :: plevel = MISSING_R8 ! pressure level (hPa)
-real(r8), dimension(MaxLevels) :: hlevel = MISSING_R8 ! height (meters)
-integer, dimension(MaxLevels) :: mlevel = MISSING_R8 ! model level (integer index)
+real(r8), dimension(MaxLevels) :: plevel = MISSING_R8 ! pressure level (hPa)
+real(r8), dimension(MaxLevels) :: hlevel = MISSING_R8 ! height (meters)
+real(r8), dimension(MaxLevels) :: mlevel = MISSING_R8 ! model level (integer index)
+real(r8), dimension(MaxLevels+1) :: plevel_edges = MISSING_R8 ! pressure level (hPa)
+real(r8), dimension(MaxLevels+1) :: hlevel_edges = MISSING_R8 ! height (meters)
+real(r8), dimension(MaxLevels+1) :: mlevel_edges = MISSING_R8 ! model levels (nondimensional)
integer :: Nregions = 0
real(r8), dimension(MaxRegions) :: lonlim1= MISSING_R8, lonlim2= MISSING_R8
-real(r8), dimension(MaxRegions) :: latlim1= MISSING_R8, latlim2= MISSING_R8
+real(r8), dimension(MaxRegions) :: latlim1= MISSING_R8, latlim2= MISSING_R8
character(len = stringlength), dimension(MaxRegions) :: reg_names = 'null'
type(location_type), dimension(MaxRegions) :: min_loc, max_loc
@@ -207,14 +213,14 @@
plevel, hlevel, mlevel, rat_cri, input_qc_threshold, &
Nregions, lonlim1, lonlim2, latlim1, latlim2, &
reg_names, print_mismatched_locs, print_obs_locations, &
- verbose, outliers_in_histogram, &
- create_rank_histogram, hlevel_edges, trusted_obs
+ create_rank_histogram, outliers_in_histogram, &
+ plevel_edges, hlevel_edges, mlevel_edges, &
+ verbose, trusted_obs
!-----------------------------------------------------------------------
! Variables used to accumulate the statistics.
!-----------------------------------------------------------------------
-
integer, parameter :: Ncopies = 22
character(len = stringlength), dimension(Ncopies) :: copy_names = &
(/ 'Nposs ', 'Nused ', 'NbigQC ', 'NbadIZ ', 'NbadUV ', &
@@ -263,17 +269,15 @@
integer, dimension(:,:,:), pointer :: NDartQC_4, NDartQC_5, NDartQC_6, NDartQC_7
end type LRV_type
-type(TLRV_type) :: analy, guess
-type( LRV_type) :: analyAVG, guessAVG
+type(TLRV_type) :: poste, prior
+type( LRV_type) :: posteAVG, priorAVG
-type(time_type), pointer, dimension(:) :: bincenter
-type(time_type), pointer, dimension(:,:) :: binedges
+type(time_type), allocatable, dimension(:) :: bincenter
+type(time_type), allocatable, dimension(:,:) :: binedges
+real(digits12), allocatable, dimension(:) :: epoch_center
+real(digits12), allocatable, dimension(:,:) :: epoch_edges
+integer, allocatable, dimension(:) :: obs_used_in_epoch
-real(digits12), pointer, dimension(:) :: epoch_center
-real(digits12), pointer, dimension(:,:) :: epoch_edges
-
-integer, pointer, dimension(:) :: obs_used_in_epoch
-
!-----------------------------------------------------------------------
! General purpose variables
!-----------------------------------------------------------------------
@@ -282,10 +286,6 @@
integer :: iregion, iepoch, ivar, ifile, num_obs_in_epoch
real(r8) :: obslon, obslat, obslevel, obsloc3(3)
-real(r8), dimension(MaxLevels+1) :: plevel_edges = MISSING_R8 ! pressure level (hPa)
-real(r8), dimension(MaxLevels+1) :: hlevel_edges = MISSING_R8 ! height (meters)
-real(r8), dimension(MaxLevels+1) :: mlevel_edges = MISSING_R8 ! model levels (nondimensional)
-
integer :: obsindex, i, iunit, ierr, io, ikind
integer :: ivert
@@ -293,7 +293,7 @@
integer :: Nlevels, ilev, num_trusted ! counters
integer :: seconds, days, Nepochs
-logical :: matched, trusted
+logical :: trusted
integer, allocatable, dimension(:) :: which_vert ! relates kind of level for each obs kind
real(r8), allocatable, dimension(:) :: scale_factor ! to convert to plotting units
@@ -303,7 +303,7 @@
! Replace calls to 'get_obs_kind_name' ---> index into 'obs_type_strings'
character(len = stringlength), pointer, dimension(:) :: obs_type_strings
-! These pairs of variables are used when we diagnose which observations
+! These pairs of variables are used when we diagnose which observations
! are far from the background.
integer, parameter :: MaxSigmaBins = 100
integer :: nsigma(0:MaxSigmaBins) = 0
@@ -313,23 +313,29 @@
type(time_type) :: TimeMin, TimeMax ! of the entire period of interest
type(time_type) :: beg_time, end_time ! of the particular bin
-type(time_type) :: binsep, binwidth, halfbinwidth
+type(time_type) :: binsep, binwidth, halfbinwidth
type(time_type) :: seqT1, seqTN ! first,last time in entire observation sequence
type(time_type) :: AllseqT1, AllseqTN ! first,last time in ALL observation sequences
type(time_type) :: obs_time, skip_time
-character(len = 129) :: ncName, string1, string2, string3
-character(len = stringlength) :: obsname, possible_obs_type
+character(len = 256) :: ncName, string1, string2, string3
+character(len = stringlength) :: obsname
integer :: Nidentity = 0 ! identity observations
+!FIXME/USEME Dont know the right syntax for this scope
+!interface CountDartQC
+! procedure CountDartQC_3D
+! procedure CountDartQC_4D
+!end interface CountDartQC
+
!=======================================================================
! Get the party started
!=======================================================================
call initialize_utilities('obs_diag')
-call register_module(source,revision,revdate)
-call static_init_obs_sequence() ! Initialize the obs sequence module
+call register_module(source,revision,revdate)
+call static_init_obs_sequence() ! Initialize the obs sequence module
!----------------------------------------------------------------------
! Define/Append the 'horizontal wind' obs_kinds to supplant the list declared
@@ -374,53 +380,8 @@
write(string3,*)'Please remove "print_obs_locations" from your namelists.'
call error_handler(E_ERR, 'obs_diag', string1, source, revision, revdate, &
text2=string2, text3=string3 )
-else
- write(string1,*)'"print_obs_locations" is no longer supported. See "obs_seq_to_netcdf.html"'
- write(string2,*)'Please remove "print_obs_locations" from your namelists.'
- call error_handler(E_WARN, 'obs_diag', string1, source, revision, revdate, text2=string2)
endif
-! Count up the number of 'trusted' observations
-! Check to make sure trusted observations desired are supported.
-num_trusted = 0
-CountTrusted : do i = 1,MaxTrusted
- if (trim(trusted_obs(i)) == 'null') exit CountTrusted
-
- matched = .false.
- VerifyTrusted : do ikind = 1,max_obs_kinds
- possible_obs_type = get_obs_kind_name(ikind)
- if (trim(trusted_obs(i)) == trim(possible_obs_type)) then
- matched = .true.
- exit VerifyTrusted
- endif
- enddo VerifyTrusted
-
- if (matched) then
- num_trusted = num_trusted + 1
- trusted_obsname(num_trusted) = trim(trusted_obs(i))
- else
- write(string1,*)'trusted_obs "',trim(trusted_obs(i)),'" is not a supported observation type.'
- call error_handler(E_WARN, 'obs_diag', string1, source, revision, revdate)
- endif
-enddo CountTrusted
-
-if (num_trusted == MaxTrusted) then
- write(string1,*)'There are ',num_trusted,' "trusted" observation types.'
- write(string2,*)'This is the maximum allowed unless you increase "MaxTrusted" and recompile.'
- call error_handler(E_WARN, 'obs_diag', string1, source, revision, revdate, text2=string2)
-endif
-
-if (num_trusted > 0 ) then
- write(string1,*)'There are ',num_trusted,' "trusted" observation types, they are:'
- call error_handler(E_MSG, 'obs_diag', string1, source, revision, revdate)
- do i = 1,num_trusted
- call error_handler(E_MSG, 'obs_diag', trusted_obsname(i), source, revision, revdate)
- enddo
-else
- write(string1,*)'There are no "trusted" observation types.'
- call error_handler(E_MSG, 'obs_diag', string1, source, revision, revdate)
-endif
-
!----------------------------------------------------------------------
! Check to see if we are including the outlier observations in the
! rank histogram calculation.
@@ -437,271 +398,36 @@
!----------------------------------------------------------------------
call set_calendar_type(GREGORIAN)
-call Convert2Time(beg_time, end_time, skip_time, binsep, binwidth, halfbinwidth)
-call ActOnNamelist( Nregions )
+call Namelist2Times() ! convert namelist times to DART times
+call DefineTimeBins() ! Set the temporal binning variables
+call SetPressureLevels(plevel, plevel_edges, Nplevels)
+call SetHeightLevels( hlevel, hlevel_edges, Nhlevels)
+call SetModelLevels( mlevel, mlevel_edges, Nmlevels)
+call DefineRegions()
+call CountTrustedObsTypes()
+call SetScaleFactors() ! for plotting purposes
-!----------------------------------------------------------------------
-! SetTime rectifies user input and the final binning sequence.
-!----------------------------------------------------------------------
+Nlevels = maxval((/ Nplevels, Nhlevels, Nmlevels /))
-call SetScaleFactors(scale_factor, logfileunit) ! for plotting purposes
-
-call SetRegionLimits(Nregions, lonlim1, lonlim2, latlim1, latlim2, &
- min_loc, max_loc)
-
-if (verbose) then
- write(*,*)
- write(*,*)'pressure levels = ',plevel( 1:Nplevels)
- write(*,*)'pressure interfaces = ',plevel_edges(1:Nplevels+1)
- write(*,*)'height levels = ',hlevel( 1:Nhlevels)
- write(*,*)'height interfaces = ',hlevel_edges(1:Nhlevels+1)
- write(*,*)'model levels = ',mlevel( 1:Nmlevels)
- write(*,*)'model interfaces = ',mlevel_edges(1:Nmlevels+1)
- write(*,*)
- do i = 1,Nregions
- write(*,'(''Region '',i02,1x,a32,'' (WESN): '',4(f10.4,1x))') i, &
- reg_names(i), lonlim1(i),lonlim2(i),latlim1(i),latlim2(i)
- enddo
-endif
-
!----------------------------------------------------------------------
-! SetTime rectifies user input and the final binning sequence.
+! allocate space, initialize the variables that will hold the statistics
!----------------------------------------------------------------------
-call SetTime(beg_time, end_time, binsep, halfbinwidth, &
- TimeMin, TimeMax, Nepochs, bincenter, binedges, epoch_center, epoch_edges, &
- obs_used_in_epoch)
-
-prior_mean(1) = 0.0_r8
-prior_spread(1) = 0.0_r8
-posterior_mean(1) = 0.0_r8
-posterior_spread(1) = 0.0_r8
-
-U_obs_loc = set_location_missing()
-
-!----------------------------------------------------------------------
-! Prepare the variables
-!----------------------------------------------------------------------
-
-Nlevels = maxval((/ Nplevels, Nhlevels, Nmlevels /))
-
-write(*,*)'level dimension can be ',Nlevels
-
allocate(obs_seq_filenames(Nepochs*400))
obs_seq_filenames = 'null'
-allocate(guess%rmse( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%bias( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%spread( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%totspread( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%observation(Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%ens_mean( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%Nposs( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%Nused( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NbigQC( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NbadIZ( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NbadUV( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NbadLV( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NbadDartQC( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NDartQC_0( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NDartQC_1( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NDartQC_2( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NDartQC_3( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NDartQC_4( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NDartQC_5( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NDartQC_6( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%NDartQC_7( Nepochs, Nlevels, Nregions, num_obs_types), &
- guess%Ntrusted( Nepochs, Nlevels, Nregions, num_obs_types) )
+call InitializeVariables( Nepochs, Nlevels, Nregions, num_obs_types)
-guess%rmse = 0.0_r8
-guess%bias = 0.0_r8
-guess%spread = 0.0_r8
-guess%totspread = 0.0_r8
-guess%observation = 0.0_r8
-guess%ens_mean = 0.0_r8
-guess%Nposs = 0
-guess%Nused = 0
-guess%NbigQC = 0
-guess%NbadIZ = 0
-guess%NbadUV = 0
-guess%NbadLV = 0
-guess%NbadDartQC = 0
-guess%NDartQC_0 = 0
-guess%NDartQC_1 = 0
-guess%NDartQC_2 = 0
-guess%NDartQC_3 = 0
-guess%NDartQC_4 = 0
-guess%NDartQC_5 = 0
-guess%NDartQC_6 = 0
-guess%NDartQC_7 = 0
-guess%Ntrusted = 0
+U_obs_loc = set_location_missing()
-guess%string = 'guess'
-guess%num_times = Nepochs
-guess%num_levels = Nlevels
-guess%num_regions = Nregions
-guess%num_variables = num_obs_types
-
-allocate(analy%rmse( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%bias( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%spread( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%totspread( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%observation(Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%ens_mean( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%Nposs( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%Nused( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NbigQC( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NbadIZ( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NbadUV( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NbadLV( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NbadDartQC( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NDartQC_0( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NDartQC_1( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NDartQC_2( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NDartQC_3( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NDartQC_4( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NDartQC_5( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NDartQC_6( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%NDartQC_7( Nepochs, Nlevels, Nregions, num_obs_types), &
- analy%Ntrusted( Nepochs, Nlevels, Nregions, num_obs_types) )
-
-analy%rmse = 0.0_r8
-analy%bias = 0.0_r8
-analy%spread = 0.0_r8
-analy%totspread = 0.0_r8
-analy%observation = 0.0_r8
-analy%ens_mean = 0.0_r8
-analy%Nposs = 0
-analy%Nused = 0
-analy%NbigQC = 0
-analy%NbadIZ = 0
-analy%NbadUV = 0
-analy%NbadLV = 0
-analy%NbadDartQC = 0
-analy%NDartQC_0 = 0
-analy%NDartQC_1 = 0
-analy%NDartQC_2 = 0
-analy%NDartQC_3 = 0
-analy%NDartQC_4 = 0
-analy%NDartQC_5 = 0
-analy%NDartQC_6 = 0
-analy%NDartQC_7 = 0
-analy%Ntrusted = 0
-
-analy%string = 'analy'
-analy%num_times = Nepochs
-analy%num_levels = Nlevels
-analy%num_regions = Nregions
-analy%num_variables = num_obs_types
-
-allocate(guessAVG%rmse( Nlevels, Nregions, num_obs_types), &
- guessAVG%bias( Nlevels, Nregions, num_obs_types), &
- guessAVG%spread( Nlevels, Nregions, num_obs_types), &
- guessAVG%totspread( Nlevels, Nregions, num_obs_types), &
- guessAVG%observation(Nlevels, Nregions, num_obs_types), &
- guessAVG%ens_mean( Nlevels, Nregions, num_obs_types), &
- guessAVG%Nposs( Nlevels, Nregions, num_obs_types), &
- guessAVG%Nused( Nlevels, Nregions, num_obs_types), &
- guessAVG%NbigQC( Nlevels, Nregions, num_obs_types), &
- guessAVG%NbadIZ( Nlevels, Nregions, num_obs_types), &
- guessAVG%NbadUV( Nlevels, Nregions, num_obs_types), &
- guessAVG%NbadLV( Nlevels, Nregions, num_obs_types), &
- guessAVG%NbadDartQC( Nlevels, Nregions, num_obs_types), &
- guessAVG%NDartQC_0( Nlevels, Nregions, num_obs_types), &
- guessAVG%NDartQC_1( Nlevels, Nregions, num_obs_types), &
- guessAVG%NDartQC_2( Nlevels, Nregions, num_obs_types), &
- guessAVG%NDartQC_3( Nlevels, Nregions, num_obs_types), &
- guessAVG%NDartQC_4( Nlevels, Nregions, num_obs_types), &
- guessAVG%NDartQC_5( Nlevels, Nregions, num_obs_types), &
- guessAVG%NDartQC_6( Nlevels, Nregions, num_obs_types), &
- guessAVG%NDartQC_7( Nlevels, Nregions, num_obs_types), &
- guessAVG%Ntrusted( Nlevels, Nregions, num_obs_types) )
-
-guessAVG%rmse = 0.0_r8
-guessAVG%bias = 0.0_r8
-guessAVG%spread = 0.0_r8
-guessAVG%totspread = 0.0_r8
-guessAVG%observation = 0.0_r8
-guessAVG%ens_mean = 0.0_r8
-guessAVG%Nposs = 0
-guessAVG%Nused = 0
-guessAVG%NbigQC = 0
-guessAVG%NbadIZ = 0
-guessAVG%NbadUV = 0
-guessAVG%NbadLV = 0
-guessAVG%NbadDartQC = 0
-guessAVG%NDartQC_0 = 0
-guessAVG%NDartQC_1 = 0
-guessAVG%NDartQC_2 = 0
-guessAVG%NDartQC_3 = 0
-guessAVG%NDartQC_4 = 0
-guessAVG%NDartQC_5 = 0
-guessAVG%NDartQC_6 = 0
-guessAVG%NDartQC_7 = 0
-guessAVG%Ntrusted = 0
-
-guessAVG%string = 'VPguess'
-guessAVG%num_levels = Nlevels
-guessAVG%num_regions = Nregions
-guessAVG%num_variables = num_obs_types
-
-allocate(analyAVG%rmse( Nlevels, Nregions, num_obs_types), &
- analyAVG%bias( Nlevels, Nregions, num_obs_types), &
- analyAVG%spread( Nlevels, Nregions, num_obs_types), &
- analyAVG%totspread( Nlevels, Nregions, num_obs_types), &
- analyAVG%observation(Nlevels, Nregions, num_obs_types), &
- analyAVG%ens_mean( Nlevels, Nregions, num_obs_types), &
- analyAVG%Nposs( Nlevels, Nregions, num_obs_types), &
- analyAVG%Nused( Nlevels, Nregions, num_obs_types), &
- analyAVG%NbigQC( Nlevels, Nregions, num_obs_types), &
- analyAVG%NbadIZ( Nlevels, Nregions, num_obs_types), &
- analyAVG%NbadUV( Nlevels, Nregions, num_obs_types), &
- analyAVG%NbadLV( Nlevels, Nregions, num_obs_types), &
- analyAVG%NbadDartQC( Nlevels, Nregions, num_obs_types), &
- analyAVG%NDartQC_0( Nlevels, Nregions, num_obs_types), &
- analyAVG%NDartQC_1( Nlevels, Nregions, num_obs_types), &
- analyAVG%NDartQC_2( Nlevels, Nregions, num_obs_types), &
- analyAVG%NDartQC_3( Nlevels, Nregions, num_obs_types), &
- analyAVG%NDartQC_4( Nlevels, Nregions, num_obs_types), &
- analyAVG%NDartQC_5( Nlevels, Nregions, num_obs_types), &
- analyAVG%NDartQC_6( Nlevels, Nregions, num_obs_types), &
- analyAVG%NDartQC_7( Nlevels, Nregions, num_obs_types), &
- analyAVG%Ntrusted( Nlevels, Nregions, num_obs_types) )
-
-analyAVG%rmse = 0.0_r8
-analyAVG%bias = 0.0_r8
-analyAVG%spread = 0.0_r8
-analyAVG%totspread = 0.0_r8
-analyAVG%observation = 0.0_r8
-analyAVG%ens_mean = 0.0_r8
-analyAVG%Nposs = 0
-analyAVG%Nused = 0
-analyAVG%NbigQC = 0
-analyAVG%NbadIZ = 0
-analyAVG%NbadUV = 0
-analyAVG%NbadLV = 0
-analyAVG%NbadDartQC = 0
-analyAVG%NDartQC_0 = 0
-analyAVG%NDartQC_1 = 0
-analyAVG%NDartQC_2 = 0
-analyAVG%NDartQC_3 = 0
-analyAVG%NDartQC_4 = 0
-analyAVG%NDartQC_5 = 0
-analyAVG%NDartQC_6 = 0
-analyAVG%NDartQC_7 = 0
-analyAVG%Ntrusted = 0
-
-analyAVG%string = 'VPanaly'
-analyAVG%num_levels = Nlevels
-analyAVG%num_regions = Nregions
-analyAVG%num_variables = num_obs_types
-
!----------------------------------------------------------------------
! Open file for histogram of innovations, as a function of standard deviation.
!----------------------------------------------------------------------
+
nsigmaUnit = open_file('LargeInnov.txt',form='formatted',action='rewind')
write(nsigmaUnit,'(a)')'Any observations flagged as bad are dumped into the last bin.'
-write(nsigmaUnit,'(a)') ' day secs lon lat level obs guess zscore key kind'
+write(nsigmaUnit,'(a)') ' day secs lon lat level obs prior zscore key kind'
+
!-----------------------------------------------------------------------
! We must assume the observation sequence files span an unknown amount
! of time. We must make some sort of assumption about the naming structure
@@ -713,23 +439,19 @@
! Directory/file names are similar to 01_03/obs_seq.final
!
!-----------------------------------------------------------------------
-! The strategy at this point is to open WAY too many files and
+! The strategy at this point is to open WAY too many files and
! check the observation sequences against ALL of the temporal bins.
! If the sequence is completely before the time period of interest, we skip.
! If the sequence is completely past the time period of interest, we stop.
-
-ObsFileLoop : do ifile=1, 1000
!-----------------------------------------------------------------------
- if (obs_sequence_list == '') then ! try to increment filename
+ObsFileLoop : do ifile=1, size(obs_seq_filenames)
+ if (obs_sequence_list == '') then
obs_seq_in_file_name = next_file(obs_sequence_name,ifile)
-
else
-
obs_seq_in_file_name = get_next_filename(obs_sequence_list,ifile)
if (obs_seq_in_file_name == '') exit ObsFileLoop
-
endif
if ( file_exist(trim(obs_seq_in_file_name)) ) then
@@ -737,21 +459,21 @@
call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
else
write(string1,*)trim(obs_seq_in_file_name), ' does not exist. Finishing.'
-
if (ifile == 1) then
call error_handler(E_ERR,'obs_diag',string1,source,revision,revdate)
else
call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
-
exit ObsFileLoop
endif
+ ! Record the filename read - for the netCDF file.
+
+ obs_seq_filenames(ifile) = trim(obs_seq_in_file_name)
+
! 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)
-
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, &
@@ -786,47 +508,11 @@
call read_obs_seq(obs_seq_in_file_name, 0, 0, 0, seq)
- !--------------------------------------------------------------------
! Determine the time encompassed in the observation sequence.
- !--------------------------------------------------------------------
+ ! also compare to first/last times of ALL sequences
- is_there_one = get_first_obs(seq, obs1)
- if ( .not. is_there_one ) then
- call error_handler(E_ERR,'obs_diag','No first observation in sequence.', &
- source,revision,revdate)
- endif
- call get_obs_def(obs1, obs_def)
- seqT1 = get_obs_def_time(obs_def)
+ call GetFirstLastObs(ifile, seq, obs1, obsN, seqT1, seqTN, AllseqT1, AllseqTN)
- is_there_one = get_last_obs(seq, obsN)
- if ( .not. is_there_one ) then
- call error_handler(E_ERR,'obs_diag','No last observation in sequence.', &
- source,revision,revdate)
- endif
- call get_obs_def(obsN, obs_def)
- seqTN = get_obs_def_time(obs_def)
-
- ! Capture a little information to assist in an error message if the
- ! namelist input does not intersect the observation sequence file.
-
- if ( ifile == 1 ) AllseqT1 = seqT1
- AllseqTN = seqTN
-
- call print_time(seqT1,'First observation time',logfileunit)
- call print_time(seqTN,'Last observation time',logfileunit)
- call print_date(seqT1,'First observation date',logfileunit)
- call print_date(seqTN,'Last observation date',logfileunit)
- if ( verbose ) then
- call print_time(seqT1,'First observation time')
- call print_time(seqTN,'Last observation time')
- call print_time(TimeMin,'TimeMin ')
- call print_time(TimeMax,'TimeMax ')
- call print_date(seqT1,'First observation date')
- call print_date(seqTN,'Last observation date')
- call print_date(TimeMin,'TimeMin ')
- call print_date(TimeMax,'TimeMax ')
- endif
-
!--------------------------------------------------------------------
! If the last observation is before the period of interest, move on.
!--------------------------------------------------------------------
@@ -858,7 +544,10 @@
!--------------------------------------------------------------------
if ( seqT1 > TimeMax ) then
- if (verbose) write(*,*)'seqT1 > TimeMax ... stopping.'
+ if (verbose) then
+ write(logfileunit,*)'seqT1 > TimeMax ... stopping.'
+ write( * ,*)'seqT1 > TimeMax ... stopping.'
+ endif
call destroy_obs(obs1)
call destroy_obs(obsN)
call destroy_obs(observation)
@@ -868,23 +557,27 @@
if (allocated(copyvals)) deallocate( copyvals )
exit ObsFileLoop
else
- if (verbose) write(*,*)'seqT1 < TimeMax ... using ',trim(obs_seq_in_file_name)
+ if (verbose) then
+ write(logfileunit,*)'seqT1 < TimeMax ... using ', &
+ trim(obs_seq_in_file_name)
+ write( * ,*)'seqT1 < TimeMax ... using ', &
+ trim(obs_seq_in_file_name)
+ endif
endif
!--------------------------------------------------------------------
- ! Find the index of obs, ensemble mean, spread ... etc.
- !--------------------------------------------------------------------
- ! Only require obs_index to be present; this allows the program
- ! to be run on obs_seq.in files which have no means or spreads.
- ! You can still plot locations, but that's it.
- !--------------------------------------------------------------------
+ ! Prepare some variables for the rank histogram.
+ ! FIXME : Make sure this observation sequence file has the same number of
+ ! ensemble members as 'the first one' (which defines the bins) ...
- ! Make sure this observation sequence file has the same number of
- ! ensemble members as 'the last one' ...
-
ens_size = GetEnsSize()
- if (ens_size > 0) then
+ if ((ens_size == 0) .and. create_rank_histogram) then
+ write(logfileunit,*) 'Cannot create rank histogram. Zero ensemble members.'
+ write( * ,*) 'Cannot create rank histogram. Zero ensemble members.'
+ create_rank_histogram = .false.
+
+ elseif ((ens_size > 0) .and. create_rank_histogram ) then
if (allocated(ens_copy_index)) then
if (size(ens_copy_index) /= ens_size) then
write(string1,'(''expecting '',i3,'' ensemble members, got '',i3)') &
@@ -893,18 +586,23 @@
endif
else
! This should happen exactly once, if at all.
- allocate(guess%hist_bin( Nepochs, Nlevels, Nregions, num_obs_types, ens_size+1))
+ allocate(prior%hist_bin( Nepochs, Nlevels, Nregions, num_obs_types, ens_size+1))
allocate(ens_copy_index(ens_size))
- guess%hist_bin = 0
+ prior%hist_bin = 0
call init_random_seq(ran_seq, seed=23)
endif
- create_rank_histogram = .true.
- if (verbose) write(*,*) 'Creating rank histogram with ',ens_size+1,' bins.'
- else
- write(*,*) 'Cannot create rank histogram.'
- create_rank_histogram = .false.
+ if ( verbose ) then
+ write(logfileunit,*) 'Creating rank histogram with ',ens_size+1,' bins.'
+ write( * ,*) 'Creating rank histogram with ',ens_size+1,' bins.'
+ endif
endif
+ ! Find the index of obs, ensemble mean, spread ... etc.
+ !--------------------------------------------------------------------
+ ! Only require obs_index to be present; this allows the program
+ ! to be run on obs_seq.in files which have no means or spreads.
+ ! You can still plot locations, but that's it.
+ !--------------------------------------------------------------------
! Each observation sequence file can have its copies in any order.
call SetIndices( obs_index, org_qc_index, dart_qc_index, &
@@ -920,12 +618,18 @@
endif
!====================================================================
- EpochLoop : do iepoch = 1, Nepochs
+ ! Loop over all potential time periods ... the observation sequence
+ ! files are not required to be in any particular order.
!====================================================================
+ EpochLoop : do iepoch = 1, Nepochs
+
beg_time = binedges(1,iepoch)
- end_time = binedges(2,iepoch)
+ end_time = binedges(2,iepoch)
+ ! Using linked list information in the observation sequence,
+ ! find the number of observations that are within this epoch.
+
call get_obs_time_range(seq, beg_time, end_time, key_bounds, &
num_obs_in_epoch, out_of_range)
@@ -940,6 +644,8 @@
write(logfileunit, *) 'num_obs_in_epoch (', iepoch, ') = ', num_obs_in_epoch
write( * , *) 'num_obs_in_epoch (', iepoch, ') = ', num_obs_in_epoch
+ ! actually get the indices (keys) to the observations of interest
+
allocate(keys(num_obs_in_epoch))
call get_time_range_keys(seq, key_bounds, num_obs_in_epoch, keys)
@@ -969,7 +675,7 @@
trusted = .false.
if ( num_trusted > 0 ) then
obsname = get_obs_kind_name(flavor)
- trusted = is_observation_trusted(obsname)
+ trusted = is_observation_trusted(obsname)
endif
! Check to see if it is an identity observation.
@@ -986,15 +692,17 @@
! cvrt to hPa
if(vert_is_pressure(obs_loc)) obslevel = 0.01_r8 * obsloc3(3)
- ! same sort of thing for the scale factors
+ ! same sort of thing for the scale factors
obs_error_variance = get_obs_def_error_variance(obs_def)
obs_err_var = obs_error_variance * &
scale_factor(flavor) * scale_factor(flavor)
!--------------------------------------------------------------
- ! Check consistency of the vertical coordinate system
+ ! Check consistency of the vertical coordinate system
+ ! Sometimes observations of the same flavor from different sources
+ ! are defined on different vertical coordinates.
!--------------------------------------------------------------
- call CheckVertical(obs_loc, flavor)
+ call CheckVertical(obs_loc, flavor)
!--------------------------------------------------------------
! Figure out which level the observation relates to ...
@@ -1002,12 +710,12 @@
level_index = ParseLevel(obs_loc, obslevel, flavor)
- if ( 1 == 1 ) then
+ if ( 1 == 2 ) then
write(8,*)'obsindx ',obsindex, keys(obsindex), obsloc3(3), level_index
endif
!--------------------------------------------------------------
- ! Convert the DART QC data to an integer and create histogram
+ ! Convert the DART QC data to an integer and create histogram
!--------------------------------------------------------------
call get_qc(observation, qc)
@@ -1074,7 +782,7 @@
write(*,*)'posterior_spread value is',posterior_spread(1)
write(*,*)'DART QC value is',qc_integer
write(*,*)'observation is trusted',trusted
- do i= 1,num_copies
+ do i= 1,num_copies
write(*,*)copyvals(i),trim(get_copy_meta_data(seq,i))
enddo
write(*,*)
@@ -1108,8 +816,8 @@
!--------------------------------------------------------------
! update the histogram of the magnitude of the innovation,
- ! where each bin is a single standard deviation.
- ! This is a one-sided histogram.
+ ! where each bin is a single standard deviation.
+ ! This is a one-sided histogram.
!--------------------------------------------------------------
pr_zscore = InnovZscore(obs(1), pr_mean, pr_sprd, obs_err_var, qc_integer, QC_MAX_PRIOR)
@@ -1119,22 +827,28 @@
nsigma(indx) = nsigma(indx) + 1
! Individual (valid) observations that are very far away get
- ! logged to a separate file.
+ ! logged to a separate file but remain in the processing stream.
+ ! In the days before a DART QC, this was the big rejection step.
if( (pr_zscore > 3.0_r8) .and. (qc_integer <= QC_MAX_PRIOR) ) then
call get_time(obs_time,seconds,days)
+ ivert = nint(query_location(obs_loc))
write(nsigmaUnit,FMT='(i7,1x,i5,1x,2f8.2,i7,1x,2f13.2,f8.1,2i7)') &
days, seconds, obslon, obslat, ivert, &
obs(1), pr_mean, pr_zscore, keys(obsindex), flavor
endif
+ !--------------------------------------------------------------
+ ! At this point, the observation has passed all checks.
+ !--------------------------------------------------------------
+
obs_used_in_epoch(iepoch) = obs_used_in_epoch(iepoch) + 1
!--------------------------------------------------------------
! If it is a U wind component, all we need to do is save it.
! It will be matched up with the subsequent V component.
- ! At some point we have to remove the dependency that the
+ ! At some point we have to remove the dependency that the
! U component MUST preceed the V component.
!--------------------------------------------------------------
@@ -1154,7 +868,7 @@
endif
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list