var/obsproc/src/module_recoverp.F90

   1 MODULE module_recoverp
   2
   3    USE module_type
   4    USE module_func
   5    USE module_mm5
   6
   7    INCLUDE 'missing.inc'
   8
   9 CONTAINS
  10
  11 !------------------------------------------------------------------------------
  12
  13  SUBROUTINE recover_pressure_from_height (max_number_of_obs , obs , &
  14                                           number_of_obs, print_hp_recover)
  15 !------------------------------------------------------------------------------
  16 !
  17 !  Since the data merging is based on the observed pressure (see
  18 !  module_obs_merge.F90/subroutine link_level), the observed pressure
  19 !  is not allowed to be missing.
  20 !
  21 !  The pressure missing from the decoded OBS sounding (SOUND, AIREP) data
  22 !  is recovered
  23 !    (1) from the observed heights available with the hydrostatic assumption
  24 !        (dP/dH)_hydr if possible,
  25 !    (2) or from  the observed heights available with the model reference
  26 !        state (dP/dH)_ref when hydrostatic interpolation is impossible.
  27 !
  28 !  There are two steps to do this:
  29 !
  30 !    (1) To guarantee the original pressure and height are correct, first,
  31 !        a sequence check, i.e.
  32 !             height at the next level >   height at the current level and
  33 !             pressure at the next level < pressure at the current level,
  34 !        is completed (subroutine hp_sequence_check).
  35 !
  36 !    (2) second, a consistency check between the observed pressure and
  37 !        height agaist the model reference state is completed
  38 !        (subroutine pressure_recover, FUNCTION hp_consist).
  39 !
  40 !    (3) To fill the missing pressure by using the observed pressure and
  41 !        heights that passed the above checks (subroutine recover_p_from_h).
  42 !
  43 !                                       Yong-Run Guo
  44 !                                       11/16/00
  45 !
  46 !------------------------------------------------------------------------------
  47
  48    IMPLICIT NONE
  49
  50    INTEGER,       INTENT ( IN )                :: max_number_of_obs
  51    TYPE (report), DIMENSION (max_number_of_obs):: obs
  52    INTEGER , INTENT ( IN )                     :: number_of_obs
  53    LOGICAL ,      INTENT (IN)                  :: print_hp_recover
  54
  55    TYPE (measurement), POINTER                 :: current
  56    INTEGER                                     :: qc_flag
  57    INTEGER                                     :: i, nlevel
  58
  59    CHARACTER (LEN=80)                          :: filename
  60    LOGICAL                                     :: connected, consistent
  61    INTEGER                                     :: iunit, io_error
  62    CHARACTER (LEN = 32)                        :: proc_name = &
  63                                                  "recover_pressure_from_height"
  64
  65    INCLUDE 'platform_interface.inc'
  66
  67 !------------------------------------------------------------------------------!
  68
  69       WRITE (0,'(A)')  &
  70 '------------------------------------------------------------------------------'
  71       WRITE (UNIT = 0, FMT = '(A,/)') 'PRESSURE RECOVERED FROM HEIGHTS:'
  72
  73       IF (print_hp_recover) THEN
  74
  75       filename = 'obs_recover_pressure.diag'
  76       iunit    = 999
  77
  78       INQUIRE (UNIT = iunit, OPENED = connected )
  79
  80       IF (connected) CLOSE (iunit)
  81
  82       OPEN (UNIT = iunit , FILE = filename , FORM = 'FORMATTED'  , &
  83             ACTION = 'WRITE' , STATUS = 'REPLACE', IOSTAT = io_error )
  84
  85       IF (io_error .NE. 0) THEN
  86           CALL error_handler (proc_name, &
  87          "Unable to open output diagnostic file. " , filename, .TRUE.)
  88       ELSE
  89           WRITE (UNIT = 0, FMT = '(A,A,/)') &
  90          "Diagnostics in file ", TRIM (filename)
  91       ENDIF
  92
  93       ENDIF
  94
  95
  96 ! 2.  ESTIMATE P
  97 ! ==============
  98
  99 ! 2.1 Print parameters
 100 !     ----------------
 101
 102       IF (print_hp_recover) THEN
 103
 104       WRITE (UNIT = IUNIT, FMT = '(A,/)') 'PRESSURE RECOVERED FROM HEIGHTS:'
 105
 106       WRITE(UNIT = iunit , FMT = '(A,/,3(A,F10.2,A),/,3(A,F10.2,A),/)') &
 107         "REFERENCE STATE: ", &
 108         " HTOP = ", htop,"m,  ",  &
 109         " PTOP = ", ptop,"Pa, ", &
 110         " PS0  = ", ps0 ,"Pa, ", &
 111         " TS0  = ", ts0 ,"K,  ", &
 112         " TLP  = ", tlp ,"K,  ",&
 113         " DIS  = ", DIS(idd),"m."
 114
 115         ENDIF
 116
 117 ! 2.3 Loop over obs
 118 !     -------------
 119
 120 loop_all_1:&
 121       DO i = 1, number_of_obs
 122
 123          IF ((obs (i) % info % discard)  .OR. .NOT. ASSOCIATED &
 124              (obs (i) % surface)) THEN
 125
 126              CYCLE loop_all_1
 127
 128          ENDIF
 129
 130 ! 3.  SINGLE LEVEL OBS
 131 ! ====================
 132
 133 levels:&
 134          IF (obs (i) % info % bogus .OR. (.NOT. obs (i) % info % is_sound)) THEN
 135
 136              ! IF pressure is missing, height should be present
 137
 138              IF (eps_equal &
 139                 (obs (i) % surface % meas % pressure % data, missing_r, 1.))THEN
 140                  obs (i) % surface % meas % pressure % data = floor (ref_pres &
 141                 (obs (i) % surface % meas % height   % data))
 142                  obs (i) % surface % meas % pressure % qc   = &
 143                                                          reference_atmosphere
 144              ELSE IF (obs (i) % surface % meas % pressure % qc == 0 .and. &
 145                       obs (i) % surface % meas % height   % qc == 0 .and. &
 146                       eps_equal &
 147                       (obs (i) % surface % meas % height % data, 0., 1.)) then
 148
 149              !hcl-note: obsproc should not check for qc on input file
 150              !hcl-note: qc should be assigned by obsproc
 151              ! When both h and p have qc = 0 and h = 0. (at sea level
 152              ! because the pressure is MSLP for surface data, check
 153              ! the consistency between h and p
 154
 155                       consistent = hp_consist (  &
 156                               obs (i) % surface % meas % height   % data, &
 157                               obs (i) % surface % meas % pressure % data, &
 158                               print_hp_recover, iunit)
 159                       if (.not.consistent) obs (i) % info % discard = .true.
 160
 161              ENDIF
 162
 163          ELSE levels
 164
 165 ! 3. MULTI LEVEL OBS
 166 ! ==================
 167
 168 ! 3.1 Height and pressure sequence check
 169 !     ----------------------------------
 170
 171           CALL hp_sequence_check (obs (i), i, nlevel, print_hp_recover, iunit)
 172
 173           CALL hp_missing_check (obs(i), i, nlevel)
 174
 175 ! 3.2 Single level sounding are removed
 176 !     ---------------------------------
 177
 178           IF ((nlevel == 1) .AND. eps_equal &
 179               (obs (i) % surface % meas % pressure% data, missing_r, 1.0)) THEN
 180
 181                obs (i) % info % discard = .TRUE.
 182
 183                IF (print_hp_recover) THEN
 184
 185                WRITE (UNIT = iunit, FMT = '(A,A,I5,A,A)') &
 186               "In recover_pressure_from_height: ","I = ", I," ID = ", &
 187                obs (i) % location % id (1:5)
 188
 189                WRITE (UNIT = iunit, FMT = '(A,I5,A,F9.0,A,F9.0)') &
 190                "nlevel = ",  nlevel, &
 191                "height = ",  obs (i) % surface % meas % height   % data, &
 192                "pressure = ",obs (i) % surface % meas % pressure % data
 193
 194                 ENDIF
 195
 196               CYCLE loop_all_1
 197
 198           ENDIF
 199
 200 ! 3.2 If the actual number of valid levels differs from expected, print it
 201 !     --------------------------------------------------------------------
 202
 203           IF (nlevel /= obs (i) % info % levels) THEN
 204
 205               IF (print_hp_recover) THEN
 206
 207                   WRITE (UNIT = iunit, FMT = '(3A,I4,A,I4)') &
 208                  "ID = ", obs(i)%location%id(1:5), ", expect ", &
 209                   obs (i) % info % levels, " levels, get ", nlevel
 210
 211               ENDIF
 212
 213           ENDIF
 214
 215 ! 3.3 Pressure is recovered based on height if the pressure is missing
 216 !     ----------------------------------------------------------------
 217
 218           CALL recover_p_from_h (obs (i), i, nlevel, print_hp_recover, iunit)
 219
 220         ENDIF levels
 221
 222       ENDDO loop_all_1
 223
 224       IF (print_hp_recover) CLOSE (iunit)
 225
 226  END SUBROUTINE recover_pressure_from_height
 227
 228 !------------------------------------------------------------------------------!
 229
 230  SUBROUTINE hp_sequence_check (obs, i, nlevel, print_hp_recover, iunit)
 231 !------------------------------------------------------------------------------!
 232
 233       INTEGER,       INTENT (in)    :: i, iunit
 234       INTEGER,       INTENT (out)   :: nlevel
 235       LOGICAL,       INTENT (in)    :: print_hp_recover
 236       TYPE (report), INTENT (inout) :: obs
 237
 238       TYPE (measurement), POINTER   :: current
 239       REAL                          :: height1, height2, &
 240                                        press1 , press2
 241       LOGICAL                       :: first_h, first_p
 242 !------------------------------------------------------------------------------
 243
 244       current => obs % surface
 245       first_h = .false.
 246       first_p = .false.
 247       nlevel  = 0
 248
 249 loop_hp_level: &
 250     DO WHILE (ASSOCIATED (current))
 251
 252        nlevel = nlevel + 1
 253
 254 ! 1.  HEIGHT CHECK
 255 ! ================
 256
 257        IF (.NOT. eps_equal (current%meas%height%data, missing_r, 1.0)) then
 258
 259            IF (.NOT. first_h) then
 260
 261                 height1 = current%meas%height%data
 262                 first_h = .true.
 263
 264            ELSE
 265
 266                 height2 = current%meas%height%data
 267
 268                 IF (height2 <= height1) then
 269
 270                     current%meas%height%qc = missing
 271                     current%meas%height%data = missing_r
 272
 273                     IF (print_hp_recover) THEN
 274
 275                         WRITE (UNIT = iunit, FMT = '(A,/,4A,A,I3,2(A,f12.2))') &
 276                         "HEIGHT VIOLATION: ",           &
 277                         " FM = ", obs % info % platform (4:6), &
 278                         " ID = ", obs % location % id,         &
 279                         " LEVEL = ", nlevel,                   &
 280                         " H2 = ", height2,                     &
 281                         " H1 = ", height1
 282
 283                     ENDIF
 284
 285                 ELSE
 286
 287                     height1 = height2
 288
 289                 ENDIF
 290
 291              ENDIF
 292
 293           ENDIF
 294
 295 ! 2.  PRESSURE CHECK
 296 ! ==================
 297
 298           IF (.NOT.eps_equal(current%meas%pressure%data, missing_r, 1.0)) then
 299
 300              IF (.NOT.first_p) THEN
 301                   press1  = current%meas%pressure%data
 302                   first_p = .true.
 303              ELSE
 304                   PRESS2 = current%meas%pressure%data
 305
 306                   IF (press2 >= press1) THEN
 307                       current%meas%pressure%data = missing_r
 308                       current%meas%pressure%qc   = missing
 309
 310                   IF (print_hp_recover) THEN
 311
 312                       WRITE (UNIT = iunit, FMT = '(A,/,4A,A,I3,2(A,f12.2))') &
 313                       "PRESSURE VIOLATION: ",           &
 314                       " FM = ", obs % info % platform (4:6), &
 315                       " ID = ", obs % location % id,         &
 316                       " LEVEL = ", nlevel,                   &
 317                       " P2 = ", press1,                      &
 318                       " P1 = ", press1
 319
 320                  ENDIF
 321
 322                ELSE
 323
 324                  press1 = press2
 325
 326                ENDIF
 327
 328              ENDIF
 329
 330           ENDIF
 331
 332           ! Next level
 333
 334           current => current%next
 335
 336       ENDDO loop_hp_level
 337
 338 END subroutine hp_sequence_check
 339 !------------------------------------------------------------------------------!
 340
 341  SUBROUTINE recover_p_from_h(obs, i, nlevel, print_hp_recover, iunit)
 342 !------------------------------------------------------------------------------!
 343
 344    IMPLICIT NONE
 345
 346    INTEGER,       INTENT (in)    :: i, iunit, nlevel
 347    LOGICAL,       INTENT (in)    :: print_hp_recover
 348    TYPE (report), INTENT (inout) :: obs
 349
 350    TYPE (measurement), POINTER                 :: current
 351    TYPE (field)      , dimension(nlevel)       :: pp, hh
 352    INTEGER                                     :: k, Lb, Le
 353    REAL                                        :: elev
 354    LOGICAL                                     :: do_it
 355 !------------------------------------------------------------------------------!
 356
 357 ! 1.  GET THE P/H PROFILE
 358 ! =======================
 359
 360       k  = 0
 361       current => obs % surface
 362
 363 loop_level:&
 364       DO WHILE (ASSOCIATED (current))
 365
 366          k = k + 1
 367          hh (k) = current % meas % height
 368          pp (k) = current % meas % pressure
 369          current => current % next
 370
 371       ENDDO loop_level
 372
 373 ! 2.  PRESSURE FILED BASED ON THE HEIGHT IF THE PRESSURE IS MISSING
 374 ! =================================================================
 375
 376 ! 2.1 Recover pressure
 377 !     ----------------
 378
 379       CALL pressure_recover (nlevel, pp, hh, iunit, LB, LE, do_it, print_hp_recover)
 380
 381 ! 2.2 Print status
 382 !     ------------
 383
 384       IF (do_it) THEN
 385
 386         IF (print_hp_recover) THEN
 387
 388             WRITE(UNIT = iunit,FMT = '(5A)', ADVANCE = 'no') &
 389             "== PRESSURE RECOVER DONE: ",    &
 390             " for FM=", obs % info % platform (4:6), &
 391             " ID=",     obs % location % id
 392
 393         ENDIF
 394
 395         IF (print_hp_recover) THEN
 396
 397              IF (LB > 1 .OR. LE < nlevel) THEN
 398                  WRITE(UNIT = iunit,FMT = '(A)') &
 399                " Reference pressure may have been used."
 400              ELSE
 401                  WRITE(UNIT = iunit,FMT = '(A)') " "
 402
 403              ENDIF
 404
 405         ENDIF
 406
 407
 408 ! 2.  SEND BACK TO OBS
 409 ! ====================
 410
 411         current => obs % surface
 412
 413         k = 0
 414
 415 back_level: &
 416         DO WHILE (ASSOCIATED (current))
 417
 418            k = k + 1
 419 !           current % meas % height % data   = CEILING (hh (k) % data)
 420            current % meas % height % data   = hh (k) % data
 421            current % meas % height % qc     =          hh (k) % qc
 422            current % meas % pressure % data = FLOOR (pp (k) % data)
 423            current % meas % pressure % qc   =        pp (k) % qc
 424            current => current % next
 425
 426         ENDDO back_level
 427
 428       ENDIF
 429
 430  END SUBROUTINE recover_p_from_h
 431 !------------------------------------------------------------------------------!
 432
 433  FUNCTION hp_consist (hin,pin,print_out,iunit) RESULT (hout)
 434 !------------------------------------------------------------------------------!
 435
 436     IMPLICIT NONE
 437
 438     REAL,    INTENT (in) :: hin, pin
 439     LOGICAL, INTENT (in) :: print_out
 440     INTEGER, INTENT (in) :: iunit
 441     LOGICAL              :: hout
 442
 443     REAL, parameter  :: hmax = 1000. ! a tolerance apart from
 444                                      ! the model reference state
 445                                      ! (50000/pin)*hmax
 446
 447     REAL             :: h_ref, hdiff
 448 !------------------------------------------------------------------------------!
 449       if (pin <= 0) then
 450         write(unit=0, fmt='(A,f12.2)') &
 451            "In function hp_consist, Pressure voilation: P=",pin
 452         hout = .false.
 453         return
 454       endif
 455
 456       hout = .true.
 457
 458       h_ref = Ref_height(pin)
 459
 460       hdiff  = ABS (hin - h_ref)
 461
 462       IF (hdiff > (50000/pin)*hmax) THEN
 463
 464           IF (print_out) THEN
 465               WRITE (UNIT = iunit, FMT = '(/,A,/,3(A,F12.3,/))') &
 466           "   Pressure / height inconsistency: ", &
 467           "   Pressure    =  ", pin, &
 468           "   height      =  ", hin, &
 469           "   ref_height  =  ", h_ref
 470           ENDIF
 471
 472           hout = .false.
 473
 474       ENDIF
 475
 476  END FUNCTION hp_consist
 477 !------------------------------------------------------------------------------!
 478 subroutine pressure_recover(level, pp, hh, iunit, LB, LE, do_it, print_hp_recover)
 479
 480    IMPLICIT NONE
 481
 482    INTEGER,                   INTENT(in)         :: level, iunit
 483    TYPE (field), dimension(level), INTENT(inout) :: hh, pp
 484    INTEGER,                   INTENT(out)        :: LB, LE
 485    LOGICAL,                   INTENT(out)        :: do_it
 486    LOGICAL,                   INTENT (in)        :: print_hp_recover
 487
 488    INTEGER               :: k, L, L1, L2, Lstart
 489    REAL                  :: height1, height2, &
 490                             press1 , press2 , bb, diff_pr
 491    LOGICAL               :: first_L, second_L, &
 492                             consist1, consist2
 493 ! -----------------------------------------------------------------
 494
 495    LB = 0
 496    LE = 0
 497    do_it = .false.
 498
 499    L1 = 0
 500    L2 = 0
 501    first_L  = .false.
 502    second_L = .false.
 503    consist1 = .true.
 504    consist2 = .true.
 505
 506 ! if there is ONLY one level OBS:
 507
 508    if (level == 1) then
 509      return
 510    endif
 511
 512    loop_level: do k = 1,level
 513
 514    if (.NOT.first_L) then
 515
 516    !  To find the first level with both height and pressure
 517
 518      if (.NOT.eps_equal (pp(k)%data, missing_r, 1.0) .and. &
 519          .NOT.eps_equal (hh(k)%data, missing_r, 1.0) ) then
 520          L1 = k
 521          LB = L1
 522          Lstart = L1
 523          first_L = .true.
 524      endif
 525    else
 526
 527    !  To find the second level with both height and pressure
 528
 529
 530      if (.NOT.eps_equal (pp(k)%data, missing_r, 1.0) .and. &
 531          .NOT.eps_equal (hh(k)%data, missing_r, 1.0) ) then
 532          L2 = k
 533          second_L = .true.
 534      endif
 535    endif
 536
 537    if (first_L) then
 538       height1 = hh(L1) % data
 539       press1  = pp(L1) % data
 540
 541       if (Lstart > 1) then
 542
 543    ! If the first level is not the level 1 and pp(1) is missing
 544           if (eps_equal(pp(1)%data, missing_r, 1.0)) then
 545           L2 = L1
 546           second_L = .true.
 547           L1 = 1
 548           Lstart = 1
 549           height1 = hh(L1) % data
 550    ! Set the correction from reference pressure (h1,h2) plus p2
 551           press1      = Ref_pres(height1)
 552           press2      = Ref_pres(hh(L2)%data)
 553           diff_pr     = press1 - press2
 554           press1      = pp(L2)%data + diff_pr
 555           pp(L1)%data = press1
 556           pp(L1)%qc   = reference_atmosphere
 557           else
 558             do L = 1, L2
 559     ! There are pressure missing below the level where both p and h available:
 560               if (eps_equal(pp(L)%data, missing_r, 1.0)) then
 561                  pp(L)%data = Ref_pres(hh(L)%data)
 562                  pp(L)%qc   = reference_atmosphere
 563               endif
 564             enddo
 565           endif
 566           do_it = .true.
 567       endif
 568
 569    ! height/pressure consistency check
 570       consist1 = hp_consist (height1,press1,print_hp_recover,iunit)
 571       if (.NOT.consist1) then
 572          WRITE(UNIT=IUNIT, FMT='(A,2F12.2)') &
 573             "FAILED IN H/P CONSISTENCY CHECK: H1,P1:",height1,press1
 574          first_L = .FALSE.
 575          hh(L1) % data = missing_r
 576          hh(L1) % qc   = missing
 577          consist1 = .true.
 578          CYCLE loop_level
 579       endif
 580    endif
 581
 582    if (second_L) then
 583
 584        height2 = hh(L2) % data
 585        press2  = pp(L2) % data
 586
 587        consist2 = hp_consist (height2,press2,print_hp_recover,iunit)
 588       if (.NOT.consist2) then
 589          WRITE(UNIT=IUNIT, FMT='(A,I4,2F12.2)') &
 590             "FAILED IN H/P CONSISTENCY CHECK: L2,H2,P2:",L2,height2,press2
 591          second_L = .FALSE.
 592          hh(L2) % data = missing_r
 593          hh(L2) % qc   = missing
 594          consist2 = .true.
 595          CYCLE loop_level
 596       endif
 597
 598    ! If L1 and L2 are consecutive
 599
 600      if (L2 <= (L1 + 1)) then
 601        L1 = L2
 602        second_L = .false.
 603      else
 604
 605    ! If L1 and L2 are not consecutive
 606
 607        BB = (height2 - height1)/ALOG(press2/press1)
 608        do L = L1+1, L2-1
 609 ! .. get the interpolated pressure from the height if pressure is missing:
 610          if (eps_equal(pp(L)%data, missing_r, 1.0)) then
 611            if (eps_equal(hh(L)%data, missing_r, 1.0)) then
 612               WRITE(0,'(/A,A)') &
 613                 "Both pressure and height are missing, ", &
 614                 "this should never be happened???"
 615               STOP 33333
 616            endif
 617            pp(L)%data = press1*EXP((hh(L)%data - height1)/BB)
 618            pp(L)%qc   = 0
 619            do_it = .true.
 620          endif
 621        enddo
 622        L1 = L2
 623        second_L = .false.
 624      endif
 625    endif
 626
 627    enddo loop_level
 628
 629    ! If L2 is not equal to level
 630
 631    LE = L2
 632    if (level > L2) then
 633      do k = L2+1, level
 634        if (eps_equal(pp(k)%data, missing_r, 1.0)) then
 635    ! Set the reference pressure
 636           if (L2 > 0) then
 637             press1  = Ref_pres(hh(L2)%data)
 638             press2  = Ref_pres(hh(k )%data)
 639             diff_pr = press1 - press2
 640             pp(k)%data = pp(L2)%data - diff_pr
 641           else
 642             pp(k)%data = Ref_pres(hh(k )%data)
 643           endif
 644           pp(k)%qc   =  reference_atmosphere
 645           do_it = .true.
 646        endif
 647      enddo
 648    endif
 649
 650 end subroutine pressure_recover
 651
 652  SUBROUTINE hp_missing_check (obs, i, nlevel)
 653 !------------------------------------------------------------------------------!
 654
 655       INTEGER,       INTENT (in)    :: i
 656       INTEGER,       INTENT (out)   :: nlevel
 657       TYPE (report), INTENT (inout) :: obs
 658
 659       TYPE (measurement), POINTER   :: current, tmp, new, pre
 660       integer                       :: nn
 661 !------------------------------------------------------------------------------
 662       current => obs % surface
 663       new => current % next
 664       nn = 0
 665
 666 loop_hp_level: &
 667     DO WHILE (ASSOCIATED (new))
 668
 669       nn = nn + 1
 670
 671       tmp => new % next
 672
 673 ! If both P and H missing at the next level?
 674       IF (eps_equal (new % meas % pressure % data, missing_r, 1.) .and. &
 675           eps_equal (new % meas % height % data, missing_r, 1.) ) THEN
 676
 677          write(0,'("==> Missing P and H:",i3,2f15.6)') nn, &
 678                 new%meas%pressure%data, new%meas%height%data
 679          nullify(current % next)
 680
 681          current % next => tmp
 682
 683 1001     continue
 684
 685          tmp => tmp%next
 686
 687 ! If both P and H missing at the next next level?
 688          if (associated (current%next) .and. &
 689              eps_equal (current%next%meas%pressure% data, missing_r, 1.) .and.&
 690              eps_equal (current%next%meas% height % data, missing_r, 1.) ) THEN
 691
 692              nullify(current % next)
 693              current % next => tmp
 694              goto 1001
 695
 696          endif
 697
 698          current => current % next
 699
 700       ELSE
 701
 702          current => current % next
 703
 704       ENDIF
 705       new => current % next
 706
 707     END DO loop_hp_level
 708
 709 ! To count the total number of levels:
 710       nn = 0
 711       current =>  obs % surface
 712
 713       do while (associated(current))
 714         nn = nn + 1
 715 !        write(0,'(I3," P,H:",2f15.5)') nn, &
 716 !           current%meas%pressure%data,  current%meas%height%data
 717         current => current%next
 718       enddo
 719
 720       nlevel = nn
 721       if (obs%info%levels /= nlevel ) then
 722         obs%info%levels = nlevel
 723         write(0, fmt='(a, i3, 2x, 2a, 2x, a, 2f10.3/)') &
 724              'After missing check: Number of levels = ', obs%info%levels, &
 725              'OBS=', obs%info%platform(1:12), &
 726              'LOC=', obs%location%latitude, obs%location%longitude
 727       endif
 728
 729 END  SUBROUTINE hp_missing_check
 730
 731 END MODULE module_recoverp