PageOutline

• [MichaelNaughton 28/05/2020] This page is a copy of Neill Bowler's ETKF Documentation wiki page from MOSRS site ETKFDocumentation.

ETKF Documentation

This page describes the Ensemble Transform Kalman Filter (ETKF) program, EnsProg_ETKF, which runs within the MOGREPS-G ensemble suites to generate the initial condition perturbations. This version of the documentation describes the ETKF valid at PS39, which corresponds to ops:browser/main/trunk/src/code@5685.

Source code layout

The ETKF source code is located at ops:browser/main/trunk/src/code. This contains the following ETKF-specific subdirectories:

• EnsProg_ETKF - contains the ETKF-specific routines which constitute the bulk of the code.
• EnsProg_TrimObstore - contains the source code for the TrimObstore program, which shares some routines with the ETKF.
• EnsMod_Obstore - contains routines for reading and writing obstore files, which are used by TrimObstore but not the ETKF.
• EnsMod_Varobs - contains routines and data used by the ETKF to read varobs files.
• EnsMod_Header - contains routines and data used by the ETKF to read and write fieldsfile headers.
• EnsMod_Utilities - contains some general purpose routines, of which only Ens_RoundUp appears to be used by the ETKF.

Building the ETKF

The ETKF build process is achieved using the standard rose-stem approach. Since the code is embedded within OPS, the standard approach is to build OPS and the ETKF and TrimObstore programs will be produced as a consequence of this. To build the ETKF check-out the appropriate version of the code, and then run:

rose stem --task=fcm_make

Running the ETKF

The ETKF runs on the supercomputer as a job submitted to the queuing system. Its operation and the files it reads are controlled by a number of environment variables:

• ETKF_CONTROL_DIR - file-path for the namelist file which provides various options to the ETKF.
• ETKF_PERT_INPUT_TEMPLATE - template for the file-path to the files to be used as input. Any occurrences of the string MEM will be replaced by the ensemble member number (three digits). So, this make look like /some/path/to/suite/share/cycle/20170724T0000Z/engl_um_MEM/englaa_pb000
• ETKF_PERT_OUTPUT_TEMPLATE - template for the file-path to the output file locations. Any occurrences of the string MEM will be replaced by the ensemble member number (three digits). So, this make look like /some/path/to/suite/share/cycle/20170724T0000Z/engl_etkf/perts/engl_pert_MEM
• ETKF_VAROBS_TEMPLATE - template for the file-path to the varobs files to be read in. Any occurrences of the string MEM will be replaced by the ensemble member number (three digits). So, this make look like /some/path/to/suite/share/cycle/20170724T0000Z/engl_processobs/ProcessMEM.varobs
• ETKF_TEXT_INPUT_DIR - the directory from which the ETKF reads inflation factor files.
• ETKF_TEXT_OUTPUT_DIR - the directory to which the ETKF writes inflation factor files and the Innovations file.
• ETKF_LOCALISATION_DIR - the directory from which the ETKF reads localisation centre information. The same files can be found at ops:browser/main/trunk/src/control/Localisation.
• ETKF_RECENTRE_MEMBERS - (optional) a subset of ensemble members to recentre around the deterministic analysis. Can be used to ensure that the subset of ensemble members which are being used to produce products have zero mean perturbation.
• ETKF_RECENTRE_SEPARATE - if using recentring, then whether to write the recentred ensemble members to separate files (in addition to producing files which has not had the separate recentring).
• ETKF_RECENTRE_TEMPLATE - if using separate recentring, then the template for the separate recentred files.
• OPS_SATRAD_Coeffs - the OPS directory from which SatRad channel selection namelists and VAR observation error estimates are read.

Other variables which may be of interest include:

• GEN_PRODUCE_HTML- set this to False to remove HTML tags from the standard output/error/stats streams, which can make them easier to read when not doing so through a web browser.
• GEN_MODE - controls the degree of diagnostic messaging (10/20/30/40 for Normal/Diagnostic/Debug/Verbose)
• GEN_TRACEUSE - controls the production of GEN execution traces.

Within rose suites, EnsProg_ETKF is run via the app engl_ens_etkf. The wrapper script which runs the code performs the following actions:

• Provide some default settings for GEN-managed diagnostics and tracing.
• Check that required environment variables have been supplied.
• Check that required input files have been supplied.
• Remove any existing perturbation files, to avoid ambiguity about whether the ETKF completes successfully or not.
• Run EnsProg_ETKF.
• Check that the expected perturbation files have been produced.

The wrapper script requires the following environment variables, in addition to those required by EnsProg_ETKF:

• ENS_MEMBERS - space-separated list of the ensemble members for which input and output files are expected.

Overview of the ETKF

The overall purpose of EnsProg_ETKF is to provide the initial conditions for the MOGREPS-G ensemble forecasts. It does this by producing global perturbation fields for wind, temperature, humidity and pressure on all model levels for NENS (currently 44) ensemble members. These are combined with the reconfigured 4DVAR analysis using the Incremental Analysis Update (IAU) mechanism within the Unified Model. The control forecast (member 0) uses the reconfigured analysis without perturbations, and hence receives no input from the ETKF.

The Ensemble Transform Kalman Filter (Wang and Bishop 2003) generates the analysis perturbations by mixing and scaling evolved perturbations (valid at the new analysis time, in our case T+6 hours) from the previous forecast cycle. EnsProg_ETKF obtains the evolved perturbations from the T+6 state of the perturbed ensemble members, provided by the pp1 fieldsfiles which it receives as input. The transform matrix which performs the mixing is calculated to produce perturbations which provide an NENS-dimensional representation of the analysis error covariance matrix of an optimal data assimilation system, given the background error covariance represented by the evolved perturbations from the previous forecast, together with the error estimates attached to each assimilated observation. This calculation requires the model equivalent of each observation for each ensemble member, to provide the estimates of background uncertainty in observation space. These 'pseudo-observations' are calculated by the Observation Processing System (OPS) as a separate step in the forecast suite, and provided to the ETKF through modelobs files beneath $ETKF_VAROBS_TEMPLATE.

The basic ETKF produces perturbations which are unrealistically small, because the small number of ensemble members leads to unrealistically large background error covariance estimates, which in turn lead to overestimation of the impact of each observation. This problem is reduced (and the relationship between spread and error as a function of latitude improved) by the use of horizontal localisation (Houtekamer and Mitchell 1998). A number of localisation centres (currently 92) are defined, equally spaced across the globe. For each centre, observations within a certain radius (now 2000 km, similar to the centre separation) are used to produce a local transform matrix. The final transform matrix for each point on the model grid is obtained by interpolation between the transform matrices for the nearest localisation centres. This prevents any observations acting through any (generally spurious) long-range correlations which may be present in the ensemble-estimated background error covariance. It also improves the rank of the analysis error covariance estimate, since there is no longer just one transform matrix being used for the whole globe.

To correct the remaining deficiency in spread, the raw transform matrix for each region is multiplied a region-specific 'inflation factor'. This inflation factor is updated by the ratio of the rms error of the ensemble mean to the rms spread, as measured against observations at T+6 hours into the previous forecast. The ensemble mean and spread in observation space are already available from the modelobs files, and the OPS is made to write the corresponding real observations (processed against the control forecast for model level position, filtering, background checks, etc) to varobs files in $ETKF_VAROBS_TEMPLATE. This mechanism assumes that the required inflation factor is reasonably stable between forecasts, since the error/spread ratio from the previous forecast is used to estimate the inflation factor which should be applied to the next forecast.

The direct mean squared difference between the ensemble mean and observations (the mean squared innovation) is the sum of the forecast and observation error variances (neglecting any correlation between forecast errors and the error of observations made after the forecast was produced). To estimate the rms error of the forecast alone, it is necessary to subtract an accurate estimate of observation error variance. Unlike data assimilation (where an exaggerated observation error estimate simply means the observation will be ignored) an inaccurate observation error estimate in this context leads directly to an overspread or underspread ensemble. In data assimilation, it is common to exaggerate error variances for correlated observations, to make up for not representing that correlation explicitly. The transform matrix calculation should use the exaggerated error, to replicate the information content assigned to the observation by 4DVAR, but ensemble spread calibration requires the true observation error variance without any exaggeration. For this reason, not all observations can be used to calibrate the spread, only those for which a 'pure' observation error variance estimate is available. Initially, only sonde data up to level 26 was used. This creates a sampling problem, with some localisation centres being calibrated by few or no observations, and the potential for very large perturbations to develop unchecked in unobserved regions, particularly in the tropics. For this reason, localisation was restricted to the extratropics, between 20 and 70 degrees latitude (LatThresholdPoles and LatThresholdTropics, respectively). Ticket ENS:#101 (internal only) introduced ATOVS data into the inflation factor calculation, to support the extension of localisation to the tropical and polar regions, as well as vertical variation in the inflation factor calculation to improve the relation between spread and skill as a function of height (ticket ENS:#205 (internal only)).

The detailed operation and implementation of EnsProg_ETKF is described in the following sections.

Initialisation

The start of [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsProg_ETKF.f90 EnsProg_ETKF.f90] performs the following initialisation steps:

• Initialise the GEN system, including the documentation URL to be associated with HTML messages.
• Open the standard statistics file and emit the 'ETKF statistics' heading.
• Print a banner to standard output showing the current date and time.
• Initialise the IO timing and GEN execution tracing mechanisms, providing default values for some settings.
• Read most (but not all) of the environment variables used by the ETKF. These are stored as local variables and passed as parameters to relevant subroutines, whereas other subroutines read their own environment variables directly.
• Read the $ENS_CONTROL_DIR/Ens_ETKFControl.nl file.

The Ens_ETKFControl.nl file contains a namelist through which various ETKF configuration options and settings can be specified. This allows them to be changed without requiring a new ETKF executable, and allows the same executable to be used for the global and regional ensemble, through the use of two different control directories. Unlike environment variables, some of whose values change from run to run, the control file settings are only changed by explicit human intervention at Parallel Suites or Emergency Changes.

The control file is read by [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_ReadETKFControl.f90 Ens_ReadETKFControl.f90] into variables contained with the EnsMod_ETKFControl module. Any settings which are not specified in the control file will retain the default values specified in [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsMod_ETKFControl.f90 EnsMod_ETKFControl.f90]. The meaning of the individual settings will be described in the following sections alongside the code which they each control.

The ticket ENS:#101 (internal only) changes introduce some variable length arrays for which Ens_ReadETKFControl performs extra processing. If no inflation categories have been defined, a single category covering sonde observations (essentially reproducing the previous behaviour) is used. The number of inflation categories is determined and stored in nInfCats. If any category has been given zero weight (so that its statistics are passively accumulated but it does not affect the final inflation factor that is actually used), a warning message is printed. The number of satwind error scaling control points is also counted and reported. For more details, see the sections on satwind error scaling and inflation factor categories. [#Verticallocalisation Vertical localisation] introduces a further array of weighting functions, and adds a bands dimension to various arrays. In the control file, the multidimensional arrays are mostly readily specified as ArrayName(:,1) # Val1, Val2, ... ArrayName(:,2) Val1, Val2, ...

Reading varobs files

The first main job of the ETKF is to assemble all relevant observation and pseudo-observation data into a form suitable for the subsequent calculations. This process is orchestrated by the [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_ProcessVarobs.f90 Ens_ProcessVarobs] routine, called by EnsProg_ETKF.

The observation and pseudo-observation data produced by the Observation Processing System (OPS) is provided to the ETKF in the directories $ETKF_VAROBS_TEMPLATE, one for each ensemble member. Observations are organised into observation groups, which represent the major divisions in how data is processed and the shape of the OPS output (number of levels, variables, etc). The observation groups currently supported by the ETKF are listed in [ops:browser/main/trunk/src/code/EnsMod_Varobs/EnsMod_VarobsNumber.f90 EnsMod_Varobs/EnsMod_VarobsNumber.f90]. For each observation group, the OPS will normally produce a GroupName.varobs file containing the processed observations, and a GroupName.modelobs file containing the corresponding pseudo-observations. The ETKF uses the pseudo-observations for the perturbed members (1-44), and takes the actual observations from the unperturbed control forecast (member 0). The control is preferred over the perturbed members as it should provide the most accurate background forecast against which to process the observations. It is preferred to the high-resolution deterministic varobs files because it has the same model levels and orography as the rest of the ensemble. This is particularly important for sonde data, which is processed onto model levels (see ticket ENS:#97 (internal only)).

One or more observation groups are occasionally missing from individual forecasts. It is important that the ETKF code correctly recognises and handles this case, setting nObs and nPresent to zero in the corresponding Varobs structure.

The varobs file format is described in OTDP16. It has the same overall structure as a fieldsfile (including a fixed length header, other headers, lookup records, and data, see UMDPF3), but the array shapes are very different. To permit parallel processing, the observations are divided into batches, each described by one lookup header. Within a given varobs file, every observation record has the same number of variables and levels, although some of these entries may be filled with the missing data indicator. The level dimension is used for channels in the case of satellite observations. Each possible variable has been allocated an entry within the Column Dependent Constants header, which indicates whether that variable is used, and if so its offset within the data records. EnsProg_ETKF uses the ColDepC entry index as a code number to identify the observation data type.

The ETKF reads each varobs file into a Varobs structure, defined in [ops:browser/main/trunk/src/code/EnsMod_Varobs/EnsMod_VarobsType.f90 EnsMod_Varobs/EnsMod_VarobsType.f90]. This mimics the observationvariablelevel arrays of the varobs files, with potential missing values, unlike the EnsObs structure to be discussed later, into which the matched observations are ultimately assembled.

The [ops:browser/main/trunk/src/code/EnsMod_Varobs/Ens_ReadVarobs.f90 Ens_ReadVarobs] routine reads a single varobs file into a Varobs structure, allocating the arrays within the structure to match the data shape found within the file. [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_ReadWholeVarobs.f90 Ens_ReadWholeVarobs] fills an array of Varobs structures from the varobs (or modelobs, as appropriate) files for each recognised observation group in a given directory. If no file is found for a given observation group, the nObs entry in the corresponding Varobs structure is set to zero and the internal arrays are not allocated. The ELEMENTAL Ens_DeallocateVarobs subroutine in [ops:browser/main/trunk/src/code/EnsMod_Varobs/EnsMod_VarobsType.f90 EnsMod_VarobsType.f90] provides convenient deallocation of arrays of varobs structures, taking account of which internal arrays were actually allocated.

After reading the observation data for all ensemble members, Ens_ProcessVarobs writes the observation counts for each member and observation group to standard output.

Reading auxiliary data

For most observation types, the observation error estimate used by VAR is included with the observation in the varobs file. For certain satellite data types, however, the observation error estimate is fixed as a function of channel and stored separately to save repeating it with every observation. The ETKF therefore has to read the error estimates from these files and merge them into the rest of the data. The pure ATOVS error estimates are provided in a similar file within $ENS_CONTROL_DIR. The construction of these error estimates has been documented in a separate intranet page. In the case of ATOVS data, there is also a channel selection that is specific to VAR rather than OPS, which is read from a further data file.

To reduce the complexity and local variable count of [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_FindNumberMatched.f90 Ens_FindNumberMatched] (discussed below), and avoid the cost and potential inconsistency of reading this data twice, it has been collected together in [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsMod_AuxData.f90 EnsMod_AuxData]. The Ens_ReadAuxData routine is called from Ens_ProcessVarobs after the varobs files have been read. This ordering arises because the observations define the satellites for which error estimates are required. If the VAR error estimate or channel selection files cannot be read, the corresponding observations cannot be used, a fact which is represented by setting nPresent = 0 in the relevant Varobs structure. If the file can be read but lacks data for a particular satellite, the internal SatInfo % RMatrix % NChans and SatInfo % nchannels will be set to zero, again leading to rejection of the observation in Ens_FindNumberMatched. This second mechanism is also used for the pure error estimates, since an observation without a pure error estimate can be used in the transform matrix, but not the inflation factor calculation.

Ticket ENS:#62 (internal only) introduced the SSMIS and IASI observation groups, which involve similar files and processing to ATOVS. It is expected that this 'SatRad' formulation will increasingly dominate satellite observation processing. To support this as cleanly and flexibly as possible, EnsMod_AuxData contains an array of SatRad_Info structures, each providing the auxiliary data for one SatRad observation group. The observation groups which fall under this umbrella are defined by the SatRadGroups array in [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsMod_AuxData.f90 EnsMod_AuxData.f90] (see [ops:browser/main/trunk/src/code/EnsMod_Varobs/EnsMod_VarobsNumber.f90 EnsMod_Varobs/EnsMod_VarobsNumber.f90] for the names associated with these code numbers).

Both Ens_ReadAuxData and Ens_FindNumberMatched use a number of specified metadata fields which are provided for particular satellite observation groups. The required search over the Varobs % ID array to identify the iPresent indices to use in the Varobs % Value array is performed in EnsMod_AuxData and the results saved within the AuxDataType structure.

Observation matching

Since each ensemble member constitutes a different background forecast, OPS filters such as the background check, cloud screening, etc will occasionally select different observations in different members. Because the control and perturbed members (and observations verses pseudo-observations) are processed differently, it is even possible for them to provide different sets of variables for the same observation group. The ETKF transform calculation requires a pseudo-observation for each ensemble member, so only observations which are available for all ensemble members can be used. The next step in Ens_ProcessVarobs is therefore to invoke Ens_MatchData for each observation group of each ensemble member, to identify whether it contains each variable and observation found in the control member, and if so at what index. The results are stored in the IDMatched and Matched arrays respectively. This process is particularly amenable to parallelisation, as indicated by the #$OMP directives in Ens_ProcessVarobs.

In performing the observation matching, Ens_MatchData requires that the latitude, longitude and time match to within a Tolerance of 0.001. For observation groups which report a unique Level with each observation (and therefore have just one level per observation), this is also required to match within the same tolerance. Otherwise, the number of levels must match, or the whole observation group is rejected. Since the IDMatched results will generally be the same for all ensemble members, messages about unmatched IDs are only output for ensemble member 1, unless VerboseMode diagnostics have been requested.

Collating observations

The final step in Ens_ProcessVarobs is to assemble the matched observations into a form suitable for calculating the ETKF transform matrix and inflation factor. This requires a flat array of scalar observations, where all levels, data types, etc, have been exploded out into a single 'observation' dimension, and all missing, invalid, unused or incompletely matched data have been removed. This format is provided by the EnsObs structure defined in [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsMod_EnsObsType.f90 EnsMod_EnsObsType.f90], which is constructed from the Varobs, AuxData and matching arrays by [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_FindNumberMatched.f90 Ens_FindNumberMatched].

To support later filtering and localisation, EnsObs retains metadata such as the observation group, latitude, longitude, level and data type, although these are now associated with each individual observation rather than a multilevel/variable block of observations as in the Varobs structure. If the number of ensemble members were known at compile time, it would be possible to define an EnsOb structure for a single observation and use an array of these to represent the complete set. However, because nEns is read from the control file, and it is undesireable to nest the EnsObs definition and all procedures which use it inside a routine which takes nEns as an argument, it is necessary to allocate at least the nEns*nObs forecast perturbation array (FcValue) separately. It then becomes natural to make each element of EnsObs a separate array, as in the Varobs structure. CONTAIN-ed allocation, deallocation and copying functions have also been defined to reduce the chance of omissions when new members are added.

In order to correctly dimension the EnsObs arrays, Ens_FindNumberMatched is actually run twice. The first pass evaluates which observations should be included and returns the total count. The second pass repeats the evaluation of which observations should be included (since there is no convenient way to store the result) and fills the supplied EnsObs structure with data for each selected observation.

To obtain the complete set of scalar observations, Ens_FindNumberMatched contains a four layer nested loop over the levels contained within each observation for each variable inside each observation group. Various tests are invoked to decide whether an individual observation or set of related observations should be rejected:

• Satellite metadata fields such as channel number, stratospheric temperature extrapolation, and SSMI total/precipitable water should not be considered by the ETKF. These variables are simply skipped before entering the loop over individual observations, although several of the metadata fields will be consulted when Ens_FindNumberMatched processes the corresponding brightness temperature observations.
• Variables which are not present for both the observations and all ensemble members are similarly skipped.
• Observations which were not matched between the observations and all ensemble members are skipped, since the ETKF requires a value for each ensemble member.
• Observations which have been set to the missing data indicator in the observation or any of the pseudo-observations, or have an unacceptably high probability of gross error (PGE), or are marked as having failed OPS quality control control tests (iPresentIR/iPresentMW where appropriate), are skipped, as are observations for which any pseudo-observation differs from the actual observation by more than 100 times the VAR observation error.
• Satellite observations for which the auxiliary data does not provide a VAR error estimate are skipped, since the ETKF requires normalisation by this value. (A lack of pure error estimate, by contrast, only excludes an observation from the inflation factor calculation).
• SatRad observation groups (ATOVS, SSMIS, IASI, AIRS, SEVIRI) consult a channel rejection file which masks out channels that are currently ignored by VAR or are used by the high-resolution system but are too near the top of the low-resolution ensemble model to be reasonably simulated in the forecast perturbations. This latter restriction stops channels the model could never hope to fit from dominating the structure of the analysis uncertainty estimate encoded in the ETKF perturbations.

In assembling the EnsObs entry for accepted observations, the forecast and observed values, observation group code, data type code, latitude, longitude and time are simply copied from the Varobs array. For observation groups which report a unique Level with each observation, that value is copied from the Varobs entry, otherwise the level index or a suitable encoding of satellite channel number (and satellite id where appropriate) is used. The VAR error estimate is copied from the Varobs entry or auxiliary data as appropriate.

Pure error estimates are assigned by the GetPureError function, which returns 0 if no suitable error estimate is available. The current implementation provides three sets of pure error estimate:

1. Sonde levels 26 and below, copying the VAR error estimate since the signals from direct error estimation methods were complex and inconclusive, and this is the set of data that was used for inflation factor calculations before ticket ENS:#101 (internal only).

2. ATOVS satellite/channel pairs for which reasonable error estimates were obtained in the observation error study. These are read from ATOVS_Rmatrix_Pure in $ENS_CONTROL_DIR (not $OPS_SATRAD_Coeffs since they have been defined by the ensemble group rather than the satellite data assimilation group). Unlike the $OPS_SATRAD_Coeffs/ATOVS_Rmatrix_Var file used for the VAR error estimates, ATOVS_Rmatrix_Pure contains zeros for satellite/channel pairs for which no suitable error estimate is available.

3. Satwind error estimates obtained in the observation error study. The ETKF receives the VAR satwind error estimates in the varobs file rather than by reading a file like ATOVS_Rmatrix_Var, and there are plans to make the OPS calculation which produces this error estimate more dynamic in future. The pure error estimate has therefore been implemented as a scaling of the supplied error rather than a replacement as was used for ATOVS. Three variables in the ETKF control file are used to define error scaling factors for one or more control pressures. The actual error scaling for each satwind observation is obtained by linear interpolation, or following the nearest control point if the observation is outside the range covered by them.

On completion of the second call to Ens_FindNumberMatched, the EnsObs structure contains all the observational information, metadata and error estimates required for subsequent calculations, so Ens_ProcessVarobs is able to release the Varobs and AuxData arrays. This helps to reduce the peak memory usage in Ens_ApplyTransform.

Miscellaneous observation processing

Having obtained the MatchedData EnsObs structure from Ens_ProcessVarobs, [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsProg_ETKF.f90 EnsProg_ETKF] performs a number of processing and diagnostic steps before passing it to the main ETKF calculation:

1. The complete contents of MatchedData is recorded in an 'Innovations' text file in $ENS_OUTPUT_DIR. This file is useful for quick examination of the set of accepted observations, the degree of agreement between ensemble members and observations, and the supplied observation error estimates. Since it collates the key observational information relevant to the ensemble, the innovations file is also very useful for statistical analysis of ensemble performance, and formed the basic input to the observation error study. Note that the raw data is now output, before debiasing or subtraction of the ensemble mean. This allows subsequent analysis to potentially examine those biases, and construct statistics such as Brier scores which require the full value of the variables they examine.

2. If there are nontrivial latitude thresholds which have not yet been applied, collate the selected observations at the start of the MatchedData arrays and adjust i4TotalPoints accordingly. This facility was added so that ensemble performance in the tropics could be analysed through the innovations file even when tropical data is excluded from the main ETKF processing.

3. Subtract the ensemble mean from the forecast and observed values for each observation. This converts them into forecast perturbations and innovations respectively, still in the units of the original observation.

4. Call [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_DebiasInnov.f90 Ens_DebiasInnov] to subtract the mean value from selected subgroups of innovations. This only affects the inflation factor calculation (since the transform matrix only uses the forecast perturbations) and currently applies to each distinct ATOVS channel. The rationale for doing this is that ensemble spread should represent random rather than systematic error, and the model can only be expected to predict and evolve errors on its own attractor, not the difference between the model and real world attractors. Ens_DebiasInnov has been written in a generic manner reminiscent of the new inflation factor calculation, where categories are defined by the assignment of a nonzero category code, and the rest of the algorithm operates generically on each unique code number produced. The GrowArr function avoids the need to hard code the maximum number of categories, instead the working array is expanded as required to fit the number of categories generated. Each expansion allocates twice the previous size, so that the net allocation cost per element is constant. Similar GrowArr routines are used in [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateInflation.f90 Ens_CalculateInflation] and [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_MonitorInflation.f90 Ens_MonitorInflation]. It has been designed to be easily cut-and-pasted, but cannot be shared because the array has a different element type in each case and Fortran lacks the template facilities found in languages such as C++.

5. The overall rms innovation, spread and count for each observation group is written to the statistics file. This permits quick examination of the overall spread/skill relationship, comparison of performance between runs, and identification of empty observation groups. Unlike the innovations file, these statistics are constructed after latitude thresholds have been applied and channel-specific biases removed, so that they more closely summarise what the transform matrix and inflation factor calculations will see.

Transform matrix calculation

The next block of EnsProg_ETKF calculates the transform matrix or matrices which will be used to derive the initial conditions for the new ensemble run. These define the analysis perturbations for each ensemble member as a linear combination of the forecast perturbations, with near-unit determinant. In all cases, this proceeds in two stages, first the calculation of a basic transform matrix, then multiplication by a scalar inflation factor, whose calculation is described in later sections.

The simplest mode of all (not currently used) is simple breeding, triggered by ETKFNotBreeding=.FALSE.. In this case, the transform matrix is a simple identity divided by FixedFactor and multiplied by the global inflation factor. The next simplest mode, where a single ETKF transform is calculated for the whole domain, is specified by ETKFNotBreeding=.TRUE. and LocalETKF=.FALSE.. This is currently used for the regional ensemble, since its domain is of a similar size to the current data selection diameter used for localisation. The localised mode, where a separate transform matrix is obtained for each local region, is described later.

The raw ETKF transform matrix for a given set of observations is found by [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateTransform.f90 Ens_CalculateTransform]. This uses the forecast perturbations divided by the VAR observation error estimate, representing the ratio of background uncertainty to observation uncertainty seen by the data assimilation system. See Wang and Bishop, 2003 for details. The eigenvector decomposition is calculated using the LAPACK DSYEV routine. In the Met Office the LAPACK routines have been downloaded to an internal repository (http://fcm1/projects/downloads/browser/lapack/branches/dev/frwd/r113_LapackForOPS/lapack-3.5.0) and these are automatically included in the build process. Collaborators will need to similarly download the LAPACK software to their site. The Wang et al, 2004 method of centring the perturbations without resorting to positive/negative pairs is used in preference to the original Wang and Bishop, 2003 transform matrix formula (this change amounts to postmultiplication by C^T^). A 'bias amelioration' option is included but not normally used.

Inflation factor calculation

As discussed in the [#OverviewoftheETKF Overview], the raw ETKF transform matrix will generally underestimate the magnitude of the analysis uncertainty, essentially because the small number of ensemble members leads to spurious long-range correlations in the background error estimate. This problem is mitigated by multiplying the transform matrix by an inflation factor, chosen on the basis of the spread/error relationship in previous runs. It should be noted that this calculation is of necessity historic, and assumes that the degree of inflation required for the new run will be similar to that required for recent previous runs.

In order to support extension of horizontal localisation to the tropical and polar regions, as well as vertical localisation, ticket ENS:#101 (internal only) introduced a generic mechanism for using more than one observation type in the inflation factor calculation, and specifically tested the addition of ATOVS data. The point of extending localisation is to improve the relationship between spread and error as a function of latitude, longitude and height, and this of necessity requires only using observations from the relevant region to calculate the corresponding inflation factor. Sondes alone do not provide good sampling in some of these regions, hence the wish to use other observation types, provided suitable pure error estimates are available.

One consequence of using multiple observation types to calculate the inflation factor is that they can and do give different answers in some cases. Without vertical localisation, they may for instance be dominated by the spread/skill relationship at different heights, or for different atmospheric variables, or in a different subset of conditions. Lumping these distinct error populations directly into a single rms error would be statistically dubious and vulnerable to run-to-run fluctuations in the relative proportions of the different observation types. For this reason, the new scheme introduces the concept of inflation factor categories, where each category is supposed to provide internally consistent sampling. Two categories are currently used, namely sonde data and ATOVS data. Neither of these is necessarily homogeneous in its internal statistics, but to the extent that each ATOVS observation includes the same set of channels and each sonde observation the same set of levels and variables, the subgroup proportions within each category should be reasonably consistent. It was also seen as desirable not to introduce a large number of distinct categories, since this would increase the number of forecasts over which it would be necessary to average to obtain reasonable results in any one category. Nonetheless, the code is very generic, and the introduction of finer category divisions could basically be achieved by making the category code assigned to each observation near the top of [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateInflation.f90 Ens_CalculateInflation] depend on more pieces of metadata.

For observations within any one category, the inflation factor that should have been applied to the previous run is estimated as

I = I,,0,,*E/S

where I,,0,, is the inflation factor that was actually used, E is the rms error of the ensemble mean for the selected observations, and S is the rms spread of the corresponding pseudo-observations. The ETKF only deals with the perturbed ensemble members, so these averages all exclude the control member. In the current ETKF, rms spread is always calculated with a division by nEns (the number of perturbed members over which the spread was calculated, i.e. 45) rather than nEns-1. The factor by which the previous run is underspread around T+6 hours, E/S, is called the "inflation scaling factor", or IScale in the standard statistics output.

The above formula implicitly assumes linear perturbation growth, so that the factor by which the initial spread should be increased is equal to the factor by which the ensemble is underspread around T+6 hours. In reality, this is probably not a good assumption, since studies have shown that spread grows rapidly up to a certain scale and then has a tendency to saturate. This would mean that the inflation scaling factor underestimates the fraction by which the inflation factor needs to change to achieve a given change in spread around T+6 hours, although this would be picked up as a continued need to adjust in subsequent runs, and would still give IScale~1 once the correct inflation factor has been achieved.

For each region, the new inflation factor calculation tracks three pieces of information between runs:

1. The final inflation factor actually used. This is multiplied by the category-specific inflation scaling factor at T+6 to obtain the inflation factor desired for that particular category and run. In many cases, the final inflation factor could be worked out from the previous run's category averages, but explicit tracking allows for better handling of changes to the calculation (since it avoids having to assume what calculation was performed on the previous run) and also covers the case where the ScalingLimit (described below) reduced the inflation factor below the standard combination of the category averages.

2. The estimated inflation factor desired by each category. To improve sampling and allow for runs without observations in a particular region, this is maintained as a running average, so that the inflation target deduced for the most recent run is averaged into a history of such results (see below for details). The average for each category must be separately tracked because, as discussed above, they potentially measure different things, so that the ATOVS average inflation factor can only be calculated from ATOVS observations.

3. A similar running average observation count for each category. This allows the weights in the running average inflation factor to be adjusted to take account of whether the current run has more or fewer observations than normal in a given category, and thus obtain a more optimal estimate of the average inflation factor. It also provides a simple way to monitor the quality of sampling in different regions. The current system takes no account of the relative accuracy of or correlation between different observations, it simply assumes information content is proportional to observation count.

The following sections describe the steps of the new inflation factor calculation in more detail.

Vertical localisation

Ticket ENS:#205 (internal only) introduces vertical localisation of the inflation factor calculation, to improve the spread/skill relationship as a function of height. The ticket includes links to design notes, code changes and verification results. For simplicity and to minimise the disruption to model balance, only the inflation factor calculation is localised - there is still only one underlying raw transform matrix per horizontal region. The localisation is achieved by dividing the atmosphere into a small number of bands, performing what is mostly a separate inflation factor calculation in each band using relevant local observations. The bands are defined in terms of model levels through the Ens_ETKFControl.nl file, and [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_ApplyTransform.f90 Ens_ApplyTransform] linearly interpolates the inflation factor between band centres as discussed [#Applyingthetransformmatrix below].

Local observations such as sonde data can be directly assigned to a single band. Integrated observations such as satellite radiances can include contributions from more than one band, and the optimal increment to the inflation factor for each band can only be determined by solving an inverse problem. The current implementation takes a simpler approach, where the contribution of each observation to each band is taken as equal to that band's contribution to the observation, as expressed through a weighting function supplied in Ens_ETKFControl.nl. This approximation tends to underestimate the required change in inflation factor. In the long term, such an approach should still converge to the correct answer, and in practice the near-surface increments tend to be restrained by the [#Scalinglimiter transform scale limiter] anyway.

Inflation factor phases

At PS29 (implemented on 28 Mar 2012) the MOGREPS-G cycle time changed from 12 to 6 hours. This brought it into line with deterministic systems, supporting downstream models by providing more frequent updates, and improving the lead-time consistency of the information provided to hybrid data assimilation.

In moving MOGREPS-G to 6-hourly cycling, there is a danger that the required inflation factor will systematically differ between 00/12Z cycles and 06/18Z cycles due to the impact of the difference in conventional observation counts on the raw transform scale. A simple 6-hour cycling would detect that the 06Z run (say) required a lower inflation factor, but apply this to the 12Z run, which actually needs a larger inflation factor. Thus, the used and required inflation factors would be out-of-phase, harming forecast spread and hence overall performance.

ENS:#482, ENS:changeset:2165 (internal only) addresses this problem by introducing the concept of inflation factor 'phases'. In the recommended two-phase configuration, separate sets of running average inflation factors and associated observation counts are maintained for 06/18Z and 00/12Z. There is still only one final inflation factor per region, since this only needs to propagate to the next cycle (where it is multiplied by the observed error/spread ratio to determine that inflation factor which should have been used in the previous cycle). Verification results demonstrate that this configuration is beneficial, although it is necessary to increase the transform scaling limit from 1.2 to something like 2.0 in order to avoid a detrimental decrease in spread. This highlights the problem of insufficient perturbation growth by the forecast model. A 6-hour cycle provides the ETKF with almost twice as many observations per day, and gives perturbations half the time to grow between cycles.

The multi-phase calculation has to interact with two sets of inflation factors on any given run. The observations on this cycle, together with the final inflation factor used on the previous cycle, can be used to update the running average inflation factors for the previous cycle. However, the new perturbations and transform scaling limit calculation should use the inflation factors appropriate to the current run, which were updated on a previous cycle. The running average data is cycled through the phases dimension, which has length nPhases+1 in memory to simplify the implementation. The input inflation factors file is read into elements (0:nPhases-1), the updated averages for the previous run are stored at (nPhases), the averages actually used for the current run are taken from index (1), and rows (1:nPhases) are written to the output file.

File interaction

[ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateGlobalInflation.f90 Ens_CalculateGlobalInflation] manages data persistence for the global inflation factor calculation. Assuming a varying inflation factor is to be calculated (VaryingInflation=.TRUE.), it calls [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_ReadInflation.f90 Ens_ReadInflation] to read data from the old inflation factor file, [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateInflation.f90 Ens_CalculateInflation] to perform the inflation factor calculation and update the persisted data, and [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_WriteInflation.f90 Ens_WriteInflation] to write that data back to the new inflation factor file. The final inflation factor is returned to the caller after deallocating the working data.

The working data for the inflation factor calculation is stored in an InflationData structure, defined in [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsMod_Inflation.f90 EnsMod_Inflation] (After ENS:#205 (internal only), this structure also contains [#Inflationsummarystatistics summary statistics] which are not persisted between runs). The category average inflation factor and observation count are both nCategories*(nPhases+1)nBandsnRegions arrays. A separate array lists the category code for each column of these arrays, and the final inflation factor for each band and region is stored in a further array. For the global inflation factor calculation and file, there is just one region, whereas multiple regions are used for the localised calculation described later.

The inflation files start with a header line giving a version number and the number of categories, phases and bands used. The number of regions is implicit in the file name (one for the global case, a formula based on nSegments for the localised case). The category codes are then listed, followed by a block for each band and region providing the associated inflation factor data. Each block starts with a line giving the final inflation factor, category average inflation factors, and category average observation counts used in the most recent run. Second and subsequent lines provide the category average inflation factors and observation counts that were used in successively older phases. Ens_ReadInflation uses list-directed input to read the required fields regardless of line breaks and formatting. In particular, this allows pre-ticket 101 inflation files to be read after the addition of a "1 0" header line. The presence or absence of band count information is automatically handled through the file version number. Ens_WriteInflation uses explicit formatting to improve layout and in particular avoid the line breaks that arise with list-directed output. The formatting of inflation factors is quite aggressively terse (only three significant figures) because any further precision is probably meaningless. However, it is crucial that the values are formatted with the 'G' specification, so that very small values (as occur in some localisation centres for reasons that are not clear) do not get rounded down to zero. When more than one band is used, a blank line is inserted between regions to improve readability.

An important feature of the new system is that categories can be added, removed, split or combined with relative ease. The set of category codes considered by the current run is defined by the InfCatCode array in the control file. The working data for the current run is always based around this category list. Ens_ReadInflation is responsible for initialising this working data based on whatever set of categories it finds in the old inflation file. First of all, it reads the file verbatim into a temporary InflationData structure. The final inflation factor for each region is copied to the output structure provided the file was read successfully, otherwise the default value supplied by the caller is used (the caller supplies the default since it is different for global and localised calculations). The input is then searched for each desired output category. Any categories which are not found are initialised from the final inflation factors array. This allows the category average inflation factor to move smoothly away from replicating the previous result as information accumulates. No similar initialisation is available for the category average observation count, so the marker value zero is written. This causes Ens_CalculateInflation to initialise this average with the first nonzero observation count it encounters, temporarily assuming this to be representative of the unknown historic observation count for the new category.

If a category were ever to be split into two or more subcategories, another possibility would be to initialise the new categories directly from the old supercategory, rather than the final inflation factor. It would be simplest if obsolete category codes were not reused, although it would in principle be possible for Ens_ReadInflation to perform a mapping between old and new code systems provided it had a way (such as a change to the Version number) to work out which system it had received as input.

If the number of vertical localisation bands requested in Ens_ETKFControl.nl differs from the number of bands recorded in the old inflation factor file, a simple linear interpolation is applied to the final and category inflation factors. This is performed in band space relative to band centres, replicating the band band centre values in the outer half of the old first and last bands. No account is taken of the relative size of the old and new bands, both to simplify the interpolation logic and avoid the need to record the old band sizes in the inflation factor file. In particular, this allows an automatic transition from the old single band system to a vertically localised inflation factor calculation, where the old inflation factors are simply replicated in each band. If a more sophisticated interpolation were required, this could always be done offline, since no interpolation is performed if the old and new band counts are identical. Where the band count changes, the category average observation counts are reset as for a newly introduced category, since there is no clear basis for apportioning the old observation count between the new bands.

If the number of phases requested in Ens_ETKFControl.nl differs from the number recorded in the old inflation factor file, the above adaptations to the requested category set and band count are applied within the phases common to both input and output. Extra phases are generated by cyclic replication if required. (This would, for example, split two to four phases initialising 06/18 from their previously common phase and 00/12 from the other one).

Updating category averages

[ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateInflation.f90 Ens_CalculateInflation] updates one row of an InflationData structure based on the observation data for the corresponding region around T+6 hours into the preceding ensemble forecast.

The first block of code loops over the supplied observations, assigning a category code to each one that has an associated pure error estimate. At present, the category code is equal to the observation group index, but information such as channel, satellite, variable and so forth could easily be included. Indeed, the number of possible codes does not even need to be known in advance, since integers a and b can be uniquely combined with the formula (a+b)*(a+b+1)/2+a, based on successive rows of triangular numbers. This is the only piece of code that has to 'know' what the category codes mean in terms of observation characteristics; everything else simply works in terms of the assigned code numbers.

Having assigned a category code, the routine searches for it in the InfCatCode array of codes which are of interest for the current run. If it is not listed, a warning message is issued and observations falling within the category are ignored. (The GrowArr technique is used to avoid duplicate warning messages, the array being persisted across calls by placing it inside [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsMod_Inflation.f90 EnsMod_Inflation]). Otherwise, contributions to the rms innovation and spread, normalised by the pure error estimate, are added to the accumulations for the identified category. In the case of vertically-integrating observations, these contributions are weighted by the appropriate WfWeights entry from Ens_ETKFControl.nl, as explained [#Verticallocalisation above]

Once all observations have been processed, the rms spread and error for each category are calculated, and from these the inflation scaling factor and single run target value. In normalised units, the rms forecast error is obtained as sqrt(AvgInnov-1), where AvgInnov is the average squared innovation (significant global mean differences having been removed at an earlier stage by Ens_DebiasInnov as discussed [#Miscellaneousobservationprocessing above]). These calculations are only possible if there was at least one observation in the specified category and the average squared innovation is greater than or equal to one; if not, then no update to the category average inflation factor can be made.

The category average observation count is obtained as a weighted average of the average from the previous run and the number of observations used for this category and region in the current run. The original implementation (used in PS18) only updated the average when an inflation scaling factor could be calculated, effectively making it a conditional average of the observation count when that count is greater than zero. This is discontinuous with the behaviour at small observation count, and leads to the inflation factor for a region with, say, 2N observations every other run adjusting much more slowly than a region with N observations every run. Changeset [576] addresses this problem by including zero observation counts in the running average. Cases where a region has observations but they produce an average squared innovation less than one are counted as having zero observations, since only observations which have actually contributed to the running average inflation factor should be counted.

Repeated weighted averages in successive runs form a geometric series, so that the contribution of the n'th previous run is (1-w)w^n^, where w (BaseHistWeight in the code) is the weight applied to the previous weighted average and (1-w) is the weight applied to the contribution from the current run. This amounts to exponential smoothing of the individual run results, with a half life given by -log(2)/log(w). This formula is used to derive BaseHistWeight from the InfHalfLife parameter in the ETKF control file, with a default value of 1.0 forecasts. This value is a compromise between speed of adjustment and currency of information on the one hand, and accumulating enough observations for accurate statistics on the other. It was reduced from the original default of 3.0 forecasts following problems in the ExtLocIF trial, noting that the turning point of the Fourier Series implied by an exponential filter occurs at a period equal to 2pi*h/ln 2 ~ 9h, where h is the half-life, with period 4 oscillations being damped by a factor of around one half at h=1. The half life has deliberately not been made a function of the average observation count in each domain, because of the wish to maintain an equal degree of currency in the different inflation factors.

To make more statistically optimal use of recent observations, the HistWeight applied to the previous average inflation factor is calculated as BaseHistWeight*OldAvgObs/NewAvgObs, where OldAvgObs and NewAvgObs are the previous and newly calculated running average observation counts. It should be noted that BaseHistWeight*OldAvgObs is the first term in NewAvgObs, so that 1-HistWeight = (1-BaseHistWeight)*nObsUsed/NewAvgObs, where nObsUsed is the number of observations used for this category and region in the latest run. The old and new estimates are thus weighted according to their allotted fraction of the new running average observation count, where historic observations are exponentially downweighted with the specified InfHalfLife. In particular, this has the limiting behaviour that HistWeight->0 as OldAvgObs/nObsUsed->0 and HistWeight->1 as nObsUsed/OldAvgObs->0.

Inflation factor rescaling

As discussed above, different inflation factor categories can and do give different results, either overall or in particular regions. The initial observation error study suggested the ATOVS signal could be quite strongly deflating, perhaps because it is dominated by the upper atmosphere, where MOGREPS-G may be overspread. Without vertical localisation, it may not be desirable to follow this signal too literally, if it is judged more important to get the spread correct near the surface where MOGREPS-G is generally underspread. To the extent that the spread measured by ATOVS and sondes are correlated, it should be possible to translate ATOVS inflation scaling factors into corresponding sonde results, and so improve sampling without overall deflation. The inflation factor rescaling option tries to achieve this by multiplying IScale for each local region and category by the ratio between the global inflation scaling factors for sondes and the category in question, with a further scaling of the departure from 1.0 by the relevant entry of InfCatRescale specified in the control file. Thus, InfCatRescale=0.0 (the default) gives no adjustment, InfCatRescale=1.0 gives full adjustment, and intermediate values give partial adjustment.

To facilitate this adjustment, Ens_CalculateInflation takes a RescaleFactors pointer argument. If NULL (as set by Ens_CalculateGlobalInflation), Ens_CalculateInflation assumes it is doing the global calculation, allocating and filling the array with the appropriate scaling factor for each category. Ens_CalculateLocalETKF propagates this result to the Ens_CalculateInflation calls for each local region, where it is applied to the category average inflation factors as they are combined to form the overall inflation factor. Note that the rescaling does not affect the category average results stored in the inflation files, these always record the actual inflation factor desired by each category, so that they seamlessly translate across changes to InfCatRescale.

In practice, the inflation factor rescaling option has not been used in any of the main trials because, at least for December 2006 with the currently specified pure observation errors, ATOVS was not found to be significantly deflating on average. With hindsight, it may not work as intended in cases where latitude thresholds have been removed, since in that case the global transform is never actually used, and thus the feedback which constrains its inflation factor calculation is removed. It is possible that the similarity of the local calculations may provide adequate constraint, but this has not been explored in detail.

When multiple inflation factor [#Inflationfactorphases phases] are used, inflation factor rescaling is the point at which the calculation changes from updating running averages appropriate to the phase of the previous cycle (stored at index nInfPhases of the phases dimension), to preparing the inflation factors that will actually be used on the current cycle (based on index 1 of the phases dimension).

Inflation factor combination

For each region, there can be only one inflation factor which is finally used to scale the transform matrix used to produce the initial condition perturbations. This is derived using a formula which minimises the weighted rms fractional discrepancy between the chosen overall inflation factor and the category average results. Since this fraction measures the ratio between the expected resulting spread and the spread desired by the given category, this amounts to minimising the weighted rms fractional overspread or underspread.

This formula is derived as follows. Let I be the chosen final inflation factor and I,,i,, be the updated (and potentially rescaled) average inflation factor for category i, with weight w,,i,,. Then I should be chosen to minimise

sum(w,,i,,*(I/I,,i,, - 1)^2^)

Expanding the square and differentiating with respect to I, this becomes

sum(w,,i,,*(2I/I,,i,,^2^ - 2/I,,i,,)) = 0

which gives

I = sum(w,,i,,/I,,i,,) / sum(w,,i,,/I,,i,,^2^)

The original implementation (used in trials up to and including ExtLocIF and hence in PS18) used the expression sum(w,,i,,*I,,i,,^2^) / sum(w,,i,,*I,,i,,), due to a confusion over whether to minimise the rms fractional discrepancy in I or 1/I. The effect is small, with a result slightly above rather than below the arithmetic mean inflation factor (and less above than the correct result is below). A simple way to see that the lower answer is correct is that a given I-I,,i,, is a larger fraction of the smaller I,,i,, than it is of the larger one. The difference is only important in cases where the different categories give significantly different target inflation factors. The formula has been corrected in changeset [569].

It should be noted that the number of observations in each category plays no role in the basic combination formula. If it did, the relative contributions of the different categories would vary in time and space with their observation counts. If the different categories have different average inflation factors, this would lead to spurious temporal and spatial fluctuations in the final inflation factor. Instead, the relative information content of runs with small and large observation counts in a particular region is taken account of in the update of the individual category averages as described [#Updatingcategoryaverages above].

However, if the long-term average observation count for a particular category and region is low, there is no basis for the inflation factor calculated from it. In the extreme case of a region with no sonde observations, the sonde inflation factor will never be changed. In such a case, it is desirable to exclude the unobserved category from the final inflation factor, relying on other categories to the extent they are present. The post-ticket 205 code therefore implements two 'fallback' provisions. For each category, RedCatWeight defaults to the fixed InfCatWeight specified in the control file, but ramped linearly to zero if the long-term average observation count falls below the MinAvgObs control file parameter (default 200). Secondly, BandObs forms an overall weighted observation count for each band across all categories. If this falls below MinAvgObs, the final inflation factor is ramped towards a fixed transform scale specified by DefaultScaling (default 0.9). In particular, this allows the definition of a special band covering the top few levels of the model, with InfCatWeight=0 and therefore a fixed transform scale of 0.9 to damp potentially damaging perturbations even when they are supported by observations. Special code eliminates these fallbacks if MinAvgObs=0, to reproduce the previous behaviour.

Scaling limiter

Having formed a best estimate final inflation factor for the region, assuming consistent behaviour in the fractional spread deficit of the raw transform matrix, Ens_CalculateInflation performs a final check on the magnitude of the transform matrix that has actually been obtained. The current implementation measures the square root of the sum of the squares of the implied transform matrix elements, divided by nEns. This reflects the expectation that the matrix will have diagonal elements just below unity, and off-diagonal elements just above zero. In most cases, ScalingFactor is expected to be less than one, reflecting the reduction in uncertainty due to the assimilation of observations. The determinant would be an alternative measure of the 'area scaling' produced by the matrix, but it should be noted that the actual magnitude of the final perturbations also depends on the incoming forecast perturbations, which need not all have the same magnitude. In addition, some of the mixing will act to reduce perturbation magnitude through cancellation, rather than increase it.

The derived ScalingFactor is compared to a ScalingLimit specified in the ETKF control file. In the general case, there are two limits which can be specified, ScalingLimit (default 1.2) for the main transforms (the local one if localisation is being used or the only one if not), and GlobalScalingLimit (default 0.9) for the fallback global transform if localisation is being used. The latter limit is smaller since in the traditional case where the global transform is primarily used in the tropical and polar regions but calculated from extratropical data, the feedback mechanism which controls its inflation factor is somewhat indirect, making it appropriate to be more cautious about perturbations near to or larger than the scale of the incoming forecast perturbations.

The ScalingLimit helps to reduce the harm caused by cases where the assumptions underlying the main inflation factor calculation are violated. A classic example of this is where a major observation type, such as ATOVS, drops out for one forecast. This will tend to increase the magnitude of the raw transform matrix by more than the increase in analysis uncertainty, so that application of the standard inflation factor would produce perturbations which are too large. Another case in which the scaling limit has proved useful is when latitude thresholds are first removed. The initial inflation factors used for the newly activated localisation centres are at best rough guesses, and the scaling limit can help to prevent excessively large perturbations until the main inflation factor calculation has time to identify the correct values.

Region-level statistics

The final section of Ens_CalculateInflation prints a table of statistics for each band, summarising results for the current region. To aid both computer and human parsing of this output, it is protected by a WriteStatsFile critical section which should be used by all parallel code that writes the standard statistics files (see ENS:#507 (internal only)). The main body of the table summarises for each category the observation statistics from the current run, the running average statistics, and the actual weight given to historical information. The first footer line gives the Overall old and new inflation factors, along with their ratio, observation count, and transform scale value. The second ('Final') footer line indicates whether these results were changed as a result of the transform scaling limit. When multiple inflation factor [#Inflationfactorphases phases] are used, the meaning of this output changes halfway through, reflecting the transition within the main calculations. Early entries (up to AvgIF) refer to the inflation factors appropriate to the previous phase, measured this cycle. Later entries (ResCatIF, and the Overall and Final lines) refer to the inflation factors to be applied in the current phase.

Subcategory monitoring

[ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_MonitorInflation.f90 Ens_MonitorInflation] is called from Ens_CalculateInflation if GeneralMode >= DebugMode, to provide a more detailed breakdown of the spread/skill behaviour underlying the relatively broad categories used in the main calculation. These will be subject to more statistical noise, but should help to indicate the extent to which behaviour is actually homogeneous within each category. The routine works in a similar way to the category accumulations in Ens_CalculateInflation, except that accumulations are made for all identified categories (using the GrowArr technique) rather than a preordained set, and an observation is allowed to belong to more than one category. In the current implementation, the category code is partitioned into three parts, representing the observation group, satellite or variable code, and level or channel number. Categories combining all satellites/variables and/or all levels/channels for a given observation group have zero in the field(s) which have been summed over. Since these category codes only have meaning within Ens_MonitorInflation and are not persisted between runs, the layout, numbers of digits used, etc can be readily changed as required. The final output shows the rms innovation and spread using both VAR and Pure error normalisation, and the inflation scaling factor calculated on the basis of the pure error results. Like the main inflation factor calculation, only observations with an assigned pure error estimate are considered.

Inflation summary statistics

Alongside vertical localisation, ENS:#205 (internal only) introduces a set of summary statistics which are collected across the full set of localised inflation factor calculations. The working data is conveniently stored inside the InflationData structure defined in [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsMod_Inflation.f90 EnsMod_Inflation], updated by each invocation of [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateInflation.f90 Ens_CalculateInflation], and finally displayed by EnsMod_Inflation.Ens_OutputInflationStats, called from [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateLocalETKF.f90 Ens_CalculateLocalETKF]. The current set of statistics include the average number of observations in each category/band, the number of localisation centres hitting various limits or special cases, and average measures of agreement between categories and the transform scale. They thus provide a quick means to assess the overall health of the localised inflation factor calculations, and could in future be extracted by system monitoring software. For more details on each individual output, see the definition comments in [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsMod_Inflation.f90 EnsMod_Inflation] and the associated data collection code in [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateInflation.f90 Ens_CalculateInflation].

Horizontal localisation

With horizontal localisation turned on (ETKFNotBreeding and LocalETKF both .TRUE.), EnsProg_ETKF prepares an array of transform matrices, one for each localisation centre. [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_ReadLocations.f90 Ens_ReadLocations] identifies the number of centres together with their latitudes and longitudes. These positions are based on subdivisions of an icosahedral grid, so that they have approximately uniform spacing across the curved surface of a sphere. The number of subdivisions is defined through the nSegments parameter in the ETKF control file, with a default value of 3, giving 92 centres. The centre locations are obtained from the corresponding UniquePoints_NN.txt file in $ENS_LOCALISATION_DIR, which is constructed by [ops:browser/main/trunk/src/control/Localisation/Icosahedron.f90]. The value of nSegments is also built into the name of the local inflation factors file (qwqgHH_Inflation_Local_NN.txt), since the inflation factors calculated for one set of regions are not directly useful for a different set.

[ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalculateLocalETKF.f90 Ens_CalculateLocalETKF] manages the process of constructing the array of transform matrices and updating the local inflation factors file. It first obtains the global transform in the same way as [ops:browser/main/trunk/src/code/EnsProg_ETKF/EnsProg_ETKF.f90 EnsProg_ETKF] in the non-localised case. This is used as a fallback in various cases described shortly. [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_ReadInflation.f90 Ens_ReadInflation] is called to read the old local inflation file, retrieving data for all centres at once. The default value procedures described [#Fileinteraction above] mean this operation always produces a complete InflationData structure, even if the old file cannot be read or lacks some of the categories required for the current run.

Ens_CalculateLocalETKF then loops over localisation centres. Where a centre falls within the regions excluded by latitude thresholds, its immediate vicinity at least will have no observations to support a local calculation, so the global transform is used (divided by two in the case of the southern polar region, to reduce the chance of developing large localised perturbations that go unchecked because they make little or no contribution to the extratropically-derived inflation factor). Otherwise, [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_FindNearObs.f90 Ens_FindNearObs] is used to obtain the subset of observations which fall within MaxRadius of the localisation centre. The new default 2000 km radius is similar to the 2000-2500 km spacing of localisation centres for nSegments=3, to avoid substantial overlap between the observations selected for different centres (Flowerdew and Bowler, 2011). An isosceles triangle construction shows that the centre point of an equilateral triangle of side 2500 km is 1767 km from each corner, so a MaxRadius of 2000 km should be sufficient to avoid any observations being excluded from all local regions. The presence of significant overlap may be responsible for the 'patchiness' observed in some runs with latitude thresholds removed, since it creates a disparity between the area controlled by a given local transform matrix and the area over which the average spread and error are calculated to calibrate its inflation factor. In particular, the local transform might perpetuate a local excess of spread if its averages include a surrounding halo of low spread, and vice-versa. Equally, a larger selection radius gives rise to more smoothly varying local transforms, which reduces the dynamical imbalance arising from localisation.

Having determined the set of 'local' observations, Ens_CalculateLocalETKF calls Ens_CalculateTransform and Ens_CalculateInflation to obtain the local transform matrix. In the case where [#Verticallocalisation vertical localisation] of the inflation factor is used (nBands > 1 in the Ens_ETKFControl.nl file), the single raw transform matrix is multiplied by the inflation factor for each band to produce a bands*centres array of transform matrices, so that vertical localisation of the raw transform matrix could be introduced without change to the Ens_ApplyTransform routine described below. The iLocal argument to Ens_CalculateInflation ensures the correct row of the InflationData structure is updated, and the RescaleFactors obtained from the global calculation are supplied. If a local transform cannot be calculated, the global transform is substituted. With the new inflation factor calculation, this could only arise due to a failure of the eigenvector decomposition in Ens_CalculateTransform, which has never been observed. Finally, Ens_CalculateLocalETKF calls Ens_WriteInflation to write the new inflation factor file and deallocates the working arrays.

Applying the transform matrix

Once the inflated transform matrix or matrices have been obtained, EnsProg_ETKF writes them to the standard statistics output and deallocates the observation data arrays. It then calls [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_ApplyTransform.f90 Ens_ApplyTransform] to construct the actual perturbation files which will be used to initialise the perturbed ensemble member forecasts. These are based on the $DATAW/operNN/qwe[gy]HH.operNN.pp1 files produced by the previous ensemble forecast, which contain full model fields at the next analysis time for the five perturbed variables (u wind, v wind, theta, specific humidity and exner pressure). Ens_ApplyTransform opens these files, uses [ops:browser/main/trunk/src/code/EnsMod_Header/Ens_ReadHeader.f90 Ens_ReadHeader.f90] to read the headers into an array of [ops:browser/main/trunk/src/code/EnsMod_Header/EnsMod_HeaderType.f90 Header] structures, and extracts information about the layout of each variable from the headers for member 1. It opens the output perturbation files and calls [ops:browser/main/trunk/src/code/EnsMod_Header/Ens_WriteHeader.f90 Ens_WriteHeader.f90] to replicate the headers it read for each ensemble member. It then reads the forecast model states and subtracts the ensemble mean to obtain the forecast perturbations.

In the non-localised case, a single transform matrix is used for all points, levels and variables, the analysis perturbations in each case being obtained by the dot product of the forecast perturbations and the output member's row of the transform matrix. In the horizontally localised case, the transform matrix changes as a function of horizontal position, which depends on the row and column indices together with the grid layout appropriate to the variable being considered. [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_CalcCurrentTransform.f90 Ens_CalcCurrentTransform.f90] defines how the transform matrix for a given latitude and longitude is obtained from the transform matrices of nearby localisation centres. It uses a simple weighted average of all centres within SearchRadius, where the weights are proportional to the great circle distance to each centre. SearchRadius is supplied by EnsProg_ETKF, which obtains it by calling [ops:browser/main/trunk/src/code/EnsProg_ETKF/Ens_FindSearchRadius.f90 Ens_FindSearchRadius]. This routine reads the $ENS_LOCALISATION_DIR/Neighbours_NN.txt (where NN is the value of nSegments) produced by [ops:browser/main/trunk/src/control/Localisation/Icosahedron.f90] and returns the maximum distance it finds between a localisation centre and its identified neighbours.

With [#Verticallocalisation vertical localisation], Ens_CalcCurrentTransform produces a horizontally interpolated transform matrix for each band. Ens_ApplyTransform vertically interpolates these in model level space about the band centres to produce a transform matrix for each level, which it then applies to the background perturbations to produce the analysis perturbations. The interpolation is performed without extrapolation, using the band 1 transform matrix for levels below the middle of band 1 and the top band transform matrix for levels above its middle. This avoids the risk of producing excessively large or negative inflation factors by extrapolation, although in reality the model underspread does tend to increase towards the surface. Whilst replicating the band centre value does not preserve the average inflation factor within each band, there is an approximate cancellation between the overspread which the interpolation produces in one band with the underspread in its neighbour. In model level space, at least, the cancellation would be exact if the transition were forced to occur over the same number of levels in each band, but the use of band centres has the advantages of simplicity and smallest gradient of change (so least disruption to model balance). ENS:#675 (internal only) introduces a SymVInterp control file option which allows the user to choose, for each transition, whether to run the interpolation to both midpoints (the smooth 'asymmetric' option) or over an equal number of levels to the midpoint of the smaller band (the 'symmetric' option). In the latter case, the middle of the larger band replicates the overall band inflation factor. This option can reduce the number of levels affected by the taper to band 4, which is used to artificially damp perturbations near the top of the model. This reduces the extent to which the artificial damping is seen by high-peaking ATOVS channels, which may otherwise lead to an artificial bump in spread lower down.

Ticket ENS:#42 (internal only) introduced a limiter on the perturbations for individual ensemble members, based on a measure loosely related to perturbation energy. This requires that all variables be processed simultaneously, which has substantially increased the memory consumption of the ETKF (it is also necessary to process all members simultaneously, due to the need to subtract the ensemble mean). It is intended to catch wildly excessive features which develop over several cycles but never encounter an observation that would cause the inflation factor mechanism to reduce their scale. Since the energy measure depends on explicitly on temperature, pressure and wind, Ens_ApplyTransform has to identify where these variables appear within the perturbation data arrays. It will also transform exner pressure to real pressure and potential temperature to real temperature if required. To help maintain balance, any scaling necessary to bring the perturbations within the specified limit is linearly tapered with a radius of ScaleSize=10 gridboxes. The calculation takes account of wrap in the global model and boundaries in the regional model, and uses the minimum indicated perturbation size in the case where scaling regions overlap.

A further change was introduced under ticket ENS:#417 (internal only). There have been issues noted regarding the large-scale perturbations by the ETKF. The data assimilation is particularly sensitive to large-scale perturbations to pressure, since this is used in the calculation of the geostrophic wind field. The code now removes the globally-averaged value of the perturbation to exner-pressure at level 1. This change to the perturbation of exner is then propagated to all levels of the model assuming hydrostatic balance applies. Additionally, we have previously had problems with perturbations to the top level of exner, since this is not truly a prognostic field. The perturbation to the top level of exner is now calculated assuming a hydrostatic balance relationship with perturbations to the level below. This replaces a previous change which set to zero any perturbations to the top level of exner.

Finally, Ens_ApplyTransform writes the transformed and potentially scale-limited perturbations and closes the data files. EnsProg_ETKF closes the standard statistics file, calls Gen_TraceReport to print execution tracing information if appropriate, io_total_timings to print IO timings, and returns.

References

Flowerdew J, Bowler NE (2011) Improving the use of observations to calibrate ensemble spread. Q. J. R. Meteorol. Soc. 137: 467-482

Houtekamer PL, Mitchell HL (1998) Data assimilation using an ensemble Kalman filter technique. Mon Wea Rev 126: 796-811

Wang X, Bishop CH (2003) A comparison of breeding and ensemble transform Kalman Filter ensemble forecast systems. J Atmos Sci 60: 1140-1158

Wang X, Bishop CH, Julier SJ (2004) Which is better, an ensemble of positive-negative pairs or a centred spherical simplex ensemble? Mon Wea Rev 132: 1590-1605

ACCESS Links

• ACCESS wiki main page
• BoM ACCESS documentation
• BoM ACCESS-GE global NWP ensemble wiki page