ufo_groundgnss_metoffice_utils_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.
!-------------------------------------------------------------------------------
 
module ufo_groundgnss_ukmo_utils_mod
 
!use iso_c_binding
use fckit_log_module, only: fckit_log
use kinds,            only: kind_real
 
! Generic routines from elsewhere in jedi
use missing_values_mod
use ufo_constants_mod, only: &
    rd,                      &   ! Gas constant for dry air
    grav,                    &   ! Gravitational field strength
    n_alpha,                 &   ! Refractivity constant a
    n_beta                       ! Refractivity constant b
 
implicit none
public                      :: ops_groundgnss_ztd
public                      :: ops_groundgnss_topcorrection
REAL(kind_real), PARAMETER  :: refrac_scale = 1.0e-6
private
 
contains
 
!-------------------------------------------------------------------------------
! Ground based GNSS Observation operator
!-------------------------------------------------------------------------------
SUBROUTINE ops_groundgnss_ztd  (nlevq,    &
                                refrac,   &
                                zb,       &
                                zStation, &
                                Model_ZTD)
 
IMPLICIT NONE
    
    INTEGER, INTENT(IN)               :: nlevq           ! no. of temperature/theta levels
    REAL(kind_real), INTENT(IN)       :: refrac(:)       ! refractivty
    REAL(kind_real), INTENT(IN)       :: zb(:)           ! The geometric height of the model theta levels
    REAL(kind_real), INTENT(IN)       :: zstation        ! Station heights
    REAL(kind_real), INTENT(INOUT)    :: model_ztd       ! Model background ZTD
    
    
    REAL(kind_real)             :: localzenithdelay      ! Zenith total delay
    INTEGER                     :: level                 ! level iterator
    REAL(kind_real)             :: stationrefrac         ! refrac at station
    REAL(kind_real)             :: c                     ! scale height
    REAL(kind_real)             :: const                 ! refrac constant
    REAL(kind_real)             :: term1                 ! refrac term1
    REAL(kind_real)             :: term2                 ! refrac term2
    INTEGER                     :: lowest_level          ! lowest height level
 
    !------------------------------------------------------------
    ! Calculate the zenith delay for each layer and add to total
    !------------------------------------------------------------
 
    stationrefrac = 0.0
 
    DO level = 1, nlevq
      IF (zb(level) > zstation) THEN
        lowest_level = level
        EXIT
      END IF
    END DO
 
    ! Start at bottom level
 
    DO level = lowest_level, nlevq
 
      localzenithdelay = 0.0
 
      IF (level == lowest_level .AND. level /= 1) THEN
 
        ! If station lies above the lowest model level, interpolate refractivity
        ! to station height
 
        c = (log(refrac(level) / refrac(level - 1))) / (zb(level - 1) - zb(level))
        stationrefrac = refrac(level - 1) * exp(-c * (zstation - zb(level - 1)))
        const = -stationrefrac / c * exp(c * zstation)
        term1 = exp(-c * (zb(level)))
        term2 = exp(-c * zstation)
        localzenithdelay = refrac_scale * const * (term1 - term2)
 
      ELSE IF (level == 1) THEN
 
        ! If station lies below model level 1 (ie. the lowest level for which refractivity is
        ! calculated, then use c from the first full layer, but integrate down to height of
        ! station
 
        c = (log(refrac(level + 1) / refrac(level))) / (zb(level) - zb(level + 1))
        const = -refrac(level) / c * exp(c * (zb(level)))
        term1 = exp(-c * (zb(level + 1)))
        term2 = exp(-c * zstation)
        localzenithdelay = refrac_scale * const * (term1 - term2)
 
      ELSE IF (level <= nlevq .AND. level > 2) THEN
 
        ! If not at top level
 
        c = (log(refrac(level) / refrac(level - 1))) / (zb(level - 1) - zb(level))
        const = -refrac(level - 1) / c * exp(c * (zb(level - 1)))
        term1 = exp(-c * (zb(level)))
        term2 = exp(-c * (zb(level - 1)))
        localzenithdelay = refrac_scale * const * (term1 - term2)
 
      END IF
 
      model_ztd = model_ztd + localzenithdelay
  END DO
END SUBROUTINE ops_groundgnss_ztd
 
 
SUBROUTINE ops_groundgnss_topcorrection(P,    &
                                        nlevq, &
                                        za,    &
                                        zb,    &
                                        TopCorrection)
                    
    IMPLICIT NONE
    
    REAL(kind_real), INTENT(IN)      :: p(:)               ! Pressure on pressure (rho) levels
    INTEGER, INTENT(IN)              :: nlevq              ! no. of temperature/theta levels
    REAL(kind_real), INTENT(IN)      :: za(:)              ! heights of pressure (rho) levels
    REAL(kind_real), INTENT(IN)      :: zb(:)              ! Heights of temperature/theta levels
    REAL(kind_real), INTENT(INOUT)   :: topcorrection      ! ZTD Top of atmos correction
 
 
    INTEGER                          :: level              ! Loop counter
    REAL(kind_real)                  :: pn(nlevq)          ! Pressure on theta levels
    REAL(kind_real)                  :: pwt1               ! Weighting variable
    REAL(kind_real)                  :: pwt2               ! Weighting variable
    REAL(kind_real)                  :: tcconstant         ! Top correction constant
    REAL(kind_real), PARAMETER       :: hpa_to_pa = 100.0  ! hPa to Pascal conversion
    
    
    DO level = 1, nlevq
    
        pwt1 = (za(level + 1) - zb(level)) / (za(level + 1) - za(level))
        pwt2 = 1.0 - pwt1
 
        pn(level) = exp(pwt1 * log(p(level)) + pwt2 * log(p(level + 1)))
 
    END DO
 
    tcconstant = (refrac_scale * n_alpha * rd)/ (hpa_to_pa * grav)
    
    topcorrection = tcconstant * pn(nlevq)
    
END SUBROUTINE ops_groundgnss_topcorrection
 
END MODULE ufo_groundgnss_ukmo_utils_mod