ufo_gnssro_bendmetoffice_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_gnssro_bendmetoffice_mod
 
use iso_c_binding
use kinds
use ufo_vars_mod
use ufo_geovals_mod
use ufo_geovals_mod_c, only: ufo_geovals_registry
use vert_interp_mod
use lag_interp_mod,    only: lag_interp_const, lag_interp_smthweights
use obsspace_mod  
use missing_values_mod
use ufo_gnssro_ukmo1d_utils_mod
use ufo_utils_refractivity_calculator, only: ufo_calculate_refractivity
use fckit_log_module,  only : fckit_log
 
implicit none
public             :: ufo_gnssro_bendmetoffice
private
 
  !> Fortran derived type for gnssro trajectory
type :: ufo_gnssro_bendmetoffice
  logical :: vert_interp_ops
  logical :: pseudo_ops
  real(kind_real) :: min_temp_grad
  contains
    procedure :: setup     => ufo_gnssro_bendmetoffice_setup
    procedure :: simobs    => ufo_gnssro_bendmetoffice_simobs
end type ufo_gnssro_bendmetoffice
 
contains
 
! ------------------------------------------------------------------------------
! Get the optional settings for the forward model, and save them in the object
! so that they can be used in the code.
! ------------------------------------------------------------------------------
subroutine ufo_gnssro_bendmetoffice_setup(self, vert_interp_ops, pseudo_ops, min_temp_grad)
 
implicit none
 
class(ufo_gnssro_bendmetoffice), intent(inout) :: self
logical(c_bool), intent(in) :: vert_interp_ops
logical(c_bool), intent(in) :: pseudo_ops
real(c_float), intent(in) :: min_temp_grad
 
self % vert_interp_ops = vert_interp_ops
self % pseudo_ops = pseudo_ops
self % min_temp_grad = min_temp_grad
 
end subroutine ufo_gnssro_bendmetoffice_setup
 
! ------------------------------------------------------------------------------
! 1-dimensional GNSS-RO forward operator for the Met Office system
! ------------------------------------------------------------------------------
subroutine ufo_gnssro_bendmetoffice_simobs(self, geovals, obss, hofx, obs_diags)
 
  implicit none
 
  ! Arguments to this routine
  class(ufo_gnssro_bendmetoffice), intent(in)    :: self      ! The object in which this operator is contained
  type(ufo_geovals),               intent(in)    :: geovals   ! The model values, interpolated to the observation locations
  real(kind_real),                 intent(inout) :: hofx(:)   ! The model forecast of the observations
  type(c_ptr), value,              intent(in)    :: obss      ! The observations, and meta-data for those observations
  type(ufo_geovals),               intent(inout) :: obs_diags ! Observations diagnostics
 
  character(len=*), parameter     :: myname_ = "ufo_gnssro_bendmetoffice_simobs"
  integer, parameter              :: max_string = 800
 
  character(max_string)              :: err_msg         ! Error message for output
  character(max_string)              :: message         ! General message for output
  integer                            :: nobs            ! Number of observations
  integer                            :: ilev            ! Loop variable, level number
  integer                            :: iobs            ! Loop variable, observation number
  type(ufo_geoval), pointer          :: q               ! Model background values of specific humidity
  type(ufo_geoval), pointer          :: prs             ! Model background values of air pressure
  type(ufo_geoval), pointer          :: theta_heights   ! Model heights of levels containing specific humidity
  type(ufo_geoval), pointer          :: rho_heights     ! Model heights of levels containing air pressure
  real(kind_real), allocatable       :: obsLat(:)             ! Latitude of the observation
  real(kind_real), allocatable       :: obsLon(:)             ! Longitude of the observation
  real(kind_real), allocatable       :: impact_param(:)       ! Impact parameter of the observation
  real(kind_real), allocatable       :: radius_curv(:)        ! Earth's radius of curvature at the observation tangent point
  real(kind_real), allocatable       :: undulation(:)         ! Undulation - height of the geoid above the ellipsoid
  logical                            :: flip_data             ! Whether to reverse the order of the model data
  logical                            :: BAErr                 ! Was there an error in the calculation?
  integer                            :: iVar                  ! Loop variable, obs diagnostics variable number
  real(kind_real), allocatable       :: refractivity(:)       ! Refractivity on various model levels
  real(kind_real), allocatable       :: model_heights(:)      ! Geopotential heights that refractivity is calculated on
 
  write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_simobs: begin"
  call fckit_log%info(err_msg)
 
  ! If output to refractivity (and heights of the refractivity levels) is needed,
  ! then use nval as a way to check whether the array has been initialised (since
  ! it is called in a loop).
  DO ivar = 1, obs_diags % nvar
    IF (obs_diags % variables(ivar) == "refractivity" .OR. &
        obs_diags % variables(ivar) == "model_heights") THEN
      write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_simobs: initialising obs_diags for " // &
        obs_diags % variables(ivar)
      call fckit_log%info(err_msg)
      obs_diags % geovals(ivar) % nval = 0
    END IF
  END DO
 
! check if nlocs is consistent in geovals & hofx
  if (geovals%nlocs /= size(hofx)) then
      write(err_msg,*) myname_, ' error: nlocs inconsistent!'
      call abor1_ftn(err_msg)
  endif
 
  write(message, *) myname_, ' Running Met Office GNSS-RO forward operator with'
  call fckit_log%info(message)
  write(message, *) 'vert_interp_ops =', self % vert_interp_ops, &
    'pseudo_ops =', self % pseudo_ops
  call fckit_log%info(message)
 
! get variables from geovals
  call ufo_geovals_get_var(geovals, var_q, q)               ! specific humidity
  call ufo_geovals_get_var(geovals, var_prsi, prs)          ! pressure
  call ufo_geovals_get_var(geovals, var_z, theta_heights)   ! Geopotential height of the normal model levels
  call ufo_geovals_get_var(geovals, var_zi, rho_heights)    ! Geopotential height of the pressure levels
 
  write(message, '(A,10I6)') 'Q: ', q%nval, q%nlocs, shape(q%vals)
  call fckit_log%info(message)
  write(message, '(A,10I6)') 'Pressure: ', prs%nval, prs%nlocs, shape(prs%vals)
  call fckit_log%info(message)
 
  nobs  = obsspace_get_nlocs(obss)
 
  flip_data = .false.
  if (prs%vals(1,1) .lt. prs%vals(prs%nval,1) ) then
    write(err_msg,'(a)') '  ufo_gnssro_bendmetoffice_simobs:'//new_line('a')//                         &
                         '  Model vertical height profile is in descending order,'//new_line('a')// &
                         '  The data will be flipped for processing'
    call fckit_log%info(err_msg)
    flip_data = .true.
  end if
 
! set obs space struture
  allocate(obslon(nobs))
  allocate(obslat(nobs))
  allocate(impact_param(nobs))
  allocate(radius_curv(nobs))
  allocate(undulation(nobs))
 
  call obsspace_get_db(obss, "MetaData", "longitude", obslon)
  call obsspace_get_db(obss, "MetaData", "latitude", obslat)
  call obsspace_get_db(obss, "MetaData", "impact_parameter", impact_param)
  call obsspace_get_db(obss, "MetaData", "earth_radius_of_curvature", radius_curv)
  call obsspace_get_db(obss, "MetaData", "geoid_height_above_reference_ellipsoid", undulation)
 
  obs_loop: do iobs = 1, nobs 
 
    if (flip_data) then
        call ops_gpsro_forwardmodel(prs % nval, &
                                    q % nval, &
                                    rho_heights % vals(rho_heights%nval:1:-1, iobs), &
                                    theta_heights % vals(theta_heights%nval:1:-1, iobs), &
                                    prs % vals(prs%nval:1:-1, iobs), &
                                    q % vals(q%nval:1:-1, iobs), &
                                    self % pseudo_ops, &
                                    self % vert_interp_ops, &
                                    self % min_temp_grad, &
                                    1, &
                                    impact_param(iobs:iobs), &
                                    radius_curv(iobs), &
                                    obslat(iobs), &
                                    undulation(iobs), &
                                    hofx(iobs:iobs), &
                                    baerr, &
                                    refractivity, &
                                    model_heights)
    else
        call ops_gpsro_forwardmodel(prs % nval, &
                                    q % nval, &
                                    rho_heights % vals(:,iobs), &
                                    theta_heights % vals(:,iobs), &
                                    prs % vals(:,iobs), &
                                    q % vals(:,iobs), &
                                    self % pseudo_ops, &
                                    self % vert_interp_ops, &
                                    self % min_temp_grad, &
                                    1, &
                                    impact_param(iobs:iobs), &
                                    radius_curv(iobs), &
                                    obslat(iobs), &
                                    undulation(iobs), &
                                    hofx(iobs:iobs), &
                                    baerr, &
                                    refractivity, &
                                    model_heights)
    end if
 
    if (baerr) then
      write(err_msg,*) "Error with observation processing ", iobs
      call fckit_log % info(err_msg)
    end if
 
    ! If output to refractivity is needed, then initialise things
    DO ivar = 1, obs_diags % nvar
        IF (obs_diags % variables(ivar) == "refractivity") THEN
            IF (iobs == 1) THEN
                obs_diags % geovals(ivar) % nval = SIZE(refractivity)
                ALLOCATE(obs_diags % geovals(ivar) % vals(SIZE(refractivity), obs_diags % nlocs))
            END IF
 
            IF (baerr) THEN
                obs_diags % geovals(ivar) % vals(:,iobs) = missing_value(obs_diags % geovals(ivar) % vals(1,1))
            ELSE
                obs_diags % geovals(ivar) % vals(:,iobs) = refractivity(:)
            END IF
        END IF
 
        IF (obs_diags % variables(ivar) == "model_heights") THEN
            IF (iobs == 1) THEN
                obs_diags % geovals(ivar) % nval = SIZE(model_heights)
                ALLOCATE(obs_diags % geovals(ivar) % vals(SIZE(model_heights), obs_diags % nlocs))
            END IF
 
            IF (baerr) THEN
                obs_diags % geovals(ivar) % vals(:,iobs) = missing_value(obs_diags % geovals(ivar) % vals(1,1))
            ELSE
                obs_diags % geovals(ivar) % vals(:,iobs) = model_heights(:)
            END IF
        END IF
    END DO
  end do obs_loop
 
  deallocate(obslat)
  deallocate(obslon)
  deallocate(impact_param)
  deallocate(radius_curv)
  deallocate(undulation)
 
  write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_simobs: completed"
  call fckit_log%info(err_msg)
 
end subroutine ufo_gnssro_bendmetoffice_simobs
! ------------------------------------------------------------------------------
 
 
SUBROUTINE ops_gpsro_forwardmodel(nlevp, &
                                  nlevq, &
                                  za, &
                                  zb, &
                                  pressure, &
                                  humidity, &
                                  GPSRO_pseudo_ops, &
                                  GPSRO_vert_interp_ops, &
                                  GPSRO_min_temp_grad, &
                                  nobs, &
                                  zobs, &
                                  RO_Rad_Curv, &
                                  Latitude, &
                                  RO_geoid_und, &
                                  ycalc, &
                                  BAErr, &
                                  refractivity, &
                                  model_heights)
 
INTEGER, INTENT(IN)            :: nlevp                  ! no. of p levels in state vec.
INTEGER, INTENT(IN)            :: nlevq                  ! no. of theta levels
REAL(kind_real), INTENT(IN)    :: za(1:nlevp)            ! heights of rho levs
REAL(kind_real), INTENT(IN)    :: zb(1:nlevq)            ! heights of theta levs
REAL(kind_real), INTENT(IN)    :: pressure(1:nlevp)      ! Model background pressure
REAL(kind_real), INTENT(IN)    :: humidity(1:nlevq)      ! Model background specific humidity
LOGICAL, INTENT(IN)            :: GPSRO_pseudo_ops       ! Option: Use pseudo-levels in vertical interpolation?
LOGICAL, INTENT(IN)            :: GPSRO_vert_interp_ops  ! Option: Use ln(p) for vertical interpolation? (rather than exner)
REAL(kind_real), INTENT(IN)    :: GPSRO_min_temp_grad    ! The minimum temperature gradient which is used
INTEGER, INTENT(IN)            :: nobs                   ! Number of observations in the profile
REAL(kind_real), INTENT(IN)    :: zobs(1:nobs)           ! Impact parameter for the obs
REAL(kind_real), INTENT(IN)    :: RO_Rad_Curv            ! Earth's radius of curvature for these observations
REAL(kind_real), INTENT(IN)    :: Latitude               ! Latitude of this profile
REAL(kind_real), INTENT(IN)    :: RO_geoid_und           ! Undulation - difference between the geoid and the ellipsoid
REAL(kind_real), INTENT(INOUT) :: ycalc(1:nobs)          ! Model forecast of the observations
LOGICAL, INTENT(OUT)           :: BAErr                  ! Was an error encountered during the calculation?
REAL(kind_real), INTENT(INOUT), ALLOCATABLE :: refractivity(:)  ! Refractivity as calculated
REAL(kind_real), INTENT(INOUT), ALLOCATABLE :: model_heights(:) ! Height of the levels for refractivity
!
! Things that may need to be output, as they are used by the TL/AD calculation
! 
INTEGER                      :: nRefLevels          ! Number of levels in refractivity calculation
REAL(kind_real), ALLOCATABLE :: nr(:)               ! Model calculation of impact parameters
! 
! Local parameters
! 
integer, parameter           :: max_string = 800  ! Length of strings
character(len=*), parameter  :: myname_ = "Ops_GPSRO_ForwardModel"
!
! Local variables
! 
INTEGER                      :: num_pseudo        ! Number of levels, including pseudo levels
REAL(kind_real)              :: x(1:nlevp+nlevq)  ! state vector
character(max_string)        :: err_msg           ! Error message to be output
 
! The model data must be on a staggered grid, with nlevp = nlevq+1
IF (nlevp /= nlevq + 1) THEN
    write(err_msg,*) myname_ // ':' // ' Data must be on a staggered grid nlevp, nlevq = ', nlevp, nlevq
    call fckit_log % warning(err_msg)
    write(err_msg,*) myname_ // ':' // ' error: number of levels inconsistent!'
    call abor1_ftn(err_msg)
END IF
 
baerr = .false.
 
CALL ufo_calculate_refractivity (nlevp,                 &
                                 nlevq,                 &
                                 za,                    &
                                 zb,                    &
                                 pressure,              &
                                 humidity,              &
                                 gpsro_pseudo_ops,      &
                                 gpsro_vert_interp_ops, &
                                 gpsro_min_temp_grad,   &
                                 baerr,                 &
                                 nreflevels,            &
                                 refractivity,          &
                                 model_heights)
 
ALLOCATE(nr(1:nreflevels))
 
! no point proceeding further if ...
IF (.NOT. baerr) THEN
    !  2.  Calculate the refractive index * radius on theta model levels (or model impact parameter)
    CALL ops_gpsrocalc_nr (nreflevels,    &           ! number of levels for refractivity calculation
                           model_heights, &           ! geopotential heights of refractivity levels
                           refractivity,  &           ! refractivity of model
                           ro_rad_curv,   &           ! radius of curvature of earth at observation
                           latitude,      &           ! latitude at observation
                           ro_geoid_und,  &           ! geoid undulation above WGS-84
                           nr)                        ! Calculated model impact parameters
 
    !  3.  Calculate model bending angle on observation impact parameters
    CALL ops_gpsrocalc_alpha (nobs,         &      ! size of ob. vector
                              nreflevels,   &      ! no. of refractivity levels
                              zobs,         &      ! obs impact parameters
                              refractivity, &      ! refractivity values
                              nr,           &      ! index * radius product
                              ycalc)               ! forward modelled bending angle
END IF
 
END SUBROUTINE ops_gpsro_forwardmodel
 
end module ufo_gnssro_bendmetoffice_mod