From a343fea2d2d088690fcd11f19766b753bb83f085 Mon Sep 17 00:00:00 2001 From: Steve Marshall Date: Fri, 21 Mar 2025 14:17:17 +0000 Subject: [PATCH] Add geopotential height conversion to GNSSRO RefMetOffice forward operator The RefMetOffice GNSSRO forward operator was previous reading geometric height from the MetaData/height variable, but treating it as a geopotential height when comparing to model data. This resulted in a significant bias, especially in the upper atmosphere, when compared to the RefNCEP FO. This code adds an explicit conversion from geometric to geopotential height, using the same function used by RefNCEP. --- .../RefMetOffice/ufo_gnssro_refmetoffice_mod.F90 | 14 +++++++++++--- .../ufo_gnssro_refmetoffice_tlad_mod.F90 | 15 ++++++++++++++- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/src/ufo/operators/gnssro/RefMetOffice/ufo_gnssro_refmetoffice_mod.F90 b/src/ufo/operators/gnssro/RefMetOffice/ufo_gnssro_refmetoffice_mod.F90 index 34db7f2eb..7ed14ebc7 100644 --- a/src/ufo/operators/gnssro/RefMetOffice/ufo_gnssro_refmetoffice_mod.F90 +++ b/src/ufo/operators/gnssro/RefMetOffice/ufo_gnssro_refmetoffice_mod.F90 @@ -28,6 +28,7 @@ module ufo_gnssro_refmetoffice_mod c_virtual, & ! Related to mw_ratio n_alpha, & ! Refractivity constant a n_beta ! Refractivity constant b +use gnssro_mod_transform, only: geometric2geop implicit none @@ -114,6 +115,7 @@ subroutine ufo_gnssro_refmetoffice_simobs(self, geovals, obss, hofx, obs_diags) 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 + real(kind_real) :: tmp_geop_height ! Tmp var used in geopotential height calculation write(err_msg,*) "TRACE: ufo_gnssro_refmetoffice_simobs: begin" call fckit_log%info(err_msg) @@ -171,6 +173,12 @@ subroutine ufo_gnssro_refmetoffice_simobs(self, geovals, obss, hofx, obs_diags) call obsspace_get_db(obss, "MetaData", "latitude", obsLat) call obsspace_get_db(obss, "MetaData", "height", obs_height) +! Convert geometric height to geopotential height + do iobs = 1, nobs + call geometric2geop(obsLat(iobs), obs_height(iobs), tmp_geop_height) + obs_height(iobs) = tmp_geop_height + end do + obs_loop: do iobs = 1, nobs call RefMetOffice_ForwardModel(prs % nval, & @@ -270,15 +278,15 @@ SUBROUTINE RefMetOffice_ForwardModel(nlevp, & 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) :: za(1:nlevp) ! geopotential heights of rho levs +REAL(kind_real), INTENT(IN) :: zb(1:nlevq) ! geopotential 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) :: zobs(1:nobs) ! Geopotential heights for the obs REAL(kind_real), INTENT(IN) :: Latitude ! Latitude of this profile REAL(kind_real), INTENT(INOUT) :: ycalc(1:nobs) ! Model forecast of the observations LOGICAL, INTENT(OUT) :: BAErr ! Was an error encountered during the calculation? diff --git a/src/ufo/operators/gnssro/RefMetOffice/ufo_gnssro_refmetoffice_tlad_mod.F90 b/src/ufo/operators/gnssro/RefMetOffice/ufo_gnssro_refmetoffice_tlad_mod.F90 index 64f55b0cd..42b690297 100644 --- a/src/ufo/operators/gnssro/RefMetOffice/ufo_gnssro_refmetoffice_tlad_mod.F90 +++ b/src/ufo/operators/gnssro/RefMetOffice/ufo_gnssro_refmetoffice_tlad_mod.F90 @@ -32,6 +32,7 @@ module ufo_gnssro_refmetoffice_tlad_mod c_virtual, & ! Related to mw_ratio n_alpha, & ! Refractivity constant a n_beta ! Refractivity constant b +use gnssro_mod_transform, only: geometric2geop private public :: ufo_gnssro_refmetoffice_tlad @@ -126,6 +127,8 @@ subroutine ufo_gnssro_refmetoffice_tlad_settraj(self, geovals, obss) integer :: iobs ! Loop variable, observation number real(kind_real), allocatable :: obs_height(:) ! Geopotential height of the observation + real(kind_real), allocatable :: obsLat(:) ! Observed latitude used in geop hgt calc + real(kind_real) :: tmp_geop_height ! Tmp var used in geopotential height calc write(err_msg,*) "TRACE: ufo_gnssro_refmetoffice_tlad_settraj: begin" call fckit_log%info(err_msg) @@ -152,7 +155,16 @@ subroutine ufo_gnssro_refmetoffice_tlad_settraj(self, geovals, obss) ! Get the meta-data from the observations allocate(obs_height(self%nlocs)) + allocate(obsLat(self%nlocs)) call obsspace_get_db(obss, "MetaData", "height", obs_height) + call obsspace_get_db(obss, "MetaData", "latitude", obsLat) + +! Convert geometric height to geopotential height. + do iobs = 1, self%nlocs + call geometric2geop(obsLat(iobs), obs_height(iobs), tmp_geop_height) + obs_height(iobs) = tmp_geop_height + end do + allocate(self % K(1:self%nlocs, 1:prs%nval + q%nval)) ! For each observation, calculate the K-matrix @@ -167,7 +179,7 @@ subroutine ufo_gnssro_refmetoffice_tlad_settraj(self, geovals, obss) self % vert_interp_ops, & ! Whether to interpolate using log(pressure) self % min_temp_grad, & ! Minimum allowed vertical temperature gradient 1, & ! Number of observations in the profile - obs_height(iobs:iobs), & ! Impact parameter for this observation + obs_height(iobs:iobs), & ! Geopotential height for this observation self % K(iobs:iobs,1:prs%nval+q%nval)) ! K-matrix (Jacobian of the observation with respect to the inputs) ! Flip the K-matrix back the right way around self % K(iobs,1:prs%nval) = self % K(iobs, prs%nval:1:-1) @@ -178,6 +190,7 @@ subroutine ufo_gnssro_refmetoffice_tlad_settraj(self, geovals, obss) self%ltraj = .true. deallocate(obs_height) + deallocate(obsLat) end subroutine ufo_gnssro_refmetoffice_tlad_settraj