Fortran module containing the routines to perform a Marquardt-Levenberg minimization. More...

Functions/Subroutines
subroutine, public	ufo_rttovonedvarcheck_minimize_ml (self, ob, r_matrix, b_matrix, b_inv, b_sigma, local_geovals, hofxdiags, rttov_simobs, profile_index, onedvar_success)
	Perform a minimization using the Marquardt-Levenberg method. More...

subroutine	ufo_rttovonedvarcheck_ml_rttov12 (self, DeltaBT, nChans, ob, H_Matrix, H_Matrix_T, nprofelements, profile_index, DeltaProfile, GuessProfile, BackProfile, b_inv, r_matrix, geovals, hofxdiags, rttov_simobs, gamma, Jold, Status)
	Updates the profile vector during the Marquardt-Levenberg minimization. More...

Detailed Description

Fortran module containing the routines to perform a Marquardt-Levenberg minimization.

Function/Subroutine Documentation

◆ ufo_rttovonedvarcheck_minimize_ml()

subroutine, public ufo_rttovonedvarcheck_minimize_ml_mod::ufo_rttovonedvarcheck_minimize_ml	(	type(ufo_rttovonedvarcheck), intent(inout)	self,
		type(ufo_rttovonedvarcheck_ob), intent(inout)	ob,
		type(ufo_rttovonedvarcheck_rsubmatrix), intent(in)	r_matrix,
		real(kind_real), dimension(:,:), intent(in)	b_matrix,
		real(kind_real), dimension(:,:), intent(in)	b_inv,
		real(kind_real), dimension(:), intent(in)	b_sigma,
		type(ufo_geovals), intent(inout)	local_geovals,
		type(ufo_geovals), intent(inout)	hofxdiags,
		type(ufo_radiancerttov), intent(inout)	rttov_simobs,
		type(ufo_rttovonedvarcheck_profindex), intent(in)	profile_index,
		logical, intent(out)	onedvar_success
	)

Perform a minimization using the Marquardt-Levenberg method.

Heritage: Ops_SatRad_MinimizeML_RTTOV12.f90

Find the most probable atmospheric state vector by minimizing a cost function through a series of iterations. if a solution exists, the iterations will converge when the iterative increments are acceptably small. A limit on the total number of iterations allowed is imposed.

Using the formulation given by Rodgers (1976) :

Delta_x = xn + (xb-xn).I' + Wn.(ym-y(xn) - H.(xb-xn))
where:
x is an atmospheric state vector, subscripted b=background,n=nth iteration
I' is a diagonal matrix with I'(J,J) = B_damped(J,J)/B_undamped(J,J) although, however, damping will no longer be used
Wn = B.Hn'.(Hn.B.Hn'+R)^-1
B is the background error covariance matrix
R is the combined forward model and ob error covariance matrix

Delta_x is checked for convergence after each iteration

The loop is exited with convergence if either of the following conditions are true, depending on whether UseJforConvergence is true or false

if UseJforConvergence is true then the change in total cost function from one iteration to the next is sufficiently small to imply convergence
if UseJforConvergence is false then the increments to the atmospheric state vector are sufficiently small to imply convergence at an acceptable solution

Either of the following two conditions will cause the 1dvar to stop and exit with an error.

The increments are sufficiently large to suppose a solution will not be found.
The maximum number of allowed iterations has been reached. In most cases, one of the above criteria will have occurred.

References:

Rodgers, Retrieval of atmospheric temperature and composition from remote measurements of thermal radiation, Rev. Geophys.Sp.Phys. 14, 1976.

Eyre, Inversion of cloudy satellite sounding radiances by nonlinear optimal estimation. I: Theory and simulation for TOVS,QJ,July 89.

Author: Met Office

Date: 09/06/2020: Created

Parameters

[in,out]	self	structure containing settings
[in,out]	ob	satellite metadata
[in]	r_matrix	observation error covariance
[in]	b_matrix	state error covariance
[in]	b_inv	inverse state error covariance
[in]	b_sigma	standard deviations of the state error covariance diagonal
[in,out]	local_geovals	model data at obs location
[in,out]	hofxdiags	model data containing the jacobian
[in,out]	rttov_simobs	rttov simulated obs object
[in]	profile_index	index array for x vector
[out]	onedvar_success	convergence flag

Definition at line 82 of file ufo_rttovonedvarcheck_minimize_ml_mod.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ufo_rttovonedvarcheck_ml_rttov12()

subroutine ufo_rttovonedvarcheck_minimize_ml_mod::ufo_rttovonedvarcheck_ml_rttov12	(	type(ufo_rttovonedvarcheck), intent(in)	self,
		real(kind_real), dimension(:), intent(in)	DeltaBT,
		integer, intent(in)	nChans,
		type(ufo_rttovonedvarcheck_ob), intent(inout)	ob,
		real(kind_real), dimension(:,:), intent(in)	H_Matrix,
		real(kind_real), dimension(:,:), intent(in)	H_Matrix_T,
		integer, intent(in)	nprofelements,
		type(ufo_rttovonedvarcheck_profindex), intent(in)	profile_index,
		real(kind_real), dimension(:), intent(inout)	DeltaProfile,
		real(kind_real), dimension(:), intent(inout)	GuessProfile,
		real(kind_real), dimension(:), intent(in)	BackProfile,
		real(kind_real), dimension(:,:), intent(in)	b_inv,
		type(ufo_rttovonedvarcheck_rsubmatrix), intent(in)	r_matrix,
		type(ufo_geovals), intent(inout)	geovals,
		type(ufo_geovals), intent(inout)	hofxdiags,
		type(ufo_radiancerttov), intent(inout)	rttov_simobs,
		real(kind_real), intent(inout)	gamma,
		real(kind_real), intent(inout)	Jold,
		integer, intent(out)	Status
	)

private

Updates the profile vector during the Marquardt-Levenberg minimization.

Heritage: Ops_SatRad_MarquardtLevenberg_RTTOV12.f90

Updates the profile vector DeltaProfile according to Rodgers (1976), Eqn. 100, extended to allow for additional cost function terms and Marquardt-Levenberg descent.

x_(n+1) = xb + U^-1.V
where
U=(B^-1 + H^T R^-1 H + J2)
V=H^T R^-1 [(ym-y(x_n))+H(x_n-xb)] - J1
and
J_extra=J0+J1.(x-xb)+(x-xb)^T.J2.(x-xb) is the additional cost function
x is an atmospheric state vector, subscripted b=background,n=nth iteration
ym is the measurement vector (i.e. observed brightness temperatures)
y(xn) is the observation vector calculated for xn
ym and y(xn) are not used individually at all, hence these are input
as a difference vector DeltaBT.
B is the background error covariance matrix
R is the combined forward model and ob error covariance matrix
H is the forward model gradient (w.r.t. xn) matrix
H' is the transpose of H

When J_extra is zero this is simply Rogers (1976), Eqn. 100.

U^-1.V is solved using Cholesky decomposition.

This routine should be used when:

The length of the observation vector is greater than the length of the state vector,
Additional cost function terms are provided and
where Newtonian minimisation is desired.