[Dart-dev] [10352] DART/trunk/diagnostics/threed_sphere/obs_diag.f90: Correctly setting the QC values for the (computed) *HORIZONTAL_WIND metrics.
nancy at ucar.edu
nancy at ucar.edu
Fri Jun 10 11:51:47 MDT 2016
Revision: 10352
Author: thoar
Date: 2016-06-10 11:51:47 -0600 (Fri, 10 Jun 2016)
Log Message:
-----------
Correctly setting the QC values for the (computed) *HORIZONTAL_WIND metrics.
Previous version were not honoring the incoming qc value and so the
number of possible observations was getting over-reported.
This version is a little closer to the style guide requirements
and reinstates some of the indenting that got hosed with commit 10348.
This version also can be used with obs_seq.in and obs_seq.out files
in a limited fashion. The number of observations in each region, the
observation value for the region, the incoming QC for the obs in the regions,
and the summary of what types are actually in the file are all useful -
even if you do not have the prior/posterior quantities.
Modified Paths:
--------------
DART/trunk/diagnostics/threed_sphere/obs_diag.f90
-------------- next part --------------
Modified: DART/trunk/diagnostics/threed_sphere/obs_diag.f90
===================================================================
--- DART/trunk/diagnostics/threed_sphere/obs_diag.f90 2016-06-10 16:41:46 UTC (rev 10351)
+++ DART/trunk/diagnostics/threed_sphere/obs_diag.f90 2016-06-10 17:51:47 UTC (rev 10352)
@@ -231,7 +231,7 @@
!-----------------------------------------------------------------------
integer, parameter :: Ncopies = 22
-character(len=stringlength), dimension(Ncopies) :: copy_names = &
+character(len=stringlength), dimension(Ncopies) :: copy_names = &
(/ 'Nposs ', 'Nused ', 'NbigQC ', 'NbadIZ ', 'NbadUV ', &
'NbadLV ', 'rmse ', 'bias ', 'spread ', 'totalspread', &
'NbadDARTQC ', 'observation', 'ens_mean ', &
@@ -607,8 +607,8 @@
! 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.
+ ! to be run on obs_seq.[in,out] files which have no means or spreads.
+ ! You can still plot obs count, incoming QC, obs values ...
!--------------------------------------------------------------------
! Each observation sequence file can have its copies in any order.
@@ -619,9 +619,11 @@
if ( any( (/ prior_mean_index, prior_spread_index, &
posterior_mean_index, posterior_spread_index /) < 0) ) then
- 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)
+ string1 = 'Observation sequence has no prior/posterior 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
!====================================================================
@@ -680,7 +682,11 @@
obslevel = obsloc3(3) ! variable-dependent
! Check to see if this is a trusted observation
- if ( num_trusted > 0 ) trusted = is_observation_trusted(obsname)
+ if ( num_trusted > 0 ) then
+ trusted = is_observation_trusted(obsname)
+ else
+ trusted = .false.
+ endif
! Check to see if it is an identity observation.
! If it is, we count them and skip them since they are better
@@ -956,7 +962,7 @@
!-----------------------------------------------------------
call CountDartQC_4D(qc_integer, iepoch, level_index, iregion, &
- flavor, prior, poste, po_mean)
+ flavor, prior, poste, posterior_mean=po_mean)
!-----------------------------------------------------------
! Count 'large' innovations
@@ -1021,6 +1027,9 @@
call IPE(poste%NbadIZ(iepoch,level_index,iregion,wflavor), 1)
endif
+ call CountDartQC_4D(qc_integer, iepoch, level_index, iregion, &
+ wflavor, prior, poste, uqc=U_qc)
+
call Bin4D(qc_integer, 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, &
@@ -1054,7 +1063,7 @@
!-----------------------------------------------------------
call CountDartQC_3D(qc_integer, level_index, iregion, flavor, &
- priorAVG, posteAVG, po_mean)
+ priorAVG, posteAVG, posterior_mean=po_mean)
! Count 'large' innovations
@@ -1076,9 +1085,6 @@
ierr = CheckMate(flavor, U_flavor, obs_loc, U_obs_loc, wflavor)
if ( ierr /= 0 ) then
- ! we have already printed this message for this observation
- ! write(string1,*)'vertical : V with no U at index ', keys(obsindex)
- ! call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
call IPE(priorAVG%NbadUV(level_index, iregion, flavor), 1)
call IPE(posteAVG%NbadUV(level_index, iregion, flavor), 1)
else
@@ -1097,6 +1103,9 @@
call IPE(posteAVG%NbadIZ(level_index,iregion,wflavor), 1)
endif
+ call CountDartQC_3D(qc_integer, level_index, &
+ iregion, wflavor, priorAVG, posteAVG, uqc=U_qc)
+
call Bin3D(qc_integer, 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, &
@@ -1249,20 +1258,18 @@
!----------------------------------------------------------------------
if ( sum(obs_used_in_epoch) == 0 ) then
-
+ write( * ,*)
+ write(logfileunit,*)
+ call print_date(AllseqT1,' First observation date',logfileunit)
+ call print_date(AllseqTN,' Last observation date',logfileunit)
+ call print_date( TimeMin,' First REQUESTED date',logfileunit)
+ call print_date( TimeMax,' Last REQUESTED date',logfileunit)
call print_date(AllseqT1,' First observation date')
call print_date(AllseqTN,' Last observation date')
- call print_date( TimeMin,' First requested date')
- call print_date( TimeMax,' Last requested date')
-
- write( * ,*)
- write(logfileunit,*)
-
- 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)
-
+ call print_date( TimeMin,' First REQUESTED date')
+ call print_date( TimeMax,' Last REQUESTED date')
+ write(string1,*)'NO OBSERVATIONS in requested time bins.'
+ call error_handler(E_ERR,'obs_diag',string1,source,revision,revdate)
endif
!----------------------------------------------------------------------
@@ -1287,8 +1294,7 @@
!======================================================================
-
-Function grok_observation_names(my_names)
+function grok_observation_names(my_names)
!----------------------------------------------------------------------
! Define/Append the 'horizontal wind' obs_kinds to supplant the list declared
! in obs_kind_mod.f90 i.e. if there is a RADIOSONDE_U_WIND_COMPONENT
@@ -1390,13 +1396,13 @@
my_names(ivar) = names(ivar)
enddo
-end Function grok_observation_names
+end function grok_observation_names
!======================================================================
-Subroutine Namelist2Times( )
+subroutine Namelist2Times( )
! Global-scope variables read in this routine:
! first_bin_center comes from namelist
@@ -1475,33 +1481,36 @@
if ( verbose ) then
call print_date( beg_time,'requested first bincenter date',logfileunit)
+ call print_date( beg_time,'requested first bincenter date')
call print_date( end_time,'requested last bincenter date',logfileunit)
+ call print_date( end_time,'requested last bincenter date')
call print_date( skip_time,'implied skip-to date',logfileunit)
+ call print_date( skip_time,'implied skip-to date')
+ write(logfileunit,*)
+ write(*,*)
call print_time( beg_time,'requested first bincenter time',logfileunit)
+ call print_time( beg_time,'requested first bincenter time')
call print_time( end_time,'requested last bincenter time',logfileunit)
+ call print_time( end_time,'requested last bincenter time')
call print_time( skip_time,'implied skip-to time',logfileunit)
+ call print_time( skip_time,'implied skip-to time')
+ write(logfileunit,*)
+ write(*,*)
call print_time( binsep,'requested bin separation',logfileunit)
+ call print_time( binsep,'requested bin separation')
call print_time( binwidth,'requested bin width',logfileunit)
+ call print_time( binwidth,'requested bin width')
call print_time(halfbinwidth,'implied halfbinwidth',logfileunit)
-
- call print_date( beg_time,'requested first bincenter date')
- call print_date( end_time,'requested last bincenter date')
- call print_date( skip_time,'implied skip-to date')
- call print_time( beg_time,'requested first bincenter time')
- call print_time( end_time,'requested last bincenter time')
- call print_time( skip_time,'implied skip-to time')
- call print_time( binsep,'requested bin separation')
- call print_time( binwidth,'requested bin width')
call print_time(halfbinwidth,'implied halfbinwidth')
endif
-end Subroutine Namelist2Times
+end subroutine Namelist2Times
!======================================================================
-Subroutine DefineTimeBins()
+subroutine DefineTimeBins()
! Determine temporal bin characteristics.
! The user input is not guaranteed to align on bin centers.
! So -- we will assume the start time is correct and take strides till we
@@ -1612,13 +1621,13 @@
TimeMin = binedges(1, 1) ! minimum time of interest
TimeMax = binedges(2,Nepochs) ! maximum time of interest
-end Subroutine DefineTimeBins
+end subroutine DefineTimeBins
!======================================================================
-Subroutine setPressureLevels(layerMiddles, layerEdges, nLayers)
+subroutine setPressureLevels(layerMiddles, layerEdges, nLayers)
! Supported scenarios:
!
@@ -1666,13 +1675,13 @@
endif
-end Subroutine setPressureLevels
+end subroutine setPressureLevels
!======================================================================
-Subroutine setHeightLevels(layerMiddles, layerEdges, nLayers)
+subroutine setHeightLevels(layerMiddles, layerEdges, nLayers)
! Supported scenarios:
!
@@ -1725,13 +1734,13 @@
endif
-end Subroutine setHeightLevels
+end subroutine setHeightLevels
!======================================================================
-Subroutine setModelLevels(layerMiddles, layerEdges, nLayers)
+subroutine setModelLevels(layerMiddles, layerEdges, nLayers)
! Supported scenarios:
!
@@ -1782,13 +1791,13 @@
endif
-end Subroutine setModelLevels
+end subroutine setModelLevels
!======================================================================
-Subroutine DefineRegions()
+subroutine DefineRegions()
!----------------------------------------------------------------------
! globally-scoped variables used
@@ -1867,13 +1876,13 @@
enddo
endif
-end Subroutine DefineRegions
+end subroutine DefineRegions
!======================================================================
-Subroutine CountTrustedObsTypes()
+subroutine CountTrustedObsTypes()
! Count up the number of 'trusted' observations types from those
! specified in the input namelist 'trusted_obs' argument.
@@ -1953,13 +1962,13 @@
call error_handler(E_MSG, 'CountTrustedObsTypes', string1, source, revision, revdate)
endif
-end Subroutine CountTrustedObsTypes
+end subroutine CountTrustedObsTypes
!======================================================================
-Subroutine SetScaleFactors()
+subroutine SetScaleFactors()
! The surface pressure in the obs_sequence is in Pa, we want to convert
! from Pa to hPa for plotting. The specific humidity is a similar thing.
@@ -1977,6 +1986,11 @@
character(len=stringlength) :: obs_string
integer :: ivar
+if (verbose) then ! whitespace
+ write( * ,*)
+ write(logfileunit,*)
+endif
+
scale_factor = 1.0_r8
! The scale_factor array is dimensioned from obs_kind_mod:num_obs_types
@@ -1994,20 +2008,17 @@
! Somehow, we should plot statistics on the dBZ scale for these ...
! scale_factor(KIND_RADAR_REFLECTIVITY) = 10log10(z)
- if (verbose) then
- write( * ,*)'scaling of ',scale_factor(ivar),obs_string
- write(logfileunit,*)'scaling of ',scale_factor(ivar),obs_string
- endif
+ ! The scaling is summarized in WriteNetCDF if the verbose option is chosen.
enddo
-end Subroutine SetScaleFactors
+end subroutine SetScaleFactors
!======================================================================
-Subroutine InitializeVariables( ntimes, nlevs, nareas, ntypes )
+subroutine InitializeVariables( ntimes, nlevs, nareas, ntypes )
! Global variables set in this routine:
! type(TLRV_type), intent(out) :: poste, prior
@@ -2224,13 +2235,13 @@
posteAVG%num_regions = nareas
posteAVG%num_variables = ntypes
-end Subroutine InitializeVariables
+end subroutine InitializeVariables
!======================================================================
-Subroutine GetFirstLastObs(my_fileid, my_sequence, my_obs1, my_obsN, my_seqT1, my_seqTN, my_AllseqT1, my_AllseqTN)
+subroutine GetFirstLastObs(my_fileid, my_sequence, my_obs1, my_obsN, my_seqT1, my_seqTN, my_AllseqT1, my_AllseqTN)
! We need to know the time of the first and last observations in the sequence,
! primarily just to see if they intersect the desired Epoch window.
! We also record these times so we can report the first/last times of all
@@ -2244,14 +2255,14 @@
type(obs_def_type) :: obs_def
if ( .not. get_first_obs(my_sequence, my_obs1) ) then
- call error_handler(E_ERR,'obs_diag','No first observation in sequence.', &
+ call error_handler(E_ERR,'GetFirstLastObs','No first observation in sequence.', &
source,revision,revdate)
endif
call get_obs_def(my_obs1, obs_def)
my_seqT1 = get_obs_def_time(obs_def)
if ( .not. get_last_obs(my_sequence, my_obsN) ) then
- call error_handler(E_ERR,'obs_diag','No last observation in sequence.', &
+ call error_handler(E_ERR,'GetFirstLastObs','No last observation in sequence.', &
source,revision,revdate)
endif
call get_obs_def(my_obsN, obs_def)
@@ -2283,13 +2294,13 @@
call print_date( TimeMax,'Last date of interest')
endif
-end Subroutine GetFirstLastObs
+end subroutine GetFirstLastObs
!======================================================================
-Function GetEnsSize()
+function GetEnsSize()
!-----------------------------------------------------------------------
! Loop over all the metadata to count the number of ensemble members
! available in the observation sequence file. We need this count to
@@ -2311,18 +2322,19 @@
enddo MetaDataLoop
write(string1,'(''There are '',i4,'' ensemble members.'')') GetEnsSize
-call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
+call error_handler(E_MSG,'GetEnsSize',string1,source,revision,revdate)
-end Function GetEnsSize
+end function GetEnsSize
!======================================================================
-Subroutine SetIndices( obs_index, org_qc_index, dart_qc_index, &
- prior_mean_index, posterior_mean_index, &
- prior_spread_index, posterior_spread_index,&
- ens_copy_index )
+subroutine SetIndices( obs_index, org_qc_index, dart_qc_index, &
+ prior_mean_index, posterior_mean_index, &
+ prior_spread_index, posterior_spread_index, &
+ ens_copy_index )
+
! There are many 'copy' indices that need to be set from the obs_sequence
! metadata. Some are required, some are optional.
@@ -2386,98 +2398,99 @@
!--------------------------------------------------------------------
if ( prior_mean_index < 0 ) then
- write(string1,*)'metadata:prior ensemble mean not found'
- call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
+ write(string1,*)'metadata:"prior ensemble mean" not found'
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if ( posterior_mean_index < 0 ) then
- write(string1,*)'metadata:posterior ensemble mean not found'
- call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
+ write(string1,*)'metadata:"posterior ensemble mean" not found'
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if ( prior_spread_index < 0 ) then
- write(string1,*)'metadata:prior ensemble spread not found'
- call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
+ write(string1,*)'metadata:"prior ensemble spread" not found'
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if ( posterior_spread_index < 0 ) then
- write(string1,*)'metadata:posterior ensemble spread not found'
- call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
+ write(string1,*)'metadata:"posterior ensemble spread" not found'
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if ( org_qc_index < 0 ) then
- write(string1,*)'metadata:Quality Control not found'
- call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
+ write(string1,*)'metadata:"Quality Control" not found'
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if ( dart_qc_index < 0 ) then
- write(string1,*)'metadata:DART quality control not found'
- call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
+ write(string1,*)'metadata:"DART quality control" not found'
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
! Only require obs_index to be present; this allows the program
-! to be run on obs_seq.in files which have no means or spread.
+! to be run on obs_seq.[in,out] files which have no means or spread.
+! Can still count number of obs, observation mean, ...
if ( obs_index < 0 ) then
if ( use_zero_error_obs ) then
- write(string1,*)'metadata:truth not found'
+ write(string1,*)'metadata:"truth" not found'
else
- write(string1,*)'metadata:observation not found'
+ write(string1,*)'metadata:"observation" not found'
endif
- call error_handler(E_ERR,'obs_diag',string1,source,revision,revdate)
+ call error_handler(E_ERR,'SetIndices',string1,source,revision,revdate)
endif
!--------------------------------------------------------------------
! Echo what we found.
!--------------------------------------------------------------------
+
if ( use_zero_error_obs ) then
- write(string1,'(''truth index '',i2,'' metadata '',a)') &
+ write(string1,'(''"truth" index '',i2,'' metadata '',a)') &
obs_index, trim(get_copy_meta_data(seq,obs_index))
else
- write(string1,'(''observation index '',i2,'' metadata '',a)') &
+ write(string1,'(''"observation" index '',i2,'' metadata '',a)') &
obs_index, trim(get_copy_meta_data(seq,obs_index))
endif
+call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
-call error_handler(E_MSG,'obs_diag',string1,source,revision,revdate)
-
if (prior_mean_index > 0 ) then
- write(string1,'(''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',string1,source,revision,revdate)
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if (posterior_mean_index > 0 ) then
- write(string1,'(''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',string1,source,revision,revdate)
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if (prior_spread_index > 0 ) then
- write(string1,'(''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',string1,source,revision,revdate)
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if (posterior_spread_index > 0 ) then
- write(string1,'(''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',string1,source,revision,revdate)
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if (org_qc_index > 0 ) then
- write(string1,'(''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',string1,source,revision,revdate)
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
if (dart_qc_index > 0 ) then
- write(string1,'(''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',string1,source,revision,revdate)
+ call error_handler(E_MSG,'SetIndices',string1,source,revision,revdate)
endif
-end Subroutine SetIndices
+end subroutine SetIndices
!======================================================================
-Function is_observation_trusted(obsname)
+function is_observation_trusted(obsname)
! Is the observation one that we 'trust'.
! If so, disregard the DART QC ==7 (outlier rejection) and use it to
@@ -2496,13 +2509,13 @@
endif
enddo rUtrusted
-end Function is_observation_trusted
+end function is_observation_trusted
!======================================================================
-Subroutine CheckVertical(obslocation, flavor)
+subroutine CheckVertical(obslocation, flavor)
! It is not allowed to have one observation flavor exist on multiple
! types of vertical coordinate systems. If the flavor has been
@@ -2562,13 +2575,13 @@
endif
-end Subroutine CheckVertical
+end subroutine CheckVertical
!======================================================================
-Function ParseLevel(obslocation, obslevel, flavor)
+function ParseLevel(obslocation, obslevel, flavor)
! returns the index of the level closest to 'obslevel'
! There are three independent 'level' arrays
type(location_type), intent(in ) :: obslocation
@@ -2629,7 +2642,7 @@
which_vert(flavor) = VERTISUNDEF
else
- call error_handler(E_ERR,'obs_diag','Vertical coordinate not recognized', &
+ call error_handler(E_ERR,'ParseLevel','Vertical coordinate not recognized', &
source,revision,revdate)
endif
@@ -2638,13 +2651,13 @@
300 format('obskey ',i7,' # ',i7,' [',3(f10.4,1x),'] mid == ',f10.4,' level ',i3,1x,a)
301 format('obskey ',i7,' # ',i7,' [',3(f10.4,1x),'] mid == ',i10 ,' level ',i3,1x,a)
-end Function ParseLevel
+end function ParseLevel
!======================================================================
-Function InnovZscore(obsval, prmean, prspred, errvar, qcval, qcmaxval)
+function InnovZscore(obsval, prmean, prspred, errvar, qcval, qcmaxval)
! This function tries to get a handle on the magnitude of the innovations.
! If the ratio of the observation to the prior mean is 'big', it is an outlier.
@@ -2683,13 +2696,13 @@
endif
-end Function InnovZscore
+end function InnovZscore
!======================================================================
-Function Rank_Histogram(copyvalues, obs_index, &
+function Rank_Histogram(copyvalues, obs_index, &
error_variance, ens_copy_index ) result(rank)
! Calculates the bin/rank
@@ -2758,13 +2771,13 @@
write(*,*)
endif
-end Function Rank_Histogram
+end function Rank_Histogram
!======================================================================
-Function CheckMate(vflavor, uflavor, obsloc1, obsloc2, flavor )
+function CheckMate(vflavor, uflavor, obsloc1, obsloc2, flavor )
! This routine ensures that the U,V components of wind
! are from the same observation location so we can convert
@@ -2864,7 +2877,7 @@
if (index(vobs_string, uobs_string(1:indx1)) < 1) then
write(string1,*) 'around OBS ', keys(obsindex), 'observation types not compatible.'
- call error_handler(E_WARN,'obs_diag',string1,source,revision,revdate, &
+ call error_handler(E_WARN,'CheckMate',string1,source,revision,revdate, &
text2=vobs_string, text3=uobs_string)
endif
@@ -2883,28 +2896,37 @@
if (CheckMate /= 0) then
write(string1,*) 'around OBS ', keys(obsindex), 'observation types not known.'
- call error_handler(E_ERR,'obs_diag',string1,source,revision,revdate,&
+ call error_handler(E_ERR,'CheckMate',string1,source,revision,revdate,&
text2=vobs_string, text3=uobs_string)
endif
-end Function CheckMate
+end function CheckMate
!======================================================================
-Subroutine CountDartQC_4D(myqc, iepoch, ilevel, iregion, itype, prior, poste, &
- posterior_mean)
+subroutine CountDartQC_4D(inqc, iepoch, ilevel, iregion, itype, prior, poste, &
+ posterior_mean, uqc)
-integer, intent(in) :: myqc
+integer, intent(in) :: inqc
integer, intent(in) :: iepoch
integer, intent(in) :: ilevel
integer, intent(in) :: iregion
integer, intent(in) :: itype
type(TLRV_type), intent(inout) :: prior
type(TLRV_type), intent(inout) :: poste
-real(r8), intent(in) :: posterior_mean
+real(r8), intent(in), optional :: posterior_mean
+integer, intent(in), optional :: uqc
+integer :: myqc
+
+if (present(uqc)) then
+ myqc = maxval((/ inqc, uqc /))
+else
+ myqc = inqc
+endif
+
if ( myqc == 0 ) then
call IPE(prior%NDartQC_0(iepoch,ilevel,iregion,itype), 1)
call IPE(poste%NDartQC_0(iepoch,ilevel,iregion,itype), 1)
@@ -2936,31 +2958,42 @@
elseif ( myqc == 7 ) then
call IPE(prior%NDartQC_7(iepoch,ilevel,iregion,itype), 1)
- if ( abs(posterior_mean - MISSING_R8) < 1.0_r8 ) then
- ! ACTUALLY A FAILED FORWARD OPERATOR - ambiguous case
- call IPE(poste%NDartQC_4(iepoch,ilevel,iregion,itype), 1)
- else
- call IPE(poste%NDartQC_7(iepoch,ilevel,iregion,itype), 1)
+ if (present(posterior_mean)) then
+ if ( abs(posterior_mean - MISSING_R8) < 1.0_r8 ) then
+ ! ACTUALLY A FAILED FORWARD OPERATOR - ambiguous case
+ call IPE(poste%NDartQC_4(iepoch,ilevel,iregion,itype), 1)
+ else
+ call IPE(poste%NDartQC_7(iepoch,ilevel,iregion,itype), 1)
+ endif
endif
endif
-end Subroutine CountDartQC_4D
+end subroutine CountDartQC_4D
!======================================================================
-Subroutine CountDartQC_3D(myqc, ilevel, iregion, itype, prior, poste, posterior_mean)
+subroutine CountDartQC_3D(inqc, ilevel, iregion, itype, prior, poste, posterior_mean, uqc)
-integer, intent(in) :: myqc
+integer, intent(in) :: inqc
integer, intent(in) :: ilevel
integer, intent(in) :: iregion
integer, intent(in) :: itype
type(LRV_type), intent(inout) :: prior
type(LRV_type), intent(inout) :: poste
-real(r8), intent(in) :: posterior_mean
+real(r8), intent(in), optional :: posterior_mean
+integer, intent(in), optional :: uqc
+integer :: myqc
+
+if (present(uqc)) then
+ myqc = maxval((/ inqc, uqc /))
+else
+ myqc = inqc
+endif
+
if ( myqc == 0 ) then
call IPE(prior%NDartQC_0(ilevel,iregion,itype), 1)
call IPE(poste%NDartQC_0(ilevel,iregion,itype), 1)
@@ -2992,22 +3025,24 @@
elseif ( myqc == 7 ) then
call IPE(prior%NDartQC_7(ilevel,iregion,itype), 1)
- if ( abs(posterior_mean - MISSING_R8) < 1.0_r8 ) then
- ! ACTUALLY A FAILED FORWARD OPERATOR - ambiguous case
- call IPE(poste%NDartQC_4(ilevel,iregion,itype), 1)
- else
- call IPE(poste%NDartQC_7(ilevel,iregion,itype), 1)
+ if (present(posterior_mean)) then
+ if ( abs(posterior_mean - MISSING_R8) < 1.0_r8 ) then
+ ! ACTUALLY A FAILED FORWARD OPERATOR - ambiguous case
+ call IPE(poste%NDartQC_4(ilevel,iregion,itype), 1)
+ else
+ call IPE(poste%NDartQC_7(ilevel,iregion,itype), 1)
+ endif
endif
endif
-end Subroutine CountDartQC_3D
+end subroutine CountDartQC_3D
!======================================================================
-Subroutine Bin4D(iqc, iepoch, ilevel, iregion, flavor, trusted, &
+subroutine Bin4D(iqc, iepoch, ilevel, iregion, flavor, trusted, &
obsval, obserrvar, prmean, prsprd, pomean, posprd, rank, &
uobs, uobserrvar, uprmean, uprsprd, upomean, uposprd, uqc)
!----------------------------------------------------------------------
@@ -3100,7 +3135,7 @@
! counted.
myrank = -99
-else if ( any(optionals) ) then
+elseif ( any(optionals) ) then
call error_handler(E_ERR,'Bin4D','wrong number of optional arguments', &
source,revision,revdate)
else
@@ -3175,13 +3210,13 @@
call IPE(poste%NbadDartQC(iepoch,ilevel,iregion,flavor), 1 )
endif
-end Subroutine Bin4D
+end subroutine Bin4D
!======================================================================
-Subroutine Bin3D(iqc, ilevel, iregion, flavor, trusted, &
+subroutine Bin3D(iqc, ilevel, iregion, flavor, trusted, &
obsval, obserrvar, prmean, prsprd, pomean, posprd, &
uobs, uobserrvar, uprmean, uprsprd, upomean, uposprd, uqc )
!----------------------------------------------------------------------
@@ -3255,7 +3290,7 @@
priorspredplus = prsprd**2 + obserrvar + uprsprd**2 + uobserrvar
postspredplus = posprd**2 + obserrvar + uposprd**2 + uobserrvar
-else if ( any(optionals) ) then
+elseif ( any(optionals) ) then
call error_handler(E_ERR,'Bin3D','wrong number of optional arguments', &
source,revision,revdate)
else
@@ -3316,20 +3351,16 @@
call IPE(posteAVG%NbadDartQC(ilevel,iregion,flavor), 1 )
endif
-end Subroutine Bin3D
+end subroutine Bin3D
!======================================================================
-Subroutine Normalize4Dvars()
+subroutine Normalize4Dvars()
-if (verbose) then
- do ivar = 1,SIZE(which_vert)
- write(logfileunit,*)'which_vert(',ivar,' of ',num_obs_types,') = ',which_vert(ivar)
- write( * ,*)'which_vert(',ivar,' of ',num_obs_types,') = ',which_vert(ivar)
- enddo
-endif
+! The vertical coordinate system definition is summarized in WriteNetCDF
+! if the verbose option is chosen.
if (verbose) then
write(logfileunit,*)'Normalizing time-level-region-variable quantities.'
@@ -3426,13 +3457,13 @@
enddo
call close_file(nsigmaUnit)
-end Subroutine Normalize4Dvars
+end subroutine Normalize4Dvars
!======================================================================
-Subroutine Normalize3Dvars()
+subroutine Normalize3Dvars()
if (verbose) then
write(logfileunit,*)'Normalize quantities for all levels.'
write( * ,*)'Normalize quantities for all levels.'
@@ -3515,13 +3546,13 @@
enddo
enddo
-end Subroutine Normalize3Dvars
+end subroutine Normalize3Dvars
!======================================================================
-Subroutine WriteNetCDF(fname)
+subroutine WriteNetCDF(fname)
character(len=*), intent(in) :: fname
integer :: ncid, i, indx1, nobs, typesdimlen
@@ -3723,61 +3754,61 @@
!----------------------------------------------------------------------------
call nc_check(nf90_def_dim(ncid=ncid, &
- name='time', len = NF90_UNLIMITED, dimid = TimeDimID), &
- 'WriteNetCDF', 'time:def_dim '//trim(fname))
+ name='time', len=NF90_UNLIMITED, dimid=TimeDimID), &
+ 'WriteNetCDF', 'time:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='bounds', len = 2, dimid = BoundsDimID), &
- 'WriteNetCDF', 'bounds:def_dim '//trim(fname))
+ name='bounds', len=2, dimid=BoundsDimID), &
+ 'WriteNetCDF', 'bounds:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='copy', len = Ncopies, dimid = CopyDimID), &
- 'WriteNetCDF', 'copy:def_dim '//trim(fname))
+ name='copy', len=Ncopies, dimid=CopyDimID), &
+ 'WriteNetCDF', 'copy:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='obstypes', len = max_obs_kinds, dimid = TypesDimID), &
- 'WriteNetCDF', 'types:def_dim '//trim(fname))
+ name='obstypes', len=max_obs_kinds, dimid=TypesDimID), &
+ 'WriteNetCDF', 'types:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='region', len = Nregions, dimid = RegionDimID), &
- 'WriteNetCDF', 'region:def_dim '//trim(fname))
+ name='region', len=Nregions, dimid=RegionDimID), &
+ 'WriteNetCDF', 'region:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='surface', len = 1, dimid = SlevelDimID), &
- 'WriteNetCDF', 'slevel:def_dim '//trim(fname))
+ name='surface', len=1, dimid=SlevelDimID), &
+ 'WriteNetCDF', 'slevel:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='undef', len = 1, dimid = UlevelDimID), &
- 'WriteNetCDF', 'ulevel:def_dim '//trim(fname))
+ name='undef', len=1, dimid=UlevelDimID), &
+ 'WriteNetCDF', 'ulevel:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='mlevel', len = Nmlevels, dimid = MlevelDimID), &
- 'WriteNetCDF', 'mlevel:def_dim '//trim(fname))
+ name='mlevel', len=Nmlevels, dimid=MlevelDimID), &
+ 'WriteNetCDF', 'mlevel:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='mlevel_edges', len = Nmlevels+1, dimid = MlevIntDimID), &
- 'WriteNetCDF', 'mlevel_edges:def_dim '//trim(fname))
+ name='mlevel_edges', len=Nmlevels+1, dimid=MlevIntDimID), &
+ 'WriteNetCDF', 'mlevel_edges:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='plevel', len = Nplevels, dimid = PlevelDimID), &
- 'WriteNetCDF', 'plevel:def_dim '//trim(fname))
+ name='plevel', len=Nplevels, dimid=PlevelDimID), &
+ 'WriteNetCDF', 'plevel:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='plevel_edges', len = Nplevels+1, dimid = PlevIntDimID), &
- 'WriteNetCDF', 'plevel_edges:def_dim '//trim(fname))
+ name='plevel_edges', len=Nplevels+1, dimid=PlevIntDimID), &
+ 'WriteNetCDF', 'plevel_edges:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='hlevel', len = Nhlevels, dimid = HlevelDimID), &
- 'WriteNetCDF', 'hlevel:def_dim '//trim(fname))
+ name='hlevel', len=Nhlevels, dimid=HlevelDimID), &
+ 'WriteNetCDF', 'hlevel:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='hlevel_edges', len = Nhlevels+1, dimid = HlevIntDimID), &
- 'WriteNetCDF', 'hlevel_edges:def_dim '//trim(fname))
+ name='hlevel_edges', len=Nhlevels+1, dimid=HlevIntDimID), &
+ 'WriteNetCDF', 'hlevel_edges:def_dim '//trim(fname))
call nc_check(nf90_def_dim(ncid=ncid, &
- name='stringlength', len = stringlength, dimid = StringDimID), &
- 'WriteNetCDF', 'stringlength:def_dim '//trim(fname))
+ name='stringlength', len=stringlength, dimid=StringDimID), &
+ 'WriteNetCDF', 'stringlength:def_dim '//trim(fname))
if (create_rank_histogram) then
call nc_check(nf90_def_dim(ncid=ncid, &
- name='rank_bins', len = ens_size+1, dimid = RankDimID), &
- 'WriteNetCDF', 'rank_bins:def_dim '//trim(fname))
+ name='rank_bins', len=ens_size+1, dimid=RankDimID), &
+ 'WriteNetCDF', 'rank_bins:def_dim '//trim(fname))
endif
!----------------------------------------------------------------------------
@@ -3882,7 +3913,7 @@
!----------------------------------------------------------------------------
call nc_check(nf90_def_var(ncid=ncid, name='hlevel', xtype=nf90_real, &
- dimids=HlevelDimID, varid=HlevelVarID), 'WriteNetCDF', 'hlevel:def_var')
+ dimids=HlevelDimID, varid=HlevelVarID), 'WriteNetCDF', 'hlevel:def_var')
call nc_check(nf90_put_att(ncid, HlevelVarID, 'long_name', 'height bin midpoints'), &
'WriteNetCDF', 'hlevel:long_name')
call nc_check(nf90_put_att(ncid, HlevelVarID, 'units', 'm'), &
@@ -3890,7 +3921,7 @@
call nc_check(nf90_put_att(ncid, HlevelVarID, 'axis', 'Z'), &
'WriteNetCDF', 'hlevel:axis')
call nc_check(nf90_put_att(ncid, HlevelVarID, 'valid_range', &
- (/ minval(hlevel(1:Nhlevels)), maxval(hlevel(1:Nhlevels)) /)), &
+ (/ minval(hlevel(1:Nhlevels)), maxval(hlevel(1:Nhlevels)) /)), &
'WriteNetCDF', 'hlevel:valid_range')
!----------------------------------------------------------------------------
@@ -3898,7 +3929,7 @@
!----------------------------------------------------------------------------
call nc_check(nf90_def_var(ncid=ncid, name='hlevel_edges', xtype=nf90_real, &
- dimids=HlevIntDimID, varid=HlevIntVarID), 'WriteNetCDF', 'hlevel_edges:def_var')
+ dimids=HlevIntDimID, varid=HlevIntVarID), 'WriteNetCDF', 'hlevel_edges:def_var')
call nc_check(nf90_put_att(ncid, HlevIntVarID, 'long_name', 'height bin edges'), &
'WriteNetCDF', 'hlevel_edges:long_name')
call nc_check(nf90_put_att(ncid, HlevIntVarID, 'units', 'm'), &
@@ -3979,7 +4010,7 @@
call nc_check(nf90_put_att(ncid, TypesMetaVarID, 'long_name', 'DART observation types'), &
'WriteNetCDF', 'typesmeta:long_name')
call nc_check(nf90_put_att(ncid, TypesMetaVarID, 'comment', &
- 'table relating integer to observation type string'), &
+ 'table relating integer to observation type string'), &
'WriteNetCDF', 'typesmeta:comment')
if ( create_rank_histogram ) then
@@ -3989,7 +4020,7 @@
call nc_check(nf90_put_att(ncid, RankVarID, 'long_name', 'rank histogram bins'), &
'WriteNetCDF', 'rank_bins:long_name')
call nc_check(nf90_put_att(ncid, RankVarID, 'comment', &
- 'position of the observation among the sorted noisy ensemble members'), &
+ 'position of the observation among the sorted noisy ensemble members'), &
'WriteNetCDF', 'rank_bins:comment')
endif
@@ -3998,7 +4029,7 @@
!----------------------------------------------------------------------------
call nc_check(nf90_set_fill(ncid, NF90_NOFILL, i), &
- 'WriteNetCDF', 'set_nofill '//trim(fname))
+ 'WriteNetCDF', 'set_nofill '//trim(fname))
!----------------------------------------------------------------------------
! Leave define mode so we can fill
@@ -4012,43 +4043,43 @@
!----------------------------------------------------------------------------
call nc_check(nf90_put_var(ncid, CopyVarId, (/ (i,i=1,Ncopies) /) ), &
- 'WriteNetCDF', 'copy:put_var')
+ 'WriteNetCDF', 'copy:put_var')
call nc_check(nf90_put_var(ncid, CopyMetaVarID, copy_names), &
- 'WriteNetCDF', 'copymeta:put_var')
+ 'WriteNetCDF', 'copymeta:put_var')
call nc_check(nf90_put_var(ncid, TypesVarId, (/ (i,i=1,max_obs_kinds) /) ), &
- 'WriteNetCDF', 'types:put_var')
+ 'WriteNetCDF', 'types:put_var')
call nc_check(nf90_put_var(ncid, TypesMetaVarID, obs_type_strings(1:max_obs_kinds)), &
- 'WriteNetCDF', 'typesmeta:put_var')
+ 'WriteNetCDF', 'typesmeta:put_var')
call nc_check(nf90_put_var(ncid, RegionVarID, (/ (i,i=1,Nregions) /) ), &
- 'WriteNetCDF', 'region:put_var')
+ 'WriteNetCDF', 'region:put_var')
call nc_check(nf90_put_var(ncid, MlevelVarID, mlevel(1:Nmlevels) ), &
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list