var/obsproc/src/obsproc.F90

   1 PROGRAM main_obsproc
   2
   3 !-----------------------------------------------------------------------------!
   4 ! Read decoded observations (output format of the MM5 "gts_decoder" and "fetch"
   5 !                           facilities, which is also the input format
   6 !                           of MM5 "rawins" and "little_r" programs).
   7 !                           Keep only records with at least height or pressure
   8 !                           in MM5 horizontal domain as defined in namelist
   9 ! Fill with pressure when missing,
  10 ! Sort observations by location and time,
  11 ! Merge space duplicate stations (same type, same location, same time),
  12 ! Remove time duplicate stations in time (same type & location, different
  13 ! time),
  14 ! Fill with height when missing,
  15 ! Check soundings (vertical consistency and super adiabatic),
  16 ! Remove observations above MM5 lid as defined by ptop in namelist,
  17 ! Estimate observational error,
  18 ! Write out for inputting into MM5 3D-VAR
  19 !
  20 !-----------------------------------------------------------------------------!
  21 !
  22 !  HISTORY:
  23 !
  24 !         F. VANDENBERGHE, March 2001
  25 !
  26 !         01/13/2003 - Updated for Profiler obs.           S. R. H. Rizvi
  27 !
  28 !         02/04/2003 - Updated for Buoy     obs.           S. R. H. Rizvi
  29 !
  30 !         02/11/2003 - Reviewed and modified for Profiler
  31 !                      and Buoy obs.                       Y.-R. Guo
  32 !
  33 !         05/23/2003 - GPS ZTD code added              L. Cucurull
  34 !         08/13/2003 - Reviewed                        Y.-R. Guo
  35 !         09/XX/2005 - Modified to output in PREPBUFR format J. Drake
  36 !         06/30/2006 - Updated for AIRS  retrievals      Syed  RH  Rizvi
  37 !         10/09/2006 - Updated for GPS Excess Phase    S.-Y. Chen
  38 !------------------------------------------------------------------------------
  39
  40    USE module_mm5
  41    USE module_map
  42    USE map_utils
  43    USE module_namelist
  44
  45    USE module_decoded
  46    USE module_type
  47    USE module_func
  48    USE module_icao
  49    USE module_sort
  50    USE module_duplicate
  51    USE module_per_type
  52    USE module_recoverp
  53    USE module_recoverh
  54    USE module_diagnostics
  55    USE module_qc
  56    USE module_err_ncep
  57    USE module_err_afwa
  58    USE module_complete
  59    USE module_thin_ob
  60    USE module_write
  61    use module_stntbl
  62 !------------------------------------------------------------------------------!
  63    IMPLICIT NONE
  64
  65    CHARACTER (LEN = 80)                 :: nml_filename
  66    CHARACTER (LEN = 80)                 :: title, caption
  67    LOGICAL                              :: exist
  68    INTEGER                              :: ii, fm_code, ntime
  69    INTEGER                              :: ins, jew
  70    INTEGER                              :: loop_index, number_of_obs
  71    INTEGER                              :: total_dups_loc, total_dups_time
  72    INTEGER                              :: total_dups
  73    INTEGER                              :: map_projection
  74    REAL                                 :: missing_flag
  75
  76    TYPE (report),      DIMENSION (:), ALLOCATABLE :: obs, obs_copy
  77    TYPE (measurement), POINTER                    :: next, current
  78    INTEGER,            DIMENSION (:), ALLOCATABLE :: index
  79
  80    real  :: lat, lon, xjd, yid, xxi, yyj, xxi1, yyj1
  81
  82    CHARACTER (LEN = 19)      :: time_min
  83    CHARACTER (LEN = 19)      :: time_fg
  84    CHARACTER (LEN = 19)      :: time_max
  85    CHARACTER (LEN = 31)      :: SLOT_TITLE
  86    LOGICAL                   :: dis0, dis1, dis2
  87 !------------------------------------------------------------------------------!
  88
  89
  90 ! 1.  READ THE NAMELIST
  91 ! ====================
  92
  93       nml_filename = 'namelist.obsproc'
  94
  95       WRITE (UNIT = 0, FMT = '(/,A,A,/,A)')        &
  96     ' READ NAMELIST FILE: ',TRIM  (nml_filename),  &
  97     ' ------------------'
  98
  99       CALL get_namelist (trim(nml_filename))
 100
 101
 102 ! 2.  SET THE MM5 MAP PROJECTION COEFICIENTS
 103 ! ==========================================
 104
 105 ! 2.1 Height at model lid
 106 !     -------------------
 107
 108       htop = 0.; h_tropo = 0.
 109 !   Tropopause height:
 110       h_tropo = Ref_height (pis0)
 111       htop = Ref_height (ptop)
 112
 113      height_max_icao = h_from_p_icao (0.0)
 114
 115 !      write(0,'(/5(A,f12.1)/)')  &
 116 !          "  htop=", htop, ", h_tropo=", h_tropo, &
 117 !          ", height_max_icao=", height_max_icao,  &
 118 !          ", ptop=", ptop, ", pis0=", pis0
 119
 120 ! 2.2 Map coeficients
 121 !     ---------------
 122
 123        CALL setup (domain_check_h, iproj, phic, xlonc, truelat1, truelat2, &
 124 #ifdef BKG
 125                    dis, xcntr, ycntr, xn, pole, psi1,  c2)
 126 #else
 127                    maxnes, nestix, nestjx, dis, numc, nesti, nestj, &
 128                    ixc, jxc, xcntr, ycntr, xn, pole, psi1,  c2, &
 129                    xim11, xjm11)
 130 #endif
 131
 132       if (fg_format == 'WRF') then
 133         lat =  map_info%lat1
 134         lon =  map_info%lon1
 135         call LLXY(lat, lon, xxi, yyj)
 136         write(0,'(/"         LLXY:lat, lon, (dot) xxi, yyj:",4f12.5 )') &
 137                                                lat, lon, xxi, yyj
 138         call latlon_to_ij(map_info, lat, lon, xxi, yyj)
 139         write(0,'( "latlon_to_ij:lat, lon, (cross)xxi, yyj:",4f12.5/)') &
 140                                                lat, lon, xxi, yyj
 141       endif
 142
 143 ! 2.3 Grid dimensions for the current domain
 144 !     --------------------------------------
 145
 146        ins = nestix (idd)
 147        jew = nestjx (idd)
 148
 149 ! 3.  LOAD GTS OBSERVATIONS
 150 ! =========================
 151
 152       icor = 0;
 153       caption (7*(icor+0)+1:7*(icor+1)) = "   READ"
 154       caption (7*(icor+1)+1:7*(icor+2)) = "  EMPTY"
 155       caption (7*(icor+2)+1:7*(icor+3)) = "OUTSIDE"
 156
 157 ! 3.1 Reset the number of ingested obs per type
 158 !     -----------------------------------------
 159
 160       nsynops = 0; nshipss = 0; nmetars = 0; npilots = 0; nsounds = 0;
 161       nsatems = 0; nsatobs = 0; naireps = 0; ngpspws = 0; namdars = 0;
 162       nssmt1s = 0; nssmt2s = 0; nssmis  = 0; ntovss  = 0; nqscats = 0;
 163       nothers = 0; nprofls = 0; nbuoyss = 0; ngpsztd = 0; ngpsref = 0;
 164       nboguss = 0; nairss  = 0; ngpseph = 0; ntamdar = 0;
 165
 166 ! 3.2 Reset the total number of ingested obs
 167 !     --------------------------------------
 168
 169       number_of_obs   = 0
 170       total_dups_time = 0
 171       total_dups_loc  = 0
 172
 173 ! 3.3 Allocate memory
 174 !     --------------
 175
 176       ALLOCATE (obs   (max_number_of_obs))
 177       ALLOCATE (index (max_number_of_obs))
 178
 179 ! 3.4 Read gts observations
 180 !     ---------------------
 181
 182       INQUIRE (FILE = obs_gts_filename, EXIST = exist )
 183
 184       IF (exist .and. LEN(TRIM(obs_gts_filename))>0) THEN
 185
 186          call read_msfc_table('msfc.tbl')
 187
 188       !  Read data from input file
 189
 190         CALL read_obs_gts (obs_gts_filename, obs, number_of_obs, &
 191           max_number_of_obs, fatal_if_exceed_max_obs, print_gts_read, &
 192           ins, jew, time_window_min, time_window_max,                 &
 193           map_projection, missing_flag)
 194
 195       !  Reset unused memory for subsequent reading
 196
 197         DO loop_index = number_of_obs+1, max_number_of_obs
 198            NULLIFY (obs (loop_index) % surface)
 199         ENDDO
 200
 201       ELSE
 202          WRITE (0,'(/,A,/)') "No decoded observation file to read."
 203          STOP
 204       ENDIF
 205
 206 ! 3.6 Check if any data have been loaded
 207 !     ----------------------------------
 208
 209       IF (number_of_obs .GT. 0) THEN
 210
 211       icor = icor + 3
 212
 213 ! 4.  OBSERVATION SORTING
 214 ! =======================
 215
 216       !  Sort the observations so that they can be merged easily.  Really
 217       !  what happens is that the index array is uniquely sorted by location
 218       !  (except for observations that are from the same "place").  This
 219       !  puts duplicate location observations next to each other.
 220       !  Then, merge the observations to (try to) remove all duplicates and
 221       !  build composite data.
 222       !  Then, sort obervations upon time, type and location
 223
 224 ! 4.1 Recover the missing pressure based on the observed heights under
 225 !     hydrostatic assumption or the model reference state.
 226 !     ---------------------------------------------------
 227       CALL recover_pressure_from_height (max_number_of_obs , &
 228                                          obs, number_of_obs, print_recoverp)
 229
 230       !hcl-note: is check_obs actually doing anything?
 231       CALL check_obs (max_number_of_obs, obs, number_of_obs, 'pressure')
 232
 233 ! 4.2 Sort station per location
 234 !     -------------------------
 235
 236       CALL sort_obs (obs ,number_of_obs , compare_loc, index )
 237
 238 ! 4.3 Merge duplicate stations (same location and same time)
 239 !     ------------------------------------------------------
 240 !     Because data merging is based on the pressure, the missing
 241 !     pressure have been recovered (Sec. 4.1)
 242
 243       caption (7*icor+1:7*(icor+1)) = "LOCDUPL"
 244
 245       CALL check_duplicate_loc (obs ,index ,number_of_obs, total_dups_loc, &
 246                                 time_analysis, print_duplicate_loc)
 247
 248       icor = icor + 1
 249
 250 ! 4.4 Remove duplicate stations (same location and different time)
 251 !     ------------------------------------------------------------
 252
 253       caption (7*icor+1:7*(icor+1)) = "TIMDUPL"
 254
 255       if (use_for /= '4DVAR' )  &
 256       CALL check_duplicate_time (obs ,index ,number_of_obs, total_dups_time, &
 257                                  time_analysis, print_duplicate_time)
 258
 259       icor = icor + 1
 260
 261 ! 4.5 Total Number of duplicate stations removed
 262 !     ------------------------------------------
 263
 264       total_dups = total_dups_loc + total_dups_time
 265
 266
 267 ! 4.7 Sort obs chronologically
 268 !     ------------------------
 269
 270       CALL sort_obs (obs ,number_of_obs , compare_tim, index )
 271
 272       ! Temporarily set back index to increasing order
 273       ! index = (/ ( loop_index , loop_index = 1 , max_number_of_obs ) /)
 274
 275 ! 5.  DIAGNOSTICS, ESTIMATE HEIGHT WHEN MISSING, OBSERVATIONAL ERROR
 276 ! ==================================================================
 277
 278 ! 5.1 Wind components, dew point relative humidity and mixing ratio
 279 !     -------------------------------------------------------------
 280 !     Gross error check on wind module (<0), direction (>360) and
 281 !     Temperature (>0), dew point (>0&<T) and RH (0< <100) are also performed
 282
 283       CALL derived_quantities (max_number_of_obs , obs , number_of_obs)
 284
 285 ! 5.2 Height
 286 !     ------
 287
 288       CALL recover_height_from_pressure (max_number_of_obs , &
 289                                          obs, number_of_obs, print_recoverh)
 290
 291       !hcl-note: is check_obs actually doing anything?
 292       CALL check_obs (max_number_of_obs, obs, number_of_obs, 'height')
 293
 294 ! 5.3 Observational error (pressure must be present)
 295 !     ----------------------------------------------
 296
 297       !  When a input file is provided, read it
 298
 299       INQUIRE (FILE = obs_err_filename, EXIST = exist )
 300
 301       IF (exist) THEN
 302
 303           CALL obs_err_afwa (obs_err_filename, &
 304                          max_number_of_obs, obs, number_of_obs)
 305
 306       ELSE
 307
 308          WRITE (0,'(/,A,/)') "No obs err input file "//trim(obs_err_filename)//" found."
 309          STOP
 310
 311          !starting from V3.7.1, disable the usage below
 312
 313          !Otherwise, use ncep value as default
 314          !write(0,'(/A/A/)') "!!! WARNING *** WARNING *** WARNING *** WARNING !!!", &
 315          !"<< NO AFWA OBS ERRORS FILE (obserr.txt) AVAILABE, USE NCEP ERROR AS DEFAULT >>"
 316          !CALL obs_err_ncep (max_number_of_obs, obs, number_of_obs)
 317
 318       ENDIF
 319
 320 ! 6.  CHECK COMPLETNESS AND QUALITY CONTROL PURELY BASED ON THE OBS DATA
 321 ! ======================================================================
 322
 323 ! 6.1 Vertical consistency check and dry convective adjustment QC
 324 !     -----------------------------------------------------------
 325
 326       CALL proc_qc1 (max_number_of_obs , obs , number_of_obs,            &
 327                      qc_test_vert_consistency , qc_test_convective_adj , &
 328                      print_qc_vert, print_qc_conv)
 329
 330 ! 6.2 Vertical domain check
 331 !     ---------------------
 332
 333       CALL proc_qc2 (max_number_of_obs , obs , number_of_obs, &
 334                      qc_test_above_lid, print_qc_lid)
 335
 336 ! 6.2 Check completness (Remove levels without data and too low or too high)
 337 !     ---------------------------------------------------------------------
 338
 339       caption (7*icor+1:7*(icor+1)) = "UNCOMPL"
 340
 341       CALL check_completness (max_number_of_obs , &
 342                               obs, number_of_obs, remove_above_lid, &
 343                               print_uncomplete)
 344
 345       icor = icor + 1
 346
 347 ! 6.3 Print statistics on obs
 348 !     -----------------------
 349
 350       caption (7*icor+1:7*(icor+1)+1) = "INGESTD"
 351       title = "INGESTED OBSERVATION AFTER CHECKS:"
 352
 353       CALL print_per_type (TRIM (title), TRIM (caption), icor)
 354
 355       icor = icor+1
 356
 357 ! 7. Reduce the QC flags and thin the satellite data
 358 ! ==================================================
 359
 360 ! 7.1 Reduce QC from 8 to 3 digits
 361 !     ----------------------------
 362
 363       CALL qc_reduction (max_number_of_obs, obs, number_of_obs)
 364
 365 ! 7.2 Thin the SATOB data before writing out
 366
 367       IF (Thining_SATOB) &
 368       CALL THIN_OB (max_number_of_obs, obs, number_of_obs, &
 369                    nestix(idd), nestjx(idd), missing_flag, &
 370                                         'SATOB', 88, 1000.0)
 371
 372       IF (Thining_SSMI) THEN
 373       CALL THIN_OB (max_number_of_obs, obs, number_of_obs, &
 374                    nestix(idd), nestjx(idd), missing_flag, &
 375                                         'SSMI_Rtvl', 125)
 376
 377       CALL THIN_OB (max_number_of_obs, obs, number_of_obs, &
 378                    nestix(idd), nestjx(idd), missing_flag, &
 379                                         'SSMI_Tb', 126)
 380       ENDIF
 381
 382       IF (Thining_QSCAT) &
 383       CALL THIN_OB (max_number_of_obs, obs, number_of_obs, &
 384                    nestix(idd), nestjx(idd), missing_flag, &
 385                                         'Qscatcat', 281)
 386
 387 ! 8.  OUTPUT DIRECT OBSERVATIONS (U,V,T,QV) IN ASCII and PREBUFR FILE
 388 ! ===================================================================
 389
 390 ! 8.1 Keep a copy of the obs
 391 !     ------------------------
 392
 393       allocate (obs_copy(1:max_number_of_obs))
 394       obs_copy = obs
 395
 396 ! 8.2 Loop over the time slots
 397 !     ------------------------
 398
 399       do ntime = 1, num_time_slots
 400
 401 ! 8.3 Set the time_min, time_fg, and time_max for the time slot: ntime
 402 !     ----------------------------------------------------------------
 403
 404         if ( num_time_slots > 1 ) then
 405
 406 ! 8.3.1 the number of time slots greater than 1 for FGAT and 4DVAR
 407 !       ..........................................................
 408
 409           if ( ntime == 1 .or. ntime == num_time_slots ) then
 410              idt = slot_len / 2 - 1
 411           else
 412              idt = slot_len - 1
 413           endif
 414
 415           if (ntime == 1) then
 416              time_min  = time_window_min
 417              time_fg  = time_window_min
 418           else
 419              call geth_newdate (time_min, time_max, 1)
 420              if (ntime == num_time_slots ) then
 421                  time_fg = time_window_max
 422              else
 423                  call geth_newdate (time_fg, time_min, slot_len/2)
 424              endif
 425           endif
 426           call geth_newdate (time_max, time_min, idt)
 427
 428         else
 429
 430 ! 8.3.2 the number of time slots equal to 1 for 3DVAR, FGAT and 4DVAR
 431 !       ..........................................................
 432
 433           time_min = time_window_min
 434           call  geth_newdate (time_max, time_window_max, -1)
 435           time_fg = time_analysis
 436
 437         endif
 438
 439 ! 8.3.3 Print the time slot information
 440 !       ...............................
 441
 442         write(0,'(//a,i2,4(2x,a))') "slot=",ntime, &
 443             "time_min, time_fg, time_max:", time_min, time_fg, time_max
 444
 445 ! 8.4 Get the original total obs back
 446 !     -------------------------------
 447
 448         deallocate (obs)
 449         allocate(obs(1:max_number_of_obs))
 450
 451         obs = obs_copy
 452
 453 ! 8.5 Loop over all the obs
 454 !     ---------------------
 455
 456         do ii = 1,number_of_obs
 457
 458 ! 8.6 determine the obs "discard" flag
 459 !     --------------------------------
 460
 461 ! 8.7.1 Initialize the "discard" switch" for a obs
 462 !       ...........................................
 463
 464         dis0 = obs(ii)%info%discard
 465         dis1 = .false.
 466         dis2 = .false.
 467
 468 ! 8.7.2 Kick out the specific type of obs based on namelist settings
 469 !       ............................................................
 470
 471         read(obs(ii)%info%platform(4:6),'(i3)') fm_code
 472         if ( (.not.write_synop .and. (fm_code==12  .or. fm_code==14)) .or. &
 473              (.not.write_ship  .and.  fm_code==13)                    .or. &
 474              (.not.write_metar .and. (fm_code==15  .or. fm_code==16)) .or. &
 475              (.not.write_buoy  .and. (fm_code==18  .or. fm_code==19)) .or. &
 476              (.not.write_pilot .and. (fm_code>=32 .and. fm_code<=34)) .or. &
 477              (.not.write_sound .and. (fm_code>=35 .and. fm_code<=38)) .or. &
 478              (.not.write_amdar .and.  fm_code==42)                    .or. &
 479              (.not.write_satem .and.  fm_code==86)                    .or. &
 480              (.not.write_satob .and.  fm_code==88)                    .or. &
 481              (.not.write_airep .and. (fm_code==96  .or. fm_code==97)) .or. &
 482              (.not.write_tamdar.and.  fm_code==101)                   .or. &
 483              (.not.write_gpspw .and.  fm_code==111)                   .or. &
 484              (.not.write_gpsztd.and.  fm_code==114)                   .or. &
 485              (.not.write_gpsref.and.  fm_code==116)                   .or. &
 486              (.not.write_gpseph.and.  fm_code==118)                   .or. &
 487              (.not.write_ssmt1 .and.  fm_code==121)                   .or. &
 488              (.not.write_ssmt2 .and.  fm_code==122)                   .or. &
 489              (.not.write_ssmi  .and. (fm_code==125 .or. fm_code==126)).or. &
 490              (.not.write_tovs  .and.  fm_code==131)                   .or. &
 491              (.not.write_qscat .and.  fm_code==281)                   .or. &
 492              (.not.write_profl .and.  fm_code==132)                   .or. &
 493              (.not.write_bogus .and.  fm_code==135)                   .or. &
 494              (.not.write_airs  .and.  fm_code==133)                   )    &
 495              dis1 = .true.
 496 ! 8.7.3 Kick out the obs outside the time slot
 497 !       ......................................
 498
 499          CALL inside_window (obs(ii)%valid_time%date_char, &
 500                               time_min, time_max, &
 501                               dis2)
 502
 503 ! 8.7.4 Reset the obs "discard" flag
 504 !       ............................
 505
 506          obs(ii)%info%discard = dis0 .or. dis1 .or.dis2
 507
 508 ! 8.8 end of the loop over obs
 509 !     ------------------------
 510
 511       enddo
 512
 513 ! 8.9 Time duplicate check within a slot for 4DVAR
 514 !     --------------------------------------------
 515
 516       if (use_for == '4DVAR' )  then
 517         CALL sort_obs (obs ,number_of_obs , compare_loc, index )
 518
 519         CALL check_duplicate_time (obs ,index ,number_of_obs, total_dups_time, &
 520                                    time_fg, print_duplicate_time)
 521       endif
 522
 523 ! 8.10 Print report per platform type and count the # levels per stations
 524 !      ------------------------------------------------------------------
 525
 526       write(SLOT_TITLE,'("OBSERVATIONS FOR OUTPUT SLOT ",I2.2)') ntime
 527       CALL sort_platform (max_number_of_obs, obs, number_of_obs, &
 528                           nsynops (icor), nshipss (icor), nmetars (icor), &
 529                           npilots (icor), nsounds (icor), nsatems (icor), &
 530                           nsatobs (icor), naireps (icor), ngpspws (icor), &
 531                           ngpsztd (icor), ngpsref (icor), ngpseph (icor), &
 532                           nssmt1s (icor), nssmt2s (icor), nssmis  (icor), &
 533                           ntovss  (icor), nothers (icor), namdars (icor), &
 534                           nqscats (icor), nprofls (icor), nbuoyss (icor), &
 535                           nboguss (icor), nairss  (icor), ntamdar (icor), &
 536                           SLOT_TITLE)
 537
 538 ! 8.11 Determine output type
 539 !      ----------------------------------------------
 540
 541       IF (output_ob_format .eq. 1 .or. output_ob_format .eq. 3) THEN
 542
 543 ! 8.11.1 Output observations in PREPBUFR
 544 !        -------------------------------
 545       CALL output_prep (max_number_of_obs, obs, number_of_obs, index, &
 546                           prepbufr_table_filename, &
 547                           prepbufr_output_filename, &
 548                           nsynops (icor), nshipss (icor), nmetars (icor), &
 549                           npilots (icor), nsounds (icor), nsatems (icor), &
 550                           nsatobs (icor), naireps (icor), ngpspws (icor), &
 551                           ngpsztd (icor), ngpsref (icor), ngpseph (icor), &
 552                           nssmt1s (icor), nssmt2s (icor), nssmis  (icor), &
 553                           ntovss  (icor), nothers (icor), namdars (icor), &
 554                           nqscats (icor), nprofls (icor), nbuoyss (icor), &
 555                           nboguss (icor), missing_flag, time_analysis)
 556
 557
 558       ENDIF
 559
 560       IF (output_ob_format .eq. 2 .or. output_ob_format .eq. 3) THEN
 561
 562 ! 8.11.2 Output observations in 3D-VAR Version 3.1 gts format
 563 !        ----------------------------------------------------
 564
 565       CALL output_gts_31 (max_number_of_obs, obs, number_of_obs, index, &
 566                           nsynops (icor), nshipss (icor), nmetars (icor), &
 567                           npilots (icor), nsounds (icor), nsatems (icor), &
 568                           nsatobs (icor), naireps (icor), ngpspws (icor), &
 569                           ngpsztd (icor), ngpsref (icor), ngpseph (icor), &
 570                           nssmt1s (icor), nssmt2s (icor), nssmis  (icor), &
 571                           ntovss  (icor), nothers (icor), namdars (icor), &
 572                           nqscats (icor), nprofls (icor), nbuoyss (icor), &
 573                           nboguss (icor), nairss  (icor), ntamdar(icor), missing_flag, time_fg)
 574
 575       CALL output_ssmi_31 (max_number_of_obs, obs, number_of_obs, index, &
 576                            nssmis  (icor), &
 577                            missing_flag, time_fg)
 578
 579       ENDIF
 580
 581
 582       enddo
 583
 584 ! 9.  END
 585 ! =======
 586       ENDIF
 587       STOP "99999"
 588
 589 END PROGRAM main_obsproc