ufo_gnssro_refmetoffice_tlad_mod.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! (C) British Crown Copyright 2021 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 refractivity Met Office's tangent linear and adjoint
 
module ufo_gnssro_refmetoffice_tlad_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 ufo_basis_tlad_mod,  only: ufo_basis_tlad
use obsspace_mod
use gnssro_mod_conf
use missing_values_mod
use fckit_log_module, only : fckit_log
use ufo_utils_refractivity_calculator, only: &
    ufo_calculate_refractivity, ufo_refractivity_partial_derivatives
use ufo_constants_mod, only: &
    rd,                      &    ! Gas constant for dry air
    cp,                      &    ! Heat capacity at constant pressure for air
    rd_over_cp,              &    ! Ratio of gas constant to heat capacity
    grav,                    &    ! Gravitational field strength
    pref,                    &    ! Reference pressure for calculating exner
    mw_ratio,                &    ! Ratio of molecular weights of water and dry air
    c_virtual,               &    ! Related to mw_ratio
    n_alpha,                 &    ! Refractivity constant a
    n_beta                        ! Refractivity constant b
 
 
integer, parameter :: max_string=800
 
!> Fortran derived type for gnssro trajectory
type, extends(ufo_basis_tlad)   ::  ufo_gnssro_refmetoffice_tlad
  private
  logical :: vert_interp_ops              !< Do vertical interpolation using ln(p) or exner?
  logical :: pseudo_ops                   !< Use pseudo-levels in the refractivity calculation
  logical :: flip_it                      !< Do the outputs need to be flipped in order
  real(kind_real) :: min_temp_grad        !< The minimum temperature gradient before a profile is considered isothermal
                                          !  Used in pseudo-levels calculation
  integer                       :: nlevp  !< The number of pressure levels
  integer                       :: nlevq  !< The number of specific humidity (or temperature) levels
  integer                       :: nlocs  !< The number of observations
  real(kind_real), allocatable  :: k(:,:) !< The K-matrix (Jacobian)
  contains
    procedure :: setup      => ufo_gnssro_refmetoffice_tlad_setup
    procedure :: delete     => ufo_gnssro_refmetoffice_tlad_delete
    procedure :: settraj    => ufo_gnssro_refmetoffice_tlad_settraj
    procedure :: simobs_tl  => ufo_gnssro_refmetoffice_simobs_tl
    procedure :: simobs_ad  => ufo_gnssro_refmetoffice_simobs_ad
end type ufo_gnssro_refmetoffice_tlad
 
contains
 
!-------------------------------------------------------------------------------
!> \brief Set up the Met Office GNSS-RO refractivity TL/AD
!!
!! \details **ufo_gnssro_refmetoffice_tlad_setup**
!! * Get the optional settings for the forward model and its linear, and save
!!   them in the object so that they can be used in the code.
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 26 May 2021
!!
!-------------------------------------------------------------------------------
subroutine ufo_gnssro_refmetoffice_tlad_setup(self, vert_interp_ops, pseudo_ops, min_temp_grad)
 
implicit none
 
class(ufo_gnssro_refmetoffice_tlad), 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_refmetoffice_tlad_setup
 
!-------------------------------------------------------------------------------
!> \brief Set up the K-matrix for (Jacobian) for the Met Office's GNSS-RO
!!        refractivity operator
!!
!! \details **ufo_gnssro_refmetoffice_tlad_settraj**
!! * It is necessary to run this routine before calling the TL or AD routines.
!! * Get the geovals specifying the state around which to linearise, flipping
!!   the vertical order if required.
!! * Call the helper function to calculate the K-matrix for each observation.
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 26 May 2021
!!
!-------------------------------------------------------------------------------
subroutine ufo_gnssro_refmetoffice_tlad_settraj(self, geovals, obss)
       
  implicit none
! Subroutine arguments
  class(ufo_gnssro_refmetoffice_tlad), intent(inout) :: self     !< The object that we use to save data in
  type(ufo_geovals),                intent(in)       :: geovals  !< The input geovals
  type(c_ptr), value,               intent(in)       :: obss     !< The input observations
 
! Local parameters
  character(len=*), parameter :: myname_="ufo_gnssro_refmetoffice_tlad_settraj"
 
! Local variables
  character(max_string)       :: err_msg                       ! Messages to be output to the user
  type(ufo_geoval), pointer   :: q                             ! The model geovals - specific humidity
  type(ufo_geoval), pointer   :: prs                           ! The model geovals - atmospheric pressure
  type(ufo_geoval), pointer   :: rho_heights                   ! The model geovals - heights of the pressure-levels
  type(ufo_geoval), pointer   :: theta_heights                 ! The model geovals - heights of the theta-levels (stores q)
  type(ufo_geoval), pointer   :: p_temp                        ! The model geovals - atmospheric pressure before flipping
  integer                     :: iobs                          ! Loop variable, observation number
 
  real(kind_real), allocatable       :: obs_height(:)          ! Geopotential height of the observation
  type(ufo_geovals)                  :: geovals_local          ! The model values, interpolated to the observation locations
                                                               ! and flipped in the vertical (if required)
 
  write(err_msg,*) "TRACE: ufo_gnssro_refmetoffice_tlad_settraj: begin"
  call fckit_log%info(err_msg)
 
! Make sure that any previous values of geovals don't get carried over
  call self%delete()
 
! make sure that the geovals are in the correct vertical order (surface first)
  call ufo_geovals_copy(geovals, geovals_local)  ! dont want to change geovals input
  call ufo_geovals_get_var(geovals_local, var_prsi, p_temp)
  if( p_temp%vals(1,1) < p_temp%vals(p_temp%nval,1) ) then 
    self%flip_it = .true.
  else
    self%flip_it = .false.
  endif
  call ufo_geovals_reorderzdir(geovals_local, var_prsi, "bottom2top")
 
! get model state variables from geovals
  call ufo_geovals_get_var(geovals_local, var_q,    q)             ! specific humidity
  call ufo_geovals_get_var(geovals_local, var_prsi, prs)           ! pressure
  call ufo_geovals_get_var(geovals_local, var_z,    theta_heights) ! Geopotential height of the normal model levels
  call ufo_geovals_get_var(geovals_local, var_zi,   rho_heights)   ! Geopotential height of the pressure levels
 
! Keep copy of dimensions
  self % nlevp = prs % nval
  self % nlevq = q % nval
  self % nlocs = obsspace_get_nlocs(obss)
  
! Get the meta-data from the observations
  allocate(obs_height(self%nlocs))
  call obsspace_get_db(obss, "MetaData", "obs_height", obs_height)
  allocate(self % K(1:self%nlocs, 1:prs%nval + q%nval))
 
! For each observation, calculate the K-matrix
  obs_loop: do iobs = 1, self % nlocs
    CALL jacobian_interface(prs % nval, &                          ! Number of pressure levels
                            q % nval, &                            ! Number of specific humidity levels
                            rho_heights % vals(:,iobs), &          ! Heights of the pressure levels
                            theta_heights % vals(:,iobs), &        ! Heights of the specific humidity levels
                            prs % vals(:,iobs), &                  ! Values of the pressure
                            q % vals(:,iobs), &                    ! Values of the specific humidity
                            self % pseudo_ops, &                   ! Whether to use pseudo-levels in the calculation
                            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
                            self % K(iobs:iobs,1:prs%nval+q%nval)) ! K-matrix (Jacobian of the observation with respect to the inputs)
  end do obs_loop
 
! Note that this routine has been run.
  self%ltraj = .true.
 
  deallocate(obs_height)
 
end subroutine ufo_gnssro_refmetoffice_tlad_settraj
 
!-------------------------------------------------------------------------------
!> \brief Given an increment to the model state, calculate an increment to the
!!        observation
!!
!! \details **ufo_gnssro_refmetoffice_simobs_tl**
!! * Check that set trajectory has been previously called.
!! * Get the geovals for the increment.
!! * For each observation apply the K-matrix to calculate the increment to the
!!   observation.
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 26 May 2021
!!
!-------------------------------------------------------------------------------
 
subroutine ufo_gnssro_refmetoffice_simobs_tl(self, geovals, hofx, obss)
 
  implicit none
 
! Subroutine arguments
  class(ufo_gnssro_refmetoffice_tlad), intent(in) :: self       !< Object which is being used to transfer information
  type(ufo_geovals),                intent(in)    :: geovals    !< Model perturbations
  real(kind_real),                  intent(inout) :: hofx(:)    !< Increment to the observations
  type(c_ptr),   value,             intent(in)    :: obss       !< Input - the observations
 
! Local parameters
  character(len=*), parameter  :: myname_="ufo_gnssro_refmetoffice_simobs_tl"
 
! Local variables
  integer                      :: iobs      ! Loop variable, observation number
  integer                      :: nlocs     ! Number of observations
  character(max_string)        :: err_msg   ! Message to be output
  type(ufo_geoval), pointer    :: q_d       ! Increment to the specific humidity
  type(ufo_geoval), pointer    :: prs_d     ! Increment to the air pressure
  real(kind_real), allocatable :: x_d(:)    ! Increment to the complete state
 
  write(err_msg,*) "TRACE: ufo_gnssro_refmetoffice_simobs_tl: begin"
  call fckit_log%info(err_msg)
 
! Check if trajectory was set
  if (.not. self%ltraj) then
     write(err_msg,*) myname_, ' trajectory wasnt set!'
     call abor1_ftn(err_msg)
  endif
 
! 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
 
! Get variables from geovals
  call ufo_geovals_get_var(geovals, var_q,     q_d)         ! specific humidity
  call ufo_geovals_get_var(geovals, var_prsi,  prs_d)       ! pressure on rho levels
 
  nlocs = self % nlocs ! number of observations
 
  allocate(x_d(1:prs_d%nval+q_d%nval))
! Loop through the obs, calculating the increment to the observation
  obs_loop: do iobs = 1, nlocs   ! order of loop doesn't matter
 
    x_d(1:prs_d%nval) = prs_d % vals(:,iobs)
    x_d(prs_d%nval+1:prs_d%nval+q_d%nval) = q_d % vals(:,iobs)
    hofx(iobs) = sum(self % K(iobs,:) * x_d)
 
  end do obs_loop
 
  deallocate(x_d)
 
  write(err_msg,*) "TRACE: ufo_gnssro_refmetoffice_simobs_tl: complete"
  call fckit_log%info(err_msg)
 
  return
    
end subroutine ufo_gnssro_refmetoffice_simobs_tl
 
!-------------------------------------------------------------------------------
!> \brief Given an increment to the observation, find the equivalent increment
!!        to the model state
!!
!! \details **ufo_gnssro_refmetoffice_simobs_ad**
!! * Check that set trajectory has previously been called.
!! * Get the geovals for the model increment, and allocate these if they have
!!   not already been allocated.
!! * For each observation calculate the observation increment using the K-matrix
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 26 May 2021
!!
!-------------------------------------------------------------------------------
 
subroutine ufo_gnssro_refmetoffice_simobs_ad(self, geovals, hofx, obss)
 
  use typesizes,     only: wp => eightbytereal
 
  implicit none
 
! Subroutine arguments
  class(ufo_gnssro_refmetoffice_tlad), intent(in) :: self       !< Object which is being used to transfer information
  type(ufo_geovals),                intent(inout) :: geovals    !< Calculated perturbations to model state
  real(kind_real),                  intent(in)    :: hofx(:)    !< Increment to the observations
  type(c_ptr),  value,              intent(in)    :: obss       !< Input - the observations
 
! Local parameters
  character(len=*), parameter     :: myname_="ufo_gnssro_refmetoffice_simobs_ad"
 
! Local variables
  real(c_double)               :: missing  ! Missing data values
  type(ufo_geoval), pointer    :: q_d      ! Pointer to the specific humidity perturbations
  type(ufo_geoval), pointer    :: prs_d    ! Pointer to the pressure perturbations
  integer                      :: iobs     ! Loop variable, observation number
  real(kind_real), allocatable :: x_d(:)   ! Perturbation to the full model state
  character(max_string)        :: err_msg  ! Message to be output
 
  write(err_msg,*) "TRACE: ufo_gnssro_refmetoffice_simobs_ad: begin"
  call fckit_log%info(err_msg)
 
! Check if trajectory was set
  if (.not. self%ltraj) then
     write(err_msg,*) myname_, ' trajectory wasnt set!'
     call abor1_ftn(err_msg)
  endif
 
! 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
     
! Get variables from geovals
  call ufo_geovals_get_var(geovals, var_q,     q_d)         ! specific humidity
  call ufo_geovals_get_var(geovals, var_prsi,  prs_d)       ! pressure
 
  missing = missing_value(missing)
  allocate(x_d(1:prs_d%nval + q_d%nval))
 
! Loop through the obs, calculating the increment to the model state
  obs_loop: do iobs = 1, self % nlocs
 
    if (hofx(iobs) /= missing) then
        x_d = self % K(iobs,:) * hofx(iobs)
        if (self%flip_it) then
            prs_d % vals(:,iobs) = prs_d % vals(:,iobs) + x_d(prs_d%nval:1:-1)
            q_d % vals(:,iobs) = q_d % vals(:,iobs) + x_d(prs_d%nval+q_d%nval:prs_d%nval+1:-1)
        else
            prs_d % vals(:,iobs) = prs_d % vals(:,iobs) + x_d(1:prs_d%nval)
            q_d % vals(:,iobs) = q_d % vals(:,iobs) + x_d(prs_d%nval+1:prs_d%nval+q_d%nval)
        end if
    end if
 
  end do obs_loop
 
  deallocate(x_d)
 
  write(err_msg,*) "TRACE: ufo_gnssro_refmetoffice_simobs_ad: complete"
  call fckit_log%info(err_msg)
 
  return
 
end subroutine ufo_gnssro_refmetoffice_simobs_ad
 
!-------------------------------------------------------------------------------
!> \brief Tidy up the variables that are used for passing information
!!
!! \details **ufo_gnssro_refmetoffice_tlad_delete**
!! * Set lengths to zero, and deallocate K-matrix
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 26 May 2021
!!
!-------------------------------------------------------------------------------
 
subroutine ufo_gnssro_refmetoffice_tlad_delete(self)
 
  implicit none
  class(ufo_gnssro_refmetoffice_tlad), intent(inout) :: self
  character(len=*), parameter :: myname_="ufo_gnssro_refmetoffice_tlad_delete"
      
  self%nlocs = 0
  self%nlevp = 0
  self%nlevq = 0
  if (allocated(self%K)) deallocate(self%K)
  self%ltraj = .false. 
 
end subroutine ufo_gnssro_refmetoffice_tlad_delete
 
!-------------------------------------------------------------------------------
!> \brief Calculate the K-matrix used in the TL/AD
!!
!! \details **jacobian_interface**
!! * Allocate temporary arrays
!! * Call partial-derivatives code which calculates various quantities on model
!!   levels
!! * Loop over each observation, doing the following
!! * If using pseudo-levels, calculate specific humidity, pressure and
!!   temperature on intermediate levels, and interpolate these appropriately.
!!   Then calculate the refractivity gradients for each observation,
!!   interpolated from the pseudo-levels.
!! * If not using pseudo-levels, assume that refractivity varies exponentially
!!   with height and calculate gradients.
!! * Calculate K-matrix from the component gradients.
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 26 May 2021
!!
!-------------------------------------------------------------------------------
 
SUBROUTINE jacobian_interface(nlevp,  &
                              nlevq,  &
                              za,     &
                              zb,     &
                              p,      &
                              q,      &
                              pseudo_ops, &
                              vert_interp_ops, &
                              min_temp_grad, &
                              nobs,   &
                              zobs,   &
                              Kmat)
 
IMPLICIT NONE
 
! Subroutine arguments:
integer, intent(in)          :: nlevp            !< Number of pressure levels
integer, intent(in)          :: nlevq            !< Number of specific humidity levels
real(kind_real), intent(in)  :: za(:)            !< Height of the pressure levels
real(kind_real), intent(in)  :: zb(:)            !< Height of the specific humidity levels
real(kind_real), intent(in)  :: p(:)             !< Input pressure
real(kind_real), intent(in)  :: q(:)             !< Input specific humidity
logical, intent(in)          :: vert_interp_ops  !< Use log(p) for vertical interpolation?
logical, intent(in)          :: pseudo_ops       !< Use pseudo-levels to reduce errors?
real(kind_real), intent(in)  :: min_temp_grad    !< Minimum value for the vertical temperature gradient
integer, intent(in)          :: nobs             !< Number of observations
real(kind_real), intent(in)  :: zobs(:)          !< Height of the observations
real(kind_real), intent(out) :: Kmat(:,:)        !< K-matrix (Jacobian)
 
! Local declarations:
integer                      :: n                          ! Loop variable, observation number
integer                      :: i                          ! Loop variable
real(kind_real)              :: Exner(nlevp)               ! Exner pressure, calculated on model pressure levels
real(kind_real)              :: Pb(nlevq)                  ! Pressure on model theta levels
real(kind_real)              :: Tv(nlevq)                  ! Virtual temperature on model theta levels
real(kind_real)              :: T(nlevq)                   ! Temperature on model theta levels
real(kind_real)              :: dExtheta_dPb(nlevq,nlevq)  ! Partial derivative of theta wrt pressure (on model theta levels)
real(kind_real), ALLOCATABLE :: drefob_dPob(:,:)           ! Partial derivative of refractivity wrt pressure (at ob location)
real(kind_real), ALLOCATABLE :: drefob_dTob(:,:)           ! Partial derivative of refractivity wrt temperature (at ob location)
real(kind_real), ALLOCATABLE :: drefob_dqob(:,:)           ! Partial derivative of refractivity wrt specific humidity (at ob location)
real(kind_real), ALLOCATABLE :: dqob_dqb(:,:)              ! Partial derivative of specific humidity at ob location wrt specific humidity on model theta levels
real(kind_real), ALLOCATABLE :: dTob_dTb(:,:)              ! Partial derivative of temperature at ob location wrt temperature on model theta levels
real(kind_real), ALLOCATABLE :: dPob_dT(:,:)               ! Partial derivative of pressure at ob location wrt temperature on model theta levels
real(kind_real), ALLOCATABLE :: dPob_dPb(:,:)              ! Partial derivative of pressure at ob location wrt pressure on model theta levels
real(kind_real), ALLOCATABLE :: dPob_dP(:,:)               ! Partial derivative of pressure at ob location wrt pressure on model pressure levels
real(kind_real), ALLOCATABLE :: dTob_dP(:,:)               ! Partial derivative of temperature at ob location wrt pressure on model pressure levels
real(kind_real), ALLOCATABLE :: dTb_dp(:,:)                ! Partial derivative of temperature on model theta levels wrt pressure on model pressure levels
real(kind_real)              :: c                          ! Continuity constant for hydrostatic pressure
real(kind_real)              :: P_ob                       ! Model pressure at the ob location
real(kind_real)              :: q_ob                       ! Model specific humidity at the ob location
real(kind_real)              :: T_ob                       ! Model temperature at the ob location
real(kind_real)              :: gamma                      ! Vertical gradient of log(q)
real(kind_real)              :: beta                       ! Vertical temperature gradient
real(kind_real)              :: Ndry                       ! Component of refractivity due to dry terms
real(kind_real)              :: Nwet                       ! Component of refractivity due to wet terms
real(kind_real)              :: refrac(nlevq)              ! Refractivity on model theta levels
real(kind_real)              :: dEx_dP(nlevp,nlevp)        ! Partial derivative of exner wrt pressure
real(kind_real)              :: dPb_dp(nlevq,nlevp)        ! Partial derivative of pressure on theta levels wrt pressure on pressure levels
real(kind_real)              :: dTv_dExtheta(nlevq,nlevq)  ! Virtual temperature divided by exner on theta levels
real(kind_real)              :: dTv_dEx(nlevq,nlevp)       ! Partial derivative of virtual temperature wrt exner
real(kind_real)              :: dT_dTv(nlevq,nlevq)        ! Partial derivative of temperature wrt virtual temperature
real(kind_real)              :: dT_dq(nlevq,nlevq)         ! Partial derivative of temperature wrt specific humidity
real(kind_real)              :: dref_dPb(nlevq,nlevq)      ! Partial derivative of refractivity wrt pressure on theta levels
real(kind_real)              :: dref_dT(nlevq,nlevq)       ! Partial derivative of refractivity wrt temperature
real(kind_real)              :: dref_dq(nlevq,nlevq)       ! Partial derivative of refractivity wrt specific humidity
real(kind_real)              :: dNref_dref(nobs,nlevq)     ! Partial derivative of refractivity at the ob loation wrt model refractivity
real(kind_real)              :: kmatP(nobs,nlevp)          ! K-matrix contribution for pressure terms
real(kind_real)              :: kmatq(nobs,nlevq)          ! K-matrix contribution for specific humidity terms
real(kind_real)              :: Nref(nobs)                 ! Model refractivity at observation locations
real(kind_real)              :: m1(nobs,nlevq)             ! Temporary matrix product used in calculation
real(kind_real)              :: m2(nobs,nlevp)             ! Temporary matrix product used in calculation
real(kind_real)              :: m3(nobs,nlevq)             ! Temporary matrix product used in calculation
real(kind_real)              :: m4(nobs,nlevq)             ! Temporary matrix product used in calculation
real(kind_real)              :: g_RB                       ! Frequently used term
real(kind_real)              :: c_ZZ                       ! Frequently used term
real(kind_real)              :: dPo_dT1                    ! dP_ob / dT_below
real(kind_real)              :: dPo_dTo                    ! dP_ob / dT_ob
real(kind_real)              :: dTo_dT1                    ! dT_ob / dT_below
real(kind_real)              :: dPo_dbeta                  ! dP_ob / dbeta
real(kind_real)              :: dbeta_dT1                  ! dbeta / dT_below
real(kind_real)              :: dTo_dbeta                  ! dT_ob / dbeta
real(kind_real)              :: dbeta_dT2                  ! dbeta / dT_above
real(kind_real)              :: dPo_dc                     ! dP_ob / dc
real(kind_real)              :: dc_dT1                     ! dc / dT_below
real(kind_real)              :: dc_dbeta                   ! dc / dbeta
real(kind_real)              :: dc_dt2                     ! dc / dT_above
real(kind_real)              :: dPo_dP1                    ! dP_ob / dP_below
real(kind_real)              :: dc_dP1                     ! dc / dP_below
real(kind_real)              :: dc_dP2                     ! dc / dP_above
real(kind_real)              :: model_height_diff          ! Difference in height between two model levels
real(kind_real)              :: obs_height_diff            ! Difference in height between the observation and the model level below it
real(kind_real)              :: height_diff_ratio          ! obs_height_diff / model_height_diff
 
!-----------------------
! 1. Allocate/initialise matrices
!-----------------------
 
IF (pseudo_ops) THEN
  ALLOCATE (drefob_dpob(nobs,nobs))
  ALLOCATE (drefob_dtob(nobs,nobs))
  ALLOCATE (drefob_dqob(nobs,nobs))
  ALLOCATE (dqob_dqb(nobs,nlevq))
  ALLOCATE (dtob_dtb(nobs,nlevq))
  ALLOCATE (dpob_dt(nobs,nlevq))
  ALLOCATE (dpob_dpb(nobs,nlevq))
  ALLOCATE (dpob_dp(nobs,nlevp))
  ALLOCATE (dtob_dp(nobs,nlevp))
  ALLOCATE (dtb_dp(nlevq,nlevp))
 
  drefob_dpob(:,:) = 0.0
  drefob_dtob(:,:) = 0.0
  drefob_dqob(:,:) = 0.0
  dqob_dqb(:,:) = 0.0
  dtob_dtb(:,:) = 0.0
  dpob_dt(:,:) = 0.0
  dpob_dpb(:,:) = 0.0
  dpob_dp(:,:) = 0.0
  dtob_dp(:,:) = 0.0
  dtb_dp(:,:) = 0.0
END IF
 
call ufo_refractivity_partial_derivatives(nlevp,           &
                                          nlevq,           &
                                          za,              &
                                          zb,              &
                                          p,               &
                                          q,               &
                                          vert_interp_ops, &
                                          dt_dtv,          &
                                          dt_dq,           &
                                          dref_dpb,        &
                                          dref_dt,         &
                                          dref_dq,         &
                                          refrac,          &
                                          t,               &
                                          pb,              &
                                          dex_dp,          &
                                          dextheta_dpb,    &
                                          dtv_dextheta,    &
                                          dpb_dp,          &
                                          dtv_dex)
 
 
dnref_dref(:,:) = 0.0
kmatp(:,:) = 0.0
kmatq(:,:) = 0.0
kmat(:,:) = 0.0
 
!-------------------------------------------------
! 2. Calculate the partial derivatives at the observation locations
!-------------------------------------------------
 
DO n = 1, nobs
  i = 1
 
  if (zobs(n) >= zb(nlevq) .or. zobs(n) == missing_value(zobs(n))) cycle
  DO
    ! Find model layer containing ob
    IF (zobs(n) < zb(i + 1)) EXIT
 
    i = i + 1
 
  END DO
 
  IF (pseudo_ops) THEN
    ! Calculate some quantities that will be re-used
    model_height_diff = zb(i + 1) - zb(i)
    obs_height_diff = zobs(n) - zb(i)
    height_diff_ratio = obs_height_diff / model_height_diff
  
    ! Interpolate P,T,q separately
    IF (min(q(i), q(i + 1)) > 0.0) THEN
      ! q varies exponentially with height
      gamma = log(q(i) / q(i + 1)) / model_height_diff
      q_ob = q(i) * exp(-gamma * obs_height_diff)
 
      dqob_dqb(n,i) = (q_ob / q(i)) * (1.0 - height_diff_ratio)
      dqob_dqb(n,i + 1) = (q_ob / q(i + 1)) * height_diff_ratio
 
    ! Assume linear variation if humidities are -ve to avoid singularity
    ELSE
      q_ob = q(i) + (q(i + 1) - q(i)) * height_diff_ratio
 
      dqob_dqb(n,i) = (zb(i + 1) - zobs(n)) / model_height_diff
      dqob_dqb(n,i + 1) = height_diff_ratio
    END IF
 
    ! T varies linearly with height
    beta = (t(i + 1) - t(i)) / model_height_diff
    t_ob = t(i) + beta * obs_height_diff
 
    dtob_dtb(n,i) = (zb(i + 1) - zobs(n)) / model_height_diff
    dtob_dtb(n,i + 1) = height_diff_ratio
 
    ! P varies to maintain hydrostatic balance
    IF (abs(t(i + 1) - t(i)) > min_temp_grad) THEN
      c = ((pb(i + 1) / pb(i)) * (t(i + 1)/t(i)) ** (grav / (rd * beta)) - 1.0) / model_height_diff
      p_ob = (pb(i) * (t_ob / t(i)) ** &
                  (-grav / (rd * beta))) * (1.0 + c * obs_height_diff)
    ELSE
      ! If layer is isothermal, assume exponential variation to
      ! avoid singularity
      p_ob = pb(i) * exp(log(pb(i + 1) / pb(i)) * height_diff_ratio)
    END IF
 
    ! Temporary variables to keep equations neater
    g_rb = grav / (rd * beta)
    c_zz = c + 1.0 / model_height_diff
 
    dpo_dt1 = (g_rb / t(i)) * p_ob
    dpo_dto = -(g_rb / t_ob) * p_ob
    dto_dt1 = dtob_dtb(n,i)
    dpo_dbeta = (g_rb / beta) * log(t_ob / t(i)) * p_ob
    dbeta_dt1 = -1.0 / model_height_diff
    dto_dbeta = obs_height_diff
    dbeta_dt2 = 1.0 / model_height_diff
 
    ! Pressure Jacobians for hydrostatic/exponential cases
    IF (abs(t(i + 1) - t(i)) > min_temp_grad) THEN
      dpo_dc = obs_height_diff * pb(i) * (t_ob / t(i)) ** (-g_rb)
      dc_dt1 = -(g_rb / (t(i))) * (c_zz)
      dc_dbeta = -(g_rb / beta) * log(t(i + 1) / t(i)) * c_zz
      dc_dt2 = (g_rb / t(i + 1)) * c_zz
      dpo_dp1 = p_ob / pb(i)
      dpob_dt(n,i) = dpo_dt1 + dpo_dto * dto_dt1 + dpo_dbeta * dbeta_dt1 + dpo_dc * &
                     (dc_dt1 + dc_dbeta * dbeta_dt1)
      dpob_dt(n,i+1) = (dpo_dbeta + dpo_dto * dto_dbeta + dpo_dc * dc_dbeta) * &
                     dbeta_dt2 + dpo_dc * dc_dt2
 
      dc_dp1 = -(1.0 / pb(i)) * c_zz
      dc_dp2 = (1.0 / pb(i + 1)) * c_zz
 
      dpob_dpb(n,i) = dpo_dp1 + dpo_dc * dc_dp1
      dpob_dpb(n,i + 1) = dpo_dc * dc_dp2
    ELSE
      dpob_dpb(n,i) = exp(log(pb(i + 1) / pb(i)) * height_diff_ratio) * &
        (1.0 - height_diff_ratio)
      dpob_dpb(n,i + 1) = (pb(i) / pb(i + 1)) * height_diff_ratio * &
        exp(log(pb(i + 1) / pb(i)) * height_diff_ratio)
    END IF
 
    ! Calculate refractivity
    ndry = n_alpha * p_ob / t_ob
    nwet = n_beta * p_ob * q_ob / (t_ob ** 2 * (mw_ratio + (1.0 - mw_ratio) * q_ob))
    nref(n) = ndry + nwet
 
    drefob_dpob(n,n) = nref(n) / p_ob
    drefob_dtob(n,n) = -(ndry + 2.0 * nwet) / t_ob
    drefob_dqob(n,n) = n_beta * p_ob * mw_ratio / (t_ob * (mw_ratio + (1.0 - mw_ratio) * q_ob)) ** 2
 
  ELSE
 
    ! Use simple assumption of exponentially varying refractivity
 
    gamma = (zb(i + 1) - zobs(n)) / (zb(i + 1) - zb(i))
 
    beta = 1.0 - gamma
 
    nref(n) = exp(gamma * log(refrac(i)) + beta * log(refrac(i + 1)))
 
    dnref_dref(n,i) = nref(n) * gamma / refrac(i)
 
    dnref_dref(n,i + 1) = nref(n) * beta / refrac(i + 1)
  END IF
 
END DO
 
!-------------------------------------------------
! 3. Evaluate the Kmatrix by matrix multiplication
!-------------------------------------------------
 
IF (pseudo_ops) THEN
 
  ! Derivatives:
  ! dPob/dP = (dPob/dPb * dPb/dP) + (dPob/dT * dT/dTv * dTv/dEx * dEx/dP)
  dpob_dp = matmul(dpob_dpb, dpb_dp) + matmul(dpob_dt, matmul(dt_dtv, matmul(dtv_dex, dex_dp)))
  ! dTob/dP = (dTob/dT * dT/dTv * dTv/dEx * dEx/dP)
  dtob_dp = matmul(dtob_dtb, matmul(dt_dtv, matmul(dtv_dex, dex_dp)))
 
  ! calc K matrix for p on rho levels
  ! dNob/dP = (dNob/dPob * dPob/dP) + (dNob/dTob *dTob/dP)
  kmatp(:,:) = matmul(drefob_dpob, dpob_dp) + matmul(drefob_dtob, dtob_dp)
 
  ! calc Kmatrix for q on theta levels
  ! dNob/dq = (dNob/dqob * dqob/dq) + (dNob/dTob * dTob/dT * dT/dq) + (dNob/dPob * dPob/dT * dT/dq)
  kmatq(:,:) = matmul(drefob_dqob, dqob_dqb) + matmul(matmul(drefob_dtob, dtob_dtb), dt_dq) + &
               matmul(matmul(drefob_dpob, dpob_dt), dt_dq)
 
ELSE
 
  ! calc K matrix for p on rho levels
  ! dNob/dP = (dNob/dN * dN/dPb * dPb/dP)   + ....
  kmatp(:,:) = matmul(dnref_dref, matmul(dref_dpb, dpb_dp))
 
  !  .... (dNob/dN * dN/dT * dT/dTv * dTv/dEx * dEx/dP) + .....
  m1(:,:) = matmul(dnref_dref, matmul(dref_dt, dt_dtv))
  m2(:,:) = matmul(m1, dtv_dex)
  kmatp(:,:) = kmatp(:,:) + matmul(m2, dex_dp)
 
  !  .... (dNob/dN * dN/dT * dT/dTv * dTv/dExtheta * dExtheta/dPb * dPb/dP)
  m3(:,:) = matmul(m1, dtv_dextheta)
  m4(:,:) = matmul(m3, dextheta_dpb)
  kmatp(:,:) = kmatp(:,:) + matmul(m4,dpb_dp)
 
  ! calc Kmatrix for q on theta levels
  !  dNob/dq = (dNob/dN * dN/dq) + (dNob/dN * dN/dT * dT/dq)
  kmatq(:,:) = kmatq(:,:) + matmul(dnref_dref, matmul(dref_dt, dt_dq))
 
END IF
 
! Calculate the full kmatrix in correct units
 
kmat(1:nobs, 1:nlevp) = 1.0e2 * kmatp(1:nobs, 1:nlevp)   ! h/Pa
 
kmat(1:nobs, nlevp+1:nlevp+nlevq) = 1.0e-3 * kmatq(1:nobs, 1:nlevq)   ! g/kg
 
!-------------------------------------------------------------------------------
! 4. Deallocate the temporary arrays
!-------------------------------------------------------------------------------
 
IF (ALLOCATED (drefob_dpob)) DEALLOCATE (drefob_dpob)
IF (ALLOCATED (drefob_dtob)) DEALLOCATE (drefob_dtob)
IF (ALLOCATED (drefob_dqob)) DEALLOCATE (drefob_dqob)
IF (ALLOCATED (dqob_dqb)) DEALLOCATE (dqob_dqb)
IF (ALLOCATED (dtob_dtb)) DEALLOCATE (dtob_dtb)
IF (ALLOCATED (dpob_dt)) DEALLOCATE (dpob_dt)
IF (ALLOCATED (dpob_dpb)) DEALLOCATE (dpob_dpb)
IF (ALLOCATED (dpob_dp)) DEALLOCATE (dpob_dp)
IF (ALLOCATED (dtob_dp)) DEALLOCATE (dtob_dp)
IF (ALLOCATED (dtb_dp)) DEALLOCATE (dtb_dp)
 
END SUBROUTINE jacobian_interface
 
end module ufo_gnssro_refmetoffice_tlad_mod