[Dart-dev] [5865] DART/branches/development/diagnostics/threed_sphere/obs_diag.f90: Beta test for three things:
nancy at ucar.edu
nancy at ucar.edu
Thu Sep 13 11:36:01 MDT 2012
Revision: 5865
Author: thoar
Date: 2012-09-13 11:36:00 -0600 (Thu, 13 Sep 2012)
Log Message:
-----------
Beta test for three things:
1) can define a list of 'trusted' observation types via a 'trusted_obs' namelist item.
2) can now define height_edges -OR- heights via the namelist, BUT NOT BOTH.
If that passes, then the methodology will be implemented for pressure and
model levels, just to be consistent.
3) can turn off the rank histogram via a 'create_rank_histogram' namelist item.
Additionally, I removed the deprecated "print observation locations" feature
that has been superceded by the 'obs_seq_to_netcdf' program. I imagine there will
be screams about that one.
Modified Paths:
--------------
DART/branches/development/diagnostics/threed_sphere/obs_diag.f90
-------------- next part --------------
Modified: DART/branches/development/diagnostics/threed_sphere/obs_diag.f90
===================================================================
--- DART/branches/development/diagnostics/threed_sphere/obs_diag.f90 2012-09-12 20:53:06 UTC (rev 5864)
+++ DART/branches/development/diagnostics/threed_sphere/obs_diag.f90 2012-09-13 17:36:00 UTC (rev 5865)
@@ -32,7 +32,7 @@
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, get_obs_name
+ get_obs_def_location, get_obs_kind
use obs_kind_mod, only : max_obs_kinds, get_obs_kind_var_type, get_obs_kind_name, &
KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT
use location_mod, only : location_type, get_location, set_location_missing, &
@@ -74,6 +74,7 @@
integer, parameter :: MaxLevels = 50
integer, parameter :: MaxRegions = 50
+integer, parameter :: MaxTrusted = 50
integer, parameter :: stringlength = 32
!---------------------------------------------------------------------
@@ -88,6 +89,7 @@
character(len = 129) :: obs_seq_in_file_name
character(len = 129), allocatable, dimension(:) :: obs_seq_filenames
+character(len = stringlength), dimension(MaxTrusted) :: trusted_obsname = 'null'
! Storage with fixed size for observation space diagnostics
real(r8), dimension(1) :: prior_mean, posterior_mean, prior_spread, posterior_spread
@@ -122,7 +124,6 @@
character(len=129) :: obs_seq_read_format
logical :: pre_I_format
-logical :: only_print_locations = .false.
integer, dimension(2) :: key_bounds
real(r8), dimension(1) :: obs
@@ -131,7 +132,7 @@
integer, allocatable, dimension(:) :: keys
integer, allocatable, dimension(:) :: ens_copy_index
-logical :: out_of_range, is_there_one, keeper, create_rank_histogram
+logical :: out_of_range, is_there_one, keeper
!---------------------------------------------------------------------
! variables associated with quality control
@@ -154,7 +155,7 @@
! 8+ reserved for future use
integer :: org_qc_index, dart_qc_index
-integer :: qc_integer, my_qc_integer
+integer :: qc_integer
integer, parameter :: QC_MAX = 8
integer, parameter :: QC_MAX_PRIOR = 3
integer, parameter :: QC_MAX_POSTERIOR = 1
@@ -163,6 +164,8 @@
real(r8), allocatable, dimension(:) :: copyvals
integer, parameter, dimension(5) :: hist_qcs = (/ 0, 1, 2, 3, 7 /)
+integer, parameter, dimension(5) :: trusted_prior_qcs = (/ 0, 1, 2, 3, 7 /)
+integer, parameter, dimension(3) :: trusted_poste_qcs = (/ 0, 1, 7 /)
integer :: numqcvals
!-----------------------------------------------------------------------
@@ -188,32 +191,36 @@
character(len = stringlength), dimension(MaxRegions) :: reg_names = 'null'
type(location_type), dimension(MaxRegions) :: min_loc, max_loc
+character(len = stringlength), dimension(MaxTrusted) :: trusted_obs = 'null'
+
real(r8):: rat_cri = 5000.0_r8 ! QC ratio
real(r8):: input_qc_threshold = 3.0_r8 ! maximum NCEP QC factor
logical :: print_mismatched_locs = .false.
-logical :: print_obs_locations = .false.
logical :: verbose = .false.
logical :: outliers_in_histogram = .false.
+logical :: create_rank_histogram = .true.
namelist /obs_diag_nml/ obs_sequence_name, obs_sequence_list, &
first_bin_center, last_bin_center, &
bin_separation, bin_width, time_to_skip, max_num_bins, &
plevel, hlevel, mlevel, rat_cri, input_qc_threshold, &
Nregions, lonlim1, lonlim2, latlim1, latlim2, &
- reg_names, print_mismatched_locs, print_obs_locations, &
- obs_sequence_list, verbose, outliers_in_histogram
+ reg_names, print_mismatched_locs, &
+ obs_sequence_list, verbose, outliers_in_histogram, &
+ create_rank_histogram, hlevel_edges, trusted_obs
!-----------------------------------------------------------------------
! Variables used to accumulate the statistics.
!-----------------------------------------------------------------------
-integer, parameter :: Ncopies = 21
+
+integer, parameter :: Ncopies = 22
character(len = stringlength), dimension(Ncopies) :: copy_names = &
(/ 'Nposs ', 'Nused ', 'NbigQC ', 'NbadIZ ', 'NbadUV ', &
'NbadLV ', 'rmse ', 'bias ', 'spread ', 'totalspread', &
'NbadDARTQC ', 'observation', 'ens_mean ', &
'N_DARTqc_0 ', 'N_DARTqc_1 ', 'N_DARTqc_2 ', 'N_DARTqc_3 ', &
- 'N_DARTqc_4 ', 'N_DARTqc_5 ', 'N_DARTqc_6 ', 'N_DARTqc_7 ' /)
+ 'N_DARTqc_4 ', 'N_DARTqc_5 ', 'N_DARTqc_6 ', 'N_DARTqc_7 ', 'N_trusted ' /)
type TLRV_type
! statistics by time-level-region-variable
@@ -223,7 +230,7 @@
integer :: variable_dim = 4
character(len=8) :: string
integer :: num_times = 0, num_levels = 0, num_regions = 0, num_variables = 0
- integer, dimension(:,:,:,:), pointer :: Nposs, Nused
+ integer, dimension(:,:,:,:), pointer :: Nposs, Nused, Ntrusted
integer, dimension(:,:,:,:), pointer :: NbigQC ! # original QC values >= input_qc_threshold
integer, dimension(:,:,:,:), pointer :: NbadIZ ! # bad (ie huge) Innovation Zscore
integer, dimension(:,:,:,:), pointer :: NbadUV ! # unmatched U/V wind pairs
@@ -243,7 +250,7 @@
integer :: variable_dim = 3
character(len=8) :: string
integer :: num_levels = 0, num_regions = 0, num_variables = 0
- integer, dimension(:,:,:), pointer :: Nposs, Nused
+ integer, dimension(:,:,:), pointer :: Nposs, Nused, Ntrusted
integer, dimension(:,:,:), pointer :: NbigQC ! # bad (original) QC values
integer, dimension(:,:,:), pointer :: NbadIZ ! # bad (ie huge) Innovation Zscore
integer, dimension(:,:,:), pointer :: NbadUV ! # unmatched U/V wind pairs
@@ -278,13 +285,15 @@
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, lunit, ierr, io
+integer :: obsindex, i, iunit, ierr, io, ikind
integer :: ivert
integer :: level_index
-integer :: Nlevels, ilev ! counters
+integer :: Nlevels, ilev, num_trusted ! counters
integer :: seconds, days, Nepochs
+logical :: matched, 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
integer, allocatable, dimension(:) :: ob_defining_vert ! obs index defining vert coord type
@@ -308,8 +317,8 @@
type(time_type) :: AllseqT1, AllseqTN ! first,last time in ALL observation sequences
type(time_type) :: obs_time, skip_time
-character(len = 129) :: ncName, locName, msgstring
-character(len = stringlength) :: str1, str2, str3
+character(len = 129) :: ncName, string1, string2, string3
+character(len = stringlength) :: obsname, possible_obs_type
integer :: Nidentity = 0 ! identity observations
@@ -351,12 +360,52 @@
if (do_nml_term()) write( * , nml=obs_diag_nml)
if ((obs_sequence_name /= '') .and. (obs_sequence_list /= '')) then
- write(msgstring,*)'specify "obs_sequence_name" or "obs_sequence_list"'
- call error_handler(E_MSG, 'obs_diag', msgstring, source, revision, revdate)
- write(msgstring,*)'set other to an empty string ... i.e. ""'
- call error_handler(E_ERR, 'obs_diag', msgstring, source, revision, revdate)
+ 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_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.
@@ -372,19 +421,10 @@
! Now that we have input, do some checking and setup
!----------------------------------------------------------------------
-call ObsLocationsExist( print_obs_locations )
call set_calendar_type(GREGORIAN)
call Convert2Time(beg_time, end_time, skip_time, binsep, binwidth, halfbinwidth)
call ActOnNamelist( Nregions )
-! layer boundaries - heights get sorted to solve single-layer case
-
-Nplevels = Rmidpoints2edges(plevel, plevel_edges)
-Nhlevels = Rmidpoints2edges(hlevel, hlevel_edges)
-Nmlevels = Imidpoints2edges(mlevel, mlevel_edges)
-
-hlevel_edges(1:Nhlevels+1) = sort(hlevel_edges(1:Nhlevels+1))
-
!----------------------------------------------------------------------
! SetTime rectifies user input and the final binning sequence.
!----------------------------------------------------------------------
@@ -395,12 +435,14 @@
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)
@@ -453,7 +495,8 @@
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%NDartQC_7( Nepochs, Nlevels, Nregions, num_obs_types), &
+ guess%Ntrusted( Nepochs, Nlevels, Nregions, num_obs_types) )
guess%rmse = 0.0_r8
guess%bias = 0.0_r8
@@ -476,6 +519,7 @@
guess%NDartQC_5 = 0
guess%NDartQC_6 = 0
guess%NDartQC_7 = 0
+guess%Ntrusted = 0
guess%string = 'guess'
guess%num_times = Nepochs
@@ -503,7 +547,8 @@
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%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
@@ -526,6 +571,7 @@
analy%NDartQC_5 = 0
analy%NDartQC_6 = 0
analy%NDartQC_7 = 0
+analy%Ntrusted = 0
analy%string = 'analy'
analy%num_times = Nepochs
@@ -553,7 +599,8 @@
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%NDartQC_7( Nlevels, Nregions, num_obs_types), &
+ guessAVG%Ntrusted( Nlevels, Nregions, num_obs_types) )
guessAVG%rmse = 0.0_r8
guessAVG%bias = 0.0_r8
@@ -576,6 +623,7 @@
guessAVG%NDartQC_5 = 0
guessAVG%NDartQC_6 = 0
guessAVG%NDartQC_7 = 0
+guessAVG%Ntrusted = 0
guessAVG%string = 'VPguess'
guessAVG%num_levels = Nlevels
@@ -602,7 +650,8 @@
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%NDartQC_7( Nlevels, Nregions, num_obs_types), &
+ analyAVG%Ntrusted( Nlevels, Nregions, num_obs_types) )
analyAVG%rmse = 0.0_r8
analyAVG%bias = 0.0_r8
@@ -625,6 +674,7 @@
analyAVG%NDartQC_5 = 0
analyAVG%NDartQC_6 = 0
analyAVG%NDartQC_7 = 0
+analyAVG%Ntrusted = 0
analyAVG%string = 'VPanaly'
analyAVG%num_levels = Nlevels
@@ -668,12 +718,17 @@
endif
if ( file_exist(trim(obs_seq_in_file_name)) ) then
- write(msgstring,*)'opening ', trim(obs_seq_in_file_name)
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,*)'opening ', trim(obs_seq_in_file_name)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
else
- write(msgstring,*)trim(obs_seq_in_file_name),&
- ' does not exist. Finishing up.'
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ 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
@@ -704,14 +759,12 @@
write(logfileunit,*)'num_obs is ',num_obs
write(logfileunit,*)'max_num_obs is ',max_num_obs
write(logfileunit,*)'obs_seq_read_format is ',trim(obs_seq_read_format)
- write(logfileunit,*)'pre_I_format is ',pre_I_format
write( * ,*)
write( * ,*)'num_copies is ',num_copies
write( * ,*)'num_qc is ',num_qc
write( * ,*)'num_obs is ',num_obs
write( * ,*)'max_num_obs is ',max_num_obs
write( * ,*)'obs_seq_read_format is ',trim(obs_seq_read_format)
- write( * ,*)'pre_I_format is ',pre_I_format
endif
! Read in the entire observation sequence
@@ -819,9 +872,9 @@
if (ens_size > 0) then
if (allocated(ens_copy_index)) then
if (size(ens_copy_index) /= ens_size) then
- write(msgstring,'(''expecting '',i3,'' ensemble members, got '',i3)') &
+ write(string1,'(''expecting '',i3,'' ensemble members, got '',i3)') &
size(ens_copy_index), ens_size
- call error_handler(E_ERR,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_ERR,'obs_diag',string1,source,revision,revdate)
endif
else
! This should happen exactly once, if at all.
@@ -846,11 +899,9 @@
if ( any( (/ prior_mean_index, prior_spread_index, &
posterior_mean_index, posterior_spread_index /) < 0) ) then
- only_print_locations = .true.
- msgstring = 'observation sequence has no prior/posterior information'
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
- else
- only_print_locations = .false.
+ string1 = 'observation sequence has no prior/posterior information'
+ string2 = 'for location information, use "obs_seq_to_netcdf"'
+ call error_handler(E_ERR,'obs_diag',string1,source,revision,revdate,text2=string2)
endif
!====================================================================
@@ -874,18 +925,6 @@
write(logfileunit, *) 'num_obs_in_epoch (', iepoch, ') = ', num_obs_in_epoch
write( * , *) 'num_obs_in_epoch (', iepoch, ') = ', num_obs_in_epoch
- ! Open each epoch file here if writing out auxiliary location files
- if (print_obs_locations) then
- ! Append epoch number to name
- write(locName,'(a,i3.3,a)') 'observation_locations.', iepoch, '.dat'
- if (file_exist(locName)) then
- lunit = open_file(locName, form='formatted', action='append')
- else
- lunit = open_file(locName, form='formatted', action='rewind')
- write(lunit, '(a)') ' lon lat lev kind key QCval'
- endif
- endif
-
allocate(keys(num_obs_in_epoch))
call get_time_range_keys(seq, key_bounds, num_obs_in_epoch, keys)
@@ -911,6 +950,13 @@
obslat = obsloc3(2) ! [-90, 90]
obslevel = obsloc3(3) ! variable-dependent
+ ! Check to see if this is a trusted observation
+ trusted = .false.
+ if ( num_trusted > 0 ) then
+ obsname = get_obs_kind_name(flavor)
+ trusted = is_observation_trusted(obsname)
+ endif
+
! Check to see if it is an identity observation.
! If it is, we count them and skip them since they are better
! explored with the model-space diagnostics.
@@ -941,7 +987,7 @@
level_index = ParseLevel(obs_loc, obslevel, flavor)
- if ( 1 == 2 ) then
+ if ( 1 == 1 ) then
write(8,*)'obsindx ',obsindex, keys(obsindex), obsloc3(3), level_index
endif
@@ -962,31 +1008,6 @@
endif
!--------------------------------------------------------------
- ! Write location of observation if namelist item is true
- !--------------------------------------------------------------
-
- if (print_obs_locations) then
-
- if (dart_qc_index > 0) then
- my_qc_integer = nint(qc(dart_qc_index))
- elseif (org_qc_index > 0) then
- my_qc_integer = -1 * nint(qc(org_qc_index))
- else
- my_qc_integer = -99
- endif
-
- write(lunit, FMT='(3(f10.2,1x),3(i7,1x))') &
- obslon, obslat, obslevel, flavor, keys(obsindex), my_qc_integer
- endif
-
- !--------------------------------------------------------------
- ! Early exit from the observation loop if the observation
- ! does not have all the required copies (attributes).
- !--------------------------------------------------------------
-
- if (only_print_locations) cycle ObservationLoop
-
- !--------------------------------------------------------------
! retrieve observation prior and posterior means and spreads
!--------------------------------------------------------------
@@ -1029,6 +1050,7 @@
if (.true.) then
write(*,*)
write(*,*)'Observation index is ',keys(obsindex),' and has:'
+ write(*,*)'observation type is',obsname
write(*,*)'flavor is',flavor
write(*,*)'obs value is',obs(1)
write(*,*)'prior_mean value is',prior_mean(1)
@@ -1036,6 +1058,7 @@
write(*,*)'prior_spread value is',prior_spread(1)
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
write(*,*)copyvals(i),trim(get_copy_meta_data(seq,i))
enddo
@@ -1213,9 +1236,9 @@
! Do all the heavy lifting
!-----------------------------------------------------------
- call Bin4D(qc_integer, iepoch, level_index, iregion, flavor, &
- obs(1), obs_err_var, pr_mean, pr_sprd, po_mean, po_sprd, &
- rank_histogram_bin)
+ call Bin4D(qc_integer, iepoch, level_index, iregion, flavor, trusted, &
+ obs(1), obs_err_var, pr_mean, pr_sprd, po_mean, po_sprd, &
+ rank_histogram_bin)
!-----------------------------------------------------------
! Additional work for horizontal wind (given U,V)
@@ -1260,7 +1283,7 @@
endif
call Bin4D(maxval( (/qc_integer, U_qc/) ), &
- iepoch, level_index, iregion, wflavor, &
+ iepoch, level_index, iregion, wflavor, .false., &
obs(1), obs_err_var, pr_mean, pr_sprd, po_mean, po_sprd, &
rank_histogram_bin, &
U_obs, U_obs_err_var, U_pr_mean, U_pr_sprd, U_po_mean, U_po_sprd )
@@ -1335,7 +1358,7 @@
call IPE(analyAVG%NbadIZ(level_index,iregion,flavor), 1)
endif
- call Bin3D(qc_integer, level_index, iregion, flavor, &
+ call Bin3D(qc_integer, level_index, iregion, flavor, trusted, &
obs(1), obs_err_var, pr_mean, pr_sprd, po_mean, po_sprd )
! Handle horizontal wind given U,V components
@@ -1360,7 +1383,7 @@
call IPE(analyAVG%NbadIZ(level_index,iregion,wflavor), 1)
endif
- call Bin3D(maxval( (/qc_integer, U_qc/) ), level_index, iregion, wflavor, &
+ call Bin3D(maxval( (/qc_integer, U_qc/) ), level_index, iregion, wflavor, .false., &
obs(1), obs_err_var, pr_mean, pr_sprd, po_mean, po_sprd, &
U_obs, U_obs_err_var, U_pr_mean, U_pr_sprd, U_po_mean, U_po_sprd )
endif
@@ -1378,8 +1401,6 @@
deallocate(keys)
- if (print_obs_locations) close(lunit)
-
enddo EpochLoop
if (verbose) then
@@ -1680,10 +1701,10 @@
write( * ,*)
write(logfileunit,*)
- write(msgstring,*)'WARNING: NO OBSERVATIONS in requested time bins.'
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,*)'WARNING: NO OBSERVATIONS in requested time bins.'
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
@@ -1703,7 +1724,7 @@
deallocate(guess%rmse, guess%bias, guess%spread, guess%totspread, &
guess%observation, guess%ens_mean, guess%Nposs, guess%Nused, &
guess%NbigQC, guess%NbadIZ, guess%NbadUV, guess%NbadLV, &
- guess%NbadDartQC, &
+ guess%NbadDartQC, guess%Ntrusted, &
guess%NDartQC_0, guess%NDartQC_1, guess%NDartQC_2, guess%NDartQC_3, &
guess%NDartQC_4, guess%NDartQC_5, guess%NDartQC_6, guess%NDartQC_7)
@@ -1714,7 +1735,7 @@
deallocate(analy%rmse, analy%bias, analy%spread, analy%totspread, &
analy%observation, analy%ens_mean, analy%Nposs, analy%Nused, &
analy%NbigQC, analy%NbadIZ, analy%NbadUV, analy%NbadLV, &
- analy%NbadDartQC, &
+ analy%NbadDartQC, analy%Ntrusted, &
analy%NDartQC_0, analy%NDartQC_1, analy%NDartQC_2, analy%NDartQC_3, &
analy%NDartQC_4, analy%NDartQC_5, analy%NDartQC_6, analy%NDartQC_7)
@@ -1724,7 +1745,8 @@
guessAVG%NbadIZ, guessAVG%NbadUV, guessAVG%NbadLV, &
guessAVG%NbadDartQC, guessAVG%NDartQC_0, guessAVG%NDartQC_1, &
guessAVG%NDartQC_2, guessAVG%NDartQC_3, guessAVG%NDartQC_4, &
- guessAVG%NDartQC_5, guessAVG%NDartQC_6, guessAVG%NDartQC_7)
+ guessAVG%NDartQC_5, guessAVG%NDartQC_6, guessAVG%NDartQC_7, &
+ guessAVG%Ntrusted )
deallocate(analyAVG%rmse, analyAVG%bias, analyAVG%spread, &
analyAVG%totspread, analyAVG%observation, analyAVG%ens_mean, &
@@ -1732,7 +1754,8 @@
analyAVG%NbadIZ, analyAVG%NbadUV, analyAVG%NbadLV, &
analyAVG%NbadDartQC, analyAVG%NDartQC_0, analyAVG%NDartQC_1, &
analyAVG%NDartQC_2, analyAVG%NDartQC_3, analyAVG%NDartQC_4, &
- analyAVG%NDartQC_5, analyAVG%NDartQC_6, analyAVG%NDartQC_7)
+ analyAVG%NDartQC_5, analyAVG%NDartQC_6, analyAVG%NDartQC_7, &
+ analyAVG%Ntrusted )
deallocate(epoch_center, epoch_edges, bincenter, obs_used_in_epoch)
@@ -1810,7 +1833,7 @@
! Global-scope variables in use in this routine:
! max_num_bins comes from namelist
! verbose comes from namelist
- ! msgstring
+ ! string1
! logfileunit
! Determine temporal bin characteristics.
@@ -1884,22 +1907,22 @@
do iepoch = 1,Nepochs
write(logfileunit,*)
write( * ,*)
- write(str1,'(''epoch '',i6,'' start'')')iepoch
- write(str2,'(''epoch '',i6,'' center'')')iepoch
- write(str3,'(''epoch '',i6,'' end'')')iepoch
- call print_time(binedges(1,iepoch),str1,logfileunit)
- call print_time(bincenter( iepoch),str2,logfileunit)
- call print_time(binedges(2,iepoch),str3,logfileunit)
- call print_time(binedges(1,iepoch),str1)
- call print_time(bincenter( iepoch),str2)
- call print_time(binedges(2,iepoch),str3)
+ write(string1,'(''epoch '',i6,'' start'')')iepoch
+ write(string2,'(''epoch '',i6,'' center'')')iepoch
+ write(string3,'(''epoch '',i6,'' end'')')iepoch
+ call print_time(binedges(1,iepoch),string1,logfileunit)
+ call print_time(bincenter( iepoch),string2,logfileunit)
+ call print_time(binedges(2,iepoch),string3,logfileunit)
+ call print_time(binedges(1,iepoch),string1)
+ call print_time(bincenter( iepoch),string2)
+ call print_time(binedges(2,iepoch),string3)
- call print_date(binedges(1,iepoch),str1,logfileunit)
- call print_date(bincenter( iepoch),str2,logfileunit)
- call print_date(binedges(2,iepoch),str3,logfileunit)
- call print_date(binedges(1,iepoch),str1)
- call print_date(bincenter( iepoch),str2)
- call print_date(binedges(2,iepoch),str3)
+ call print_date(binedges(1,iepoch),string1,logfileunit)
+ call print_date(bincenter( iepoch),string2,logfileunit)
+ call print_date(binedges(2,iepoch),string3,logfileunit)
+ call print_date(binedges(1,iepoch),string1)
+ call print_date(bincenter( iepoch),string2)
+ call print_date(binedges(2,iepoch),string3)
enddo
write(logfileunit,*)
write( * ,*)
@@ -1938,23 +1961,23 @@
! do some error-checking first
if ( (bin_separation(1) /= 0) .or. (bin_separation(2) /= 0) ) then
- write(msgstring,*)'bin_separation:year,month must both be zero, they are ', &
+ write(string1,*)'bin_separation:year,month must both be zero, they are ', &
bin_separation(1),bin_separation(2)
- call error_handler(E_WARN,'obs_diag:Convert2Time',msgstring,source,revision,revdate)
+ call error_handler(E_WARN,'obs_diag:Convert2Time',string1,source,revision,revdate)
error_out = .true.
endif
if ( (bin_width(1) /= 0) .or. (bin_width(2) /= 0) ) then
- write(msgstring,*)'bin_width:year,month must both be zero, they are ', &
+ write(string1,*)'bin_width:year,month must both be zero, they are ', &
bin_width(1),bin_width(2)
- call error_handler(E_WARN,'obs_diag:Convert2Time',msgstring,source,revision,revdate)
+ call error_handler(E_WARN,'obs_diag:Convert2Time',string1,source,revision,revdate)
error_out = .true.
endif
if ( (time_to_skip(1) /= 0) .or. (time_to_skip(2) /= 0) ) then
- write(msgstring,*)'time_to_skip:year,month must both be zero, they are ', &
+ write(string1,*)'time_to_skip:year,month must both be zero, they are ', &
time_to_skip(1),time_to_skip(2)
- call error_handler(E_WARN,'obs_diag:Convert2Time',msgstring,source,revision,revdate)
+ call error_handler(E_WARN,'obs_diag:Convert2Time',string1,source,revision,revdate)
error_out = .true.
endif
@@ -2065,19 +2088,31 @@
if (all(plevel < 1.0)) then
plevel(1:11) = (/ 1000, 925, 850, 700, 500, 400, &
300, 250, 200, 150, 100 /)
+ Nplevels = Rlevels2edges(plevel, plevel_edges)
endif
! Height levels
- if (all(hlevel < 1.0)) then
+ if( all(hlevel_edges == MISSING_R8) .and. &
+ all(hlevel == MISSING_R8)) then
hlevel(1:11) = (/ 1000, 2000, 3000, 4000, 5000, 6000, &
7000, 8000, 9000, 10000, 11000 /)
+ Nhlevels = Rlevels2edges(hlevel, hlevel_edges)
+ hlevel_edges(1:Nhlevels+1) = sort(hlevel_edges(1:Nhlevels+1))
+ elseif ( any(hlevel_edges /= MISSING_R8) ) then
+ ! user specified edges and we have to find the midpoints
+ Nhlevels = Redges2levels(hlevel_edges, hlevel)
+ else
+ ! user specified midpoints and we have to find the edges
+ Nhlevels = Rlevels2edges(hlevel, hlevel_edges)
+ hlevel_edges(1:Nhlevels+1) = sort(hlevel_edges(1:Nhlevels+1))
endif
! Model levels
if (all(mlevel < 1)) then
mlevel(1:11) = (/ (i,i=1,11) /) ! set model levels to indices
+ Nmlevels = Ilevels2edges(mlevel, mlevel_edges)
endif
end Subroutine ActOnNamelist
@@ -2105,8 +2140,8 @@
GetEnsSize = GetEnsSize + 1
enddo MetaDataLoop
- write(msgstring,'(''There are '',i4,'' ensemble members.'')') GetEnsSize
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,'(''There are '',i4,'' ensemble members.'')') GetEnsSize
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
end Function GetEnsSize
@@ -2172,28 +2207,28 @@
!--------------------------------------------------------------------
if ( prior_mean_index < 0 ) then
- write(msgstring,*)'metadata:prior ensemble mean not found'
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,*)'metadata:prior ensemble mean not found'
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if ( posterior_mean_index < 0 ) then
- write(msgstring,*)'metadata:posterior ensemble mean not found'
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,*)'metadata:posterior ensemble mean not found'
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if ( prior_spread_index < 0 ) then
- write(msgstring,*)'metadata:prior ensemble spread not found'
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,*)'metadata:prior ensemble spread not found'
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if ( posterior_spread_index < 0 ) then
- write(msgstring,*)'metadata:posterior ensemble spread not found'
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,*)'metadata:posterior ensemble spread not found'
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if ( org_qc_index < 0 ) then
- write(msgstring,*)'metadata:Quality Control not found'
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,*)'metadata:Quality Control not found'
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if ( dart_qc_index < 0 ) then
- write(msgstring,*)'metadata:DART quality control not found'
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,*)'metadata:DART quality control not found'
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
! Only require obs_index to be present; this allows the program
@@ -2201,63 +2236,63 @@
! less info from them, but you can still plot locations, etc.
if ( obs_index < 0 ) then
- write(msgstring,*)'metadata:observation not found'
- call error_handler(E_ERR,'obs_diag',msgstring,source,revision,revdate)
+ write(string1,*)'metadata:observation not found'
+ call error_handler(E_ERR,'obs_diag',string1,source,revision,revdate)
endif
!--------------------------------------------------------------------
! Echo what we found.
!--------------------------------------------------------------------
- write(msgstring,'(''observation index '',i2,'' metadata '',a)') &
+ write(string1,'(''observation index '',i2,'' metadata '',a)') &
obs_index, trim(get_copy_meta_data(seq,obs_index))
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
if (prior_mean_index > 0 ) then
- write(msgstring,'(''prior mean index '',i2,'' metadata '',a)') &
+ write(string1,'(''prior mean index '',i2,'' metadata '',a)') &
prior_mean_index, trim(get_copy_meta_data(seq,prior_mean_index))
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if (posterior_mean_index > 0 ) then
- write(msgstring,'(''posterior mean index '',i2,'' metadata '',a)') &
+ write(string1,'(''posterior mean index '',i2,'' metadata '',a)') &
posterior_mean_index, trim(get_copy_meta_data(seq,posterior_mean_index))
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if (prior_spread_index > 0 ) then
- write(msgstring,'(''prior spread index '',i2,'' metadata '',a)') &
+ write(string1,'(''prior spread index '',i2,'' metadata '',a)') &
prior_spread_index, trim(get_copy_meta_data(seq,prior_spread_index))
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if (posterior_spread_index > 0 ) then
- write(msgstring,'(''posterior spread index '',i2,'' metadata '',a)') &
+ write(string1,'(''posterior spread index '',i2,'' metadata '',a)') &
posterior_spread_index, trim(get_copy_meta_data(seq,posterior_spread_index))
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if (org_qc_index > 0 ) then
- write(msgstring,'(''Quality Control index '',i2,'' metadata '',a)') &
+ write(string1,'(''Quality Control index '',i2,'' metadata '',a)') &
org_qc_index, trim(get_qc_meta_data(seq,org_qc_index))
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
if (dart_qc_index > 0 ) then
- write(msgstring,'(''DART quality control index '',i2,'' metadata '',a)') &
+ write(string1,'(''DART quality control index '',i2,'' metadata '',a)') &
dart_qc_index, trim(get_qc_meta_data(seq,dart_qc_index))
- call error_handler(E_MSG,'obs_diag',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
endif
end Subroutine SetIndices
- Function Rmidpoints2edges(levels, midpoints)
- ! determine layer midpoints ... ascending vs. descending complications
- integer :: Rmidpoints2edges
+ Function Rlevels2edges(levels, leveledges)
+ ! determine layer edges ... ascending vs. descending complications
+ integer :: Rlevels2edges
real(r8), dimension(:), intent(inout) :: levels
- real(r8), dimension(:), intent( out) :: midpoints
+ real(r8), dimension(:), intent( out) :: leveledges
logical :: increasing
integer :: n, i
@@ -2265,10 +2300,10 @@
integer, dimension(SIZE(levels)) :: indxs
logical, dimension(SIZE(levels)) :: logicarray
- Rmidpoints2edges = SIZE(levels,1)
- if (Rmidpoints2edges > MaxLevels) then
- write(msgstring,*)Rmidpoints2edges,' is too many levels - max is ',Maxlevels,' (Maxlevels)'
- call error_handler(E_ERR,'obs_diag:Rmidpoints2edges', msgstring,source,revision,revdate)
+ Rlevels2edges = SIZE(levels,1)
+ if (Rlevels2edges > MaxLevels) then
+ write(string1,*)Rlevels2edges,' is too many levels - max is ',Maxlevels,' (Maxlevels)'
+ call error_handler(E_ERR,'obs_diag:Rlevels2edges', string1,source,revision,revdate)
endif
! find length of useful portion of levels array
@@ -2277,12 +2312,12 @@
dummyindxs = maxloc(indxs, logicarray)
n = dummyindxs(1)
- Rmidpoints2edges = n
+ Rlevels2edges = n
! single-level case works for pressures, must change sign for heights
if ( n == 1 ) then
- midpoints(1) = levels(1) + 0.1 * levels(1)
- midpoints(2) = levels(1) - 0.1 * levels(1)
+ leveledges(1) = levels(1) + 0.1 * levels(1)
+ leveledges(2) = levels(1) - 0.1 * levels(1)
return
endif
@@ -2291,33 +2326,33 @@
increasing = .true.
levels(1:n) = sort(levels(1:n))
else
- increasing = .false.
- midpoints(1:n) = sort(levels(1:n))
- levels(1:n) = midpoints(n:1:-1)
+ increasing = .false.
+ leveledges(1:n) = sort(levels(1:n))
+ levels(1:n) = leveledges(n:1:-1)
endif
- midpoints( 1) = levels(1) - (levels(2) - levels( 1))/2.0_r8
- midpoints(n+1) = levels(n) + (levels(n) - levels(n-1))/2.0_r8
+ leveledges( 1) = levels(1) - (levels(2) - levels( 1))/2.0_r8
+ leveledges(n+1) = levels(n) + (levels(n) - levels(n-1))/2.0_r8
do i = 2,n
- midpoints(i) = levels(i-1) + (levels(i) - levels(i-1))/2.0_r8
+ leveledges(i) = levels(i-1) + (levels(i) - levels(i-1))/2.0_r8
enddo
if ( increasing ) then ! must be heights
- midpoints( 1) = max( 0.0_r8, midpoints( 1))
+ leveledges( 1) = max( 0.0_r8, leveledges( 1))
else ! must be pressures
- midpoints( 1) = min(1025.0_r8, midpoints( 1))
- midpoints(n+1) = max( 0.0_r8, midpoints(n+1))
+ leveledges( 1) = min(1025.0_r8, leveledges( 1))
+ leveledges(n+1) = max( 0.0_r8, leveledges(n+1))
endif
- end Function Rmidpoints2edges
+ end Function Rlevels2edges
- Function Imidpoints2edges(levelmids, leveledges)
+ Function Ilevels2edges(levelmids, leveledges)
! determine layer edges ... needed for plotting purposes
! the midpoints are integers
- integer :: Imidpoints2edges
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list