[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