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.

Author: steve-marshall-piq

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?

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