[Dart-dev] [10353] DART: Bringing up-to-date with respect to the trunk:
nancy at ucar.edu
nancy at ucar.edu
Fri Jun 10 14:37:47 MDT 2016
Revision: 10353
Author: thoar
Date: 2016-06-10 14:37:47 -0600 (Fri, 10 Jun 2016)
Log Message:
-----------
Bringing up-to-date with respect to the trunk:
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.
Removed some debug code that was causing problems (i.e. a fatal error) when people
input the _edges arrays as opposed to the midpoints to define the vertical binning.
When obs_diag supported the use of the _edges, I had code in to check to make sure
the new algorithm generated the correct answer under all circumstances. The old
algorithm did not work when the user specified the _edges. However, the test that
was inadvertently left in was to check to make sure the old level index matched
the new level index (a good thing if you specify the midpoints,
a bad thing if you specify the edges).
The new algorithm has been working for years now, time to remove the debug code.
Modified Paths:
--------------
DART/branches/FEOM/diagnostics/threed_sphere/obs_diag.f90
DART/branches/UCOAM/diagnostics/threed_sphere/obs_diag.f90
DART/branches/rma_openggcm/diagnostics/threed_sphere/obs_diag.f90
DART/branches/rma_roms/diagnostics/threed_sphere/obs_diag.f90
DART/branches/rma_trunk/diagnostics/threed_cartesian/obs_diag.f90
DART/branches/rma_trunk/diagnostics/threed_sphere/obs_diag.f90
DART/branches/tiegcm/diagnostics/threed_sphere/obs_diag.f90
DART/releases/Lanai/diagnostics/threed_sphere/obs_diag.f90
DART/trunk/diagnostics/threed_sphere/obs_diag.f90
-------------- next part --------------
Modified: DART/branches/FEOM/diagnostics/threed_sphere/obs_diag.f90
===================================================================
--- DART/branches/FEOM/diagnostics/threed_sphere/obs_diag.f90 2016-06-10 17:51:47 UTC (rev 10352)
+++ DART/branches/FEOM/diagnostics/threed_sphere/obs_diag.f90 2016-06-10 20:37:47 UTC (rev 10353)
@@ -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 ', &
@@ -573,7 +573,7 @@
endif
!--------------------------------------------------------------------
- ! Prepare some variables for the rank histogram.
+ ! 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) ...
@@ -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
@@ -799,15 +805,15 @@
endif
!--------------------------------------------------------------
- ! There is a ambiguous case wherein the prior is rejected (DART QC ==7)
- ! and the posterior forward operator fails (DART QC ==4). In this case,
- ! the DART_QC only reflects the fact the prior was rejected - HOWEVER -
- ! the posterior mean,spread are set to MISSING.
+ ! There is a ambiguous case wherein the prior is rejected (DART QC ==7)
+ ! and the posterior forward operator fails (DART QC ==4). In this case,
+ ! the DART_QC only reflects the fact the prior was rejected - HOWEVER -
+ ! the posterior mean,spread are set to MISSING.
!
! If it is your intent to compare identical prior and posterior TRUSTED
! observations, then you should enable the following few lines of code.
! and realize that the number of observations rejected because of the
- ! outlier threshold will be wrong.
+ ! outlier threshold will be wrong.
!
! This is the only block of code you should need to change.
!--------------------------------------------------------------
@@ -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
-
- 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(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(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
@@ -1404,7 +1410,7 @@
! bin_separation comes from namelist
! bin_width comes from namelist
! time_to_skip comes from namelist
-!
+!
! Global-scope variables set in this routine:
! type(time_type), intent(out) :: beg_time ! first_bin_center
! type(time_type), intent(out) :: end_time ! last_bin_center
@@ -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( binwidth,'requested bin width',logfileunit)
- 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',logfileunit)
call print_time( binwidth,'requested bin width')
+ call print_time(halfbinwidth,'implied halfbinwidth',logfileunit)
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.
@@ -1994,20 +2003,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 +2230,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 +2250,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 +2289,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 +2317,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 +2393,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 +2504,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 +2570,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
@@ -2588,7 +2596,7 @@
! at one point, negative (model) levels were used to indicate
! surface observations - ancient history ...
if (obslevel < 1.0_r8 ) obslevel = 1.0_r8 ! TJH negative model level
- ParseLevel = ClosestLevelIndex(obslevel, VERTISLEVEL)
+ ParseLevel = find_my_level(obslevel, VERTISLEVEL)
which_vert(flavor) = VERTISLEVEL
if ( (1 == 2) .and. (ParseLevel > 0) .and. (ParseLevel < MaxLevels)) &
@@ -2597,7 +2605,7 @@
mlevel(ParseLevel), ParseLevel, 'level'
elseif(vert_is_pressure(obslocation)) then
- ParseLevel = ClosestLevelIndex(obslevel, VERTISPRESSURE)
+ ParseLevel = find_my_level(obslevel, VERTISPRESSURE)
which_vert(flavor) = VERTISPRESSURE
if ( (1 == 2) .and. (ParseLevel > 0) .and. (ParseLevel < MaxLevels)) &
@@ -2606,7 +2614,7 @@
plevel(ParseLevel), ParseLevel, 'pressure'
elseif(vert_is_height(obslocation)) then
- ParseLevel = ClosestLevelIndex(obslevel, VERTISHEIGHT)
+ ParseLevel = find_my_level(obslevel, VERTISHEIGHT)
which_vert(flavor) = VERTISHEIGHT
if ( (1 == 2) .and. (ParseLevel > 0) .and. (ParseLevel < MaxLevels)) &
@@ -2615,7 +2623,7 @@
hlevel(ParseLevel), ParseLevel, 'height'
elseif(vert_is_scale_height(obslocation)) then
- ParseLevel = ClosestLevelIndex(obslevel, VERTISSCALEHEIGHT)
+ ParseLevel = find_my_level(obslevel, VERTISSCALEHEIGHT)
which_vert(flavor) = VERTISSCALEHEIGHT
if ( (1 == 2) .and. (ParseLevel > 0) .and. (ParseLevel < MaxLevels)) &
@@ -2629,7 +2637,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 +2646,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.
@@ -2654,11 +2662,11 @@
!
! In Jan of 2014 I ran into a special case. We are performing a perfect model
! experiment - perturbing a single state and then assimilating. We wanted
-! to compare against the true observation value (errvar = 0.0) and for the
+! to compare against the true observation value (errvar = 0.0) and for the
! first time step, there is no ensemble spread either. In this case the denom
! was really zero and the calculation blew up. Since we only use this to track
! how far apart these things are, we can restrict the distance to the worst-case
-! scenario ... the last bin.
+! scenario ... the last bin.
real(r8) :: InnovZscore
real(r8), intent(in) :: obsval, prmean, prspred, errvar
@@ -2683,13 +2691,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 +2766,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 +2872,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 +2891,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 +2953,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 +3020,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)
!----------------------------------------------------------------------
@@ -3054,9 +3084,9 @@
posterior_qc = iqc
! There is a ambiguous case wherein the prior is rejected (DART QC ==7)
-! and the posterior forward operator fails (DART QC ==4). In this case,
-! the DART_QC reflects the fact the prior was rejected - HOWEVER -
-! the posterior mean,spread are set to MISSING.
+! and the posterior forward operator fails (DART QC ==4). In this case,
+! the DART_QC reflects the fact the prior was rejected - HOWEVER -
+! the posterior mean,spread are set to MISSING.
if ((prior_qc == 7) .and. (abs(pomean - MISSING_R8) > 1.0_r8)) posterior_qc = 4
@@ -3100,7 +3130,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
@@ -3139,7 +3169,7 @@
call IPE(poste%Nposs(iepoch,ilevel,iregion,flavor), 1)
!----------------------------------------------------------------------
-! Select which set of qcs are valid and accrue everything
+! Select which set of qcs are valid and accrue everything
!----------------------------------------------------------------------
if ( trusted ) then
@@ -3175,13 +3205,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 )
!----------------------------------------------------------------------
@@ -3217,9 +3247,9 @@
posterior_qc = iqc
! There is a ambiguous case wherein the prior is rejected (DART QC ==7)
-! and the posterior forward operator fails (DART QC ==4). In this case,
-! the DART_QC reflects the fact the prior was rejected - HOWEVER -
-! the posterior mean,spread are set to MISSING.
+! and the posterior forward operator fails (DART QC ==4). In this case,
+! the DART_QC reflects the fact the prior was rejected - HOWEVER -
+! the posterior mean,spread are set to MISSING.
if ((prior_qc == 7) .and. (abs(pomean - MISSING_R8) > 1.0_r8)) posterior_qc = 4
@@ -3255,7 +3285,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
@@ -3280,7 +3310,7 @@
call IPE(posteAVG%Nposs(ilevel,iregion,flavor), 1)
!----------------------------------------------------------------------
-! Select which set of qcs are valid and accrue everything
+! Select which set of qcs are valid and accrue everything
!----------------------------------------------------------------------
if ( trusted ) then
@@ -3301,7 +3331,7 @@
else
call IPE(priorAVG%NbadDartQC(ilevel,iregion,flavor), 1 )
endif
-
+
! Accrue the POSTERIOR quantities
if (( trusted .and. any(trusted_poste_qcs == posterior_qc)) .or. &
(.not. trusted .and. any( good_poste_qcs == posterior_qc))) then
@@ -3316,20 +3346,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 +3452,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 +3541,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
@@ -3702,10 +3728,10 @@
string1 = obs_type_strings(ivar)
endif
- call nc_check(nf90_put_att(ncid, NF90_GLOBAL, string1, ivar), &
+ call nc_check(nf90_put_att(ncid, NF90_GLOBAL, string1, ivar), &
'WriteNetCDF', 'put_att:obs_type_string '//trim(obs_type_strings(ivar)))
endif
-
+
enddo
if (typesdimlen < 1) then
@@ -3723,61 +3749,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), &
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list