RefractivityCalculator.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.
!-------------------------------------------------------------------------------
 
!-------------------------------------------------------------------------------
!> \brief Module for containing a general refractivity forward operator and its K-matrix
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 6 May 2021
!!
!-------------------------------------------------------------------------------
 
module ufo_utils_refractivity_calculator
 
use fckit_log_module, only: fckit_log
use missing_values_mod
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
    pref,                    &    ! Reference pressure for calculating exner
    pi,                      &    ! Something to do with circles...
    grav,                    &    ! Gravitational field strength
    ecc,                     &    ! eccentricity
    k_somig,                 &    ! Somigliana's constant
    g_equat,                 &    ! equatorial gravity (ms-2)
    a_earth,                 &    ! semi-major axis of earth (m)
    flatt,                   &    ! flattening
    m_ratio,                 &    ! gravity ratio
    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
use kinds,            only: kind_real
 
implicit none
 
public :: ufo_calculate_refractivity
public :: ufo_refractivity_kmat
 
contains
 
!-------------------------------------------------------------------------------
!> \brief Calculation of the refractivity from the pressure and specific humidity
!!
!! \details **ufo_calculate_refractivity**
!! * Checks the inputs that the pressure is not missing, is monotonically
!!   decreasing and not negative.
!! * Perform the calculation on the model temperature levels.  It first calculates
!!   the exner function and virtual temperature from the inputs.  Then using
!!   these calculates the wet and dry components of the refractivity.
!! * If pseudo-level processing is being used, then calculate the temperature,
!!   pressure and specific humidity on the pseudo-levels by interpolation.  Then
!!   calculate the refractivity on all levels.
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 6 May 2021
!!
!-------------------------------------------------------------------------------
 
SUBROUTINE ufo_calculate_refractivity (nlevP,           &
                                       nlevq,           &
                                       za,              &
                                       zb,              &
                                       P,               &
                                       q,               &
                                       vert_interp_ops, &
                                       pseudo_ops,      &
                                       min_temp_grad,   &
                                       refracerr,       &
                                       nRefLevels,      &
                                       refractivity,    &
                                       model_heights,   &
                                       temperature,     &
                                       interp_pressure)
 
IMPLICIT NONE
 
! Subroutine arguments:
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(nlevp)               !< heights of rho levs
REAL(kind_real), INTENT(IN)               :: zb(nlevq)               !< heights of theta levs
REAL(kind_real), INTENT(IN)               :: p(nlevp)                !< state vector
REAL(kind_real), INTENT(IN)               :: q(nlevq)                !< state vector
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
LOGICAL, INTENT(OUT)                      :: refracerr               !< refractivity error
INTEGER, INTENT(OUT)                      :: nreflevels              !< no. of pseudo levs
REAL(kind_real), ALLOCATABLE, INTENT(OUT) :: refractivity(:)         !< Ref. on pseudo levs
REAL(kind_real), ALLOCATABLE, INTENT(OUT) :: model_heights(:)        !< height of pseudo levs
REAL(kind_real), OPTIONAL, INTENT(OUT)    :: temperature(nlevq)      !< Calculated temperature on model levels
REAL(kind_real), OPTIONAL, INTENT(OUT)    :: interp_pressure(nlevq)  !< Model pressure, interpolated to temperature levels
 
! Local declarations:
CHARACTER(len=*), PARAMETER               :: routinename = "ufo_calculate_refractivity"
INTEGER                                   :: i
INTEGER                                   :: counter
REAL(kind_real)                           :: exner(nlevp)
REAL(kind_real)                           :: pb(nlevq)
REAL(kind_real)                           :: t_local(nlevq)      ! Temp. on theta levs
REAL(kind_real)                           :: t_virtual
REAL(kind_real)                           :: ex_theta
REAL(kind_real)                           :: pwt1
REAL(kind_real)                           :: pwt2
REAL(kind_real)                           :: ndry
REAL(kind_real)                           :: nwet
REAL(kind_real), ALLOCATABLE              :: p_pseudo(:)
REAL(kind_real), ALLOCATABLE              :: q_pseudo(:)
REAL(kind_real), ALLOCATABLE              :: t_pseudo(:)
REAL(kind_real)                           :: refracmodel(1:nlevq)
REAL(kind_real)                           :: gamma
REAL(kind_real)                           :: beta
REAL(kind_real)                           :: c ! continuity constant for hydrostatic pressure
CHARACTER(LEN=200)                        :: message   ! Message to be output to user
 
! Allocate arrays for pseudo-level processing and output
IF (pseudo_ops) THEN
  nreflevels = 2 * nlevq - 1
  ALLOCATE (p_pseudo(nreflevels))
  ALLOCATE (q_pseudo(nreflevels))
  ALLOCATE (t_pseudo(nreflevels))
ELSE
  nreflevels = nlevq
END IF
ALLOCATE (model_heights(nreflevels))
ALLOCATE (refractivity(nreflevels))
 
! Set up the P and q vectors from x
 
! Initialise refractivity arrays to missing Data
refracmodel(:) = missing_value(refracmodel(1))
refractivity(:) = missing_value(refractivity(1))
t_local(:) = missing_value(t_local(1))
refracerr = .false.
 
DO i = 1, nlevq
  IF (p(i) == missing_value(p(i))) THEN  ! pressure missing
    refracerr = .true.
    WRITE(message, *) routinename, " Input pressure missing", i
    CALL fckit_log % warning(message)
    EXIT
  END IF
 
  IF (p(i) - p(i + 1) < 0.0) THEN  ! or non-monotonic pressure
    refracerr = .true.
    WRITE(message, *) routinename, " Input pressure non-monotonic", i, p(i), p(i+1)
    CALL fckit_log % warning(message)
    EXIT
  END IF
END DO
 
IF (any(p(:) <= 0.0)) THEN        ! pressure zero or negative
  refracerr = .true.
  WRITE(message, *) routinename, " Input pressure not physical"
  CALL fckit_log % warning(message)
END IF
 
! only proceed if pressure is valid
IF (.NOT. refracerr) THEN
 
  ! Calculate exner on rho levels.
 
  exner(:) = (p(:) / pref) ** rd_over_cp
 
  ! Calculate the refractivity on the same model levels as specific humidity
  DO i = 1, nlevq
    ! Calc. pressure pb
    pwt1 = (za(i + 1) - zb(i)) / (za(i + 1) - za(i))
 
    pwt2 = 1.0 - pwt1
 
    ! Calculate the pressure on the theta level.
    IF (vert_interp_ops) THEN
      ! Assume ln(P) linear with height
      pb(i) = exp(pwt1 * log(p(i)) + pwt2 * log(p(i + 1)))
    ELSE
      ! Assume Exner varies linearly with height
      pb(i) = pref * (pwt1 * (p(i) / pref) ** rd_over_cp + pwt2 * (p(i + 1) / pref) ** rd_over_cp) ** (1.0 / rd_over_cp)
    END IF
 
    ! Calculate the Exner on the theta level.
    ex_theta = (pb(i) / pref) ** rd_over_cp
 
    ! Calculate mean layer T_virtual on staggered vertical levels
 
    t_virtual = grav * (za(i + 1) - za(i)) * ex_theta / &
        (cp * (exner(i) - exner(i + 1)))
 
    t_local(i) = t_virtual / (1.0 + c_virtual * q(i))
 
    ! Wet component
    nwet = n_beta * pb(i) * q(i) / (t_local(i) ** 2 * (mw_ratio + (1.0 - mw_ratio) * q(i)))
 
    ndry = n_alpha * pb(i) / t_local(i)
 
    refracmodel(i) = ndry + nwet
 
  END DO
 
  ! Do pseudo-level processing
  IF (pseudo_ops) THEN
    counter = 1
    DO i = 1, nreflevels
 
      ! Odd 'i' (i.e. copies of actual model level values)
      IF (mod(i, 2) > 0) THEN
        model_heights(i) = zb(counter)
        p_pseudo(i) = pb(counter)
        q_pseudo(i) = q(counter)
        t_pseudo(i) = t_local(counter)
        counter = counter + 1
 
      ! Even 'i' (i.e. intermediate pseudo-levels)
      ELSE
        model_heights(i) = (zb(counter - 1) + zb(counter)) / 2.0
 
        ! Assume exponential variation when humidities are positive
        IF (min(q(counter - 1), q(counter)) > 0.0) THEN
          gamma = log(q(counter - 1) / q(counter)) / (zb(counter) - zb(counter - 1))
          q_pseudo(i) = q(counter - 1) * exp(-gamma * (model_heights(i) - model_heights(i - 1)))
 
        ! Assume linear variation if humidities are -ve
        ELSE
          q_pseudo(i) = q(counter - 1) + (q(counter) - q(counter - 1)) / (zb(counter) - &
                        zb(counter - 1)) * (model_heights(i) - zb(counter - 1))
        END IF
 
        ! T varies linearly with height
        beta = (t_local(counter) - t_local(counter - 1)) / (zb(counter) - zb(counter - 1))
        t_pseudo(i) = t_local(counter - 1) + beta * (model_heights(i) - zb(counter - 1))
 
        ! Pressure varies to maintain hydrostatic balance
        IF (abs(t_local(counter) - t_local(counter - 1)) > min_temp_grad) THEN
          c = ((pb(counter) / pb(counter - 1)) * (t_local(counter) / t_local(counter - 1)) ** (grav / (rd * beta)) - &
              1.0) / (zb(counter) - zb(counter - 1))
          p_pseudo(i) = (pb(counter - 1) * (t_pseudo(i) / t_local(counter - 1)) ** &
                      (-grav / (rd * beta))) * (1.0 + c * (model_heights(i) - zb(counter - 1)))
        ELSE
          ! If layer is isothermal, explicitly force P to vary exponentially to avoid singularity
          p_pseudo(i) = pb(counter - 1) * exp(log(pb(counter) / pb(counter - 1)) * &
                      ((model_heights(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1))))
        END IF
      END IF
    END DO
 
    refractivity = n_alpha * p_pseudo / t_pseudo + n_beta * p_pseudo * q_pseudo / &
               (t_pseudo ** 2 * (mw_ratio + (1.0 - mw_ratio) * q_pseudo))
  ELSE
    refractivity = refracmodel
    model_heights = zb
  END IF
 
END IF
 
IF (PRESENT(temperature)) THEN
  temperature(:) = t_local(:)
END IF
 
IF (PRESENT(interp_pressure)) THEN
  interp_pressure(:) = pb(:)
END IF
 
IF (pseudo_ops) THEN
  IF (ALLOCATED (p_pseudo)) DEALLOCATE (p_pseudo)
  IF (ALLOCATED (q_pseudo)) DEALLOCATE (q_pseudo)
  IF (ALLOCATED (t_pseudo)) DEALLOCATE (t_pseudo)
END IF
 
END SUBROUTINE ufo_calculate_refractivity
 
!-------------------------------------------------------------------------------
!> \brief Calculate general refractivity K matrix.
!!
!! \details **ufo_refractivity_kmat**
!! * Allocate intermediate matrices used in calculation.
!! * Calculate the refractivity and basic gradients on the temperature/theta
!!   levels.
!! * If using pseudo-levels, then calculate the partial matrices on these
!!   levels, and use these to calculate the final matrices (dref_dp and dref_dq).
!! * If not using pseudo-levels, then calculate the final matrices.
!! * Deallocate the temporary matrices.
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 6 May 2021
!!
!-------------------------------------------------------------------------------
 
!-------------------------------------------------------------------------------
!> Calculate general refractivity K matrix.
!-------------------------------------------------------------------------------
 
SUBROUTINE ufo_refractivity_kmat(nlevP,           &
                                 nlevq,           &
                                 nRefLevels,      &
                                 za,              &
                                 zb,              &
                                 P,               &
                                 q,               &
                                 pseudo_ops,      &
                                 vert_interp_ops, &
                                 min_temp_grad,   &
                                 dref_dP,         &
                                 dref_dq,         &
                                 dPb_dP,          &
                                 refractivity)
 
IMPLICIT NONE
 
! Subroutine arguments:
INTEGER, INTENT(IN)                       :: nlevp                     !< no. of pressure levels
INTEGER, INTENT(IN)                       :: nlevq                     !< no. of specific humidity levels
INTEGER, INTENT(IN)                       :: nreflevels                !< no. of refractivity levels
REAL(kind_real), INTENT(IN)               :: za(nlevp)                 !< Height of the pressure levels
REAL(kind_real), INTENT(IN)               :: zb(nlevq)                 !< Height of the specific humidity levels
REAL(kind_real), INTENT(IN)               :: p(nlevp)                  !< Pressure
REAL(kind_real), INTENT(IN)               :: q(nlevq)                  !< Specific humidity
LOGICAL, INTENT(IN)                       :: pseudo_ops                !< Whether to use pseudo-levels in calculation
LOGICAL, INTENT(IN)                       :: vert_interp_ops           !< Whether to interpolate vertically using exner or ln(p)
REAL(kind_real), INTENT(IN)               :: min_temp_grad             !< Minimum value for the vertical temperature gradient
REAL(kind_real), ALLOCATABLE, INTENT(OUT) :: dref_dp(:,:)              !< kmatrix for p
REAL(kind_real), ALLOCATABLE, INTENT(OUT) :: dref_dq(:,:)              !< kmatrix for q
REAL(kind_real), OPTIONAL, INTENT(OUT)    :: dpb_dp(nlevq,nlevp)       !< Gradient of pressure on theta levels wrt pressure on pressure levels
REAL(kind_real), OPTIONAL, INTENT(OUT), ALLOCATABLE :: refractivity(:) !< Calculated refractivity
 
! Local declarations:
INTEGER                                   :: i
INTEGER                                   :: counter
REAL(kind_real)                           :: pb(nlevq)
REAL(kind_real)                           :: t_virtual(nlevq)
REAL(kind_real)                           :: t(nlevq)
REAL(kind_real)                           :: refrac(nlevq)                   ! Refractivity on model levels
REAL(kind_real)                           :: extheta
REAL(kind_real)                           :: pwt1
REAL(kind_real)                           :: pwt2
REAL(kind_real)                           :: ndry
REAL(kind_real)                           :: nwet
REAL(kind_real)                           :: dex_dp(nlevp,nlevp)
REAL(kind_real)                           :: dpb_dp_local(nlevq,nlevp)
REAL(kind_real)                           :: dextheta_dpb(nlevq,nlevq)
REAL(kind_real)                           :: dtv_dextheta(nlevq,nlevq)
REAL(kind_real)                           :: dtv_dex(nlevq,nlevp)
REAL(kind_real)                           :: dt_dtv(nlevq,nlevq)
REAL(kind_real)                           :: dt_dq(nlevq,nlevq)
REAL(kind_real)                           :: dref_dpb(nlevq,nlevq)
REAL(kind_real)                           :: dref_dt(nlevq,nlevq)
REAL(kind_real)                           :: dref_dq_model(nlevq,nlevq)      ! Gradient of refractivity wrt specific humidity on model levels
REAL(kind_real)                           :: m1(nreflevels,nlevq)
REAL(kind_real)                           :: m2(nreflevels,nlevp)
REAL(kind_real)                           :: m3(nreflevels,nlevq)
REAL(kind_real)                           :: m4(nreflevels,nlevq)
REAL(kind_real), ALLOCATABLE              :: p_pseudo(:)
REAL(kind_real), ALLOCATABLE              :: q_pseudo(:)
REAL(kind_real), ALLOCATABLE              :: t_pseudo(:)
REAL(kind_real), ALLOCATABLE              :: model_heights(:)
REAL(kind_real)                           :: gamma
REAL(kind_real)                           :: beta
REAL(kind_real)                           :: c           ! continuity constant for hydrostatic pressure
REAL(kind_real)                           :: g_rb        ! Frequently used term
REAL(kind_real)                           :: c_zz        ! Frequently used term
REAL(kind_real)                           :: dpp_dt1     ! dP_pseudo / dT_below
REAL(kind_real)                           :: dpp_dtp     ! dP_pseudo / dT_pseudo
REAL(kind_real)                           :: dtp_dt1     ! dT_pseudo / dT_below
REAL(kind_real)                           :: dpp_dbeta   ! dP_pseudo / dbeta
REAL(kind_real)                           :: dbeta_dt1   ! dbeta / dT_below
REAL(kind_real)                           :: dtp_dbeta   ! dT_pseudo / dbeta
REAL(kind_real)                           :: dbeta_dt2   ! dbeta / dT_above
REAL(kind_real)                           :: dpp_dc      ! dP_pseudo / 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)                           :: dpp_dp1     ! dP_pseudo / dP_below
REAL(kind_real)                           :: dc_dp1      ! dc / dP_below
REAL(kind_real)                           :: dc_dp2      ! dc / dP_above
REAL(kind_real), ALLOCATABLE              :: dref_dppseudo(:,:)
REAL(kind_real), ALLOCATABLE              :: dref_dtpseudo(:,:)
REAL(kind_real), ALLOCATABLE              :: dref_dqpseudo(:,:)
REAL(kind_real), ALLOCATABLE              :: dppseudo_dpb(:,:)
REAL(kind_real), ALLOCATABLE              :: dtpseudo_dtb(:,:)
REAL(kind_real), ALLOCATABLE              :: dqpseudo_dqb(:,:)
REAL(kind_real), ALLOCATABLE              :: dppseudo_dt(:,:)
REAL(kind_real), ALLOCATABLE              :: dtb_dp(:,:)
REAL(kind_real), ALLOCATABLE              :: dppseudo_dp(:,:)
REAL(kind_real), ALLOCATABLE              :: dtpseudo_dp(:,:)
 
ALLOCATE (dref_dp(nreflevels,nlevp))
ALLOCATE (dref_dq(nreflevels,nlevq))
 
IF (pseudo_ops) THEN
  ALLOCATE (p_pseudo(nreflevels))
  ALLOCATE (q_pseudo(nreflevels))
  ALLOCATE (t_pseudo(nreflevels))
  ALLOCATE (model_heights(nreflevels))
  IF (PRESENT(refractivity)) ALLOCATE (refractivity(nreflevels))
 
  ALLOCATE (dref_dppseudo(nreflevels,nreflevels))
  ALLOCATE (dref_dtpseudo(nreflevels,nreflevels))
  ALLOCATE (dref_dqpseudo(nreflevels,nreflevels))
  ALLOCATE (dppseudo_dpb(nreflevels,nlevq))
  ALLOCATE (dppseudo_dt(nreflevels,nlevq))
  ALLOCATE (dtpseudo_dtb(nreflevels,nlevq))
  ALLOCATE (dqpseudo_dqb(nreflevels,nlevq))
 
  ALLOCATE (dtb_dp(nlevq,nlevp))
  ALLOCATE (dppseudo_dp(nreflevels,nlevp))
  ALLOCATE (dtpseudo_dp(nreflevels,nlevp))
 
  dref_dppseudo(:,:) = 0.0
  dref_dtpseudo(:,:) = 0.0
  dref_dqpseudo(:,:) = 0.0
  dppseudo_dpb(:,:) = 0.0
  dppseudo_dt(:,:) = 0.0
  dtpseudo_dtb(:,:) = 0.0
  dqpseudo_dqb(:,:) = 0.0
  dppseudo_dp(:,:) = 0.0
  dtpseudo_dp(:,:) = 0.0
END IF
 
!-----------------------
! 1. Initialise matrices
!-----------------------
 
dref_dp(:,:) = 0.0
 
call ufo_refractivity_partial_derivatives(nlevp,           &
                                          nlevq,           &
                                          za,              &
                                          zb,              &
                                          p,               &
                                          q,               &
                                          vert_interp_ops, &
                                          dt_dtv,          &
                                          dt_dq,           &
                                          dref_dpb,        &
                                          dref_dt,         &
                                          dref_dq_model,   &
                                          refrac,          &
                                          t,               &
                                          pb,              &
                                          dex_dp,          &
                                          dextheta_dpb,    &
                                          dtv_dextheta,    &
                                          dpb_dp_local,    &
                                          dtv_dex)
 
IF (pseudo_ops) THEN
  !----------------------------------!
  !- Add intermediate pseudo-levels -!
  !----------------------------------!
  dtb_dp = matmul(dt_dtv, matmul(matmul(dtv_dextheta, dextheta_dpb), dpb_dp_local) + matmul(dtv_dex, dex_dp))
  counter = 1
  DO i = 1, nreflevels
    ! Odd 'i' (i.e. copies of actual model level values)
    IF (mod(i, 2) > 0) THEN
      model_heights(i) = zb(counter)
      p_pseudo(i) = pb(counter)
      q_pseudo(i) = q(counter)
      t_pseudo(i) = t(counter)
      IF (PRESENT(refractivity)) refractivity(i) = n_alpha * p_pseudo(i) / t_pseudo(i) + &
          n_beta * p_pseudo(i) * q_pseudo(i) / &
          (t_pseudo(i) ** 2 * (mw_ratio + (1.0 - mw_ratio) * q_pseudo(i)))
 
      dref_dppseudo(i,i) = dref_dpb(counter,counter)
      dref_dtpseudo(i,i) = dref_dt(counter,counter)
      dref_dqpseudo(i,i) = dref_dq_model(counter,counter)
 
      dppseudo_dpb(i,counter) = 1.0
      dtpseudo_dtb(i,counter) = 1.0
      dqpseudo_dqb(i,counter) = 1.0
 
      counter = counter + 1
 
    ! Even 'i' (i.e. intermediate pseudo-levels)
    ELSE
      model_heights(i) = (zb(counter - 1) + zb(counter)) / 2.0
      IF (min(q(counter - 1), q(counter)) > 0.0) THEN
        gamma = log(q(counter - 1) / q(counter)) / (zb(counter) - zb(counter - 1))
        q_pseudo(i) = q(counter - 1) * exp(-gamma * (model_heights(i) - model_heights(i - 1)))
      ELSE
        q_pseudo(i) = q(counter - 1) + (q(counter) - q(counter - 1)) / &
                     (zb(counter) - zb(counter - 1)) * (model_heights(i) - zb(counter - 1))
      END IF
 
      beta = (t(counter) - t(counter - 1)) / (zb(counter) - zb(counter - 1))
      t_pseudo(i) = t(counter - 1) + beta * (model_heights(i) - zb(counter - 1))
      IF (abs(t(counter) - t(counter - 1)) > min_temp_grad) THEN
        c = ((pb(counter) / pb(counter - 1)) * (t(counter) / t(counter - 1)) ** (grav / (rd * beta)) - 1.0) / &
             (zb(counter) - zb(counter - 1))
        p_pseudo(i) = (pb(counter - 1) * (t_pseudo(i) / t(counter - 1)) ** &
                    (-grav / (rd * beta))) * (1.0 + c * (model_heights(i) - zb(counter - 1)))
      ELSE
        p_pseudo(i) = pb(counter - 1) * exp(log(pb(counter) / pb(counter - 1)) * &
                    ((model_heights(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1))))
      END IF
 
      ndry = n_alpha * p_pseudo(i) / t_pseudo(i)
      nwet = n_beta * p_pseudo(i) * q_pseudo(i) / (t_pseudo(i) ** 2 &
                    *(mw_ratio + (1.0 - mw_ratio) * q_pseudo(i)))
 
      IF (PRESENT(refractivity)) refractivity(i) = ndry + nwet
 
      dref_dppseudo(i,i) = (ndry + nwet) / p_pseudo(i)
      dref_dtpseudo(i,i) = -(ndry + 2.0 * nwet) / t_pseudo(i)
      dref_dqpseudo(i,i) = n_beta * p_pseudo(i) * mw_ratio / (t_pseudo(i) * (mw_ratio + &
                           (1.0 - mw_ratio) * q_pseudo(i))) ** 2
 
      ! Transform P, q and T to pseudo-levels
      IF (min(q(counter - 1), q(counter)) > 0.0) THEN
        dqpseudo_dqb(i,counter - 1) = (q_pseudo(i) / q(counter - 1)) * (1.0 - &
             (model_heights(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1)))
 
        dqpseudo_dqb(i,counter) = (q_pseudo(i) / q(counter)) * ((model_heights(i) - zb(counter - 1)) / &
                                  (zb(counter) - zb(counter - 1)))
      ELSE
        dqpseudo_dqb(i,counter - 1) = (zb(counter) - model_heights(i)) / (zb(counter) - zb(counter - 1))
 
        dqpseudo_dqb(i,counter) = (model_heights(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1))
      END IF
 
      dtpseudo_dtb(i,counter - 1) = (zb(counter) - model_heights(i)) / (zb(counter) - zb(counter - 1))
 
      dtpseudo_dtb(i,counter) = (model_heights(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1))
 
      ! Items that need changing for dPpseudo/dTb:
 
      g_rb = grav / (rd * beta)
      c_zz = c + 1.0 / (zb(counter) - zb(counter - 1))
 
      dpp_dt1 = (g_rb / t(counter - 1)) * p_pseudo(i)
      dpp_dtp = -(g_rb / t_pseudo(i)) * p_pseudo(i)
      dtp_dt1 = dtpseudo_dtb(i,counter - 1)
      dpp_dbeta = (g_rb / beta) * log(t_pseudo(i) / t(counter - 1)) * p_pseudo(i)
      dbeta_dt1 = -1.0 / (zb(counter) - zb(counter - 1))
      dtp_dbeta = model_heights(i) - zb(counter - 1)
      dbeta_dt2 = 1.0 / (zb(counter) - zb(counter - 1))
 
      IF (abs(t(counter) - t(counter - 1)) > min_temp_grad) THEN
        ! Incomplete computation of dPpseudo_dPb. Temperature derivatives done below
        dpp_dc = (model_heights(i) - zb(counter - 1)) * pb(counter - 1) * (t_pseudo(i) / t(counter - 1)) ** (-g_rb)
        dc_dt1 = -(g_rb / (t(counter - 1))) * c_zz
        dc_dbeta = -(g_rb / beta) * log(t(counter) / t(counter - 1)) * c_zz
        dc_dt2 = (g_rb / t(counter)) * c_zz
        dpp_dp1 = p_pseudo(i) / pb(counter - 1)
        dppseudo_dt(i,counter - 1) = dpp_dt1 + dpp_dtp * dtp_dt1 + dpp_dbeta * dbeta_dt1 + dpp_dc * &
                                   (dc_dt1 + dc_dbeta * dbeta_dt1)
        dppseudo_dt(i,counter) = (dpp_dbeta + dpp_dtp * dtp_dbeta + dpp_dc * dc_dbeta) * &
                                  dbeta_dt2 + dpp_dc * dc_dt2
 
        dc_dp1 = -(1.0 / pb(counter - 1)) * c_zz
        dc_dp2 =  (1.0 / pb(counter)) * c_zz
 
        dppseudo_dpb(i,counter - 1) = dpp_dp1 + dpp_dc * dc_dp1
        dppseudo_dpb(i,counter) = dpp_dc * dc_dp2
      ELSE
        dppseudo_dpb(i,counter - 1) = exp(log(pb(counter) / pb(counter - 1)) * ((model_heights(i) - zb(counter - 1)) / &
          (zb(counter) - zb(counter - 1)))) * (1.0 - ((model_heights(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1))))
        dppseudo_dpb(i,counter) = (pb(counter - 1) / pb(counter)) * ((model_heights(i) - zb(counter - 1)) / &
          (zb(counter) - zb(counter - 1))) * exp(log(pb(counter) / pb(counter - 1)) * ((model_heights(i) - zb(counter - 1)) / &
          (zb(counter) - zb(counter - 1))))
      END IF
 
    END IF
 
  END DO
 
  ! temperature derivatives:
  dppseudo_dp = matmul(dppseudo_dpb, dpb_dp_local) + matmul(dppseudo_dt, dtb_dp)
  dtpseudo_dp = matmul(dtpseudo_dtb, dtb_dp)
 
  !-------------------------------------------------
  ! 3. Evaluate the Kmatrix by matrix multiplication
  !-------------------------------------------------
 
  ! calc K matrix for P on rho levels
  dref_dp(:,:) = matmul(dref_dppseudo, dppseudo_dp) + matmul(dref_dtpseudo, dtpseudo_dp)
 
  ! calc Kmatrix for q on theta levels
  dref_dq(:,:) = matmul(dref_dqpseudo, dqpseudo_dqb) + matmul(matmul(dref_dtpseudo, dtpseudo_dtb), dt_dq) + &
                 matmul(matmul(dref_dppseudo, dppseudo_dt), dt_dq)
 
! Normal model levels
ELSE
 
  if (PRESENT(refractivity)) refractivity = refrac
 
  !-------------------------------------------------
  ! 3. Evaluate the Kmatrix by matrix multiplication
  !-------------------------------------------------
 
  ! calc K matrix for P on rho levels
  !  dNmod/dP = (dNmod/dPb * dPb/dP)   + ....
  dref_dp(:,:) = matmul(dref_dpb, dpb_dp_local)
 
  !  .... (dNmod/dT * dT/dTv * dTv/dEx *dEx/dP) + .....
  m1(:,:) = matmul(dref_dt, dt_dtv)
  m2(:,:) = matmul(m1, dtv_dex)
  dref_dp(:,:) = dref_dp(:,:) + matmul(m2, dex_dp)
 
  !  .... (dNmod/dT * dT/dTv * dTv/dExtheta *dExtheta/dPb*dPb/dP)
  m3(:,:) = matmul(m1, dtv_dextheta)
  m4(:,:) = matmul(m3, dextheta_dpb)
  dref_dp(:,:) = dref_dp(:,:) + matmul(m4, dpb_dp_local)
 
  ! calc Kmatrix for q on theta levels
  !  dNmod/dq = (dNmod/dq) + (dNmod/dT*dT/dq)
  dref_dq(:,:) = dref_dq_model(:,:) + matmul(dref_dt, dt_dq)
 
END IF
 
if (present(dpb_dp)) dpb_dp(:,:) = dpb_dp_local(:,:)
 
IF (ALLOCATED (p_pseudo)) DEALLOCATE (p_pseudo)
IF (ALLOCATED (q_pseudo)) DEALLOCATE (q_pseudo)
IF (ALLOCATED (t_pseudo)) DEALLOCATE (t_pseudo)
IF (ALLOCATED (model_heights)) DEALLOCATE (model_heights)
IF (ALLOCATED (dref_dppseudo)) DEALLOCATE (dref_dppseudo)
IF (ALLOCATED (dref_dtpseudo)) DEALLOCATE (dref_dtpseudo)
IF (ALLOCATED (dref_dqpseudo)) DEALLOCATE (dref_dqpseudo)
IF (ALLOCATED (dppseudo_dpb)) DEALLOCATE (dppseudo_dpb)
IF (ALLOCATED (dtpseudo_dtb)) DEALLOCATE (dtpseudo_dtb)
IF (ALLOCATED (dqpseudo_dqb)) DEALLOCATE (dqpseudo_dqb)
IF (ALLOCATED (dppseudo_dt)) DEALLOCATE (dppseudo_dt)
IF (ALLOCATED (dtb_dp)) DEALLOCATE (dtb_dp)
IF (ALLOCATED (dppseudo_dp)) DEALLOCATE (dppseudo_dp)
IF (ALLOCATED (dtpseudo_dp)) DEALLOCATE (dtpseudo_dp)
 
END SUBROUTINE ufo_refractivity_kmat
 
!-------------------------------------------------------------------------------
!> \brief Calculate some partial derivatives of refractivity on model levels
!!
!! \details **ufo_refractivity_partial_derivatives**
!! * Calculate the pressure on model theta levels
!! * Calculate exner on model theta level
!! * Calculate mean layer T_virtual, and then various partial derivatives
!!
!! \author Neill Bowler (Met Office)
!!
!! \date 26 May 2021
!!
!-------------------------------------------------------------------------------
 
subroutine ufo_refractivity_partial_derivatives(nlevP,           &
                                                nlevq,           &
                                                za,              &
                                                zb,              &
                                                P,               &
                                                q,               &
                                                vert_interp_ops, &
                                                dT_dTv,          &
                                                dT_dq,           &
                                                dref_dPb,        &
                                                dref_dT,         &
                                                dref_dq,         &
                                                refractivity,    &
                                                T,               &
                                                Pb,              &
                                                dEx_dP,          &
                                                dExtheta_dPb,    &
                                                dTv_dExtheta,    &
                                                dPb_dP_local,    &
                                                dTv_dEx)
 
IMPLICIT NONE
 
! Subroutine arguments:
INTEGER, INTENT(IN)                       :: nlevP                     !< no. of pressure levels
INTEGER, INTENT(IN)                       :: nlevq                     !< no. of specific humidity levels
REAL(kind_real), INTENT(IN)               :: za(nlevp)                 !< Height of the pressure levels
REAL(kind_real), INTENT(IN)               :: zb(nlevq)                 !< Height of the specific humidity levels
REAL(kind_real), INTENT(IN)               :: P(nlevp)                  !< Pressure
REAL(kind_real), INTENT(IN)               :: q(nlevq)                  !< Specific humidity
LOGICAL, INTENT(IN)                       :: vert_interp_ops           !< Whether to interpolate vertically using exner or ln(p)
REAL(kind_real), INTENT(OUT)              :: dT_dTv(nlevq,nlevq)       !< Partial derivative of temperature wrt virtual temperature
REAL(kind_real), INTENT(OUT)              :: dT_dq(nlevq,nlevq)        !< Partial derivative of temperature wrt specific humidity
REAL(kind_real), INTENT(OUT)              :: dref_dPb(nlevq,nlevq)     !< Partial derivative of refractivity wrt pressure on model theta levels
REAL(kind_real), INTENT(OUT)              :: dref_dT(nlevq,nlevq)      !< Partial derivative of refractivity wrt temperature
REAL(kind_real), INTENT(OUT)              :: dref_dq(nlevq,nlevq)      !< Partial derivative of refractivity wrt specific humidity
REAL(kind_real), INTENT(OUT)              :: refractivity(nlevq)       !< Calculated refractivity
REAL(kind_real), INTENT(OUT)              :: T(nlevq)                  !< Calculated temperature
REAL(kind_real), INTENT(OUT)              :: Pb(nlevq)                 !< Pressure on model theta levels
REAL(kind_real), INTENT(OUT)              :: dEx_dP(nlevP,nlevP)       !< Partial derivative of exner wrt pressure
REAL(kind_real), INTENT(OUT)              :: dExtheta_dPb(nlevq,nlevq) !< Partial derivative of refractivity wrt pressure (at ob location)
REAL(kind_real), INTENT(OUT)              :: dTv_dExtheta(nlevq,nlevq) !< Virtual temperature divided by exner on theta levels
REAL(kind_real), INTENT(OUT)              :: dPb_dP_local(nlevq,nlevP) !< Partial derivative of pressure on theta levels wrt pressure on pressure levels
REAL(kind_real), INTENT(OUT)              :: dTv_dEx(nlevq,nlevP)      !< Partial derivative of virtual temperature wrt exner
 
! Local declarations:
INTEGER                                   :: i                         ! Loop variable
REAL(kind_real)                           :: Exner(nlevP)              ! Exner on model pressure levels
REAL(kind_real)                           :: T_virtual(nlevq)          ! Virtual temperature
REAL(kind_real)                           :: Extheta                   ! Exner on model theta levels
REAL(kind_real)                           :: pwt1                      ! Weight given to the lower model level in interpolation
REAL(kind_real)                           :: pwt2                      ! Weight given to the upper model level in interpolation
REAL(kind_real)                           :: Ndry                      ! Contribution to refractivity from dry terms
REAL(kind_real)                           :: Nwet                      ! Contribution to refractivity from wet terms
 
!-----------------------
! 1. Initialise matrices
!-----------------------
 
dpb_dp_local(:,:) = 0.0
dextheta_dpb(:,:) = 0.0
dtv_dextheta(:,:) = 0.0
dtv_dex(:,:) = 0.0
dt_dtv(:,:) = 0.0
dt_dq(:,:) = 0.0
dref_dpb(:,:) = 0.0
dref_dt(:,:) = 0.0
dref_dq(:,:) = 0.0
dex_dp(:,:) = 0.0
 
! Calculate exner on rho levels.
exner(:) = (p(:) / pref) ** rd_over_cp
DO i = 1,nlevp
  dex_dp(i,i) = rd_over_cp / pref * (p(i) / pref) ** (rd_over_cp - 1.0)
END DO
 
!----------------------------------------------
! 2. Calculate the refractivity on the temperature/theta levels
!----------------------------------------------
 
DO i = 1, nlevq
 
  pwt1 = (za(i + 1) - zb(i)) / (za(i + 1) - za(i))
 
  pwt2 = 1.0 - pwt1
 
  ! calculate the pressure on the theta level.
  IF (vert_interp_ops) THEN
    pb(i) = exp(pwt1 * log(p(i)) + pwt2 * log(p(i + 1)))
 
    dpb_dp_local(i,i) = pb(i) * pwt1 / p(i)
    dpb_dp_local(i,i + 1) = pb(i) * pwt2 / p(i + 1)
  ELSE
    ! Assume Exner varies linearly with height
    pb(i) = pref * (pwt1 * (p(i) / pref) ** rd_over_cp + pwt2 * (p(i + 1) / pref) ** rd_over_cp) ** (1.0 / rd_over_cp)
 
    dpb_dp_local(i,i) = pwt1 * (pwt1 * (p(i) / pref) ** rd_over_cp + pwt2 * &
     (p(i + 1) / pref) ** rd_over_cp) ** (1.0 / rd_over_cp - 1.0) * (p(i) / pref) ** (rd_over_cp - 1.0)
    dpb_dp_local(i,i + 1) = pwt2 * (pwt1 * (p(i) / pref) ** rd_over_cp + pwt2 * &
     (p(i + 1) / pref) ** rd_over_cp) ** (1.0 / rd_over_cp - 1.0) * (p(i + 1) / pref) ** (rd_over_cp-1.0)
  END IF
 
  ! calculate Exner on the theta level.
 
  extheta = (pb(i) / pref) ** rd_over_cp
 
  dextheta_dpb(i,i) = rd_over_cp * (pb(i) ** (rd_over_cp - 1.0)) / (pref ** rd_over_cp)
 
  ! Calculate mean layer T_virtual on staggered vertical levels
 
  t_virtual(i) = grav * (za(i + 1) - za(i)) * extheta / (cp * (exner(i) - exner(i + 1)))
 
  dtv_dextheta(i,i) = t_virtual(i) / extheta
 
  dtv_dex(i,i) = -t_virtual(i) / (exner(i) - exner(i + 1))
 
  dtv_dex(i,i + 1) = t_virtual(i) / (exner(i) - exner(i + 1))
 
  t(i) = t_virtual(i) / (1.0 + c_virtual * q(i))
 
  dt_dtv(i,i) = 1.0 / (1.0 + c_virtual * q(i))
 
  dt_dq(i,i) = -c_virtual * t(i) / (1.0 + c_virtual * q(i))
 
  ! wet component
 
  nwet = n_beta * pb(i) * q(i) / (t(i) ** 2 * (mw_ratio + (1.0 - mw_ratio) * q(i)))
 
  dref_dq(i,i) = n_beta * pb(i) * mw_ratio / (t(i) * (mw_ratio + (1.0 - mw_ratio) * q(i))) ** 2
 
  ndry = n_alpha * pb(i) / t(i)
 
  refractivity(i) = ndry + nwet
 
  dref_dpb(i,i) = (ndry + nwet) / pb(i)
 
  dref_dt(i,i) = -(ndry + 2.0 * nwet) / t(i)
 
END DO
 
end subroutine ufo_refractivity_partial_derivatives
 
 
end module ufo_utils_refractivity_calculator