ufo_rttovonedvarcheck_minimize_ml_mod.f90

Go to the documentation of this file.

! (C) Copyright 2020 Met Office
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 
!> Fortran module containing the routines to perform a Marquardt-Levenberg 
!! minimization.
 
module ufo_rttovonedvarcheck_minimize_ml_mod
 
use kinds
use ufo_constants_mod, only: zero, ten
use fckit_log_module, only : fckit_log
use ufo_geovals_mod
use ufo_radiancerttov_mod
use ufo_rttovonedvarcheck_constants_mod
use ufo_rttovonedvarcheck_minimize_jacobian_mod
use ufo_rttovonedvarcheck_minimize_utils_mod
use ufo_rttovonedvarcheck_profindex_mod
use ufo_rttovonedvarcheck_rsubmatrix_mod
use ufo_rttovonedvarcheck_ob_mod
use ufo_rttovonedvarcheck_utils_mod
use ufo_utils_mod, only: ops_cholesky
 
implicit none
private
 
public ufo_rttovonedvarcheck_minimize_ml
 
contains
 
!-------------------------------------------------------------------------------
!> Perform a minimization using the Marquardt-Levenberg method.
!!
!! \details 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)) <br>
!!   where: <br>
!!          x is an atmospheric state vector, subscripted b=background,n=nth
!!           iteration <br>
!!          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 <br>
!!          Wn = B.Hn'.(Hn.B.Hn'+R)^-1 <br>
!!          B is the background error covariance matrix <br>
!!          R is the combined forward model and ob error covariance matrix <br>
!!
!!   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
!!
subroutine ufo_rttovonedvarcheck_minimize_ml(self,      &
                                         ob,            &
                                         r_matrix,      &
                                         b_matrix,      &
                                         b_inv,         &
                                         b_sigma,       &
                                         local_geovals, &
                                         hofxdiags,     &
                                         rttov_simobs,  &
                                         profile_index, &
                                         onedvar_success)
 
implicit none
 
type(ufo_rttovonedvarcheck), intent(inout) :: self   !< structure containing settings
type(ufo_rttovonedvarcheck_ob), intent(inout) :: ob  !< satellite metadata
type(ufo_rttovonedvarcheck_rsubmatrix), intent(in) :: r_matrix !< observation error covariance
real(kind_real), intent(in)       :: b_matrix(:,:)   !< state error covariance
real(kind_real), intent(in)       :: b_inv(:,:)      !< inverse state error covariance
real(kind_real), intent(in)       :: b_sigma(:)      !< standard deviations of the state error covariance diagonal
type(ufo_geovals), intent(inout)  :: local_geovals   !< model data at obs location
type(ufo_geovals), intent(inout)  :: hofxdiags       !< model data containing the jacobian
type(ufo_radiancerttov), intent(inout) :: rttov_simobs !< rttov simulated obs object
type(ufo_rttovonedvarcheck_profindex), intent(in) :: profile_index !< index array for x vector
logical, intent(out)              :: onedvar_success !< convergence flag
 
! Local declarations:
character(len=*), parameter     :: routinename = "ufo_rttovonedvarcheck_minimize_ml"
integer                         :: inversionstatus ! status of minimisation
logical                         :: outofrange      ! out of range flag for check iteration
logical                         :: converged       ! converged flag
logical                         :: error           ! error flag
integer                         :: iter            ! iteration counter
integer                         :: rterrorcode     ! error code for RTTOV
integer                         :: nchans          ! number of satellite channels used
integer                         :: ii              ! counter
integer                         :: nprofelements   ! number of profile elements
real(kind_real)                 :: gamma           ! steepness of descent parameter
real(kind_real)                 :: jcost           ! current cost value
real(kind_real)                 :: jcostold        ! previous iteration cost value
real(kind_real)                 :: jcostorig       ! initial cost value
real(kind_real)                 :: deltaj          ! change in cost during iteration
real(kind_real)                 :: deltajo         ! change in cost from original
 
real(kind_real), allocatable    :: oldprofile(:)
real(kind_real), allocatable    :: guessprofile(:)
real(kind_real), allocatable    :: guessprofilebefore(:)
real(kind_real), allocatable    :: backprofile(:)
real(kind_real), allocatable    :: h_matrix(:,:)
real(kind_real), allocatable    :: diffprofile(:)
real(kind_real), allocatable    :: absdiffprofile(:)
real(kind_real), allocatable    :: ydiff(:)
real(kind_real), allocatable    :: y(:)
real(kind_real), allocatable    :: y0(:)
real(kind_real), allocatable    :: out_h_matrix(:,:)
real(kind_real), allocatable    :: out_y(:)
real(kind_real)                 :: jout(3)
type(ufo_geovals)               :: geovals
 
!------------
! Setup
!------------
 
converged = .false.
onedvar_success = .false.
error = .false.
nchans = size(ob % channels_used)
gamma = 1.0e-4_kind_real
jcost = 1.0e4_kind_real
jcostold = 1.0e4_kind_real
nprofelements = profile_index % nprofelements
allocate(oldprofile(nprofelements))
allocate(guessprofile(nprofelements))
allocate(guessprofilebefore(nprofelements))
allocate(backprofile(nprofelements))
allocate(h_matrix(nchans,nprofelements))
allocate(diffprofile(nprofelements))
allocate(absdiffprofile(nprofelements))
allocate(ydiff(nchans))
allocate(y(nchans))
allocate(y0(nchans))
geovals = local_geovals
 
call fckit_log % debug("Using ML solver")
 
! Map GeovaLs to 1D-var profile using B matrix profile structure
call ufo_rttovonedvarcheck_geovals2profvec(geovals, profile_index, ob, guessprofile(:))
 
iterations: do iter = 1, self % max1DVarIterations
 
  !-------------------------
  ! 1. Generate new profile
  !-------------------------
 
  ! Save current profile
  oldprofile(:) = guessprofile(:)
 
  ! Get jacobian and new hofx
  call ufo_rttovonedvarcheck_get_jacobian(self, geovals, ob, ob % channels_used, &
                                       profile_index, guessprofile(:), &
                                       hofxdiags, rttov_simobs, y(:), h_matrix)
 
  if (iter == 1) then
    rterrorcode = 0
    diffprofile(:) = zero
    backprofile(:) = guessprofile(:)
    y0(:) = y(:)
  end if
 
  ! exit on error
  if (rterrorcode /= 0) then
    write(*,*) "Radiative transfer error"
    exit iterations
  end if
 
  ! Calculate Obs diff and initial cost
  ydiff(:) = ob % yobs(:) - y(:)
  if (iter == 1) then
    diffprofile(:) = guessprofile(:) - backprofile(:)
    call ufo_rttovonedvarcheck_costfunction(diffprofile, b_inv, ydiff, r_matrix, jout)
    jcost = jout(1)
    jcostorig = jcost
  end if
 
  !-----------------------------------------------------
  ! 1a. if UseJForConvergence is true
  !     then compute costfunction for latest profile
  !     and determine convergence using change in cost fn
  !-----------------------------------------------------
 
  if (self % UseJForConvergence) then
 
    ! exit on error
    !if (inversionStatus /= 0) exit Iterations
 
    ! check for convergence
    if (iter > 1) then
 
      if (self % JConvergenceOption == 1) then
 
        ! percentage change tested between iterations
        deltaj = abs((jcost - jcostold) / max(jcost, tiny(zero)))
 
        ! default test for checking that overall cost is getting smaller
        deltajo = -1.0_kind_real
 
      else
 
        ! absolute change tested between iterations
        deltaj = abs(jcost - jcostold)
 
        ! change between current cost and initial
        deltajo = jcost - jcostorig
 
      end if
 
      if (deltaj < self % cost_convergencefactor .and. &
          deltajo < zero)  then ! overall is cost getting smaller?
        converged = .true.
        exit iterations
      end if
 
    end if
 
  end if ! end of specific code for cost test convergence
 
  ! store cost function from previous cycle
  jcostold = jcost
 
  ! Iterate (Guess) profile vector
  call ufo_rttovonedvarcheck_ml_rttov12 ( self,            &
                                      ydiff,               &
                                      nchans,              &
                                      ob,                  &
                                      h_matrix,            &
                                      transpose(h_matrix), &
                                      nprofelements,       &
                                      profile_index,       &
                                      diffprofile,         &
                                      guessprofile,        &
                                      backprofile,         &
                                      b_inv,               &
                                      r_matrix,            &
                                      geovals,             &
                                      hofxdiags,           &
                                      rttov_simobs,        &
                                      gamma,               &
                                      jcost,               &
                                      inversionstatus)
 
  if (inversionstatus /= 0) then
    inversionstatus = 1
    write(*,*) "inversion failed"
    exit iterations
  end if
 
  !---------------------------------------------------------
  ! 2. Check new profile and transfer to forward model format
  !---------------------------------------------------------
 
  ! Check profile and constrain humidity variables
  call ufo_rttovonedvarcheck_checkiteration (self, & ! in
                                  geovals,         & ! in
                                  profile_index,   & ! in
                                  guessprofile(:), & ! inout
                                  outofrange)        ! out
 
  ! Update geovals to be the same as guess profile
  call ufo_rttovonedvarcheck_profvec2geovals(geovals, profile_index, ob, &
                                             guessprofile, self % UseQtSplitRain)
  
  ! if qtotal in retrieval vector check cloud
  ! variables for current iteration
 
  if ((.not. outofrange) .and. profile_index % qt(1) > 0) then
 
    if (iter >= self % IterNumForLWPCheck) then
 
        call ufo_rttovonedvarcheck_checkcloudyiteration( geovals, & ! in
                                              profile_index,      & ! in
                                              self % nlevels,     & ! in
                                              outofrange )          ! out
 
    end if                                                                                  
 
  end if
 
  !-------------------------------------------------
  ! 3. Check for convergence using change in profile
  !    This is performed if UseJforConvergence is false
  !-------------------------------------------------
 
  absdiffprofile(:) = zero
 
  if ((.not. outofrange) .and. (.not. self % UseJForConvergence)) then
    absdiffprofile(:) = abs(guessprofile(:) - oldprofile(:))
    if (all(absdiffprofile(:) <= b_sigma(:) * self % ConvergenceFactor)) then
      write(*,*) "Profile used for convergence"
      converged = .true.
    end if
  end if
 
  !---------------------
  ! 4. output diagnostics
  !---------------------
 
  if (self % FullDiagnostics) then
    write (*, '(A,I0)') 'Iteration', iter
    write (*, '(A)') '------------'
    write (*, '(A,L1)') 'Status: converged = ', converged
    if (outofrange) write (*, '(A)') 'exiting with bad increments'
    write (*, '(A)') 'New profile:'
    call ufo_geovals_print(geovals, 1)
    write (*, '(A)')
  end if
 
  ! exit conditions
 
  if (converged .OR. outofrange) exit iterations
 
end do iterations
 
! Pass convergence flag out
onedvar_success = converged
 
! Recalculate final cost - to make sure output when profile has not converged
call ufo_rttovonedvarcheck_costfunction(diffprofile, b_inv, ydiff, r_matrix, jout)
ob % final_cost = jout(1)
ob % niter = iter
ob % final_bt_diff = ydiff
 
! Pass output profile, final BTs and final cost out
if (converged) then
  ob % output_profile(:) = guessprofile(:)
 
  ! If lwp output required then recalculate
  if (self % Store1DVarLWP) then
    call ufo_rttovonedvarcheck_checkcloudyiteration( geovals, & ! in
                                            profile_index,    & ! in
                                            self % nlevels,   & ! in
                                            outofrange,       & ! out
                                            outlwp = ob % LWP ) ! out
  end if
 
  ! Recalculate final BTs for all channels
  allocate(out_h_matrix(size(ob % channels_all),nprofelements))
  allocate(out_y(size(ob % channels_all)))
  call ufo_rttovonedvarcheck_get_jacobian(self, geovals, ob, ob % channels_all, &
                                          profile_index, guessprofile(:), &
                                          hofxdiags, rttov_simobs, out_y(:), out_h_matrix)
  ob % output_BT(:) = out_y(:)
  deallocate(out_y)
  deallocate(out_h_matrix)
end if
 
!----------------------
! 4. output diagnostics
!----------------------
 
if (self % UseJForConvergence .and. self % FullDiagnostics) then
  write(*,'(A45,3F10.3,I5,L5)') "ML J initial, final, lowest, iter, converged = ", &
                                 jcostorig, jcost,  jcost, iter, onedvar_success
end if
 
! ----------
! Tidy up
! ----------
if (allocated(oldprofile))         deallocate(oldprofile)
if (allocated(guessprofile))       deallocate(guessprofile)
if (allocated(guessprofilebefore)) deallocate(guessprofilebefore)
if (allocated(backprofile))        deallocate(backprofile)
if (allocated(h_matrix))           deallocate(h_matrix)
if (allocated(diffprofile))        deallocate(diffprofile)
if (allocated(absdiffprofile))     deallocate(absdiffprofile)
if (allocated(ydiff))              deallocate(ydiff)
if (allocated(y))                  deallocate(y)
if (allocated(y0))                 deallocate(y0)
 
end subroutine ufo_rttovonedvarcheck_minimize_ml
 
!---------------------------------------------------------------------
!> Updates the profile vector during the Marquardt-Levenberg 
!! minimization.
!!
!! \details 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 <br>
!!   where <br>
!!         U=(B^-1 + H^T R^-1 H + J2) <br>
!!         V=H^T R^-1 [(ym-y(x_n))+H(x_n-xb)] - J1 <br>
!!   and   <br>
!!         J_extra=J0+J1.(x-xb)+(x-xb)^T.J2.(x-xb) is the additional cost
!!         function <br>
!!         x is an atmospheric state vector, subscripted b=background,n=nth
!!         iteration <br>
!!         ym is the measurement vector (i.e. observed brightness temperatures) <br>
!!         y(xn) is the observation vector calculated for xn <br>
!!         ym and y(xn) are not used individually at all, hence these are input <br>
!!         as a difference vector DeltaBT. <br>
!!         B is the background error covariance matrix <br>
!!         R is the combined forward model and ob error covariance matrix <br>
!!         H is the forward model gradient (w.r.t. xn) matrix <br>
!!         H' is the transpose of H <br>
!!
!!   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.
!!
!! References:
!!
!!   Rodgers, Retrieval of atmospheric temperature and composition from
!!   remote measurements of thermal radiation, Rev. Geophys.Sp.Phys. 14, 1976.
!!
!!   Rodgers, Inverse Methods for Atmospheres: Theory and Practice.  World
!!            Scientific Publishing, 2000.
!!
!! \author Met Office
!!
!! \date 09/06/2020: Created
!!
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)
 
implicit none
 
! subroutine arguments:
type(ufo_rttovonedvarcheck), intent(in) :: self            !< structure containing settings
real(kind_real), intent(in)             :: DeltaBT(:)      !< y-y(x)
integer, intent(in)                     :: nChans          !< number of channels
type(ufo_rttovonedvarcheck_ob), intent(inout) :: ob        !< satellite metadata
real(kind_real), intent(in)             :: H_Matrix(:,:)   !< Jacobian
real(kind_real), intent(in)             :: H_Matrix_T(:,:) !< transpose  of the Jacobian
integer, intent(in)                     :: nprofelements   !< number of profile elements
type(ufo_rttovonedvarcheck_profindex), intent(in) :: profile_index !< index array for x vector
real(kind_real), intent(inout)          :: DeltaProfile(:) !< x_(n+1) - xb
real(kind_real), intent(inout)          :: GuessProfile(:) !< x_(n+1)
real(kind_real), intent(in)             :: BackProfile(:)  !< xb
real(kind_real), intent(in)             :: b_inv(:,:)      !< inverse of the state error covariance
type(ufo_rttovonedvarcheck_rsubmatrix), intent(in) :: r_matrix !< observation error covariance
type(ufo_geovals), intent(inout)        :: geovals         !< model data at obs location
type(ufo_geovals), intent(inout)        :: hofxdiags       !< model data containing the jacobian
type(ufo_radiancerttov), intent(inout)  :: rttov_simobs
real(kind_real), intent(inout)          :: gamma           !< steepness of descent parameter
real(kind_real), intent(inout)          :: Jold            !< previous steps cost which gets update by good step
integer, intent(out)                    :: Status          !< code to capture failed Cholesky decomposition
 
! Local declarations:
character(len=*), parameter         :: RoutineName = 'ufo_rttovonedvarcheck_ML_RTTOV12'
integer                             :: i
integer                             :: Iter_ML             ! M-L iteration count
real(kind_real)                     :: JOut(3)             ! J, Jb, Jo
real(kind_real)                     :: JCost               ! Cost function,
logical                             :: OutOfRange
real(kind_real)                     :: HTR(nprofelements, nChans)              ! Scratch vector
real(kind_real)                     :: U(nprofelements, nprofelements)         ! U = H.B.H^T + R
real(kind_real)                     :: V(nprofelements)                        ! V = (y-y(x_n))-H^T(xb-x_n)
real(kind_real)                     :: USave(nprofelements, nprofelements)     ! U = H.B.H^T + R
real(kind_real)                     :: VSave(nprofelements)                    ! V = (y-y(x_n))-H^T(xb-x_n)
real(kind_real)                     :: New_DeltaProfile(nprofelements)
real(kind_real)                     :: Tau_Surf(nchans)                ! Transmittance from surface
real(kind_real)                     :: BriTemp(nchans)                 ! Forward modelled brightness temperatures
real(kind_real)                     :: Ydiff(nchans)
real(kind_real)                     :: Emiss(nchans)
real(kind_real)                     :: H_matrix_tmp(nchans,nprofelements)
logical                             :: CalcEmiss(nchans)
integer                             :: StatusOK = 0
integer                             :: RTStatus = 0
 
status = statusok
 
!---------------------------------------------------------------------------
! 1. Calculate HTR-1 for the three allowed forms of R matrix. if required, the
!    matrix is tested to determine whether it is stored as an inverse and
!    inverted if not.
!---------------------------------------------------------------------------
 
!HTR(-1) = matmul(H_matrix_T, R_inverse)
call r_matrix % multiply_inverse_matrix(h_matrix_t, htr)
 
!---------------------------------------------------------------------------
! 2. Calculate U and V
!---------------------------------------------------------------------------
 
u = matmul(htr, h_matrix)
v = matmul(htr, deltabt)
v = v + matmul(u, deltaprofile)
 
!---------------------------------------------------------------------------
! 3. Add on inverse of B-matrix to U.
!---------------------------------------------------------------------------
 
u = u + b_inv
 
!---------------------------------------------------------------------------
! 5. Start Marquardt-Levenberg loop.
!    Search for a value of gamma that reduces the cost function.
!---------------------------------------------------------------------------
 
jcost = 1.0e10_kind_real
gamma = gamma / ten
usave = u
vsave = v
iter_ml = 0
 
descentloop : do while (jcost > jold .and.              &
                        iter_ml < self % MaxMLIterations .and. &
                        gamma <= 1.0e25_kind_real)
 
  iter_ml = iter_ml + 1
 
  !------------------------------------------------------------------------
  ! 5.1 Add on Marquardt-Levenberg terms to U and V.
  !------------------------------------------------------------------------
 
  gamma = gamma * ten
 
  do i = 1, nprofelements
    u(i,i) = usave(i,i) + gamma
  end do
  v = vsave + gamma * deltaprofile
 
  !------------------------------------------------------------------------
  ! 5.2 Calculate new profile increment.
  !------------------------------------------------------------------------
 
  call ops_cholesky (u,                &
                     v,                &
                     nprofelements,    &
                     new_deltaprofile, &
                     status)
  if (status /= 0) then
     write(*,*) 'Error in Cholesky decomposition'
  end if
 
  !------------------------------------------------------------------------
  ! 5.3. Calculate the radiances for the new profile.
  !------------------------------------------------------------------------
 
  guessprofile(:) = backprofile(:) + new_deltaprofile(:)
 
  call ufo_rttovonedvarcheck_checkiteration (self, & ! in
                                  geovals,         & ! in
                                  profile_index,   & ! in
                                  guessprofile(:), & ! inout
                                  outofrange)        ! out
 
  if (.not. outofrange) then
    call ufo_rttovonedvarcheck_profvec2geovals(geovals, profile_index, ob, &
                                               guessprofile, self % UseQtSplitRain)
 
    !Emiss(1:nchans) = RTprof_Guess % Emissivity(Channels(1:nchans))
    !CalcEmiss(1:nchans) = RTprof_Guess % CalcEmiss(Channels(1:nchans))
 
   ! Get Jabobian and new hofx
   call ufo_rttovonedvarcheck_get_jacobian(self, geovals, ob, ob % channels_used, &
                                           profile_index, guessprofile(:), &
                                           hofxdiags, rttov_simobs, britemp(:), h_matrix_tmp)
 
    !------------------------------------------------------------------------
    ! 5.6. Calculate the new cost function.
    !------------------------------------------------------------------------
 
    if (rtstatus > 0) then
 
      ! In case of RT error, set cost to large value and
      ! continue with next iteration
 
      jcost = 1.0e10_kind_real
    else
      deltaprofile(:) = guessprofile(:) - backprofile(:)
      ydiff(:) = ob % yobs(:) - britemp(:)
      call ufo_rttovonedvarcheck_costfunction(deltaprofile, b_inv, ydiff, r_matrix, jout)
      jcost = jout(1)
      if (status /= statusok) goto 9999
    end if
 
  else
 
    ! if profile out of range, set cost to large value and
    ! continue with next iteration
 
    jcost = 1.0e10_kind_real
 
  end if
 
end do descentloop
 
deltaprofile = new_deltaprofile
gamma = gamma / ten
jold = jcost
 
9999 continue
 
end subroutine ufo_rttovonedvarcheck_ml_rttov12
 
end module ufo_rttovonedvarcheck_minimize_ml_mod