| 1 | ! $Id: ropp_io_thin_fixed.f90 1959 2008-11-13 12:15:18Z frhl $
|
|---|
| 2 |
|
|---|
| 3 | !****s* Thin/ropp_io_thin_fixed *
|
|---|
| 4 | !
|
|---|
| 5 | ! NAME
|
|---|
| 6 | ! ropp_io_thin_fixed - Thin data by interpolating onto fixed altitude
|
|---|
| 7 | ! levels.
|
|---|
| 8 | !
|
|---|
| 9 | ! SYNOPSIS
|
|---|
| 10 | ! call ropp_io_thin_fixed ( nLev, FullLev, FullVal, &
|
|---|
| 11 | ! nThinLev, ThinLev, ThinVal, &
|
|---|
| 12 | ! Method, DEBUG, sigma )
|
|---|
| 13 | !
|
|---|
| 14 | !
|
|---|
| 15 | ! INPUTS
|
|---|
| 16 | ! nLev int Number of original levels
|
|---|
| 17 | ! FullLev dflt Original level values of the data
|
|---|
| 18 | ! FullVal dflt Original data values
|
|---|
| 19 | ! nThinLev int Number of thinned levels
|
|---|
| 20 | ! ThinLev dflt Thinned level values
|
|---|
| 21 | ! Method chr Required interpolation method. One of:
|
|---|
| 22 | ! SGLOG : Savitzky-Golay smoothing + log interp.
|
|---|
| 23 | ! SGLIN : Savitzky-Golay smoothing + lin interp.
|
|---|
| 24 | ! ASGLOG : Adaptive S-G smoothing + log interp.
|
|---|
| 25 | ! ASGLIN : Adaptive S-G smoothing + lin interp.
|
|---|
| 26 | ! LOG : logarithmic interpolation
|
|---|
| 27 | ! LIN : linear interpolation
|
|---|
| 28 | ! Default: LIN
|
|---|
| 29 | ! DEBUG log .T. to write additional diagnostics info to stdout
|
|---|
| 30 | ! sigma log .T. indicates error scaling. Applicable to ASGLIN
|
|---|
| 31 | ! and ASGLOG only (optional, default: .F.)
|
|---|
| 32 | !
|
|---|
| 33 | ! OUTPUTS
|
|---|
| 34 | ! ThinVal dlft Thinned data values on ThinLev levels
|
|---|
| 35 | !
|
|---|
| 36 | ! CALLS
|
|---|
| 37 | ! ropp_io_thin_sg
|
|---|
| 38 | !
|
|---|
| 39 | ! CALLED BY
|
|---|
| 40 | ! ropp_io_thin_select
|
|---|
| 41 | !
|
|---|
| 42 | ! USES
|
|---|
| 43 | ! typesizes
|
|---|
| 44 | ! ropp_io
|
|---|
| 45 | ! ropp_utils
|
|---|
| 46 | !
|
|---|
| 47 | ! DESCRIPTION
|
|---|
| 48 | ! This subroutine thins input data arrays to fixed number of levels.
|
|---|
| 49 | ! The levels are assumed to be altitudes with data values on those
|
|---|
| 50 | ! levels, but the algorithm is generic.
|
|---|
| 51 | ! Several interpolation methods are supported, including:
|
|---|
| 52 | ! a) Savitzky-Golay pre-smoothing & logarithmic interplation
|
|---|
| 53 | ! with optional adaptive filter parameters
|
|---|
| 54 | ! b) Logarithmic interpolation (no pre-smoothing)
|
|---|
| 55 | ! c) Linear interpolation (no pre-smoothing)
|
|---|
| 56 | ! (For LOG interpolations, an offset of the minimum data value is
|
|---|
| 57 | ! applied to the array of data values being interpolated (and removed
|
|---|
| 58 | ! after interpolation) to ensure no attempt to take the log of a negative
|
|---|
| 59 | ! value.)
|
|---|
| 60 | ! The input data values are interpolated to the given fixed (thinned)
|
|---|
| 61 | ! levels in the returned thinned array.
|
|---|
| 62 | ! Notes:
|
|---|
| 63 | ! 1) The pairs of input (full) levels/values arrays and output (thinned)
|
|---|
| 64 | ! levels/values arrays must be of the same order, rank and size.
|
|---|
| 65 | ! 2) The Adaptive SG takes account of the relative height 'resolutions'
|
|---|
| 66 | ! (maximum level separations) of the original and thinned levels, and
|
|---|
| 67 | ! modifies the width of the filter accordingly and the fitting polynomial
|
|---|
| 68 | ! order is set to 2 (1 for standard SG). Thus, thining of hi-res data
|
|---|
| 69 | ! applies more smoothing than for low-res data while still minimising
|
|---|
| 70 | ! vertical correlations or loss of information.
|
|---|
| 71 | ! 3) If the full (input) height range does not fill the range of thinned
|
|---|
| 72 | ! (fixed) heights, the out-of-range data values will be set missing.
|
|---|
| 73 | ! 4) If any input data values are flagged with 'missing' values
|
|---|
| 74 | ! (e.g. -99999) these will be included in the smoothing and/or
|
|---|
| 75 | ! interpolation which may not preserve any special missing data value for
|
|---|
| 76 | ! the thinned data values (but ought to remain invalid for that
|
|---|
| 77 | ! parameter in practice).
|
|---|
| 78 | ! 5) Because of the offsetting for LOG interpolation, data with
|
|---|
| 79 | ! large negative values (e.g. missing data) may reduce the accuracy
|
|---|
| 80 | ! of the interpolation at all levels.
|
|---|
| 81 | !
|
|---|
| 82 | ! TODO
|
|---|
| 83 | ! Possible future enhancements:
|
|---|
| 84 | ! 1) Replace (or add as option) the LOG interpolation post-SG
|
|---|
| 85 | ! with a spline interpolation (recommendation in Ref.)
|
|---|
| 86 | ! 2) Add (unsmoothed) spline interpolation as a new option
|
|---|
| 87 | ! 3) Other methods we haven't yet thought of.
|
|---|
| 88 | !
|
|---|
| 89 | ! REFERENCE
|
|---|
| 90 | ! Mono-dimensional data thinning for GPS radio occulation data
|
|---|
| 91 | ! SAF/GRAS/METO/ALG/ROPP/001
|
|---|
| 92 | !
|
|---|
| 93 | ! AUTHOR
|
|---|
| 94 | ! Met Office, Exeter, UK.
|
|---|
| 95 | ! Any comments on this software should be given via the GRAS SAF
|
|---|
| 96 | ! Helpdesk at http://www.grassaf.org
|
|---|
| 97 | !
|
|---|
| 98 | ! COPYRIGHT
|
|---|
| 99 | ! (c) EUMETSAT. All rights reserved.
|
|---|
| 100 | ! For further details please refer to the file COPYRIGHT
|
|---|
| 101 | ! which you should have received as part of this distribution.
|
|---|
| 102 | !
|
|---|
| 103 | !****
|
|---|
| 104 |
|
|---|
| 105 | subroutine ropp_io_thin_fixed ( nLev, FullLev, FullVal, &
|
|---|
| 106 | nThinLev, ThinLev, ThinVal, &
|
|---|
| 107 | Method, DEBUG, sigma )
|
|---|
| 108 |
|
|---|
| 109 | ! Modules
|
|---|
| 110 |
|
|---|
| 111 | USE typeSizes, ONLY: wp => EightByteReal
|
|---|
| 112 | USE ropp_io, ONLY: ropp_io_thin_sg
|
|---|
| 113 | USE ropp_io_types, ONLY: ropp_io_mdfv, &
|
|---|
| 114 | ropp_io_mdtv
|
|---|
| 115 | USE ropp_utils, ONLY: sort
|
|---|
| 116 |
|
|---|
| 117 |
|
|---|
| 118 | IMPLICIT NONE
|
|---|
| 119 |
|
|---|
| 120 | ! Fixed values
|
|---|
| 121 |
|
|---|
| 122 | INTEGER, PARAMETER :: nl = 1 ! SG: No. of points left of center
|
|---|
| 123 | INTEGER, PARAMETER :: nr = 1 ! SG: No. of points right of center
|
|---|
| 124 | INTEGER, PARAMETER :: m = 1 ! SG: Order of polynomial to be fitted
|
|---|
| 125 | INTEGER, PARAMETER :: lder = 0 ! SG: Order of derivative
|
|---|
| 126 |
|
|---|
| 127 | ! Argument list parameters
|
|---|
| 128 |
|
|---|
| 129 | CHARACTER (LEN=*), INTENT(IN) :: Method
|
|---|
| 130 | INTEGER, INTENT(IN) :: nLev
|
|---|
| 131 | INTEGER, INTENT(IN) :: nThinLev
|
|---|
| 132 | REAL(wp), INTENT(IN) :: FullLev(1:nLev)
|
|---|
| 133 | REAL(wp), INTENT(IN) :: FullVal(1:nLev)
|
|---|
| 134 | REAL(wp), INTENT(IN) :: ThinLev(1:nThinLev)
|
|---|
| 135 | REAL(wp), INTENT(OUT) :: ThinVal(1:nThinLev)
|
|---|
| 136 | LOGICAL, INTENT(IN) :: DEBUG
|
|---|
| 137 | LOGICAL, OPTIONAL, INTENT(IN) :: sigma
|
|---|
| 138 |
|
|---|
| 139 | ! Local variables
|
|---|
| 140 |
|
|---|
| 141 | REAL(wp) :: Val(1:nLev)
|
|---|
| 142 | REAL(wp) :: Lev(1:nLev)
|
|---|
| 143 | REAL(wp) :: SortVal(1:nLev)
|
|---|
| 144 | REAL(wp) :: SortLev(1:nLev)
|
|---|
| 145 | REAL(wp) :: ThinInt
|
|---|
| 146 | REAL(wp) :: minLev, maxLev, fullres, thinres
|
|---|
| 147 | REAL(wp) :: sigmafac, Val1, Val2
|
|---|
| 148 | INTEGER :: i, j, k, nsmooth
|
|---|
| 149 | LOGICAL :: ascending, lsigma, logint
|
|---|
| 150 | INTEGER :: idx(1:nLev)
|
|---|
| 151 |
|
|---|
| 152 | !----------------------------------------------------------------------------
|
|---|
| 153 | ! 0. Initialise
|
|---|
| 154 | !----------------------------------------------------------------------------
|
|---|
| 155 |
|
|---|
| 156 | ThinVal(:) = ropp_io_mdfv
|
|---|
| 157 | Val(:) = FullVal(:)
|
|---|
| 158 | Lev(:) = FullLev(:)
|
|---|
| 159 |
|
|---|
| 160 | IF ( .NOT. PRESENT(sigma) ) THEN
|
|---|
| 161 | lsigma = .FALSE.
|
|---|
| 162 | ELSE
|
|---|
| 163 | lsigma = sigma
|
|---|
| 164 | END IF
|
|---|
| 165 | sigmafac = 1.0_wp
|
|---|
| 166 |
|
|---|
| 167 | !----------------------------------------------------------------------------
|
|---|
| 168 | ! 1. Sovitzky-Golay pre-smoothing (optional)
|
|---|
| 169 | ! - with or without adaptive smoothing control
|
|---|
| 170 | !----------------------------------------------------------------------------
|
|---|
| 171 |
|
|---|
| 172 | IF ( INDEX ( Method, "SGLOG" ) == 1 .OR. &
|
|---|
| 173 | INDEX ( Method, "SGLIN" ) == 1 ) THEN
|
|---|
| 174 |
|
|---|
| 175 | CALL ropp_io_thin_sg ( nLev, Val )
|
|---|
| 176 |
|
|---|
| 177 | ELSE IF ( INDEX ( Method, "ASGLOG" ) == 1 .OR. &
|
|---|
| 178 | INDEX ( Method, "ASGLIN" ) == 1 ) THEN
|
|---|
| 179 |
|
|---|
| 180 | ! Get resolutions & smoothing distance (half-width)
|
|---|
| 181 |
|
|---|
| 182 | fullres = ABS ( MAXVAL ( Lev(2:nLev) - Lev(1:nLev-1) ) )
|
|---|
| 183 | thinres = ABS ( MINVAL ( ThinLev(2:nThinlev) - ThinLev(1:nThinlev-1) ) )
|
|---|
| 184 | nsmooth = MAX ( 2, FLOOR ( thinres / MAX(fullres,0.1_wp) ) )
|
|---|
| 185 | nsmooth = 2 * ( nsmooth / 2 )
|
|---|
| 186 |
|
|---|
| 187 | IF ( lsigma ) sigmafac = SQRT(nsmooth+1.0_wp)
|
|---|
| 188 |
|
|---|
| 189 | IF ( DEBUG ) THEN
|
|---|
| 190 | WRITE ( *, FMT="(A,I3,A)" ) " Adaptive SG smoothing over ", &
|
|---|
| 191 | nsmooth+1, " points"
|
|---|
| 192 | IF ( lsigma ) &
|
|---|
| 193 | WRITE ( *, FMT="(A,F5.2)" ) " Uncorrelated error reduction by factor ", &
|
|---|
| 194 | sigmafac
|
|---|
| 195 | END IF
|
|---|
| 196 |
|
|---|
| 197 | CALL ropp_io_thin_sg ( nLev, Val, npoints=INT(nsmooth/2), order=2 )
|
|---|
| 198 | END IF
|
|---|
| 199 |
|
|---|
| 200 |
|
|---|
| 201 | !----------------------------------------------------------------------------
|
|---|
| 202 | ! 2. Sort impact parameters. The input impact parameter data will be a
|
|---|
| 203 | ! monotonic function of height
|
|---|
| 204 | !----------------------------------------------------------------------------
|
|---|
| 205 |
|
|---|
| 206 | idx = sort(Lev)
|
|---|
| 207 | SortVal = Val(idx)
|
|---|
| 208 | SortLev = Lev(idx)
|
|---|
| 209 |
|
|---|
| 210 | Val = SortVal
|
|---|
| 211 | Lev = SortLev
|
|---|
| 212 |
|
|---|
| 213 | !----------------------------------------------------------------------------
|
|---|
| 214 | ! 3. Interpolation (Log or Lin). Allow for an ascending or descending
|
|---|
| 215 | ! input profile. The output (thinned) profile will be in the order
|
|---|
| 216 | ! of the given thinned levels. Do not extrapolate!
|
|---|
| 217 | !----------------------------------------------------------------------------
|
|---|
| 218 |
|
|---|
| 219 | IF ( Lev(1) < Lev(nLev) ) THEN
|
|---|
| 220 | ascending = .TRUE.
|
|---|
| 221 | ELSE
|
|---|
| 222 | ascending = .FALSE.
|
|---|
| 223 | END IF
|
|---|
| 224 |
|
|---|
| 225 | minLev = MINVAL(Lev)
|
|---|
| 226 | maxLev = MAXVAL(Lev)
|
|---|
| 227 |
|
|---|
| 228 | !----------------------------------------------------------------------------
|
|---|
| 229 | ! 3.1 Linear interpolation
|
|---|
| 230 | !----------------------------------------------------------------------------
|
|---|
| 231 |
|
|---|
| 232 | IF ( INDEX ( Method, "ASGLIN" ) == 1 .OR. &
|
|---|
| 233 | INDEX ( Method, "SGLIN" ) == 1 .OR. &
|
|---|
| 234 | INDEX ( Method, "LIN" ) == 1 ) THEN
|
|---|
| 235 |
|
|---|
| 236 | DO k = 1, nThinLev
|
|---|
| 237 |
|
|---|
| 238 | ! Skip upper/lower parts of target profile not covered by source profile
|
|---|
| 239 |
|
|---|
| 240 | IF ( ThinLev(k) < minLev .OR. &
|
|---|
| 241 | ThinLev(k) > maxLev) CYCLE
|
|---|
| 242 |
|
|---|
| 243 | ! Search for the pair of source levels just above & below the
|
|---|
| 244 | ! target level, allowing for ascending or descending profile.
|
|---|
| 245 |
|
|---|
| 246 | j = 2
|
|---|
| 247 | IF ( ascending ) THEN
|
|---|
| 248 | DO WHILE ( j < nLev .AND. Lev(j) < ThinLev(k) )
|
|---|
| 249 | j = j + 1
|
|---|
| 250 | END DO
|
|---|
| 251 | i = j - 1
|
|---|
| 252 | ELSE
|
|---|
| 253 | DO WHILE ( j < nLev-1 .AND. Lev(j) > ThinLev(k) )
|
|---|
| 254 | j = j + 1
|
|---|
| 255 | END DO
|
|---|
| 256 | i = j
|
|---|
| 257 | j = j - 1
|
|---|
| 258 | END IF
|
|---|
| 259 |
|
|---|
| 260 | ! If either source value is missing, either level value is missing, both source levels
|
|---|
| 261 | ! are the same, or the gap between levels is larger than interval between thinned
|
|---|
| 262 | ! levels, skip this target level.
|
|---|
| 263 |
|
|---|
| 264 | IF ( k == nThinLev ) THEN
|
|---|
| 265 | ThinInt = ABS ( ThinLev(k) - ThinLev(k-1) )
|
|---|
| 266 | ELSE
|
|---|
| 267 | ThinInt = ABS ( ThinLev(k+1) - ThinLev(k) )
|
|---|
| 268 | ENDIF
|
|---|
| 269 |
|
|---|
| 270 | IF ( Val(i) < ropp_io_mdtv .OR. Val(j) < ropp_io_mdtv .OR. &
|
|---|
| 271 | Lev(i) < ropp_io_mdtv .OR. Lev(j) < ropp_io_mdtv .OR. &
|
|---|
| 272 | Lev(i) == Lev(j) .OR. &
|
|---|
| 273 | ABS(Lev(i)-Lev(j)) > ThinInt ) CYCLE
|
|---|
| 274 |
|
|---|
| 275 | ! Simple linear interpolation
|
|---|
| 276 |
|
|---|
| 277 | ThinVal(k) = ( Val(j) + ( (ThinLev(k) - Lev(j)) &
|
|---|
| 278 | / (Lev(i) - Lev(j)) * (Val(i) - Val(j)) ) ) &
|
|---|
| 279 | * sigmafac
|
|---|
| 280 |
|
|---|
| 281 | END DO
|
|---|
| 282 |
|
|---|
| 283 | !----------------------------------------------------------------------------
|
|---|
| 284 | ! 3.2 Logarithic interpolation.
|
|---|
| 285 | !----------------------------------------------------------------------------
|
|---|
| 286 |
|
|---|
| 287 | ELSE IF ( INDEX ( Method, "ASGLOG" ) == 1 .OR. &
|
|---|
| 288 | INDEX ( Method, "SGLOG" ) == 1 .OR. &
|
|---|
| 289 | INDEX ( Method, "LOG" ) == 1 ) THEN
|
|---|
| 290 |
|
|---|
| 291 | DO k = 1, nThinLev
|
|---|
| 292 |
|
|---|
| 293 | ! Skip upper/lower parts of target profile not covered by source profile
|
|---|
| 294 |
|
|---|
| 295 | IF ( ThinLev(k) < minLev .OR. &
|
|---|
| 296 | ThinLev(k) > maxLev) CYCLE
|
|---|
| 297 |
|
|---|
| 298 | ! Search for the pair of source levels just above & below the
|
|---|
| 299 | ! target level, allowing for ascending or descending profile.
|
|---|
| 300 |
|
|---|
| 301 | j = 2
|
|---|
| 302 | IF ( ascending ) THEN
|
|---|
| 303 | DO WHILE ( j < nLev .AND. Lev(j) < ThinLev(k) )
|
|---|
| 304 | j = j + 1
|
|---|
| 305 | END DO
|
|---|
| 306 | i = j - 1
|
|---|
| 307 | ELSE
|
|---|
| 308 | DO WHILE ( j < nLev-1 .AND. Lev(j) > ThinLev(k) )
|
|---|
| 309 | j = j + 1
|
|---|
| 310 | END DO
|
|---|
| 311 | i = j
|
|---|
| 312 | j = j - 1
|
|---|
| 313 | END IF
|
|---|
| 314 |
|
|---|
| 315 | ! If either source value is missing, either level value is missing, both source levels
|
|---|
| 316 | ! are the same, or the gap between levels is larger than interval between thinned
|
|---|
| 317 | ! levels, skip this target level.
|
|---|
| 318 |
|
|---|
| 319 | IF ( k == nThinLev ) THEN
|
|---|
| 320 | ThinInt = ABS ( ThinLev(k) - ThinLev(k-1) )
|
|---|
| 321 | ELSE
|
|---|
| 322 | ThinInt = ABS ( ThinLev(k+1) - ThinLev(k) )
|
|---|
| 323 | ENDIF
|
|---|
| 324 |
|
|---|
| 325 | IF ( Val(i) < ropp_io_mdtv .OR. Val(j) < ropp_io_mdtv .OR. &
|
|---|
| 326 | Lev(i) < ropp_io_mdtv .OR. Lev(j) < ropp_io_mdtv .OR. &
|
|---|
| 327 | Lev(i) == Lev(j) .OR. &
|
|---|
| 328 | ABS(Lev(i)-Lev(j)) > ThinInt ) CYCLE
|
|---|
| 329 |
|
|---|
| 330 | ! If both source values are larger than zero do interpolation in log
|
|---|
| 331 | ! space (and then invert target value), else fall back on linear
|
|---|
| 332 | ! interpolation
|
|---|
| 333 |
|
|---|
| 334 | IF ( Val(i) > 0.0_wp .AND. &
|
|---|
| 335 | Val(j) > 0.0_wp ) THEN
|
|---|
| 336 | Val1 = LOG(Val(i))
|
|---|
| 337 | Val2 = LOG(Val(j))
|
|---|
| 338 | logint = .TRUE.
|
|---|
| 339 | ELSE
|
|---|
| 340 | Val1 = Val(i)
|
|---|
| 341 | Val2 = Val(j)
|
|---|
| 342 | logint = .FALSE.
|
|---|
| 343 | END IF
|
|---|
| 344 |
|
|---|
| 345 | ThinVal(k) = ( Val2 + ( (ThinLev(k) - Lev(j)) &
|
|---|
| 346 | / (Lev(i) - Lev(j)) * (Val1 - Val2) ) )
|
|---|
| 347 |
|
|---|
| 348 | IF ( logint ) ThinVal(k) = EXP(ThinVal(k))
|
|---|
| 349 |
|
|---|
| 350 | ThinVal(k) = ThinVal(k) * sigmafac
|
|---|
| 351 | END DO
|
|---|
| 352 |
|
|---|
| 353 | !----------------------------------------------------------------------------
|
|---|
| 354 | ! 3.3 Unknown option
|
|---|
| 355 | !----------------------------------------------------------------------------
|
|---|
| 356 |
|
|---|
| 357 | ELSE
|
|---|
| 358 |
|
|---|
| 359 | IF ( DEBUG ) THEN
|
|---|
| 360 | WRITE ( *, FMT="(A/A/A)" ) "*** Unknown thinning option", &
|
|---|
| 361 | " The available options for the fixed level " // &
|
|---|
| 362 | "thinning are:", &
|
|---|
| 363 | " ASGLOG, ASGLIN, SGLOG, SGLIN, LIN and LOG"
|
|---|
| 364 | END IF
|
|---|
| 365 |
|
|---|
| 366 |
|
|---|
| 367 | END IF
|
|---|
| 368 |
|
|---|
| 369 | END subroutine ropp_io_thin_fixed
|
|---|
