gnssro_mod_obserror.F90

Go to the documentation of this file.

!==========================================================================
module gnssro_mod_obserror
!==========================================================================
 
use kinds
use gnssro_mod_constants
use ufo_roobserror_utils_mod
use fckit_log_module, only: fckit_log
 
contains
subroutine bending_angle_obserr_ecmwf(obsImpH, obsValue, nobs,  obsErr, QCflags, missing)
implicit none
integer,                         intent(in)  :: nobs
real(kind_real), dimension(nobs),intent(in)  :: obsImpH, obsValue
integer(c_int),  dimension(nobs),intent(in)  :: QCflags(:)
real(kind_real), dimension(nobs),intent(out) :: obsErr
real(kind_real)                  :: H_km, missing
integer :: i
 
obserr = missing
 
do i = 1, nobs
   
if (qcflags(i) .eq. 0) then
 
    h_km  = obsimph(i)/1000.0_kind_real
 
    if ( h_km <= 10.0 ) then
       obserr(i) = (h_km*1.25 + (10-h_km)*20)/10.0
       obserr(i) = obserr(i)/100.0*obsvalue(i)
    else if ( h_km > 10.0 .and. h_km <= 32.0 ) then
      obserr(i) = 1.25/100.0*obsvalue(i) 
    else
      obserr(i) = 3.0*1e-6
    end if
end if
 
end do
 
end subroutine bending_angle_obserr_ecmwf
!----------------------------------------
 
subroutine bending_angle_obserr_nrl(obsLat, obsImpH, obsValue, nobs,  obsErr, QCflags, missing)
implicit none
integer,                         intent(in)  :: nobs
real(kind_real), dimension(nobs),intent(in)  :: obsImpH, obsValue, obsLat
integer(c_int),  dimension(nobs),intent(in)  :: QCflags(:)
real(kind_real), dimension(nobs),intent(out) :: obsErr
real(kind_real):: H_m, missing
real(kind_real):: lat_in_rad,trop_proxy,damping_factor,errfac
real(kind_real):: max_sfc_error, min_ba_error
 
integer :: i
 
obserr = missing 
max_sfc_error = 20.0   ! %
min_ba_error  = 1.25   ! %
 
do i = 1, nobs
 
if (qcflags(i) .eq. 0) then
   
   h_m  = obsimph(i)
   lat_in_rad = deg2rad * obslat(i)
   trop_proxy = 8666.66 + 3333.33*cos(2.0*lat_in_rad)
   damping_factor = 0.66 + cos(lat_in_rad)/3.0
   errfac = max_sfc_error*damping_factor*(trop_proxy - h_m)/trop_proxy
   errfac = max(min_ba_error, errfac)
   obserr(i) =  max(obsvalue(i)*errfac/100.0, 3.0*1e-6) ! noise floor at top of RO profile
 
end if
 
end do
 
end subroutine bending_angle_obserr_nrl
!--------------------------------------
 
subroutine  bending_angle_obserr_nbam(obsLat, obsImpH, obsSaid, nobs, obsErr, QCflags, missing)
implicit none
integer,                         intent(in)  :: nobs
real(kind_real), dimension(nobs),intent(in)  :: obsImpH, obsLat
integer(c_int),  dimension(nobs),intent(in)  :: obsSaid, QCflags(:)
real(kind_real), dimension(nobs),intent(out) ::  obsErr
real(kind_real)                 :: H_km, missing
 
integer :: i
 
obserr = missing
 
do i = 1, nobs
 
if (qcflags(i) .eq. 0) then
 
   h_km  = obsimph(i)/1000.0_kind_real
   if( (obssaid(i)==41).or.(obssaid(i)==722).or.(obssaid(i)==723).or.   &
       (obssaid(i)==4).or.(obssaid(i)==42).or.(obssaid(i)==3).or.       &
       (obssaid(i)==5).or.(obssaid(i)==821.or.(obssaid(i)==421)).or.    &
       (obssaid(i)==440).or.(obssaid(i)==43)) then
       if( abs(obslat(i))>= 40.00 ) then
         if(h_km>12.0) then
           obserr(i)=0.19032 +0.287535 *h_km-0.00260813*h_km**2
         else
           obserr(i)=-3.20978 +1.26964 *h_km-0.0622538 *h_km**2
         endif
       else
         if(h_km>18.) then
           obserr(i)=-1.87788 +0.354718 *h_km-0.00313189 *h_km**2
         else
           obserr(i)=-2.41024 +0.806594 *h_km-0.027257 *h_km**2
         endif
       endif
 
   else !!!! CDAAC processing
     if( abs(obslat(i))>= 40.00 ) then
       if ( h_km>12.00 ) then
          obserr(i)=-0.685627 +0.377174 *h_km-0.00421934 *h_km**2
       else
          obserr(i)=-3.27737 +1.20003 *h_km-0.0558024 *h_km**2
       endif
     else
       if( h_km>18.00 ) then
          obserr(i)=-2.73867 +0.447663 *h_km-0.00475603 *h_km**2
       else
          obserr(i)=-3.45303 +0.908216 *h_km-0.0293331 *h_km**2
       endif
     endif
     obserr(i) = 0.001 /abs(exp(obserr(i)))
 
  endif
 
end if
 
end do
 
end subroutine bending_angle_obserr_nbam
!---------------------------------------
 
subroutine refractivity_obserr_nbam(obsLat, obsZ, nobs, obsErr, QCflags,missing)
implicit none
integer,                         intent(in)  :: nobs
real(kind_real), dimension(nobs),intent(in)  :: obsLat, obsZ
real(kind_real), dimension(nobs),intent(out) :: obsErr
integer(c_int),  dimension(nobs),intent(in)  :: QCflags(:)
real(kind_real)                   :: H_km, missing
 
integer :: i
 
obserr = missing
 
do i = 1, nobs
 
  if (qcflags(i) .eq. 0) then
     h_km  = obsz(i)/1000.0_kind_real
     if( abs(obslat(i))>= 20.0 ) then
         obserr(i)=-1.321+0.341*h_km-0.005*h_km**2
     else
       if(h_km > 10.0) then
          obserr(i)=2.013-0.060*h_km+0.0045*h_km**2
       else
          obserr(i)=-1.18+0.058*h_km+0.025*h_km**2
       endif
     endif
     obserr(i) = 1.0_kind_real/abs(exp(obserr(i)))
  end if
end do
 
end subroutine refractivity_obserr_nbam
 
 
subroutine gnssro_obserr_avtemp(nobs, n_horiz, rmatrix_filename, obsSatid, obsOrigC, nlevs, &
                                air_temperature, geopotential_height, obsZ, obsValue, obsErr, &
                                QCflags, missing)
 
implicit none
 
! Subroutine arguments
integer, intent(in)          :: nobs                     ! Number of observations
integer, intent(in)          :: n_horiz                  ! Number of geovals per observation
character(len=*), intent(in) :: rmatrix_filename         ! Name of the R-matrix file
integer, intent(in)          :: obsSatid(:)              ! Satellite identifier
integer, intent(in)          :: obsOrigC(:)              ! Originating centre number
integer, intent(in)          :: nlevs                    ! Number of model levels
real, intent(in)             :: air_temperature(:,:)     ! Temperature of the model background
real, intent(in)             :: geopotential_height(:,:) ! Geopotential height of the model levels
real(kind_real), intent(in)  :: obsZ(:)                  ! Height of the observation
real(kind_real), intent(in)  :: obsValue(:)              ! The observed value
real(kind_real), intent(out) :: obsErr(:)                ! The calculated observation error (uncertainty)
integer(c_int),  intent(in)  :: QCflags(:)               ! Quality control flags for the observations
real(kind_real), intent(in)  :: missing                  ! Missing value indicator
 
! Local parameters
integer, parameter :: Rmax_num = 1000                    ! Max number of R matrices to be read in
 
! Local variables
type(rmatrix_type), allocatable :: Rmatrix_list(:) ! List of all the R matrices to use
type(rmatrix_type) :: Rmatrix                     ! The chosen R matrix
real(kind_real) :: frac_err                       ! Fractional observation error
real :: av_temp
integer :: npoints
integer :: ilev
integer :: R_num_sats                             ! Actual number of R-matrices read in
character(len=200) :: Message                     ! Message to be output
integer :: iob                                    ! Loop variable, observation number
integer :: igeoval                                ! Loop variable, geoval number
integer :: iheight                                ! Loop variable, height in profile
 
! Read in R matrix data
CALL ufo_roobserror_getrmatrix(rmax_num,         &  ! Max number of R matrices to read in
                               rmatrix_filename, &  ! The name of the file to be read in
                               rmatrix_list,     &  ! List of all R matrices to use
                               r_num_sats)          ! Number of R matrices read in
!--------------------------------------------------------
! Choose the R-matrix values.  We use different code if we are passed
! a matrix which depends on latitude or average temperature
!--------------------------------------------------------
 
do iob = 1, nobs
  if (qcflags(iob) .eq. 0) then
 
    IF (rmatrix_list(1) % av_temp > 0) THEN
      !--------------------------------------------------------
      ! Choose R matrix depending on satid, origctr and the average
      ! background temperature between the surface and 20km
      !--------------------------------------------------------
 
      ! Calculate the average troposphere temperature for this profile
 
      av_temp = 0
      npoints = 0
      igeoval = (iob-1) * n_horiz + (n_horiz + 1) / 2
      DO ilev = 1, nlevs
        IF (geopotential_height(igeoval, ilev) < rmatrix_list(1) % max_height) THEN
          av_temp = av_temp + air_temperature(igeoval, ilev)
          npoints = npoints + 1
        END IF
      END DO
 
      IF (npoints > 0) THEN
        av_temp = av_temp / npoints
      ELSE
        av_temp = missing
      END IF
 
      ! Find the observation error matrix which best matches the average
      ! temperature we found
 
      CALL ufo_roobserror_interpolate_rmatrix(obssatid(iob),   &
                                              obsorigc(iob),   &
                                              av_temp,         &
                                              r_num_sats,      &
                                              rmatrix_list,    &
                                              rmatrix)
 
    ELSE
      WRITE (message, '(2A)') "RMatrices must have positive average ", &
                                   "temperature"
      CALL abor1_ftn(message)
    END IF
 
    do iheight = 1, rmatrix % num_heights - 1
      if (obsz(iob) < rmatrix % height(iheight + 1)) then
        exit
      end if
    end do
 
    ! Fractional error
    frac_err = rmatrix % frac_err(iheight) + &
               (rmatrix % frac_err(iheight + 1) - rmatrix % frac_err(iheight)) * &
               (obsz(iob) - rmatrix % height(iheight)) / &
               (rmatrix % height(iheight + 1) - rmatrix % height(iheight))
 
    WRITE(message,'(A,I8,2F16.4,2E26.8)') 'Result', iob, obsz(iob), frac_err, &
        obserr(iob), max(frac_err * obsvalue(iob), rmatrix % min_error)
    CALL fckit_log % debug(message)
 
    ! Standard deviation
    obserr(iob) = max(frac_err * obsvalue(iob), rmatrix % min_error)
 
  else
    obserr(iob) = missing
  end if
end do
 
end subroutine gnssro_obserr_avtemp
 
 
subroutine gnssro_obserr_latitude(nobs, rmatrix_filename, obsSatid, obsOrigC, obsLat, obsZ, obsValue, obsErr, QCflags, missing)
 
implicit none
 
! Subroutine arguments
integer, intent(in)          :: nobs              ! Number of observations
character(len=*), intent(in) :: rmatrix_filename  ! Name of the R-matrix file
integer, intent(in)          :: obsSatid(:)       ! Satellite identifier
integer, intent(in)          :: obsOrigC(:)       ! Originating centre number
real(kind_real), intent(in)  :: obsLat(:)         ! Latitude of the observation
real(kind_real), intent(in)  :: obsZ(:)           ! Height of the observation
real(kind_real), intent(in)  :: obsValue(:)       ! The observed value
real(kind_real), intent(out) :: obsErr(:)         ! The calculated observation error (uncertainty)
integer(c_int),  intent(in)  :: QCflags(:)        ! Quality control flags for the observations
real(kind_real), intent(in)  :: missing           ! Missing value indicator
 
! Local parameters
integer, parameter :: Rmax_num = 1000             ! Max number of R matrices to be read in
 
! Local variables
type(rmatrix_type), allocatable :: Rmatrix_list(:) ! List of all the R matrices to use
type(rmatrix_type) :: Rmatrix                     ! The chosen R matrix
real(kind_real) :: frac_err                       ! Fractional observation error
integer :: R_num_sats                             ! Actual number of R-matrices read in
character(len=200) :: Message                     ! Message to be output
integer :: iob                                    ! Loop variable, observation number
integer :: iheight                                ! Loop variable, height in profile
 
! Read in R matrix data
CALL ufo_roobserror_getrmatrix(rmax_num,         &  ! Max number of R matrices to read in
                               rmatrix_filename, &  ! The name of the file to be read in
                               rmatrix_list,     &  ! List of all R matrices to use
                               r_num_sats)          ! Number of R matrices read in
 
!--------------------------------------------------------
! Choose R matrix depending on satid, origctr and latitude
! No interpolation between matrices to match old code
!--------------------------------------------------------
 
do iob = 1, nobs
  if (qcflags(iob) .eq. 0) then
    IF (rmatrix_list(1) % latitude > 0) THEN
      CALL ufo_roobserror_findnearest_rmatrix(obssatid(iob),  &
                                              obsorigc(iob),  &
                                              obslat(iob),    &
                                              r_num_sats,     &
                                              rmatrix_list,   &
                                              rmatrix)
    ELSE
      WRITE (message, '(2A)') "RMatrices must have positive average ", &
                                   "temperature or latitude"
      CALL abor1_ftn(message)
    END IF
 
    do iheight = 1, rmatrix % num_heights
      if (obsz(iob) < rmatrix % height(iheight)) then
        exit
      end if
    end do
 
    ! Calculate fractional error
    if (iheight == 1) then
      frac_err = rmatrix % frac_err(iheight)
    else if (iheight > rmatrix % num_heights) then
      frac_err = rmatrix % frac_err(rmatrix % num_heights)
    else
      frac_err = rmatrix % frac_err(iheight - 1) + &
                 (rmatrix % frac_err(iheight) - rmatrix % frac_err(iheight - 1)) * &
                 (obsz(iob) - rmatrix % height(iheight - 1)) / &
                 (rmatrix % height(iheight) - rmatrix % height(iheight - 1))
    end if
 
    WRITE(message,'(A,I8,2F16.4,2E21.8,F12.4)') 'Result', iob, obsz(iob), &
        frac_err, obserr(iob), max(frac_err * obsvalue(iob), rmatrix % min_error), &
        obslat(iob)
    CALL fckit_log % debug(message)
 
    ! Standard deviation
    obserr(iob) = max(frac_err * obsvalue(iob), rmatrix % min_error)
 
  else
    WRITE(message,'(A,I8,2F16.4)') 'Missing', iob, obsz(iob), obslat(iob)
    CALL fckit_log % debug(message)
    obserr(iob) = missing
  end if
end do
 
end subroutine gnssro_obserr_latitude
 
end module gnssro_mod_obserror