! $Id: ropp_apps_utils.f90 4750 2016-03-22 10:08:19Z idculv $

MODULE ropp_apps_utils

!****m* Modules/ropp_apps_utils *
!
! NAME
!    ropp_apps_utils - Utility routines for ROPP applications module
!
! SYNOPSIS
!    USE ropp_apps_utils
!
! DESCRIPTION
!    This module provides signal and matrix processing routines used by
!    the ROPP pre-processor ionospheric correction package.
!
! AUTHOR
!   Met Office, Exeter, UK and ECMWF, Reading, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

CONTAINS

! ==============================================================================

!****s* APPSUtils/ropp_apps_calc_geop *
!
! NAME
!    ropp_apps_calc_geop - calculate geopotential for ECMWF-like bgr fields
!                          from {T, q, Ak, Bk}
!
! SYNOPSIS
!    CALL ropp_apps_calc_geop(ro_data, geop)
!
!**** 

  SUBROUTINE ropp_apps_calc_geop(ro_data, geop)
!
! Calculate geopotential for ECMWF-like bgr files from {T, q, Ak, Bk}.
! Slightly modified code lifted from relevant section of
! ropp_fm/common/ropp_fm_roprof2state.f90,
! to which the reader is referred for explanatory comments.

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_utils
    USE ropp_io_types, ONLY: ROprof
    USE ropp_apps_constants, ONLY: R_dry, g_wmo
!    USE ropp_apps_utils, not_this => ropp_apps_calc_geop

    IMPLICIT NONE

! In/out
    TYPE(ROprof), INTENT(IN)              :: ro_data ! RO profile containing lev2b data
    REAL(wp), DIMENSION(:), INTENT(INOUT) :: geop

! Local variables
    REAL(wp), DIMENSION(:), ALLOCATABLE   :: p_hlv, geop_hlv
    REAL(wp), DIMENSION(:), ALLOCATABLE   :: Tvflv, ln_prflv
    REAL(wp), DIMENSION(:), ALLOCATABLE   :: del_p, del_geop, alpha

    INTEGER                               :: n_hlv, n_flv, lvl

    CHARACTER(len=256)                    :: level_type
    CHARACTER(len=256)                    :: routine

! 1.1 Initialisation
! ------------------

    CALL message_get_routine(routine)

    CALL message_set_routine('ropp_apps_calc_geop')

! 1.2 Set up
! ----------

    level_type = ro_data%Lev2d%level_type
    IF ( INDEX(level_type,'UNKNOWN') > 0) level_type = ro_data%bg%source

    CALL To_Upper(level_type)

    IF ( (INDEX(level_type,'HYBRID') > 0) .OR. &
         (INDEX(level_type,'ECMWF')  > 0) ) THEN  ! ECMWF background

      n_hlv = SIZE(ro_data%lev2d%level_coeff_a)
      n_flv = SIZE(ro_data%lev2b%press)

      IF ( n_hlv /= n_flv+1 ) &
        CALL message(msg_fatal, "Inconsistent number of pressure levels " // &
                                "and Ak coefficients \n")

      ALLOCATE (p_hlv(n_hlv), geop_hlv(n_hlv))
      ALLOCATE (Tvflv(n_flv), ln_prflv(n_flv))
      ALLOCATE (del_p(n_flv), del_geop(n_flv), alpha(n_flv))

      p_hlv = ro_data%lev2d%level_coeff_a + &
              ro_data%lev2d%level_coeff_b * ro_data%lev2c%press_sfc

      WHERE (ro_data%lev2b%shum > ropp_MDTV) ! very occasionally they're not
        Tvflv  = (1.0_wp + 0.61e-3_wp * ro_data%lev2b%shum) * ro_data%lev2b%temp !NB: shum in g/kg
      ELSEWHERE
        Tvflv  = (1.0_wp + 0.61e-3_wp *             0.0_wp) * ro_data%lev2b%temp !NB: shum in g/kg
      END WHERE

      del_p(1:n_flv) = p_hlv(1:n_hlv-1) - p_hlv(2:n_hlv)

      ln_prflv(1:n_flv-1) = LOG(p_hlv(1:n_hlv-2)/p_hlv(2:n_hlv-1))
      ln_prflv(n_flv)     = 0.0_wp ! For completeness (never used)

      alpha(1:n_flv-1) = 1.0_wp - p_hlv(2:n_hlv-1)/del_p(1:n_flv-1) * ln_prflv(1:n_flv-1)
      alpha(n_flv)     = LOG(2.0_wp)

      del_geop = R_dry * Tvflv * ln_prflv / g_wmo

! 1.3 Integration
! ---------------

      geop_hlv(1) = ro_data%lev2c%geop_sfc
      DO lvl = 2, n_hlv-1
        geop_hlv(lvl) = geop_hlv(lvl-1) + del_geop(lvl-1)
      END DO

! 1.4 Interpolation
! -----------------

      geop = geop_hlv(1:n_hlv-1) + alpha * R_dry * Tvflv / g_wmo

! 1.5 Cleaning up
! ---------------

      DEALLOCATE (del_p, del_geop, alpha)
      DEALLOCATE (Tvflv, ln_prflv)
      DEALLOCATE (p_hlv, geop_hlv)

    ELSE

! 1.6 Computer says no
! --------------------

      CALL message(msg_info, "Background profile not in ECMWF format ... " // &
                             "cannot calculate geopotentials \n")
      geop = ropp_MDFV

    END IF

! 1.7 Clear up
! ------------

    CALL message_set_routine(routine)

  END SUBROUTINE ropp_apps_calc_geop


! ==============================================================================

!****s* APPSUtils/ropp_apps_calc_tdry *
!
! NAME
!    ropp_apps_calc_tdry - calculate dry temperature from refractivity profile
!
! SYNOPSIS
!    CALL ropp_apps_calc_tdry(ro_data, tdry)
!
!**** 

  SUBROUTINE ropp_apps_calc_tdry(ro_data, tdry)
!
! Calculate geopotential dry temperature from refractivity profile. 
! We integrate dp/dz = (-g(z)/k_1 R) N using a very crude climatological 
! estimate of dT/dz at the top, based on z(top).
! Euler forward, so 1st order accurate.  Probably good enough here - we just 
! need something with the same gradients as Tdry in the lower troposphere.

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_utils
    USE ropp_io, ONLY: ropp_io_ascend
    USE ropp_io_types, ONLY: ROprof
    USE ropp_apps_constants, ONLY: R_dry, g_wmo, kappa1

    IMPLICIT NONE

! In/out
    TYPE(ROprof), INTENT(IN)              :: ro_data ! RO profile containing lev2a data
    REAL(wp), DIMENSION(:), INTENT(INOUT) :: tdry

! Local variables
    REAL(wp), DIMENSION(:), ALLOCATABLE   :: z, ref, t, lnp
    REAL(wp)                              :: const ! g_wmo/R/kappa1
    REAL(wp)                              :: p_top, z_top, dlognbydz_top, dtbydz_top

! 'Climatological' lapse rate parameters (actually dT/dz = - usual lapse rate)
    REAL(wp), PARAMETER                   :: z_surf=0.0_wp            ! gpm
    REAL(wp), PARAMETER                   :: z_trop_lower=11.0E3_wp   ! gpm
    REAL(wp), PARAMETER                   :: z_trop_upper=13.0E3_wp   ! gpm
    REAL(wp), PARAMETER                   :: z_strat_lower=47.5E3_wp  ! gpm
    REAL(wp), PARAMETER                   :: z_strat_upper=52.5E3_wp  ! gpm
    REAL(wp), PARAMETER                   :: z_mes_lower=87.5E3_wp    ! gpm
    REAL(wp), PARAMETER                   :: z_mes_upper=92.5E3_wp    ! gpm

    REAL(wp), PARAMETER                   :: lrt_trop=-8.5E-3_wp      ! K/m
    REAL(wp), PARAMETER                   :: lrt_strat=1.7E-3_wp      ! K/m
    REAL(wp), PARAMETER                   :: lrt_mes=-2.5E-3_wp       ! K/m
    REAL(wp), PARAMETER                   :: lrt_therm=4.8E-3_wp      ! K/m

    INTEGER                               :: i, nlev

    CHARACTER(len=256)                    :: routine
    CHARACTER(len=10)                     :: s_z_top, s_dtbydz_top

    LOGICAL, DIMENSION(:), ALLOCATABLE    :: lvalid

! 1.1 Initialisation
! ------------------

    CALL message_get_routine(routine)

    CALL message_set_routine('ropp_apps_calc_tdry')

! 1.2 Set up
! ----------

    ALLOCATE ( lvalid(ro_data%lev2a%npoints) )

    lvalid = .FALSE.

    WHERE ( (ro_data%lev2a%refrac     > ropp_ZERO) .AND. &
            (ro_data%lev2a%alt_refrac > ropp_MDTV) ) lvalid = .TRUE.

    nlev = COUNT(lvalid)

    IF ( nlev < 3 ) THEN

      CALL message(msg_info, "Not enough valid refractivities and heights ... " // &
                             "cannot calculate dry temperatures \n")

      tdry = ropp_MDFV

      CALL message_set_routine(routine)

      RETURN

    END IF

! 1.3 Make sure heights are increasing
! ------------------------------------

    CALL ropp_io_ascend(ro_data)

! 1.4 Define constants and variables
! ----------------------------------

    ALLOCATE ( z(nlev), ref(nlev), t(nlev), lnp(nlev) )

    z = PACK(ro_data%lev2a%alt_refrac, lvalid)

    z = geometric2geopotential(ro_data%georef%lat, z)  ! Because we're assuming a constant g_wmo

    ref = PACK(ro_data%lev2a%refrac, lvalid)

    const = g_wmo / r_dry / kappa1

! 1.5 Boundary condition on p
! ---------------------------

! Estimate dlogN/dz(nlev-1).

    dlognbydz_top = LOG(ref(nlev)/ref(nlev-2)) / (z(nlev) - z(nlev-2))

    IF ( dlognbydz_top > -ropp_ZDTV ) THEN

      CALL message(msg_warn, "dlogN/dz non-negative at top of profile ... " // &
                             "cannot calculate dry temperature.\n")

      tdry = ropp_MDFV

      CALL message_set_routine(routine)

      RETURN

    END IF

! Set dT/dz(nlev-1).  Do not rearrange the order of these statements!

    z_top = z(nlev-1)

    IF (z_top > z_surf)        dtbydz_top = lrt_trop
    IF (z_top > z_trop_lower)  dtbydz_top = 0.0_wp
    IF (z_top > z_trop_upper)  dtbydz_top = lrt_strat
    IF (z_top > z_strat_lower) dtbydz_top = 0.0_wp
    IF (z_top > z_strat_upper) dtbydz_top = lrt_mes
    IF (z_top > z_mes_lower)   dtbydz_top = 0.0_wp
    IF (z_top > z_mes_upper)   dtbydz_top = lrt_therm

    p_top = -ref(nlev-1) * ((g_wmo/R_dry) + dtbydz_top) / &
                            kappa1 / dlognbydz_top

    WRITE (s_z_top,      FMT='(F10.5)')      z_top / 1.E3_wp
    WRITE (s_dtbydz_top, FMT='(F10.5)') dtbydz_top / 1.E-3_wp

    CALL message(msg_diag, "Setting dT/dz = " // s_dtbydz_top // " K/km," // &
                           " at z_top = " // s_z_top // " km.\n")

    lnp(nlev-1) = LOG(p_top)

! 1.6 Integrate downwards from z(nlev-1) to z(1)
! ----------------------------------------------

    DO i=nlev-1,2,-1

      lnp(i-1) = lnp(i) - &
                 0.5_wp * const * ref(i) * (z(i-1) - z(i)) * EXP(-lnp(i))  ! Estimate of lnp(i-1/2)

      lnp(i-1) = lnp(i) - &
                 const * SQRT(ref(i)*ref(i-1)) * (z(i-1) - z(i)) * EXP(-lnp(i-1))  ! Using previous estimate of lnp(i-1/2) on RHS

    END DO

! 1.7 Calculate lnp at top level
! ------------------------------

    lnp(nlev) = lnp(nlev-2) - const * ref(nlev-1) * (z(nlev) - z(nlev-2)) * EXP(-lnp(nlev-1))

! 1.8 Convert p to t
! ------------------

    t = kappa1 * EXP(lnp) / ref

! 1.9 Unpack t to return tdry
! ---------------------------

    tdry = UNPACK(t, lvalid, ropp_MDFV)

! 1.10 Clear up
! -------------

    DEALLOCATE (lnp, t, ref, z, lvalid)

    CALL message_set_routine(routine)

  END SUBROUTINE ropp_apps_calc_tdry

! ==============================================================================

!****s* APPSUtils/ropp_apps_impact2geom *
!
! NAME
!    ropp_apps_impact2geom - calculate geometric height from impact parameter
!
! SYNOPSIS
!    CALL ropp_apps_impact2geom(ro_data, geom, pblh_qc_flag)
!
!**** 

  SUBROUTINE ropp_apps_impact2geom(ro_data, geom, pblh_qc_flag)
!
!   Calculate geometric height h corresponding to a given impact parameter a
!   by solving (iteratively) a = r n(r), where 
!   n = 1 + refrac*10^-6 and
!   r = h + RoC + und
!   Refractivity on (lev1b) r values is found by interpolating (lev2a)
!   logs of the refracs.

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_utils
    USE ropp_io
    USE ropp_io_types, ONLY: ROprof
    USE ropp_apps_types
!    USE ropp_apps_utils, not_this => ropp_apps_impact2geom

    IMPLICIT NONE

! In/out
    TYPE(ROprof), INTENT(IN)              :: ro_data ! RO profile containing lev2b data
    REAL(wp), DIMENSION(:), INTENT(INOUT) :: geom
    INTEGER, INTENT(INOUT)                :: pblh_qc_flag

! Local parameters
    REAL(wp), PARAMETER                   :: refrac_surf_default=300.0_wp
    REAL(wp), PARAMETER                   :: tol=1.0e-3_wp ! m
    INTEGER, PARAMETER                    :: iter_max=250

! Local variables
    REAL(wp)                              :: roc, und
    REAL(wp)                              :: refrac_surf
    REAL(wp), DIMENSION(MAX(1, ro_data%lev1b%npoints)) &
                                          :: old_rad, new_rad, refrac

    INTEGER,  DIMENSION(:), POINTER       :: idx => NULL()   ! Array indices
    INTEGER                               :: nidx

    INTEGER                               :: iter
    CHARACTER(LEN=3)                      :: siter_max
    CHARACTER(len=256)                    :: routine

! 1.1 Initialisation
! ------------------

    CALL message_get_routine(routine)

    CALL message_set_routine('ropp_apps_impact2geom')

! 1.2 Definition
! --------------

    IF ( ro_data%lev1b%npoints < 1 ) THEN
      CALL message(msg_fatal, "No level 1b data in profile \n")
      CALL message_set_routine(routine)
      RETURN
    END IF

    IF ( ro_data%GEOref%roc < ropp_MDTV ) THEN
      CALL message(msg_error, "Missing radius of curvature ... " // &
                              "cannot calculate geometric height \n")
      pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
      CALL message_set_routine(routine)
      RETURN
    ELSE
      roc = ro_data%GEOref%roc
    END IF

    IF ( ro_data%GEOref%undulation < ropp_MDTV ) THEN
      CALL message(msg_warn, "Missing undulation ... assuming it to be zero \n")
      und = 0.0_wp
    ELSE
      und = ro_data%GEOref%undulation
    END IF

    idx => WHERE (ro_data%lev2a%refrac > ropp_ZERO, nidx)
    IF (nidx == 0) THEN
      CALL message(msg_error, "No positive refractivities ... " // &
                              "cannot calculate heights from impact parameters \n")
      pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
      CALL message_set_routine(routine)
      RETURN
    END IF

    IF ( (ro_data%lev2a%refrac(idx(1)) <   10.0_wp) .OR. &
         (ro_data%lev2a%refrac(idx(1)) > 1000.0_wp) ) THEN
      CALL message(msg_info, "First refractivity < 10 N-units or > 1000 N-units ... " // &
                             "using 300 N-units on first iteration \n")
      refrac_surf = refrac_surf_default
    ELSE
      refrac_surf = ro_data%lev2a%refrac(idx(1))
    END IF

! 1.3 Iteration
! -------------

    iter = 0

    old_rad = ro_data%lev1b%impact / (1.0_wp + 1.0e-6_wp*refrac_surf)

    DO

! Check monotonicity

      IF ( ANY( (old_rad(2:ro_data%lev1b%npoints  ) - &
                 old_rad(1:ro_data%lev1b%npoints-1)) < 0) ) THEN
        CALL message(msg_error, 'Calculation generated non-monotonic refractivity altitudes')
        pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
        EXIT
      END IF

      refrac = EXP( ropp_apps_interpol_1d(ro_data%lev2a%alt_refrac(idx),  &
                                          LOG(ro_data%lev2a%refrac(idx)), &
                                          old_rad - roc - und, .TRUE.) )

      new_rad = ro_data%lev1b%impact / (1.0_wp + 1.0e-6_wp*refrac)

      IF ( MAXVAL( ABS( old_rad - new_rad ) ) <= tol ) EXIT

      iter = iter + 1

      IF ( iter > iter_max ) THEN
        WRITE ( siter_max, '(i3)' ) iter_max
        CALL message(msg_error, 'Iteration failed to converge within ' // &
                                 siter_max // ' iterations')
        pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
        EXIT
      END IF

      old_rad = new_rad

    END DO

! 1.4 Calculation
! ---------------

  IF ( .NOT. BTEST(pblh_qc_flag, PBLH_QC_data_invalid) ) &
    geom = new_rad - roc - und

! 1.5 Clear up
! ------------

    CALL message_set_routine(routine)


  END SUBROUTINE ropp_apps_impact2geom

! ==============================================================================

!****s* PlanetaryBoundaryLayerHeight/ropp_apps_pblh_locate *
!
! NAME
!   ropp_apps_pblh_locate
!
! SYNOPSIS
!   Locate planetary boundary layer height diagnostic by interpolating dlapse/dz
!
!   CALL ropp_apps_pblh_locate(pblh_index, var_x, var_y, var_grad, pblh, pblx)
!
! DESCRIPTION
!   Hone PBLH estimate by find zero of lapse rate in vicinity of max/min.
!
! INPUTS
!  REAL(wp), DIMENSION(:), INTENT(IN)   :: var_x      ! Variable
!  REAL(wp), DIMENSION(:), INTENT(IN)   :: var_y      ! Height
!  REAL(wp), DIMENSION(:), INTENT(IN)   :: var_grad   ! Gradient
!  INTEGER, INTENT(IN)                  :: pblh_index ! Location of max/min
!
! OUTPUTS
!  REAL(wp), INTENT(OUT)                :: pblh         ! Height of PBL
!  REAL(wp), INTENT(OUT)                :: pblx         ! Variable at PBLH
!  INTEGER, INTENT(INOUT)               :: pblh_qc_flag ! QC flag
!
! AUTHOR
!   Met Office, Exeter, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

  SUBROUTINE ropp_apps_pblh_locate(pblh_index, var_x, var_y, var_grad, pblh, pblx, pblh_qc_flag)

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_utils
    USE ropp_io
    USE ropp_apps_types

    IMPLICIT NONE

! I/O
    REAL(wp), INTENT(OUT)                :: pblh ! Height of PBL
    REAL(wp), INTENT(OUT)                :: pblx ! Variable at PBLH

    REAL(wp), DIMENSION(:), INTENT(IN)   :: var_x
    REAL(wp), DIMENSION(:), INTENT(IN)   :: var_y
    REAL(wp), DIMENSION(:), INTENT(IN)   :: var_grad

    INTEGER, INTENT(IN)                  :: pblh_index
    INTEGER, INTENT(INOUT)               :: pblh_qc_flag

! Local
    REAL(wp)                             :: dgamma_by_dz_m,dgamma_by_dz_p ! Lapse rates
    REAL(wp)                             :: alpha, beta, dy_p, dy_m
    INTEGER                              :: n_points1
    CHARACTER(LEN=256)                   :: routine

! 1.1 Initialisation
! ------------------

    CALL message_get_routine(routine)

    CALL message_set_routine('ropp_apps_pblh_locate')

! 1.2 Estimate location of PBLH by interpolating dlapse/dz
! --------------------------------------------------------

    n_points1 = SIZE(var_x)

    IF (.NOT. BTEST(pblh_qc_flag, PBLH_QC_data_invalid)) THEN

      IF ( (pblh_index <= 2) .OR. (pblh_index >= (n_points1-1)) ) THEN

        CALL message(msg_diag, "Maximum gradient at top or bottom of profile \n")

        pblh = (var_y(MAX(pblh_index, 1)) + var_y(MIN(pblh_index+1, n_points1))) / 2.0_wp

        pblx = var_x(MAX(pblh_index, 1)) + &  ! using var_grad between pblh_index and pblh_index+1
               var_grad(MAX(pblh_index, 1)) * (pblh - var_y(MAX(pblh_index, 1)))

      ELSE

        CALL message(msg_diag, "Calculating PBLH by interpolating lapse rate gradient")

        dy_m = 0.5_wp * (var_y(pblh_index+1) - var_y(pblh_index-1)) ! at i=pblh_index
        dy_p = 0.5_wp * (var_y(pblh_index+2) - var_y(pblh_index  )) ! at i=pblh_index+1

        dgamma_by_dz_m = (var_grad(pblh_index  ) - var_grad(pblh_index-1)) / dy_m ! at i=pblh_index

        dgamma_by_dz_p = (var_grad(pblh_index+1) - var_grad(pblh_index  )) / dy_p ! at i=pblh_index+1

        IF ( (dgamma_by_dz_m*dgamma_by_dz_p > ropp_ZERO) .OR. &
             (ABS(dgamma_by_dz_m-dgamma_by_dz_p) < ropp_ZDTV) ) THEN

          CALL message(msg_info, "Could not find PBLH by interpolating ... " // &
                                 "leaving as missing data \n")
          pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)

        ELSE

! Fit gamma = gamma(i*) + alpha*(y-y(i*)) + beta*(y-y(i*))^2 through (i*-1, i*, i*+1) where i*=pblh_index

          alpha = (dgamma_by_dz_p*dy_m + dgamma_by_dz_m*dy_p) / (dy_m + dy_p)

          beta  = (dgamma_by_dz_p      - dgamma_by_dz_m     ) / (dy_m + dy_p)

          pblh = 0.5_wp * (var_y(pblh_index) + var_y(pblh_index+1)) - &
                 (0.5_wp*alpha/beta)

          pblx = 0.5_wp * (var_x(pblh_index) + var_x(pblh_index+1)) + &
                 (0.5_wp*alpha/beta)*(-var_grad(pblh_index) + (alpha**2/6.0_wp/beta))

        END IF

      END IF

    END IF

! 1.3 Clear up
! ------------

    CALL message_set_routine(routine)

  END SUBROUTINE ropp_apps_pblh_locate

! ==============================================================================

!****s* PlanetaryBoundaryLayerHeight/ropp_apps_pblh_check_min_height *
!
! NAME
!   ropp_apps_pblh_check_min_height
!
! SYNOPSIS
!   Check that the derived PBLH(s) is/are not less than pblh_min.
!
!   CALL ropp_apps_pblh_check_min_height(n_pblh, pblh_min, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)
!
! DESCRIPTION
!   Check that (up to) 2 PBLHs are >= pblh_min
!
! INPUTS
!  REAL(wp), INTENT(IN)                 :: pblh_min ! Minimum allowed PBLH
!  INTEGER, INTENT(IN)                  :: n_pblh ! No. of PBLHs
!
!  REAL(wp), INTENT(INOUT)              :: pblh1, pblh2 ! Heights of PBL
!  REAL(wp), INTENT(INOUT)              :: pblx1, pblx2 ! Variable at PBLH
!  INTEGER, INTENT(INOUT)               :: pblh_qc_flag  ! PBLH QC flag
!
! OUTPUTS
!  REAL(wp), INTENT(INOUT)              :: pblh1, pblh2 ! Heights of PBL
!  REAL(wp), INTENT(INOUT)              :: pblx1, pblx2 ! Variable at PBLH
!  INTEGER, INTENT(INOUT)               :: pblh_qc_flag  ! PBLH QC flag
!
! AUTHOR
!   Met Office, Exeter, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

  SUBROUTINE ropp_apps_pblh_check_min_height(n_pblh, pblh_min, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_utils
    USE ropp_io
    USE ropp_apps_types

    IMPLICIT NONE

! I/O
    REAL(wp), INTENT(IN)                 :: pblh_min ! Minimum allowed PBLH
    INTEGER, INTENT(IN)                  :: n_pblh ! No. of PBLHs

    REAL(wp), INTENT(INOUT)              :: pblh1, pblh2 ! Heights of PBL
    REAL(wp), INTENT(INOUT)              :: pblx1, pblx2 ! Variable at PBLH
    INTEGER, INTENT(INOUT)               :: pblh_qc_flag  ! PBLH QC flag

! Local
    INTEGER                              :: n_pblh_temp ! No. of PBLHs
    CHARACTER(LEN=7)                     :: str_pblh,str_pblh_min
    CHARACTER(LEN=256)                   :: routine

! 1.1 Initialisation
! ------------------

    CALL message_get_routine(routine)

    CALL message_set_routine('ropp_apps_pblh_check_min_height')

    CALL message(msg_diag, "Checking that the PBLH(s) is/are not too low")

! 1.2 Check PBLH(s) >= pblh_min
! -----------------------------

    n_pblh_temp = n_pblh

    IF (n_pblh > 1) THEN

      IF ((pblh1 < pblh_min) .AND. (pblh1 > ropp_MDTV)) THEN
        WRITE(str_pblh    , FMT='(F7.3)') 1.0e-3_wp*pblh1
        WRITE(str_pblh_min, FMT='(F7.3)') 1.0e-3_wp*pblh_min
        CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
                               " is below minimum of " // str_pblh_min // " km")
        pblh1 = ropp_MDFV
        pblx1 = ropp_MDFV
        n_pblh_temp = n_pblh_temp - 1
      END IF

      IF ((pblh2 < pblh_min) .AND. (pblh2 > ropp_MDTV)) THEN
        WRITE(str_pblh    , FMT='(F7.3)') 1.0e-3_wp*pblh2
        WRITE(str_pblh_min, FMT='(F7.3)') 1.0e-3_wp*pblh_min
        CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
                               " is below minimum of " // str_pblh_min // " km")
        pblh2 = ropp_MDFV
        pblx2 = ropp_MDFV
        n_pblh_temp = n_pblh_temp - 1
      END IF

      IF ((pblh1 < ropp_MDTV) .AND. (pblh2 < ropp_MDTV)) THEN
        pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_too_low)
      ELSE
        IF (pblh1 < ropp_MDTV) THEN ! swap them around
          pblh1 = pblh2  ;  pblh2 = ropp_MDFV
          pblx1 = pblx2  ;  pblx2 = ropp_MDFV
        END IF
      END IF

    ELSE  ! n_pblh = 1

      IF ((pblh1 < pblh_min) .AND. (pblh1 > ropp_MDTV)) THEN
        WRITE(str_pblh    , FMT='(F7.3)') 1.0e-3_wp*pblh1
        WRITE(str_pblh_min, FMT='(F7.3)') 1.0e-3_wp*pblh_min
        CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
                               " is below minimum of " // str_pblh_min // " km")
        pblh1 = ropp_MDFV
        pblx1 = ropp_MDFV
        pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_too_low)
        n_pblh_temp = n_pblh_temp - 1
      END IF

    END IF

! 1.3 Reset PBLH_QC_flag if necessary
! -----------------------------------

    SELECT CASE (n_pblh_temp)

      CASE(:1)
        pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_double_pblh)
        pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_multiple_pblh)

      CASE(2)
        pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_double_pblh)
        pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_multiple_pblh)

      CASE(3:)
        pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_double_pblh)
        pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_multiple_pblh)

    END SELECT

! 1.4 Clear up
! ------------

    CALL message_set_routine(routine)

  END SUBROUTINE ropp_apps_pblh_check_min_height

! ==============================================================================

!****s* PlanetaryBoundaryLayerHeight/ropp_apps_pblh_check_max_height *
!
! NAME
!   ropp_apps_pblh_check_max_height
!
! SYNOPSIS
!   Check that the derived PBLH(s) is/are not greater than pblh_max.
!
!   CALL ropp_apps_pblh_check_max_height(n_pblh, pblh_max, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)
!
! DESCRIPTION
!   Check that (up to) 2 PBLHs are <= pblh_max
!
! INPUTS
!  REAL(wp), INTENT(IN)                 :: pblh_max ! Maximum allowed PBLH
!  INTEGER, INTENT(IN)                  :: n_pblh ! No. of PBLHs
!
!  REAL(wp), INTENT(INOUT)              :: pblh1, pblh2 ! Heights of PBL
!  REAL(wp), INTENT(INOUT)              :: pblx1, pblx2 ! Variable at PBLH
!  INTEGER, INTENT(INOUT)               :: pblh_qc_flag  ! PBLH QC flag
!
! OUTPUTS
!  REAL(wp), INTENT(INOUT)              :: pblh1, pblh2 ! Heights of PBL
!  REAL(wp), INTENT(INOUT)              :: pblx1, pblx2 ! Variable at PBLH
!  INTEGER, INTENT(INOUT)               :: pblh_qc_flag  ! PBLH QC flag
!
! AUTHOR
!   Met Office, Exeter, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

  SUBROUTINE ropp_apps_pblh_check_max_height(n_pblh, pblh_max, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_utils
    USE ropp_io
    USE ropp_apps_types

    IMPLICIT NONE

! I/O
    REAL(wp), INTENT(IN)                 :: pblh_max ! Maximum allowed PBLH
    INTEGER, INTENT(IN)                  :: n_pblh ! No. of PBLHs

    REAL(wp), INTENT(INOUT)              :: pblh1, pblh2 ! Heights of PBL
    REAL(wp), INTENT(INOUT)              :: pblx1, pblx2 ! Variable at PBLH
    INTEGER, INTENT(INOUT)               :: pblh_qc_flag  ! PBLH QC flag

! Local
    INTEGER                              :: n_pblh_temp ! No. of PBLHs
    CHARACTER(LEN=7)                     :: str_pblh,str_pblh_max
    CHARACTER(LEN=256)                   :: routine

! 1.1 Initialisation
! ------------------

    CALL message_get_routine(routine)

    CALL message_set_routine('ropp_apps_pblh_check_max_height')

    CALL message(msg_diag, "Checking that the PBLH(s) is/are not too high")

! 1.2 Check PBLH(s) <= pblh_max
! -----------------------------

    n_pblh_temp = n_pblh

    IF (n_pblh > 1) THEN

      IF ((pblh1 > pblh_max) .AND. (pblh1 > ropp_MDTV)) THEN
        WRITE(str_pblh    , FMT='(F7.3)') 1.0e-3_wp*pblh1
        WRITE(str_pblh_max, FMT='(F7.3)') 1.0e-3_wp*pblh_max
        CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
                               " is above maximum of " // str_pblh_max // " km")
        pblh1 = ropp_MDFV
        pblx1 = ropp_MDFV
        n_pblh_temp = n_pblh_temp - 1
      END IF

      IF ((pblh2 > pblh_max) .AND. (pblh2 > ropp_MDTV)) THEN
        WRITE(str_pblh    , FMT='(F7.3)') 1.0e-3_wp*pblh2
        WRITE(str_pblh_max, FMT='(F7.3)') 1.0e-3_wp*pblh_max
        CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
                               " is above maximum of " // str_pblh_max // " km")
        pblh2 = ropp_MDFV
        pblx2 = ropp_MDFV
        n_pblh_temp = n_pblh_temp - 1
      END IF

      IF ((pblh1 < ropp_MDTV) .AND. (pblh2 < ropp_MDTV)) THEN
        pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_too_high)
      ELSE
        IF (pblh1 < ropp_MDTV) THEN ! swap them around
          pblh1 = pblh2  ;  pblh2 = ropp_MDFV
          pblx1 = pblx2  ;  pblx2 = ropp_MDFV
        END IF
      END IF

    ELSE  ! n_pblh = 1

      IF ((pblh1 > pblh_max) .AND. (pblh1 > ropp_MDTV)) THEN
        WRITE(str_pblh    , FMT='(F7.3)') 1.0e-3_wp*pblh1
        WRITE(str_pblh_max, FMT='(F7.3)') 1.0e-3_wp*pblh_max
        CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
                               " is above maximum of " // str_pblh_max // " km")
        pblh1 = ropp_MDFV
        pblx1 = ropp_MDFV
        pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_too_high)
        n_pblh_temp = n_pblh_temp - 1
      END IF

    END IF

! 1.3 Reset PBLH_QC_flag if necessary
! -----------------------------------

    SELECT CASE (n_pblh_temp)

      CASE(:1)
        pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_double_pblh)
        pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_multiple_pblh)

      CASE(2)
        pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_double_pblh)
        pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_multiple_pblh)

      CASE(3:)
        pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_double_pblh)
        pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_multiple_pblh)

    END SELECT

! 1.4 Clear up
! ------------

    CALL message_set_routine(routine)

  END SUBROUTINE ropp_apps_pblh_check_max_height

! ==============================================================================

!****f* APPSUtils/ropp_apps_interpol_1d *
!
! NAME
!    ropp_apps_interpol_1d - interpolate in 1D
!
! SYNOPSIS
!    y_out = ropp_apps_interpol_1d(x_in, y_in, x_out, flat)
!
!**** 
  FUNCTION ropp_apps_interpol_1d(x_in, y_in, x_out, flat) RESULT (y_out)

!   Linearly interpolates 1d arrays, with constant or linear extrapolation.
!   Assumes x_in and x_out are monotonically increasing.

    USE typesizes, ONLY: wp => EightByteReal
!    USE ropp_apps_utils, not_this => ropp_apps_interpol

    REAL(wp), DIMENSION(:), INTENT(IN)  :: x_in, y_in, x_out
    REAL(wp), DIMENSION(SIZE(x_out))    :: y_out
    LOGICAL, OPTIONAL                   :: flat ! If you want flat extrapolation

    REAL(wp)                            :: grad
    INTEGER                             :: i_in, i_out, n_in, n_out
    LOGICAL                             :: l_flat

    n_in  = SIZE(x_in)
    n_out = SIZE(x_out)

    l_flat = .FALSE.
    IF ( PRESENT(flat) ) l_flat = flat

    DO i_out=1,n_out

      grad = 0.0_wp

      IF ( x_out(i_out) <= x_in(1) ) THEN
        i_in = 1
        IF ( .NOT. l_flat ) grad = (y_in(i_in+1) - y_in(i_in)) / (x_in(i_in+1) - x_in(i_in))
      ELSE IF ( x_out(i_out) >= x_in(n_in) ) THEN
        i_in = n_in
        IF ( .NOT. l_flat ) grad = (y_in(i_in) - y_in(i_in-1)) / (x_in(i_in) - x_in(i_in-1))
      ELSE
        i_in = SUM ( MINLOC ( x_out(i_out)-x_in, MASK=(x_out(i_out)-x_in >= 0) ) )
        grad = (y_in(i_in+1) - y_in(i_in)) / (x_in(i_in+1) - x_in(i_in))
      END IF

      y_out(i_out) = y_in(i_in) + grad * (x_out(i_out) - x_in(i_in))  ! Linear (inter/extra)polation

    END DO

  END FUNCTION ropp_apps_interpol_1d


END MODULE ropp_apps_utils