var/obsproc/src/module_recoverh.F90

   1 MODULE module_recoverh
   2
   3 !------------------------------------------------------------------------------!
   4 ! Recover observation height, when missing, based on pressure.
   5 !
   6 ! Y.-R. GUO, September 2000
   7 !------------------------------------------------------------------------------!
   8
   9    USE module_type
  10    USE module_func
  11    USE module_per_type
  12    USE module_mm5
  13
  14    INCLUDE 'missing.inc'
  15
  16 CONTAINS
  17
  18 !------------------------------------------------------------------------------
  19 SUBROUTINE recover_height_from_pressure(max_number_of_obs , &
  20                              obs , number_of_obs, print_hp_recover)
  21
  22 !  This routine recovers the missing heights based on the pressure
  23 !  under the hydrostatic assumption, or the model reference state.
  24 !  for the multi-level OBS data (SOUND, AIREP, etc.).
  25 !
  26
  27    IMPLICIT NONE
  28
  29    INTEGER,       INTENT ( IN )                :: max_number_of_obs
  30    TYPE (report), DIMENSION (max_number_of_obs):: obs
  31    INTEGER , INTENT ( IN )                     :: number_of_obs
  32    LOGICAL , INTENT ( IN )                     :: print_hp_recover
  33
  34    TYPE (measurement), POINTER                 :: current
  35    INTEGER                                     :: iunit
  36    INTEGER                                     :: qc_flag
  37    INTEGER                                     :: i, j, nlevel, k, &
  38                                                   k_start, k_top
  39    CHARACTER (LEN = 80)                        :: filename
  40    CHARACTER (LEN = 80)                        :: proc_name = &
  41                                                  "recover_height_from_pressure"
  42    LOGICAL                                     :: connected, correct, failed
  43    INTEGER                                     :: io_error
  44
  45
  46    TYPE (field)  , dimension(9000)              :: hh
  47    REAL          , dimension(9000)              :: pp, tt, qq
  48
  49    INCLUDE 'platform_interface.inc'
  50
  51 !------------------------------------------------------------------------------!
  52              WRITE (0,'(A)')  &
  53 '------------------------------------------------------------------------------'
  54       WRITE (UNIT = 0, FMT = '(A,/)') 'HEIGHT RECOVERED FROM P, T, Q,..:'
  55 !
  56
  57       !  Open diagnostic file
  58
  59       IF (print_hp_recover) THEN
  60
  61       filename = 'obs_recover_height.diag'
  62       iunit    = 999
  63
  64       INQUIRE (UNIT = iunit, OPENED = connected )
  65
  66       IF (connected) CLOSE (iunit)
  67
  68       OPEN (UNIT = iunit , FILE = filename , FORM = 'FORMATTED'  , &
  69             ACTION = 'WRITE' , STATUS = 'REPLACE', IOSTAT = io_error )
  70
  71       IF (io_error .NE. 0) THEN
  72           CALL error_handler (proc_name, &
  73          "Unable to open output diagnostic file. " , filename, .TRUE.)
  74       ELSE
  75           WRITE (UNIT = 0, FMT = '(A,A,/)') &
  76          "Diagnostics in file ", TRIM (filename)
  77       ENDIF
  78
  79       ENDIF
  80
  81       IF (print_hp_recover) &
  82       WRITE (UNIT = IUNIT, FMT = '(/A/)') &
  83         'HEIGHT RECOVERED FROM PRESSURE FOR MULTI-LEVEL OBS DATA:'
  84
  85       failed = .false.
  86
  87 ! 1.  ESTIMATE H
  88 ! ==============
  89
  90       j = 0
  91
  92 ! 1.1 Loop over obs
  93 !     -------------
  94
  95 loop_all: &
  96       DO i = 1, number_of_obs
  97
  98          IF ((obs (i) % info % discard)  .OR. .NOT. ASSOCIATED &
  99              (obs (i) % surface)) THEN
 100
 101              CYCLE loop_all
 102
 103          ENDIF
 104
 105
 106 ! 2.  SINGLE LEVEL OBS
 107 ! ====================
 108
 109 surface:&
 110          IF ((ASSOCIATED (obs (i) % surface)) .AND. &
 111         (.NOT.ASSOCIATED (obs (i) % surface % next))) THEN
 112
 113              ! IF height is missing, pressure should be present
 114
 115              IF (eps_equal (&
 116                  obs (i) % surface % meas % height % data, missing_r, 1.)) THEN
 117
 118                  obs (i) % surface % meas % height   % data = ref_height &
 119                 (obs (i) % surface % meas % pressure % data)
 120                  obs (i) % surface % meas % height   % qc   = Reference_atmosphere
 121
 122                  obs (i) % surface % meas % height % data = NINT &
 123                 (obs (i) % surface % meas % height % data + .5)
 124                  obs (i) % surface % meas % height % data = MAX &
 125                 (obs (i) % surface % meas % height % data, 0.)
 126
 127
 128                  IF (print_hp_recover) THEN
 129
 130                      WRITE (UNIT = iunit,FMT = '(/,A,A5,1X,A23,2F9.3)')        &
 131                     "Recover 1 level  station id = ",                          &
 132                      TRIM  (obs (i) % location % id ) ,                        &
 133                      TRIM  (obs (i) % location % name),                        &
 134                             obs (i) % location % latitude,                     &
 135                             obs (i) % location % longitude
 136                      WRITE (UNIT = iunit, FMT = '(2(A,I5),A)')                 &
 137                     "Use reference state to infer height (",                   &
 138                      INT (obs (i) % surface % meas % height % data),           &
 139                      "m) from pressure (",&
 140                      INT (obs (i) % surface % meas % pressure % data/100.),"hPa)."
 141                  ENDIF
 142
 143              ENDIF
 144
 145              !  Recover station elevation for single level obs
 146
 147 !            IF (eps_equal (&
 148 !                obs (i) % info % elevation, missing_r, 1.)) THEN
 149 !                obs (i) % info % elevation =  &
 150 !                obs (i) % surface % meas % height % data
 151 !            ENDIF
 152
 153              CYCLE loop_all
 154
 155          ENDIF surface
 156
 157 ! 3. MULTI LEVEL OBS
 158 ! ==================
 159 ! Y.-R. Guo (10/25/2005):.......
 160          call reorder(obs(i), i, 'pressure', failed)
 161          if (failed) then
 162             obs(i) % info % discard = .true.
 163             cycle loop_all
 164          endif
 165
 166 ! ..........................................
 167 ! 3.1 Get the OBS profile and count the number of levels
 168 !     --------------------------------------------------
 169
 170          nlevel  = 0
 171          correct = .FALSE.
 172          current => obs(i)%surface
 173
 174 count_level_1:&
 175          DO WHILE (ASSOCIATED (current))
 176
 177             nlevel = nlevel + 1
 178
 179             hh (nlevel) = current%meas%height
 180             pp (nlevel) = current%meas%pressure%data
 181             tt (nlevel) = current%meas%temperature%data
 182             qq (nlevel) = current%meas%qv%data
 183
 184             IF (eps_equal(current%meas%height%data, missing_r, 1.)) THEN
 185                 correct = .TRUE.
 186             ENDIF
 187
 188             current => current%next
 189
 190          ENDDO count_level_1
 191
 192 ! 3.2 If all levels have height, go to next station
 193 !     ---------------------------------------------
 194
 195          IF (.not.correct) CYCLE loop_all
 196
 197
 198 ! 3.2 Otherwise recover missing height for upper-level
 199 !     ------------------------------------------------
 200
 201 levels:&
 202          IF (nlevel <= 1) THEN
 203
 204              IF (print_hp_recover) THEN
 205                  WRITE (UNIT = iunit , FMT = '(A,A5,1X,A23,2F9.3)')     &
 206                 "No level found for sound id= " ,                       &
 207                  TRIM  (obs (i) % location % id ) ,                     &
 208                  TRIM  (obs (i) % location % name),                     &
 209                         obs (i) % location % latitude,                  &
 210                         obs (i) % location % longitude
 211              ENDIF
 212
 213              STOP 'in recover_height.F90'
 214
 215          ELSE IF (nlevel > 1) THEN levels
 216
 217              CALL recover_h_from_ptq (pp, tt, qq, hh, nlevel,k_start,k_top)
 218
 219              IF (k_start >= 1 .or. k_top <= nlevel) THEN
 220
 221                  IF (print_hp_recover) &
 222                      WRITE (UNIT = iunit, FMT = '(/,A,A5,1X,A23,2F9.3)') &
 223                     "Recover upperair station id = ",                    &
 224                      TRIM  (obs (i) % location % id ) ,                  &
 225                      TRIM  (obs (i) % location % name),                  &
 226                             obs (i) % location % latitude,               &
 227                             obs (i) % location % longitude
 228
 229                  ENDIF
 230
 231          ENDIF levels
 232
 233 ! 3.  CORRECT OBS
 234 ! ===============
 235
 236          k = 0
 237          current => obs(i)%surface
 238
 239 correct_levels: &
 240          DO WHILE (ASSOCIATED (current))
 241
 242             k = k + 1
 243
 244 !            IF (eps_equal(current%meas%height%data, missing_r, 1.)) THEN
 245
 246                  IF (print_hp_recover) THEN
 247                     WRITE (UNIT = iunit, FMT = '(2(A,I5),A)')    &
 248                    "Height missing set to ",INT (hh (k) % data), &
 249                    "m, pressure = ",INT (pp(k)/100.),"hpa."
 250                  ENDIF
 251
 252                  current%meas%height % data = CEILING (hh(k)%data)
 253                  current%meas%height % qc   =          hh(k)%qc
 254
 255 !            ENDIF
 256
 257             current => current%next
 258
 259
 260          ENDDO correct_levels
 261
 262          call reorder(obs(i), i, 'pressure', failed)
 263          if (failed) then
 264             obs(i) % info % discard = .true.
 265             cycle loop_all
 266          endif
 267
 268 ! 3.4 Go to next station
 269 !     ------------------
 270
 271       ENDDO loop_all
 272
 273      IF (print_hp_recover) CLOSE (IUNIT)
 274
 275 END SUBROUTINE recover_height_from_pressure
 276
 277 !-----------------------------------------------------------------------
 278
 279  SUBROUTINE recover_h_from_ptq(P,T,Q,HGT,KXS,K_START,K_TOP)
 280 !-----------------------------------------------------------------------
 281 ! To recover the missing heights for the MULTI-LEVEL (KXS>1) OBS data
 282 !
 283 ! (1) To compute the heights based on pressure(P), and available
 284 !     temperature(T) and specific humidity(Q) under the hydrostatic
 285 !     assumption.
 286 !
 287 !     Using the available observed heights to calibrate the computed
 288 !     heights, then these calibrated and computed heights are used
 289 !     where the observed heights are missing.
 290 !
 291 ! (2) When the hydrostatic assumption can not be applied at the bottom
 292 !     (K_START > 1) or top (K_TOP < KXS) part of the atmosphere, the
 293 !     model reference state is used to derive the missing heights based
 294 !     on the observed pressure with the calibration from the available
 295 !     observed heights.
 296 !
 297 ! Note:
 298 !
 299 !   In case (1), the derived heights have very high accuracy since only
 300 !     a high accurate "hydrostatic assumption" used here.
 301 !
 302 !   In case (2), the model reference state parameters are used. These
 303 !     parameters are consistent with the model BKG, and independent of
 304 !     the model domain settings. Of course, the model reference state
 305 !     atmosphere is not a real atmosphere, but it's consistent with the model.
 306 !     With the observed heights calibrated, these derived heights
 307 !     also included some information from OBS (pressure and height).
 308 !
 309 !   This subroutine can not be applied to the SIGLE-LEVEL(KXS=1) data.
 310 !   For those OBS, the missing heights can only be derived from the
 311 !   observed pressure and the model BKG and reference state (subroutine
 312 !   recover_hpt).
 313 !
 314 !                                    Yong-Run Guo
 315 !                                      11/18/00
 316 !
 317 !-------------------------------------------------------------------------
 318
 319   IMPLICIT NONE
 320
 321   INTEGER,                      INTENT(in)    :: KXS
 322   REAL        , DIMENSION(KXS), INTENT(in)    :: T, Q, P
 323   TYPE (field), DIMENSION(KXS), INTENT(inout) :: HGT
 324   INTEGER,                      INTENT(out)   :: k_start,k_top
 325
 326   INTEGER                :: k, kk, kwk, L, K0, K1, K2
 327   REAL                   :: height0, height1, diff_hp, &
 328                             TM1, TM2, ALNP, DDH, DDP, DHH, &
 329                             dh1, dh2, dp1, dp2, aa, bb, cc
 330   INTEGER,dimension(9000) :: KP
 331   REAL   ,DIMENSION(9000) :: TWK, QWK, PWK, HWK, DHGT
 332   LOGICAL                :: Vert_ok
 333
 334   TYPE (field), DIMENSION(KXS) :: HGT0, HGT1
 335
 336   include 'constants.inc'
 337 ! -------------------------------------------------------------------
 338
 339      HGT0 = HGT
 340
 341 ! .. get data at the first levels where both pressure and height
 342 !    available
 343 !
 344
 345 !    Find the first level with height, pressure and temperature
 346      K_top   = kxs+1
 347      K_start = 0
 348
 349 first_level: &
 350      DO k = 1, kxs
 351
 352         IF (.NOT. eps_equal(P  (k)     , missing_r, 1.) .AND. &
 353             .NOT. eps_equal(HGT(k)%data, missing_r, 1.) .AND. &
 354             .NOT. eps_equal(T  (k)     , missing_r, 1.)) THEN
 355              K_start = k
 356              EXIT first_level
 357         ENDIF
 358
 359      ENDDO first_level
 360
 361     !
 362     ! Obs without temperature and/or height are processed here
 363     !
 364      IF (K_start == 0) THEN
 365
 366         DO k = 1, kxs
 367
 368            IF (.NOT. eps_equal(P  (k)     , missing_r, 1.) .AND. &
 369                      eps_equal(HGT(k)%data, missing_r, 1.)) THEN
 370
 371                                HGT(k)%data = ref_height (P  (k) )
 372                                HGT(k)%qc   = reference_atmosphere
 373            ENDIF
 374
 375         ENDDO
 376
 377     ! To avoid the duplicated heights between the computed and
 378     ! observed heights:
 379
 380 1999    Vert_ok = .True.
 381
 382         DO k = 2, kxs
 383           if ((P(k) < P(K-1)) .and. (HGT(k)%data > HGT(k-1)%data)) then
 384             cycle
 385           else
 386             Vert_ok = .False.
 387             if (HGT(k)%qc <= 0) then
 388             ! HGT at level k-1 is computed:
 389                if (k > 2) then
 390                   AA = P(k  )-P(k-2)
 391                   BB = P(k-1)-P(k-2)
 392                   CC = AA - BB
 393                   HGT(k-1)%data = (HGT(k)%data*BB + HGT(k-2)%data * CC) / AA
 394 ! Y.-R. Guo (10/25/2005), Y.-R. Guo (01/16/2006):
 395                   if (HGT(k-1)%qc>0) HGT(k-1)%qc   = - HGT(k-1)%qc
 396 !     print '(" Computed k=",i3,"  pp,hh,qc:",2f11.2,i8)', k-1,p(k-1),HGT(k-1)%data, HGT(k-1)%qc
 397
 398                else
 399                   if ( k <= kxs-1 ) then
 400                     AA = P(k+1) - P(k)
 401                     BB = P(k-1) - P(k)
 402                     CC = AA - BB
 403                     HGT(k-1)%data = (HGT(k+1)%data*BB + HGT(k)%data * CC) / AA
 404                   else
 405 ! Y.-R. Guo (1/31/2008): must be processed separately when k >= kxs
 406                     AA = HGT(k)%data - ref_height (P  (k) )
 407                     HGT(k-1)%data = HGT(k-1)%data + AA
 408                   endif
 409 ! Y.-R. Guo (10/25/2005), Y.-R. Guo (01/16/2006):
 410                   if (HGT(k-1)%qc>0) HGT(k-1)%qc   = - HGT(k-1)%qc
 411                endif
 412             else
 413             ! HGT at level k is computed
 414                if ( k <= kxs-1 ) then
 415                  AA = P(k+1) - P(k-1)
 416                  BB = P(k)   - P(k-1)
 417                  CC = AA - BB
 418                  HGT(k)%data = (HGT(k+1)%data*BB + HGT(k-1)%data * CC) / AA
 419                else
 420 ! Y.-R. Guo (1/31/2008): must be processed separately when k >= kxs
 421                  AA = HGT(k-1)%data - ref_height (P  (k-1) )
 422                  HGT(k)%data = HGT(k)%data + AA
 423                endif
 424 ! Y.-R. Guo (10/25/2005), Y.-R. Guo (01/16/2006):
 425                if (HGT(k-1)%qc>0) HGT(k-1)%qc   = - HGT(k-1)%qc
 426             endif
 427           endif
 428         ENDDO
 429         if (.Not. Vert_ok) goto 1999
 430
 431         K_start = 1
 432         K_top   = kxs
 433
 434 ! Y.-R. Guo (10/25/2005):
 435         DO k = 1, kxs
 436            HGT(k)%qc   = abs(HGT(k)%qc)
 437 !           print '(i3,"  p,h,qc:",2f12.2,i8)', k, p(k), HGT(k)%data, HGT(k)%qc
 438         ENDDO
 439
 440         RETURN
 441
 442      ENDIf
 443
 444 !  .. if k_start > 1, using the model reference state to get
 445 !     the height at any level k below the level k_start,
 446 !     keep the height difference between level k_start and level k
 447 !     same as between the observed height and the model reference height
 448 !     at level k_start
 449
 450      IF (k_start > 1) THEN
 451
 452          !  Model reference height
 453
 454          height1 = Ref_height (p(k_start))
 455
 456          DO k = k_start-1, 1, -1
 457
 458            !  Missing height
 459
 460            IF (eps_equal (hgt(k)%data, missing_r, 1.)) THEN
 461
 462                IF (p(k)-p(k_start) > 20000.) THEN
 463
 464                ! too far below (200mb) the level where the OBS height available
 465
 466                   HGT(k)%data  = Ref_height(P(k))
 467                   HGT(k)%qc    = reference_atmosphere
 468
 469                ELSE
 470
 471                   HGT(k)%data  = HGT(k_start)%data - height1 &
 472                                + Ref_height(P(k))
 473                   HGT(k)%qc    = reference_OBS_scaled
 474
 475                ENDIF
 476
 477           ENDIF
 478
 479        ENDDO
 480
 481      ENDIF
 482
 483      !  Use the hydrostatic equation to correct the heigt
 484
 485      kwk = 0
 486
 487 temp_search: &
 488       DO k = k_start, KXS
 489
 490        IF (.NOT.eps_equal(T(k), missing_r, 1.)) THEN
 491             kwk = kwk+1
 492             pwk (kwk) = P(k)
 493             twk (kwk) = T(k)
 494             qwk (kwk) = q(k)
 495       ENDIF
 496
 497      ENDDO temp_search
 498
 499      HWK(1) = HGT(k_start)%data * G
 500
 501 !     ... INTEGRATE HYDROSTATIC EQN
 502 !
 503 hydro_int:  &
 504       DO K=2,KWK
 505
 506          ALNP = ALOG (PWK(K-1) /PWK(K)) * gasr
 507
 508          IF (.NOT.eps_equal(QWK(k), missing_r, 1.)) THEN
 509               TM2 = TWK(K  )*(1.+0.608*QWK(K  ))
 510          ELSE
 511               TM2 = TWK(K  )
 512          ENDIF
 513
 514          IF (.NOT.eps_equal(QWK(k-1), missing_r, 1.)) THEN
 515               TM1 = TWK(K-1)*(1.+0.608*QWK(K-1))
 516          ELSE
 517               TM1 = TWK(K-1)
 518          ENDIF
 519               HWK(K) = HWK(K-1) + .5*(TM1+TM2) * ALNP
 520       ENDDO hydro_int
 521
 522 !
 523 ! .. CALIBRATION OF THE HWK BASED ON THE AVAILABLE HGT:
 524
 525       K0 = 1
 526       KP(1) = 1
 527       DHGT(1) = 0
 528
 529 calibration: &
 530       DO K = 1,KXS
 531          IF (eps_equal(HGT(K)%data, missing_r, 1.)) then
 532           ! nothing
 533         ELSE
 534           DO KK = 1,KWK
 535           IF (P(K).EQ.PWK(KK)) THEN
 536             K0 = K0+1
 537             KP(K0) = KK
 538             DHGT(K0) = HWK(KK)/G - HGT(K)%data
 539             CYCLE calibration
 540           ENDIF
 541           ENDDO
 542         ENDIF
 543
 544      ENDDO calibration
 545
 546 !  .. levels between k_start(K0=1) and K0-1
 547 !
 548      DO L = 1,K0-1
 549         K1 = KP(L)
 550         K2 = KP(L+1)
 551         DO KK = K1,K2-1
 552           DDH = DHGT(L+1) - DHGT(L)
 553           DDP = ALOG(PWK(K2)/PWK(K1))
 554           DHH = DHGT(L) + ALOG(PWK(KK)/PWK(K1))*DDH/DDP
 555           HWK(KK) = HWK(KK)/G - DHH
 556         END DO
 557      END DO
 558 !
 559 !  .. Above the level KP(K0):
 560
 561      DO K = KP(K0),KWK
 562          HWK(K) = HWK(K)/G - DHGT(K0)
 563      END DO
 564 !     WRITE(0,'(I3,2X,4E15.5)') &
 565 !        (L,PWK(L),TWK(L),QWK(L),HWK(L),L=1,KWK)
 566
 567 !
 568 ! .. FILL BACK THE HEIGHTS AT THE LEVELS WHERE TEMPERATURE AVAILABLE:
 569
 570      height1 = Ref_height(Pwk(kwk))
 571
 572 above_k_start: DO K = k_start,KXS
 573
 574      if (abs(P(K) - PWK(KWK)) < 0.01) k_top = k
 575
 576 !     if (eps_equal(HGT(k)%data, missing_r, 1.)) then
 577
 578       IF (P(K) >= PWK(KWK)) then
 579
 580         DO KK = 1,KWK
 581           IF (P(K).EQ.PWK(KK)) THEN
 582              HGT(K)%data = HWK(KK)
 583           ELSE IF (KK.LT.KWK .AND. &
 584                  P(K).LT.PWK(KK) .AND. P(K).GT.PWK(KK+1)) THEN
 585 !
 586 ! .. Get the interpolated heights at the levels temperature unavailable:
 587
 588             ALNP = ALOG (P(K)/PWK(KK)) / ALOG(PWK(KK+1)/PWK(KK))
 589             HGT(K)%data = HWK(KK) + ALNP*(HWK(KK+1)-HWK(KK))
 590           ENDIF
 591         ENDDO
 592         HGT(k)%qc    =  Hydrostatic_recover
 593       ELSE
 594         if ((PWK(KWK)-P(K)) > 10000.) then
 595
 596      ! too far above (100mb) the level where the OBS height available
 597
 598            HGT(k)%data  = Ref_height(P(k))
 599            HGT(k)%qc    = reference_atmosphere
 600         else
 601            HGT(k)%data  = Hwk(kwk) - height1 &
 602                         + Ref_height(P(k))
 603            HGT(k)%qc    = reference_OBS_scaled
 604         endif
 605       ENDIF
 606
 607 !     endif
 608
 609     ENDDO above_k_start
 610
 611 !
 612 ! Non-hydrastatic adjustment:
 613 !
 614 ! In case of only P and H observed without T, TD, the hydrostatic
 615 ! height may be different from the observed H. We must keep this
 616 ! observed H, and adjust the calculated heights at the adjacent
 617 ! levels to avoid inconsistency between P and H when having
 618 ! the dense observed levels.
 619 !
 620 !                                 Yong-Run Guo  12/06/2001
 621 ! ------------------------------------------------------------
 622 ! .. Find the starting level with the observed height:
 623 !                                     Fixed the bug on 07/08/2004
 624     k0 = 1
 625     do k = 1,kxs
 626       if (HGT0(k) % qc == 0) then
 627         k0 = k
 628         exit
 629       endif
 630     enddo
 631 ! ------------------------------------------------------------
 632     k1 = 0
 633 ! .. Keep the calculated height:
 634     HGT1 = HGT
 635     do k = k0, kxs
 636
 637 ! .. Find the starting level (k1+1) without the observed height:
 638       if (HGT0(k) % qc /= 0 .and. k1 == 0) then
 639          k1 = k-1
 640       else if (HGT0(k) % qc == 0 .and. &
 641                abs(HGT0(k)%data - HGT(k)%data) <= 0.10*HGT(k)%data ) then
 642          HGT1(k) = HGT0(k)
 643       endif
 644
 645 ! .. Find the ending level (k2-1) without the observed height:
 646       if (HGT0(k) % qc == 0 .and. k1 > 0 .and. k > k1+1) then
 647          k2 = k
 648
 649          if (abs(HGT0(k2)%data - HGT(k2)%data) > 0.10*HGT(k2)%data) then
 650
 651 ! ....When the difference between the observed height and retrieved height
 652 !     greater than 10% of Hydro. Retrv. height, ignore the observed height,
 653 !     and no adjustment done:
 654
 655            HGT0(k) = HGT1(k)
 656
 657          else
 658
 659            HGT1(k) = HGT0(k)
 660
 661 ! .. Adjust the calculated hydristatic heights from level k1+1 to k2-1:
 662          dp2 = p(k2) - p(k1)
 663          dh1 = HGT0(k1)%data - HGT(k1)%data
 664          dh2 = HGT0(k2)%data - HGT(k2)%data
 665          do kk = k1+1, k2-1
 666            dp1 = p(kk) - p(k1)
 667            HGT1(kk)%data = HGT(kk)%data + dh1*(1-dp1/dp2) + dh2*dp1/dp2
 668            HGT1(kk)%qc   = HGT(kk)%qc
 669            HGT1(kk)%error= HGT(kk)%error
 670          enddo
 671         endif
 672          k1 = 0
 673          k2 = 0
 674       endif
 675
 676     enddo
 677
 678 ! .. Fill back to HGT with HGT1:
 679     HGT = HGT1
 680
 681  END subroutine recover_h_from_ptq
 682 !------------------------------------------------------------------------------!
 683
 684  SUBROUTINE reorder(obs, i, order_component, failed)
 685 !------------------------------------------------------------------------------!
 686
 687   IMPLICIT NONE
 688
 689   INTEGER,       INTENT (in)    :: i
 690   CHARACTER*(*), INTENT (in)    :: order_component
 691   TYPE (report), INTENT (inout) :: obs
 692   logical                       :: failed
 693
 694   INCLUDE 'missing.inc'
 695
 696   TYPE (measurement), pointer :: current, new, tmp, pre
 697 !
 698   integer :: num, ii, num_loop, nlevels
 699   logical :: need_check
 700
 701   failed = .false.
 702
 703   ii = 0
 704   current => obs % surface
 705   count_levels:  do while (associated(current))
 706      ii = ii + 1
 707      current => current%next
 708   enddo count_levels
 709   nlevels = ii
 710
 711   current => obs%surface
 712
 713 !--Missing Check:
 714
 715    num = 0
 716    ii  = 0
 717
 718    missing_check: do while (associated(current))
 719       num = num + 1
 720       if(eps_equal(comp_field(current,order_component), &
 721                                        missing_r, 1.0)) then
 722 ! remove the missing data link:
 723          num=num-1
 724       end if
 725       current => current%next
 726    end do missing_check
 727
 728    if (num /= nlevels) then
 729      write(0,'(/I5,A,A,A/3X,A,A,A,A,A,2I5,A,f10.3,A,f10.3)') I, &
 730         ' There are ',order_component,&
 731         ' at several levels missing. Reordering can not be done.', &
 732         ' FM=',obs%info%platform(4:5),' ID=',obs%location%id(1:5), &
 733         ' nlevels, num:', nlevels, num, &
 734         ' lat=', obs%location%latitude, ' long=', obs%location%longitude
 735 !     STOP 'in_reorder'
 736      failed = .true.
 737    endif
 738
 739 !  Re-set the number of levels to OBS
 740
 741    obs % info % levels = nlevels
 742
 743    if (nlevels <= 1) return
 744
 745    ii = 0
 746    num_loop = 0
 747
 748    need_check = .true.
 749
 750    put2order: do while(need_check)
 751
 752    num_loop = num_loop + 1
 753
 754       head_check: do while(need_check)
 755
 756          current => obs%surface
 757
 758 !--Head Check
 759
 760          if(comp_field(current,order_component) < &
 761             comp_field(current%next,order_component)) then
 762 !--------Cut the first on off and put it in pointer tmp:
 763             tmp => current
 764             current => current%next
 765             nullify(tmp%next)
 766
 767             new  => current
 768
 769             obs%surface => new
 770
 771 !--------if first pressure is not good, through tmp away
 772
 773             if(eps_equal(comp_field(current,order_component), &
 774                                             missing_r, 1.0)) then
 775                num=num-1
 776           ! Don't allow happened
 777                STOP 'in reorder_missing data'
 778 !               cycle head_check
 779             end if
 780
 781 !--------Now we need to put pointer tmp to a place.
 782
 783             current => obs%surface
 784
 785             new => obs%surface%next
 786
 787             nullify(current%next)
 788
 789             allocate(current%next)
 790             current%next => tmp
 791
 792             tmp%next => new
 793          end if
 794
 795          need_check = .false.
 796       end do head_check
 797
 798       if(num < 3) exit put2order
 799
 800 !     write(0,'(/)')
 801
 802 !     current => obs%surface
 803
 804 !     ii=0
 805
 806 !     do while (associated(current))
 807 !        ii = ii + 1
 808 !        write(0,'(A,I3,A,I3,3f12.2)') 'num_loop=',num_loop, &
 809 !              '  I,P,H,T:', ii, &
 810 !              current%meas%pressure%data,  current%meas%height%data, &
 811 !              current%meas%temperature%data
 812 !        current => current%next
 813 !     end do
 814
 815 !-----Put others in order
 816
 817       pre     => obs%surface
 818       current => obs%surface%next
 819       new     => obs%surface%next%next
 820
 821       do while (associated(new))
 822          if(comp_field(current,order_component) < &
 823             comp_field(current%next,order_component)) then
 824
 825             tmp => new%next
 826
 827             nullify(pre%next)
 828             nullify(current%next)
 829             nullify(new%next)
 830
 831             current%next => tmp
 832             new%next => current
 833             pre%next => new
 834
 835             need_check = .true.
 836
 837             exit
 838          end if
 839
 840          pre     =>     pre%next
 841          new     =>     new%next
 842          current => current%next
 843       end do
 844
 845 !     write(0,'(/)')
 846
 847 !     current => obs%surface
 848
 849 !     ii=0
 850
 851 !     do while (associated(current))
 852 !        ii = ii + 1
 853 !        write(0,'(A,I3,A,I3,3f12.2)') 'num_loop=',num_loop, &
 854 !              '  I,P,H,T:', ii, &
 855 !              current%meas%pressure%data,  current%meas%height%data, &
 856 !              current%meas%temperature%data
 857 !        current => current%next
 858 !     end do
 859    end do put2order
 860
 861 !   write(0,'(//)')
 862
 863 !   current => obs%surface
 864
 865 !   ii=0
 866 !   do while (associated(current))
 867 !      ii = ii + 1
 868 !      write(0,'(A,I3,3f12.2)') 'I,P,H,T:', ii, &
 869 !            current%meas%pressure%data,  current%meas%height%data, &
 870 !            current%meas%temperature%data
 871 !       current => current%next
 872 !   end do
 873
 874  END subroutine reorder
 875 !------------------------------------------------------------------------------!
 876
 877  FUNCTION comp_field(current, order_component) result(xxx)
 878 !------------------------------------------------------------------------------!
 879
 880 type (measurement), pointer :: current
 881 character*(*),   intent(in) :: order_component
 882 real                        :: xxx
 883
 884      SELECT CASE (order_component)
 885
 886        CASE ('pressure')
 887
 888          xxx = current%meas%pressure%data
 889
 890        CASE ('height')
 891
 892          xxx =  current%meas%height%data
 893
 894        CASE DEFAULT
 895
 896          WRITE(0,'(A,A,A)') 'order_component=',order_component, &
 897                             ' is not defined correctly'
 898          STOP 'in_reorder'
 899
 900      END SELECT
 901
 902  END function comp_field
 903 !------------------------------------------------------------------------------!
 904
 905 END MODULE module_recoverh