[Dart-dev] DART/branches Revision: 13012

[dart at ucar.edu]

thoar at ucar.edu
2019-03-12 16:00:59 -0600 (Tue, 12 Mar 2019)
77
Supports with or without DART QC8 and the availability of posterior values.




Modified: DART/branches/recam/assimilation_code/programs/obs_diag/threed_sphere/obs_diag.f90
===================================================================
--- DART/branches/recam/assimilation_code/programs/obs_diag/threed_sphere/obs_diag.f90	2019-03-12 17:07:01 UTC (rev 13011)
+++ DART/branches/recam/assimilation_code/programs/obs_diag/threed_sphere/obs_diag.f90	2019-03-12 22:00:59 UTC (rev 13012)
@@ -127,6 +127,7 @@
 integer,  allocatable, dimension(:) :: ens_copy_index
 
 logical :: out_of_range, keeper
+logical :: has_posteriors = .true.
 
 !---------------------------------------------------------------------
 ! variables associated with quality control
@@ -181,7 +182,7 @@
 
 !-----------------------------------------------------------------------
 ! Namelist with default values
-!
+
 character(len=256) :: obs_sequence_name(max_num_input_files) = ''
 character(len=256) :: obs_sequence_list = ''
 integer, dimension(6) :: first_bin_center = (/ 2003, 1, 1, 0, 0, 0 /)
@@ -395,8 +396,9 @@
 call SetScaleFactors() ! for plotting purposes
 
 Nlevels = maxval((/ Nplevels, Nhlevels, Nmlevels /))
-call InitializeVariables( Nepochs, Nlevels, Nregions, num_obs_types)
 
+call InitializeALLVariables( Nepochs, Nlevels, Nregions, num_obs_types)
+
 U_obs_loc = set_location_missing()
 
 ! Open file for histogram of innovations, as a function of standard deviation.
@@ -512,6 +514,23 @@
             prior_spread_index, posterior_spread_index, &
             ens_copy_index )
 
+   if ( any( (/ prior_mean_index, prior_spread_index/) < 0) ) then
+      string1 = 'Observation sequence has no prior information.'
+      string2 = 'You will still get a count, maybe observation value, incoming qc, ...'
+      string3 = 'For simple information, you may want to use "obs_seq_to_netcdf" instead.'
+      call error_handler(E_MSG, 'obs_diag', string1, &
+                 source, revision, revdate, text2=string2, text3=string3)
+   endif
+
+   has_posteriors = .true.
+   if ( any( (/ posterior_mean_index, posterior_spread_index /) < 0) ) then
+      has_posteriors = .false.
+      string1 = 'Observation sequence has no posterior information,'
+      string2 = 'therefore - posterior diagnostics are not possible.'
+      call error_handler(E_WARN, 'obs_diag', string1, &
+                 source, revision, revdate, text2=string2)
+   endif
+
    ! Loop over all potential time periods ... the observation sequence
    ! files are not required to be in any particular order.
 
@@ -612,6 +631,8 @@
          posterior_mean(1)   = 0.0_r8
          prior_spread(1)     = 0.0_r8
          posterior_spread(1) = 0.0_r8
+         pr_zscore           = 0.0_r8
+         po_zscore           = 0.0_r8
 
             call get_obs_values(observation,              obs,              obs_index)
          if (prior_mean_index > 0) &
@@ -723,6 +744,9 @@
          po_zscore = InnovZscore(obs(1), po_mean, po_sprd, obs_error_variance, &
                                  qc_value, QC_MAX_POSTERIOR)
 
+         if (has_posteriors) po_zscore = InnovZscore(obs(1), po_mean, po_sprd, &
+                                    obs_error_variance, qc_value, QC_MAX_POSTERIOR)
+
          indx         = min(int(pr_zscore), MaxSigmaBins)
          nsigma(indx) = nsigma(indx) + 1
 
@@ -788,8 +812,10 @@
             if ((level_index < 1) .or. (level_index > Nlevels)) then
                prior%NbadLV(iepoch,:,iregion,flavor) = &
                prior%NbadLV(iepoch,:,iregion,flavor) + 1
-               poste%NbadLV(iepoch,:,iregion,flavor) = &
-               poste%NbadLV(iepoch,:,iregion,flavor) + 1
+               if (has_posteriors) then
+                  poste%NbadLV(iepoch,:,iregion,flavor) = &
+                  poste%NbadLV(iepoch,:,iregion,flavor) + 1
+               endif
                cycle Areas
             endif
 
@@ -798,8 +824,9 @@
 
             if (      org_qc_index  > 0 ) then
                if (qc(org_qc_index) > input_qc_threshold ) then
-               call IPE(prior%NbigQC(iepoch,level_index,iregion,flavor), 1)
-               call IPE(poste%NbigQC(iepoch,level_index,iregion,flavor), 1)
+                  call IPE(prior%NbigQC(iepoch,level_index,iregion,flavor), 1)
+                  if (has_posteriors) &
+                     call IPE(poste%NbigQC(iepoch,level_index,iregion,flavor), 1)
                endif
             endif
 
@@ -815,7 +842,7 @@