NWPSAF Met Office 1D-Var User Manual

Fig. A1. NWPSAF_1DVar Calling Tree

Fig. A1 shows the calling tree for the NWPSAF_1DVar code. The interfaces between these subroutines are explicitly specified through .interface files (explicit interfaces are required when using optional arguments in Fortran90 and desirable otherwise as a useful check during code development).

In addition to these subroutines, nine modules are used:

NWPSAFMod_Channellist.f90          NWPSAFMod_ObsInfo.f90
NWPSAFMod_Constants.f90            NWPSAFMod_Params.f90
NWPSAFMod_CovarianceMatrices.f90   NWPSAFMod_RTmodel.f90

Structure	Module	Description
SatIDInfo_Type	NWPSAFMod_RTmodel	This is used to define the component instruments in composite instruments such as ATOVS
RTParams_Type	NWPSAFMod_RTModel	This contains profile, instrument choice and geometry information that is required by the RT model plus the radiances and Jacobians returned from the model.
ModelOB_Type	NWPSAFMod_ObsInfo	This contains background information for each observations including modelled radiances
OB_Type	NWPSAFMod_ObsInfo	This contains the observed radiances plus auxiliary information (e.g., latitude, longitude, surface elevation etc.). It is also used to store the retrieved profile and information on the retrieval process (e.g., final cost function).
R_Matrix_Type	NWPSAFMod_Covariance_Matrices	This stores the observational error covariance matrices (in a variety of forms).
Channel_Selection_Type	NWPSAFMod_Channellist	Stores the channel selection choices.

There are three main stages in the processing chain: setting up the input structures, initialising the retrieval and the retrieval itself. In this section, the program is examined on a subroutine-by-subroutine basi s (the main program is listed first after which the programs are listed in alphabetical order).

NWPSAF_1DVar_Driver This is the main driver program
Called by: < NONE >	Calls (in order of first call): NWPSAF_Read_ControlData NWPSAF_Read_Observations NWPSAF_Read_Background NWPSAF_SetUpRetrievals NWPSAF_1DVar NWPSAF_DeAllocate

Gen_Error_Report A generic error output subroutine
Called by: Many subroutines	Calls (in order of first call): Gen_MessageReport (3 times)

Gen_FreeUnit Frees a fortran I/O unit number previously assigned by Gen_GetUnit.
Called by: NWPSAF_Initialise	Calls (in order of first call): Gen_Message_Report

Gen_GetUnit Supplies an unused fortran I/O unit number.
Called by: NWPSAF_OpenFile	Calls (in order of first call): Gen_Message_Report

Gen_MessageReport A generic routine for outputting informational messages.
Called by: Gen_Error_Report NWPSAF_1DVar NWPSAF_InitBmatrix NWPSAF_InitRmatrix NWPSAF_Initialise NWPSAF_OpenFile NWPSAF_ProcessData NWPSAF_Gastro_Interface NWPSAF_RTIASI_Interface NWPSAF_RTTOV7_Interface	Calls (in order of first call): < NONE >

NWPSAF_1DVar The main calling routine for initialising and minimising the data (the main input routines have already been called).
Called by: NWPSAF_1DVar_Driver	Calls (in order of first call): NWPSAF_Initialise Gen_MessageReport (2 times) NWPSAF_ProcessData NWPSAF_SetUpRetrievals

NWPSAF_AdditionalCost_Cloud When performing cloudy retrievals an additional cost is added to the cost function if the cloud top pressure or cloud fraction fall outside the allowed bounds (200hPa-Surface Pressure and 0-1 respectively). The additional cost increases as the cube of the distance from the bound.
Called by: NWPSAF_Calculate_Cost_Function NWPSAF_Minimize_100ML	Calls (in order of first call): <NONE>

NWPSAF_BandInverse Invert a band diagonal, symmetric, positive definite matrix. This is used with the observational error covariance matrix.
Called by: NWPSAF_Calculate_Cost_Function NWPSAF_Minimize_100 NWPSAF_Minimize_100ML	Calls (in order of first call): Gen_Error_Report

NWPSAF_BandMultiply Multiply a matrix by a band diagonal matrix.
Called by: NWPSAF_Calculate_Cost_Function NWPSAF_Minimize_100 NWPSAF_Minimize_100ML	Calls (in order of first call): < NONE >

NWPSAF_Calculate_Cost_Function Evaluates the cost function and, optionally, its gradient (returned as a vector with length equal to the length of the state vector). A subroutine that calculates an additional cost term for cloud retrievals may be called when required. Note that if the Rodgers (1976) Eqn. 101 form of the minimisation is being used, the cost function is evaluated using the diagonal elements of the observational error covariance matrix only.
Called by: NWPSAF_Minimize NWPSAF_Minimize_100ML	Calls (in order of first call): NWPSAF_SatMatInv NWPSAF_BandInverse NWPSAF_BandMultiply NWPSAF_AdditionalCost_Cloud Gen_Error_Report (5 times)

NWPSAF_Channellist Reads in channels to be used in various retrieval situations. The use of each channel that is used at all is specified in the ChannelChoice.dat file through a binary code (described in the header to the file and in this subroutine). From the same file, this subroutine reads in which channels are to be used in monitoring (stored in the BackChans structure) and in cloud detection (stored in DetectCloudChans).
Called by: NWPSAF_Initialise	Calls (in order of first call): Gen_Error_Report (6 times) NWPSAF_IntegerSort (3 times)

NWPSAF_CheckIteration Resets the profile variables to within pre-defined soft-limits (as contained in the Soft_Limits structure) and sets the Profile_Variables_Reset to true in this case. If the profile is beyond the wider predefined gross-limits, Out_of_Range is then set to true.
Called by: NWPSAF_Minimize NWPSAF_Minimize_100ML	Calls (in order of first call): QSAT (2 times)

NWPSAF_CheckTemperatures A FUNCTION that checks for the input temperatures being inside gross limits.
Called by: NWPSAF_SetUpBackground	Calls (in order of first call): < NONE >

NWPSAF_Cholesky Solves the linear equation Uq=V for q where U is a symmetric positive definite matrix and v and q are vectors of equal length. See Golub and Van Loan (1996). If U is not positive definite this will be detected by the program and flagged as an error. U is assumed to be symmetric as only the upper triangle is in fact used.
Called by: NWPSAF_CloudCost NWPSAF_Minimize_100 NWPSAF_Minimize_100ML NWPSAF_Minimize_101	Calls (in order of first call): Gen_Error_Report

NWPSAF_CloudCost Calculates the cloudy cost function as defined in Eqn. 6 of English, Eyre and Smith (1999).
Called by: NWPSAF_CloudyOrNot	Calls (in order of first call): NWPSAF_Fastmodel_Interface NWPSAF_RMatrix_ChanSelect Gen_Error_Report (2 times) NWPSAF_Cholesky

NWPSAF_CloudyOrNot This is the cloud detection module and uses the method described by English, Eyre and Smith (1999).
Called by: NWPSAF_ProcessData	Calls (in order of first call): NWPSAF_CloudCost

NWPSAF_CO2Slice The CO2-slicing technique is used to obtain a better first guess of the cloud top height and fraction before the 1DVar is used to refine this guess. The cloud detection channels are used for this purpose.
Called by: NWPSAF_ProcessData	Calls (in order of first call): Gen_ErrorReport (3 times)

NWPSAF_DeAllocate Deallocates any remaining allocated arrays before exiting the program.
Called by: NWPSAF_1DVar_Driver	Calls (in order of first call): < NONE >

NWPSAF_Fastmodel_Interface This is the interface routine between the main 1DVar code and the fastmodels (see NWPSAF_RTTOV7_Interface and NWPSAF_RTIASI_Interface for details.
Called by: NWPSAF_CloudCost NWPSAF_Initialise NWPSAF_Minimize NWPSAF_Minimize_100ML NWPSAF_ProcessData	Calls (in order of first call): NWPSAF_Gastro_Interface NWPSAF_RTTOV7_Interface NWPSAF_RTTOV8_Interface NWPSAF_RTIASI_Interface Gen_Error_Report

NWPSAF_InitBmatrix Reads in the background error covariance matrix, cuts out any unwanted elements and calculates the inverse.
Called by: NWPSAF_Initialise	Calls (in order of first call): Gen_Error_Report (2 times) NWPSAF_SatMatInv Gen_Message_Report

NWPSAF_InitRmatrix Reads in the observational error covariance matrices for all required instruments. There are three possible formats: 1) Full matrix representation 2) Band-diagonal representations (diagonal being a sub-set of this) 3) Eigenvalue-Eigenvector representation The error covariances are given in terms of brightness temperature and implicitly include forward modelling error.
Called by: NWPSAF_Initialise	Calls (in order of first call): Gen_Error_Report (3 times) Gen_Message_Report (4 times)

NWPSAF_Initialise This is the main initialisation routine. It starts by opening the files that are required by the program Next the NWPSAF_Channellist routine is called which sets up the channel selection choices for various cloud and surface type scenarios. The fastmodels are initialised next throught the NWPSAF_Fastmodel_Interface subroutine. The user may then initialise their bias correction scheme before the background and observation error covariance matrices are read in.
Called by: NWPSAF_1DVar	Calls (in order of first call): Gen_Message_Report (4 times) NWPSAF_OpenFile (7 times) NWPSAF_Channellist NWPSAF_Fastmodel_Interface Gen_Error_Report (2 times) NWPSAF_InitRmatrix NWPSAF_InitBmatrix Gen_FreeUnit

NWPSAF_Minimize This subroutine contains the main iterative loop for minimisation. One of three minimisation algorithms are used depending on the dimensions and non-linearity of the problem (see below). The minimisation process is iterated until the convergence criteria are met. These are: 1) The last iteration resulted in the cost function being reduced by l ess than a pre-defined fraction (defined by the DeltaJ variable) 2) The cost function gradient is small 3) Gamma is not decreasing (Marquardt-Levenberg minimisation only)
Called by: NWPSAF_ProcessData	Calls (in order of first call): NWPSAF_CheckIteration (2 times) NWPSAF_RMatrix_ChanSelect NWPSAF_Fastmodel_Interface (2 times) NWPSAF_Calculate_Cost_Function (2 times) NWPSAF_Minimize_101 NWPSAF_Minimize_100ML NWPSAF_Minimize_100 NWPSAF_SatMatInv (2 times)

NWPSAF_Minimize_100 Updates the input profile vector (Delta_Profile) according to Rodgers (1976), Eqn. 100. This is the optimal equation to use for linear problems where the size of the observation vector is larger than the size of the state vector.
Called by: NWPSAF_Minimize	Calls (in order of first call): NWPSAF_SatMatInv NWPSAF_BandInverse NWPSAF_BandMultiply Gen_Error_Report (2 times) NWPSAF_Cholesky

NWPSAF_Minimize_100ML Updates the input profile vector (Delta_Profile) according to Rodgers (1976), Eqn. 100, modified from the Gauss-Newton form to the Marquardt-Levenberg form (see Rodgers, 2000). This min imisation scheme should be used for non-linear problems where Gauss-Newton descent is inefficient.
Called by: NWPSAF_Minimize	Calls (in order of first call): NWPSAF_SatMatInv NWPSAF_BandInverse NWPSAF_BandMultiply Gen_Error_Report (2 times) NWPSAF_Cholesky NWPSAF_CheckIteration (2 times) NWPSAF_Fastmodel_Interface (2 times) NWPSAF_Calculate_Cost_Function

NWPSAF_Minimize_101 Updates the input profile vector (Delta_Profile) according to Rodgers (1976), Eqn. 101. This is the optimal equation to use for linear problems where the size of the observation vector is smaller than the size of the state vector (e.g., ATOVS).
Called by: NWPSAF_Minimize	Calls (in order of first call): Gen_Error_Report (2 times) NWPSAF_Cholesky

NWPSAF_OpenFile Open a file for input/output using the next free unit number.
Called by: NWPSAF_Initialise NWPSAF_Read_Background NWPSAF_Read_ControlData NWPSAF_Read_Observations NWPSAF_SetUpRetrievals RTIASI_Subroutines	Calls (in order of first call): Gen_Error_Report (2 times) Gen_GetUnit Gen_MessageReport (2 times)

NWPSAF_ProcessData This is the main processing routine where each observation in turn is looped through and processed accordingly. The observation to be processed is first read in (an earlier version of this code had all the observations read in before this routine is called but memory considerations required this stage to be moved here - all the background data is still read in before this subroutine) and then the observation and background data are converted into the correct forms for the minimisation code to use. If required, the cloud detection step is performed next after which the channels to use are determined based on the cloud conditions and surface type. After the retrieval step itself (NWPSAF_Minimize), the results output as required using the NWPSAF_TranslateDataOut routine.
Called by: NWPSAF_1DVar	Calls (in order of first call): NWPSAF_Read_Observations Gen_Error_Report NWPSAF_TranslateDataIn NWPSAF_SetupBackground NWPSAF_Fastmodel_Interface (2 times) NWPSAF_CO2Slice NWPSAF_CloudyOrNot NWPSAF_Minimize NWPSAF_TranslateDataOut Gen_MessageReport

NWPSAF_RMatrix_ChanSelect Picks out the selected channels from the observational error covariance matrix. (This is a little complicated for the band-diagonal representation of this matrix!)
Called by: NWPSAF_CloudCost NWPSAF_Minimize	Calls (in order of first call): Gen_Error_Report (2 times)

NWPSAF_RTIASI_Interface This is the routine that interfaces with the RTIASI fast radiative transfer model (Matricardi and Saunders, 1999). It may be called in three different modes - initialise, direct and gradient.
Called by: NWPSAF_Fastmodel_Interface	Calls (in order of first call): Gen_MessageReport RTIASI_RTTVI Gen_Error_Report (2 times) RTIASI_Direct RTIASI_K

NWPSAF_RTTOV_Interface This is the routine that interfaces with the RTTOV7 fast radiative transfer model. It may be called in three different modes - initialise, direct and gradient.
Called by: NWPSAF_Fastmodel_Interface	Calls (in order of first call): Gen_MessageReport (2 times) Gen_Error_Report (6 times) RTTVI CLEANUP RTTOV RTTOVK

NWPSAF_RTTOV8_Interface This is the routine that interfaces with the RTTOV8 fast radiative transfer model. It may be called in three different modes - initialise, direct and gradient.
Called by: NWPSAF_Fastmodel_Interface	Calls (in order of first call): Gen_MessageReport (2 times) Gen_Error_Report (6 times) RTTOV_SETUP RTTOV_SETUPCHAN RTTOV_SETUPINDEX RTTOV_DIRECT RTTOV_K

NWPSAF_RTTOV9_Interface This is the routine that interfaces with the RTTOV9 fast radiative transfer model. It may be called in three different modes - initialise, direct and gradient.
Called by: NWPSAF_Fastmodel_Interface	Calls (in order of first call): Gen_MessageReport (2 times) Gen_Error_Report (6 times) RTTOV_SETUP RTTOV_SETUPCHAN RTTOV_SETUPINDEX RTTOV_DIRECT RTTOV_K

NWPSAF_RTTOV10_Interface This is the routine that interfaces with the RTTOV10 fast radiative transfer model. It may be called in three different modes - initialise, direct and gradient.
Called by: NWPSAF_Fastmodel_Interface	Calls (in order of first call): Gen_MessageReport (2 times) Gen_Error_Report (6 times) RTTOV_SETUP RTTOV_SETUPCHAN RTTOV_SETUPINDEX RTTOV_DIRECT RTTOV_K

NWPSAF_RTTOV11_Interface This is the routine that interfaces with the RTTOV11 fast radiative transfer model. It may be called in three different modes - initialise, direct and gradient.
Called by: NWPSAF_Fastmodel_Interface	Calls (in order of first call): Gen_MessageReport (2 times) Gen_Error_Report (6 times) RTTOV_SETUP RTTOV_SETUPCHAN RTTOV_SETUPINDEX RTTOV_DIRECT RTTOV_K

NWPSAF_ReadHeaders Read in colon-separated fields in observation and background file headers.
Called by: NWPSAF_Read_Background NWPSAF_Read_Observations	Calls (in order of first call): <NONE>

NWPSAF_Read_Background An example input routine for reading in the background profiles used in retrievals and placing the data in the BackGrModelOb structure (of ModelOb_Type derived type). This version reads from the file BACKGROUND.dat which is a formatted ASCII file. This routine expects the number of profiles in the input file to be either equal to the number of observations or to be one - in which case the same profile background is used for each observation.
Called by: NWPSAF_1DVar_Driver	Calls (in order of first call): NWPSAF_OpenFile (2 times) Gen_Error_Report (12 times) NWPSAF_ReadHeaders (9 times)

NWPSAF_Read_ControlData Read in control data (from ControlData.NL). Details are in Section 3a.
Called by: NWPSAF_1DVar_Driver	Calls (in order of first call): NWPSAF_OpenFile Gen_Error_Report (7 times)

NWPSAF_Read_Observations An example input routine for reading in the observed radiances to be processed, plus auxiliary data, and placing the data in the Observations structure (of Ob_Type derived type). This version reads from the file ObsFile.dat which is a formatted ASCII file. The file also contains information on which of the channels belongs to which sub-instrument when a composite instrument (such as ATOVS) is being used. This version only reads one observation at a time (or none if all that is required is for the file to be opened and variables initialised) as the memory requirements are less than if all observations are read in at once. See Section 3c for more details.
Called by: NWPSAF_1DVar_Driver NWPSAF_ProcessData	Calls (in order of first call): NWPSAF_OpenFile Gen_Error_Report (17 times) NWPSAF_ReadHeaders (7 times)

NWPSAF_SatMatInv Uses Cholesky decomposition (Golub and Van Loan, 1996) to return the inverse of a symmetric positive definite matrix, A. If the the optional parameter Matrix is present, it is replaced by (Matrix).A-1 on exit and A is left u nchanged. If U is not positive definite this will be detected by the program and flagged as an error. U is assumed to be symmetric as only the upper triangle is in fact used.
Called by: NWPSAF_CloudCost NWPSAF_InitBmatrix NWPSAF_Minimize NWPSAF_Minimize_100 NWPSAF_Minimize_100ML	Calls (in order of first call): Gen_Error_Report (2 times)

NWPSAF_SetUpBackground Converts the background profile into a useable for for the minimisation code. Values are quality-controlled and the profiles may be flagged as invalid, modified (e.g., if supersaturated), or extrapolated (above model top, or to fill in data near the surface). Also the water vapour units are converted to those used in minimisation (Ln(kg/kg))
Called by: NWPSAF_ProcessData	Calls (in order of first call): NWPSAF_StratosExtrap QSAT (2 times)

NWPSAF_SetUpRetrievals This reads in the variables to be retrieved from the Retrieval.NL namelist file as explained in Section 3b.
Called by: NWPSAF_1DVar_Driver	Calls (in order of first call): NWPSAF_OpenFile Gen_Error_Report (6 times)

NWPSAF_StratosExtrap Code used to extrapolate the temperature above the model top (a very simple version is given here). This is only called when DoTExtrapolation is set to TRUE. in NWPSAF_ProcessData .
Called by: NWPSAF_SetUpBackground	Calls (in order of first call): < NONE >

NWPSAF_TranslateDataIn Transform the input observation data into a form that can be better processed by the Fastmodel_Interface and the minimisation routines. If bias correction of the brightness temperatures is to be performed it should be done here.
Called by: NWPSAF_ProcessData	Calls (in order of first call): (User supplied bias correction scheme)

NWPSAF_TranslateDataOut Output results from processing. In an operational system, this subroutine would be used for converting the format of the output data into a form suitable for use in the next processing step.
Called by: NWPSAF_ProcessData	Calls (in order of first c all): QSAT

NWPSAF_IntegerSort A sorting routine.
Called by: NWPSAF_Channellist	Calls (in order of first call): < NONE >

QSAT A routine taken from the Met Office's Unified Model to return the water vapour saturation mixing ratio given pressure and temperature.
Called by: NWPSAF_CheckIteration NWPSAF_SetUpBackground NWPSAF_TranslateDataOut	Calls (in order of first call): < NONE >

[ Return to Top ]

7. References.

English, S.J., J.R. Eyre and J.A. Smith (1999). A cloud-detection scheme for use with satellite sounding radiances in the context of data assimilation for numerical weather prediction. Q.J.R. Meteorol. Soc., 125, 2359--2378.

Golub, G.H. and C.F. Van Loan (1996). Matrix Computations. The Johns Hopkins University Press.

Matricardi, M., and R.W. Saunders (1999). A fast radiative transfer model for simulation of IASI radiances. Applied Optics, 38, 5679-5691.

Rodgers, C. D. (1976). Retrieval of atmospheric temperature and composition from remote measurements of thermal radiation. R. Geophys. Space Phys., 14, 609-624.

Rodgers, C. D. (2000). Inverse Methods for Atmospheres: Theory and Practice. World Scientific Publishing.
 

[ Return to Top ]

[ Return to Main Manual ]