[Dart-dev] DART/branches Revision: 12991

[dart at ucar.edu]

thoar at ucar.edu
2019-02-27 17:23:40 -0700 (Wed, 27 Feb 2019)
147
Fully supporting the lack of posterior values. If there are no posterior values,
the 'analy' variables will not appear in the output netCDF file.




Modified: DART/branches/cesm_clm/assimilation_code/programs/obs_diag/threed_sphere/obs_diag.f90
===================================================================
--- DART/branches/cesm_clm/assimilation_code/programs/obs_diag/threed_sphere/obs_diag.f90	2019-02-27 21:31:38 UTC (rev 12990)
+++ DART/branches/cesm_clm/assimilation_code/programs/obs_diag/threed_sphere/obs_diag.f90	2019-02-28 00:23:40 UTC (rev 12991)
@@ -128,6 +128,7 @@
 integer,  allocatable, dimension(:) :: ens_copy_index
 
 logical :: out_of_range, keeper
+logical :: has_posteriors = .true.
 
 !---------------------------------------------------------------------
 ! variables associated with quality control
@@ -398,7 +399,7 @@
 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()
 
@@ -631,7 +632,7 @@
          if ( dart_qc_index > 0 ) then
             qc_value = qc(dart_qc_index)
          else
-            ! If there is no dart_qc, this must be a case where we 
+            ! If there is no dart_qc, this must be a case where we
             ! are interested only in getting the location information.
             qc_value = 0
          endif
@@ -723,9 +724,10 @@
 
          pr_zscore = InnovZscore(obs(1), pr_mean, pr_sprd, obs_error_variance, &
                                  qc_value, QC_MAX_PRIOR)
-         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
 
@@ -791,8 +793,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
 
@@ -802,7 +806,8 @@
             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)
+               if (has_posteriors) &
+                  call IPE(poste%NbigQC(iepoch,level_index,iregion,flavor), 1)
                endif
             endif
 
@@ -818,7 +823,7 @@
                call IPE(prior%NbadIZ(iepoch,level_index,iregion,flavor), 1)
             endif
 
-            if( po_zscore > rat_cri ) then
+            if( has_posteriors .and. po_zscore > rat_cri ) then
                call IPE(poste%NbadIZ(iepoch,level_index,iregion,flavor), 1)
             endif
 
@@ -830,7 +835,8 @@
 
             ! Additional work for horizontal wind (given U,V)
 
-            ObsIsWindCheck: if ( get_quantity_for_type_of_obs(flavor) == QTY_V_WIND_COMPONENT ) then
+            ObsIsWindCheck: if ( get_quantity_for_type_of_obs(flavor) == &
+                                               QTY_V_WIND_COMPONENT ) then
 
                ! The big assumption is that the U wind component has
                ! immediately preceeded the V component and has been saved.
@@ -843,7 +849,8 @@
 
                if ( ierr /= 0 ) then
                   call IPE(prior%NbadUV(iepoch, level_index, iregion, flavor), 1)
-                  call IPE(poste%NbadUV(iepoch, level_index, iregion, flavor), 1)
+                  if (has_posteriors) &
+                     call IPE(poste%NbadUV(iepoch, level_index, iregion, flavor), 1)
                else
 
                   ! The next big assumption is that the 'horizontal wind' flavors
@@ -859,13 +866,13 @@
 
                   zscoreU = InnovZscore(U_obs, U_pr_mean, U_pr_sprd, U_obs_err_var, &
                                         U_qc, QC_MAX_PRIOR)
-                  if( (pr_zscore > rat_cri) .or. (zscoreU > rat_cri) )  then
+                  if(pr_zscore > rat_cri .or. zscoreU > rat_cri) then