ufo_refractivity_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_refractivity_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
    cp,                     &    ! Heat capacity at constant pressure for air
    rd_over_cp,             &    ! Ratio of gas constant to heat capacity
    pref,                   &    ! Reference pressure for calculating exner
    grav,                   &    ! Gravitational field strength
    n_alpha,                &    ! Refractivity constant a
    n_beta,                  &    ! Refractivity constant b
    mw_ratio,               &    ! Ratio of molecular weights of water and dry air
    c_virtual                    ! Related to mw ratio
 
implicit none
 
public :: ufo_refractivity
private
 
contains
 
 
SUBROUTINE ufo_refractivity(nlevq,     &
                            nlevP,     &
                            za,        &
                            zb,        &
                            x,         &
                            pN,        &
                            refracerr, &
                            refrac)
                
IMPLICIT NONE
 
! Subroutine arguments:
 
INTEGER, INTENT(IN)                    :: nlevq        ! no. of levels of wet refractivity required
INTEGER, INTENT(IN)                    :: nlevp        ! no. of levels of dry refractivity required
REAL(kind_real), INTENT(IN)            :: za(:)    ! heights of pressure levels
 
REAL(kind_real), INTENT(IN)            :: zb(:)             ! Heights of theta levels
REAL(kind_real), INTENT(IN)            :: x(:) ! state vector
REAL(kind_real), INTENT(INOUT)         :: pn(:)        ! pressure on refractivity (theta) levels
LOGICAL, INTENT(OUT)                   :: refracerr    ! errors in refractivity calculation
REAL(kind_real), INTENT(INOUT)         :: refrac(:)    ! refrac on refractivity (theta) levels
 
 
! Local declarations:
CHARACTER(len=*), PARAMETER :: routinename = "ufo_refractivity"
INTEGER                     :: i
INTEGER                     :: level
REAL, ALLOCATABLE           :: exnern(:)      ! Exner on refractivity (theta) levels
REAL, ALLOCATABLE           :: exner(:)  ! Exner on pressure (rho) levels
REAL(kind_real)             :: t(nlevq)       ! Temp. on theta levs
REAL                        :: tv          ! Tv on refractivity (theta) level
REAL                        :: ndry        ! Dry refractivity
REAL                        :: nwet        ! Wet refractivity
 
INTEGER                     :: numplevs    ! Number of levels of pInter
INTEGER                     :: numqlevs    ! Numbers of levels of qN
LOGICAL                     :: nonmon      ! non-monotonic pressure warning
LOGICAL                     :: unphys      ! zero or negative pressure warning
LOGICAL                     :: levelerr    ! nlevq greater than nlevp error
integer, parameter          :: max_string = 800
CHARACTER(len=max_string)   :: message
 
REAL(kind_real)             :: p(nlevp)
REAL(kind_real)             :: q(nlevq)
REAL(kind_real)             :: pwt1
REAL(kind_real)             :: pwt2
INTEGER                     :: nstate
 
 
! Get constants into right units
 
nstate = nlevp + nlevq
p(:) = x(1:nlevp)
q(:) = x(nlevp + 1:nstate)
 
ALLOCATE (exnern(nlevq))
ALLOCATE (exner(nlevp))
 
refrac(:) = missing_value(refrac(1))
t(:) = missing_value(t(1))
nonmon = .false.
unphys = .false.
levelerr = .false.
refracerr = .false.
 
DO i = 1, nlevp
  IF (p(i) == missing_value(p(i))) THEN  !pressure missing
    refracerr = .true.
    WRITE(message, *) routinename, "Missing value P", i
    CALL fckit_log % warning(message)
    EXIT
  END IF
END DO
 
 
DO i = 1, nlevp-1
  IF (p(i) - p(i + 1) < 0.0) THEN       !or non-monotonic pressure
    refracerr = .true.
    nonmon = .true.
    WRITE(message,*) "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.
  unphys = .true.
END IF
 
IF (nlevq >= nlevp) THEN  ! nlevq must be < than nlevP
  refracerr = .true.
  levelerr = .true.
END IF
 
! only proceed if pressure is valid
IF (refracerr) THEN
  IF (nonmon) THEN
    CALL fckit_log%warning (routinename // ": Pressure non-monotonic")
  ELSE IF (unphys) THEN
    CALL fckit_log%warning (routinename // ": Pressure <= zero")
  ELSE
    CALL fckit_log%warning (routinename // ": Pressure missing")
  END IF
  IF (levelerr) THEN
    CALL fckit_log%warning (routinename // ": Too many wet levels")
  END IF
 
ELSE
 
  ! Calculate exner on rho(pressure) levels.
 
  exner(:) = (p(:) / pref) ** rd_over_cp
 
  ! Calculate the refractivity on the b levels
 
 
  DO level = 1, nlevp - 1
    
    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
 
  DO i = 1, nlevq
 
    ! Calculate Exner on the refractivity level.
    exnern(i) = (pn(i) / pref) ** rd_over_cp
 
    ! Calculate mean layer Tv using ND definition
 
    tv = grav * (za(i + 1) - za(i)) * exnern(i) / &
        (cp * (exner(i) - exner(i + 1)))
 
    IF (i > nlevq) THEN
 
      t(i) = tv
 
      ! No wet component
 
      nwet = 0.0
 
    ELSE
 
      t(i) = tv / (1.0 + c_virtual * q(i))
 
      ! Wet component
 
      nwet = n_beta * pn(i) * q(i) / (t(i) ** 2 * (mw_ratio + (1.0 - mw_ratio) * q(i)))
 
 
    END IF
 
    IF (i > nlevp ) THEN
      ! No dry component
 
      ndry = 0.0
 
    ELSE
      ! Dry component
 
      ndry = n_alpha * pn(i) / t(i)
 
    END IF
 
    refrac(i) = ndry + nwet
 
  END DO
 
END IF
 
END SUBROUTINE ufo_refractivity
 
END MODULE ufo_refractivity_utils_mod