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 theta levels
    REAL(kind_real), INTENT(IN)       :: refrac(:)
    REAL(kind_real), INTENT(IN)       :: zb(:)
    REAL(kind_real), INTENT(IN)       :: zstation
    REAL(kind_real), INTENT(INOUT)    :: model_ztd
    
    
    REAL(kind_real)             :: localzenithdelay
    INTEGER                     :: level
    REAL(kind_real)             :: stationrefrac
    REAL(kind_real)             :: c
    REAL(kind_real)             :: const
    REAL(kind_real)             :: term1
    REAL(kind_real)             :: term2
    INTEGER                     :: lowest_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(pN,    &
                                        nlevq, &
                                        TopCorrection)
                    
    IMPLICIT NONE
    
    REAL(kind_real), INTENT(IN)      :: pn(:)
    INTEGER, INTENT(IN)              :: nlevq
    REAL(kind_real), INTENT(INOUT)   :: topcorrection
 
    REAL(kind_real)                  :: tcconstant
    REAL(kind_real), PARAMETER       :: hpa_to_pa = 100.0
 
    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