ufo_gnssroonedvarcheck_do1dvar_mod.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! (C) British Crown 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 for gnssro bending angle Met Office forward operator
 
module ufo_gnssroonedvarcheck_do1dvar_mod
 
use kinds, only: kind_real
use missing_values_mod, only: missing_value
use fckit_log_module, only: fckit_log
 
private
public :: ops_gpsro_do1dvar_ba
 
contains
 
!-------------------------------------------------------------------------------
! Find a solution to the satellite sounding inverse problem
!-------------------------------------------------------------------------------
SUBROUTINE ops_gpsro_do1dvar_ba (nlevp,                  &
                                 nlevq,                  &
                                 BM1,                    &
                                 Bsig,                   &
                                 Back,                   &
                                 Ob,                     &
                                 GPSRO_pseudo_ops,       &
                                 GPSRO_vert_interp_ops,  &
                                 GPSRO_min_temp_grad,    &
                                 GPSRO_cost_funct_test,  &   ! Threshold value for the cost function convergence test
                                 GPSRO_y_test,           &   ! Threshold value for the yobs-ysol tes
                                 GPSRO_n_iteration_test, &   ! Maximum number of iterations
                                 GPSRO_Delta_factor,     &   ! Delta
                                 GPSRO_Delta_ct2,        &   ! Delta observations
                                 GPSRO_OB_test,          &   ! Threshold value for the O-B test
                                 capsupersat,            &
                                 BAerr,                  &
                                 Tb,                     &
                                 Ts,                     &
                                 O_Bdiff,                &
                                 DFS)
 
use ufo_gnssroonedvarcheck_utils_mod, only: &
    singlebg_type,             &
    singleob_type
 
use ufo_gnssroonedvarcheck_rootsolv_mod, only: &
    ops_gpsro_rootsolv_ba
 
use ufo_utils_refractivity_calculator, only: &
    ufo_calculate_refractivity
 
IMPLICIT NONE
 
! Subroutine arguments:
INTEGER, INTENT(IN)                 :: nlevp
INTEGER, INTENT(IN)                 :: nlevq
REAL(kind_real), INTENT(IN)         :: bm1(:,:)
REAL(kind_real), INTENT(IN)         :: bsig(:)
TYPE (singlebg_type), INTENT(INOUT) :: back
TYPE (singleob_type), INTENT(INOUT) :: ob
LOGICAL, INTENT(IN)                 :: gpsro_pseudo_ops
LOGICAL, INTENT(IN)                 :: gpsro_vert_interp_ops
REAL(kind_real), INTENT(IN)         :: gpsro_min_temp_grad
REAL(kind_real), INTENT(IN)         :: gpsro_cost_funct_test
REAL(kind_real), INTENT(IN)         :: gpsro_y_test
INTEGER, INTENT(IN)                 :: gpsro_n_iteration_test
REAL(kind_real), INTENT(IN)         :: gpsro_delta_ct2
REAL(kind_real), INTENT(IN)         :: gpsro_delta_factor
REAL(kind_real), INTENT(IN)         :: gpsro_ob_test
LOGICAL, INTENT(IN)                 :: capsupersat
LOGICAL, INTENT(OUT)                :: baerr
REAL(kind_real), INTENT(INOUT)      :: tb(nlevq)
REAL(kind_real), INTENT(INOUT)      :: ts(nlevq)
REAL(kind_real), INTENT(INOUT)      :: o_bdiff  ! measure of O-B for whole profile
REAL(kind_real), INTENT(INOUT)      :: dfs      ! measure of degrees of freedom of signal for whole profile
 
! Local parameters
CHARACTER(len=*), PARAMETER         :: routinename = "Ops_GPSRO_Do1DVar_BA"
INTEGER, PARAMETER                  :: max_string = 800
 
! Local variables
INTEGER                             :: nobs        ! size of the 1DVar observation vector
INTEGER                             :: i
INTEGER                             :: j
INTEGER                             :: it
LOGICAL                             :: om1_error
LOGICAL                             :: converged
LOGICAL                             :: do1dvar_error
LOGICAL                             :: ran_iteration
REAL(kind_real)                     :: j_pen
REAL(kind_real)                     :: xb(nlevp+nlevq)                 ! background profile used in the 1D-Var
REAL(kind_real)                     :: x(nlevp+nlevq)                  ! 1Dvar solution profile
REAL(kind_real)                     :: amat(nlevp+nlevq,nlevp+nlevq)        ! solultion error cov matrix
REAL(kind_real), ALLOCATABLE        :: zobs(:)
REAL(kind_real), ALLOCATABLE        :: yobs(:)
REAL(kind_real), ALLOCATABLE        :: yb(:)
REAL(kind_real), ALLOCATABLE        :: ycalc(:)
REAL(kind_real), ALLOCATABLE        :: om1(:,:)
INTEGER, ALLOCATABLE                :: index_packed(:)
REAL(kind_real), ALLOCATABLE        :: model_heights(:)           ! Heights of model and pseudo-levels
REAL(kind_real), ALLOCATABLE        :: refractivity(:)            ! Refractivity on model and pseudo_levels
INTEGER                             :: nreflevels                 ! Number of levs to calculate ref on
CHARACTER(LEN=max_string)           :: message
 
REAL(kind_real) :: temp_rad_curv     ! Temporary store of the earth's radius of curvature
REAL(kind_real) :: temp_latitude     ! Temporary store of the observation's latitude
REAL(kind_real) :: temp_undulation   ! Temporary store of the undulation
 
!--------------
! 1. Initialise
!--------------
 
do1dvar_error = .false.
baerr = .false.
ran_iteration = .false.
om1_error = .false.
 
! Set all the PGE values to gross error
 
ob % BendingAngle(:) % PGEFinal = 1.0
 
! Calculate refractivity on theta levels, to find appropriate
! impact height vertical range
 
CALL ufo_calculate_refractivity (nlevp,                  &
                                 nlevq,                  &
                                 back % za,              &
                                 back % zb,              &
                                 back % p,               &
                                 back % q,               &
                                 gpsro_pseudo_ops,       &
                                 gpsro_vert_interp_ops,  &
                                 gpsro_min_temp_grad,    &
                                 baerr,                  &
                                 nreflevels,             &
                                 refractivity,           &
                                 model_heights)
 
! Set the background vector
xb(1:nlevp) = 1.0e-2 * back % p(:)          ! in hPa
xb(nlevp + 1:nlevp+nlevq) = 1.0e3 * back % q(:)    ! in g/kg
 
! Set size of obs vector used in 1D- Var
 
nobs = count(ob % BendingAngle(:) % value /= missing_value(ob % BendingAngle(1) % value) .AND. & ! not missing bending angle
              ob % ImpactParam(:) % value /= missing_value(ob % ImpactParam(1) % value)   .AND. & ! not missing impact parameter
              ob % qc_flags(:) == 0)
 
WRITE (message, '(A,I0)') 'size of input obs vector ', SIZE (ob % BendingAngle(:) % value)
CALL fckit_log % info(message)
WRITE (message, '(A,I0)') 'size of packed obs vector ', nobs
CALL fckit_log % info(message)
 
! Only continue if we have some observations to process
IF (nobs > 0) THEN
 
  ! calculate an array of indices of the packed elements
 
  ALLOCATE (index_packed(nobs))                           ! allocate the packed index vector
  index_packed = missing_value(index_packed(1))           ! initialise
 
  ! Allocate arrays used in 1D-Var after test to stop allocating size nobs=0
  ALLOCATE (om1(nobs,nobs))
  ALLOCATE (yobs(nobs))
  ALLOCATE (zobs(nobs))
  ALLOCATE (yb(nobs))
  ALLOCATE (ycalc(nobs))
 
  ! Pack observation arrays for valid values
  ! Note: This hard-codes the R-matrix to be diagonal, since that is all that
  ! is currently available in JEDI.  This will need to be revisted once the
  ! full capability is available.
 
  om1 = 0
  j = 1
  DO i = 1, SIZE (ob % BendingAngle(:) % Value)
    IF (ob % BendingAngle(i) % Value /= missing_value(ob % BendingAngle(i) % Value) .AND. &
        ob % ImpactParam(i) % value /= missing_value(ob % ImpactParam(i) % value)   .AND. &
        ob % qc_flags(i) == 0) THEN
      index_packed(j) = i
      zobs(j) = ob % ImpactParam(i) % value
      yobs(j) = ob % BendingAngle(i) % Value
      om1(j,j) = (ob % BendingAngle(i) % oberr)**(-2)
      j = j + 1
    END IF
  END DO
 
  !-----------------------------------------------
  ! 2. If no errors so far, call the 1DVar routine
  !-----------------------------------------------
 
  IF (all(zobs(:) /= missing_value(zobs(1))) .AND.  &
      all(yobs(:) /= missing_value(yobs(1))) .AND.  &
      .NOT. om1_error) THEN
 
    temp_rad_curv = ob % RO_Rad_Curv % Value
    temp_latitude = ob % Latitude
    temp_undulation = ob % RO_geoid_und % value
    CALL ops_gpsro_rootsolv_ba (nlevp+nlevq,               &    ! size of state vector
                                nlevp,                     &    ! no. of press. levels
                                nlevq,                     &    ! no. of theta levels
                                nlevq,                     &    ! no. of q levels
                                nobs,                      &    ! no of obs
                                back % za,                 &    ! height of rho levels
                                back % zb,                 &    ! height of theta levels
                                xb,                        &    ! background vector
                                yobs,                      &    ! ob. vector
                                zobs,                      &    ! ob. impact parameters
                                bsig,                      &    ! standard dev. of B errors
                                bm1,                       &    ! Inverse Back. cov matrix
                                om1,                       &    ! Inverse Ob cov matrix
                                it,                        &    ! no of iterations
                                x,                         &    ! solution vector
                                yb,                        &    ! obs at first guess
                                ycalc,                     &    ! obs at solution
                                j_pen,                     &    ! penalty value
                                amat,                      &    ! solution cov. matrix
                                converged,                 &    ! convergence flag
                                baerr,                     &    ! error flag
                                do1dvar_error,             &    ! error flag
                                gpsro_n_iteration_test,    &    !
                                gpsro_delta_factor,        &    !
                                gpsro_delta_ct2,           &    !
                                gpsro_pseudo_ops,          &
                                gpsro_vert_interp_ops,     &
                                gpsro_min_temp_grad,       &
                                capsupersat,               &
                                o_bdiff,                   &    ! observed -background BA value
                                temp_rad_curv,             &    ! Radius of curvature of ellipsoid
                                temp_latitude,             &    ! Latitude of occ
                                temp_undulation,           &    ! geoid undulation
                                tb,                        &
                                ts,                        &
                                dfs)
    ran_iteration = .true.
  ELSE
 
    do1dvar_error = .true.
 
    ob % BendingAngle(:) % PGEFinal = 0.99
 
  END IF
 
  IF (.NOT. do1dvar_error) THEN
 
    ! store iteration and cost
 
    ob % Niter = it
 
    ob % Jcost = 2.0_kind_real * j_pen / real(nobs)
 
    ! map the retrieval information back into the ob structures
 
    ob % p(:) % Value = 1.0e2 * x(1:nlevp)
 
    DO i = 1, nlevp
 
      ob % p(i) % ObErr = 1.0e2 * sqrt(amat(i,i))   ! error. est. from cov
 
    END DO
 
    ob % q(:) % Value = 1.0e-3 * x(nlevp + 1:nlevp+nlevq)
 
    DO i = 1, nlevq
 
      j = i + nlevp
 
      ob % q(i) % ObErr = 1.0e-3 *  sqrt(amat(j,j))   ! error est. from cov
 
    END DO
 
    ! Bending angle calculated with solution
 
    ob % SolutBendingAngle(index_packed) = ycalc(1:nobs)
 
!    ! store the bending angle calculated from the background
 
!    Back % BendingAngle(index_packed) = yb(1:nobs)
 
    ! PROBABILTY OF GROSS ERROR. Use the cost function value.
 
    IF (j_pen > gpsro_cost_funct_test * real(nobs, kind=kind_real)) THEN
 
      ! the cost function for profile is too high - GROSS ERROR -
      ! set all pge's in profile to 0.8
 
      ob % BendingAngle(:) % PGEFinal = 0.8
 
    ELSE
 
      ! For each value in profile, estimate probability of gross error
      ! from difference between solution and observed value
 
      DO i = 1,nobs
 
        IF (abs(ob % BendingAngle(index_packed(i)) % Value - ob % SolutBendingAngle(index_packed(i))) < &
              (gpsro_y_test * ob % BendingAngle(index_packed(i)) % ObErr)) THEN
          ob % BendingAngle(index_packed(i)) % PGEFinal = 0.1
        ELSE
          ob % BendingAngle(index_packed(i)) % PGEFinal = 0.7
        END IF
 
      END DO
 
    END IF
 
    ! if the initial 2J/m value exceeds read-in value then flag
    IF  (o_bdiff > gpsro_ob_test) THEN
      ob % BendingAngle(:) % PGEFinal = 0.6
    END IF
 
    ! check for the BAerr being set
    IF (baerr) THEN
      ob % BendingAngle(:) % PGEFinal = 0.58     ! flag BAerr
    END IF
 
  ELSE  ! Do 1D-Var_error
 
    IF (ran_iteration) THEN
      ob % BendingAngle(:) % PGEFinal = 0.9     ! flag lack of convergence
 
      ! Bending angle calculated with solution
      ob % SolutBendingAngle(index_packed) = ycalc(1:nobs)
 
      ! map the retrieval information back into the ob structures
      ob % p(:) % Value = 1.0e2 * x(1:nlevp)
      ob % q(:) % Value = 1.0e-3 * x(nlevp + 1:nlevp+nlevq)
 
      ! store iterations and cost
      ob % Niter = it
      ob % Jcost = 2.0_kind_real * j_pen / nobs
 
    ELSE IF (om1_error) THEN
 
      ob % BendingAngle(:) % PGEFinal = 0.85    ! Flag error in getting
                                                ! observation error inverse
    END IF
 
  END IF
 
ELSE
  IF (nobs <= 10) THEN
    WRITE (message, '(A)') 'nobs is less than 10: exit Ops_GPSRO_Do1DVar_BA'
    CALL fckit_log % info(message)
    ob % BendingAngle(:) % PGEFinal = 0.55     ! flag lack of observation data
  END IF
 
  IF (baerr) THEN
    WRITE (message, '(A)') 'Error in Ops_Refractivity: exit Ops_GPSRO_Do1DVar_BA'
    CALL fckit_log % info(message)
    ob % BendingAngle(:) % PGEFinal = 0.58     ! flag BAerr
  END IF
END IF
 
END SUBROUTINE ops_gpsro_do1dvar_ba
 
end module ufo_gnssroonedvarcheck_do1dvar_mod