var/da/da_obs_io/da_scan_obs_bufr.inc

   1 subroutine da_scan_obs_bufr (iv, filename)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Read BUFR observation file for input to wrfvar
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (iv_type),             intent(inout) :: iv
  10    character(len=*), optional, intent(in)    :: filename
  11
  12 #ifdef BUFR
  13
  14    logical                      :: match, end_of_file
  15    character(len=8)             :: subst2, csid, csid2
  16    integer                      :: idate2, nlevels2, lv1, lv2
  17    real                         :: hdr2(7), hdr_save(7), r8sid, r8sid2
  18    real                         :: pob1, pob2
  19    real                         :: obs2(8,255), obs_save(8,255)
  20    equivalence                     (r8sid, csid), (r8sid2, csid2)
  21    ! for thinning
  22    real                         :: tdiff                       ! DHR
  23    real                         :: dlat_earth,dlon_earth,crit
  24    integer                      :: itt,itx,iout
  25    logical                      :: iuse
  26
  27
  28    type (multi_level_type)      :: platform
  29    logical                      :: outside, outside_all, outside_time
  30
  31    character(len=40)     :: obstr,hdstr
  32    character(len=8)      :: subset
  33    character(len=14)     :: cdate, dmn, obs_date
  34    real                  :: hdr(7)
  35    real                  :: obs(8,255)
  36    real                  :: obs_time
  37    integer               :: iyear, imonth, iday, ihour, imin
  38
  39    integer               :: iost, ndup, n, i, j, k, surface_level, num_report
  40    integer               :: iret, idate, kx, nlevels, t29
  41    integer               :: iunit, fm, obs_index
  42    integer               :: num_outside_all, num_outside_time, num_thinned
  43
  44    if (trace_use) call da_trace_entry("da_scan_obs_bufr")
  45
  46    ! open file
  47    !  ---------
  48    call da_get_unit(iunit)
  49    if (present(filename)) then
  50       call closbf(iunit)
  51       open(unit   = iunit, FILE   = trim(filename), &
  52          iostat =  iost, form = 'unformatted', STATUS = 'OLD')
  53       if (iost /= 0) then
  54          write(unit=message(1),fmt='(A,I5,A)') &
  55             "Error",iost," opening PREPBUFR obs file "//trim(filename)
  56          call da_warning(__FILE__,__LINE__,message(1:1))
  57          call da_free_unit(iunit)
  58          if (trace_use) call da_trace_exit("da_scan_obs_bufr")
  59          return
  60       end if
  61    end if
  62
  63    hdstr='SID XOB YOB DHR TYP ELV T29     '
  64    obstr='POB QOB TOB ZOB UOB VOB PWO CAT '
  65
  66    !--------------------------------
  67    ! open bufr file then check date
  68    !--------------------------------
  69
  70    call openbf(iunit,'IN',iunit)
  71    call datelen(10)
  72    call readns(iunit,subset,idate,iret)  ! read in the next subset
  73    if ( iret /= 0 ) then
  74       write(unit=message(1),fmt='(A,I5,A)') &
  75          "Error",iret," reading PREPBUFR obs file "//trim(filename)
  76       call da_warning(__FILE__,__LINE__,message(1:1))
  77       call closbf(iunit)
  78       call da_free_unit(iunit)
  79       if (trace_use) call da_trace_exit("da_scan_obs_bufr")
  80       return
  81    end if
  82    write(unit=message(1),fmt='(a,i10)') 'BUFR file date is: ', idate
  83    call da_message(message(1:1))
  84
  85    num_report       = 0
  86    num_outside_all  = 0
  87    num_outside_time = 0
  88    num_thinned      = 0
  89
  90    match        = .false.
  91    end_of_file  = .false.
  92    outside_all  = .false.
  93    outside_time = .false.
  94
  95    reports: do while ( .not. end_of_file )
  96
  97       if ( match .or. outside_all .or. outside_time ) then
  98          call readns(iunit,subset,idate,iret)  ! read in the next subset
  99          if ( iret /= 0 ) then
 100             write(unit=message(1),fmt='(A,I3,A,I3)') &
 101                "return code from readns",iret,       &
 102                "reach the end of PREPBUFR obs unit",iunit
 103             !call da_warning(__FILE__,__LINE__,message(1:1))
 104             exit reports
 105          end if
 106       end if
 107
 108       num_report = num_report+1
 109
 110       call ufbint(iunit,hdr,7,1,iret,hdstr)
 111       call ufbint(iunit,obs,8,255,nlevels,obstr)
 112
 113       r8sid = hdr(1)
 114       platform % info % name(1:8) = subset
 115       platform % info % id        = csid(1:5)
 116       platform % info % dhr       = hdr(4)     ! difference in hour
 117       platform % info % elv       = hdr(6)
 118       platform % info % lon       = hdr(2)
 119       platform % info % lat       = hdr(3)
 120
 121       ! Put a check on Lon and Lat
 122       if ( platform%info%lon >= 180.0) platform%info%lon = platform%info%lon -360.0
 123       ! Fix funny wind direction at Poles
 124       !if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
 125       !   platform%info%lon = 0.0
 126       !end if
 127       platform%info%lat = max(platform%info%lat, -89.95)
 128       platform%info%lat = min(platform%info%lat,  89.95)
 129
 130       ! Restrict to a range of reports, useful for debugging
 131
 132       if (num_report < report_start) cycle reports
 133       if (num_report > report_end) exit reports
 134
 135       call da_llxy (platform%info, platform%loc, outside, outside_all)
 136
 137       if (outside_all) then
 138          num_outside_all = num_outside_all + 1
 139          if ( print_detail_obs ) then
 140             write(unit=stdout,fmt='(a,1x,a,2(1x,f8.3),a)')  &
 141                platform%info%name(1:8),platform%info%id(1:5), &
 142                platform%info%lat, platform%info%lon, '  -> outside_domain'
 143          end if
 144          cycle reports
 145       end if
 146
 147       ! check date
 148       write(cdate,'(i10)') idate
 149       write(dmn,'(i4,a1)') int(platform%info%dhr*60.0), 'm'
 150       call da_advance_time (cdate(1:10), trim(dmn), obs_date)
 151       if ( obs_date(13:14) /= '00' ) then
 152          write(0,*) 'wrong date: ', trim(cdate), trim(dmn), trim(obs_date)
 153          call da_error(__FILE__,__LINE__,(/"Wrong date"/))
 154       else
 155          read (obs_date(1:12),'(i4,4i2)') iyear, imonth, iday, ihour, imin
 156       end if
 157       call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 158       if (obs_time < time_slots(0) .or.  &
 159          obs_time >= time_slots(num_fgat_time)) then
 160          outside_time = .true.
 161          num_outside_time = num_outside_time + 1
 162          if ( print_detail_obs ) then
 163             write(unit=stdout,fmt='(a,1x,a,1x,a,a)')  &
 164                platform%info%name(1:8),platform%info%id(1:5), &
 165                trim(obs_date), '  -> outside_time'
 166          end if
 167          cycle reports
 168       else
 169          outside_time = .false.
 170       end if
 171
 172       t29 = int(0.1 + hdr(7))
 173       kx=int(0.1+hdr(5))
 174
 175       call readns(iunit,subst2,idate2,iret)
 176
 177       if ( iret /= 0 ) then
 178          end_of_file = .true.
 179       else
 180          match_check: do
 181             call ufbint(iunit,hdr2,7,1,iret,hdstr)
 182             ! check if this subset and the previous one are matching mass and wind
 183             match = .true.
 184             if ( subset /= subst2 ) then
 185                match = .false.
 186                exit match_check
 187             end if
 188             r8sid2 = hdr2(1)
 189             if ( csid /= csid2 ) then   ! check SID
 190                match = .false.
 191                exit match_check
 192             end if
 193             do i = 2, 4   ! check XOB, YOB, DHR
 194                if ( hdr(i) /= hdr2(i) ) then
 195                   match = .false.
 196                   exit match_check
 197                end if
 198             end do
 199             if ( hdr(6) /= hdr2(6) ) then   ! check ELV
 200                match = .false.
 201                exit match_check
 202             end if
 203             ! match found, merge two reports and find out the merged nlevel
 204             call ufbint(iunit,obs2,8,255,nlevels2,obstr)
 205             ! If this is a surface report, the wind subset precedes the
 206             ! mass subset - switch the subsets around in order to combine
 207             ! the surface pressure properly
 208             if ( kx == 280 .or. kx == 281 .or. kx == 284  .or.  &
 209                  kx == 287 .or. kx == 288  ) then
 210                obs_save = obs2
 211                obs2 = obs
 212                obs = obs_save
 213                hdr_save = hdr2
 214                hdr2 = hdr
 215                hdr = hdr_save
 216             end if
 217             lev_loop: do lv2 = 1, nlevels2
 218                do lv1 = 1, nlevels
 219                   pob1 = obs(1,lv1)
 220                   pob2 = obs2(1,lv2)
 221                   if ( pob1 == pob2 ) then
 222                      cycle lev_loop
 223                   else if ( (pob2 > pob1) .or. (lv1 .eq. nlevels) ) then
 224                      nlevels = nlevels + 1
 225                      cycle lev_loop
 226                   end if
 227                end do
 228             end do lev_loop
 229             exit match_check
 230          end do match_check
 231
 232          if ( .not. match ) then
 233             subset = subst2
 234             idate  = idate2
 235          end if
 236
 237       end if
 238
 239       ! skip some types
 240       !  61: Satellite soundings/retrievals/radiances
 241       !  66: SSM/I rain rate product
 242       !  72: NEXTRAD VAD winds
 243       if ( t29 == 61 .or. t29 == 66 .or. t29 == 72 ) cycle reports
 244
 245       if (nlevels > max_ob_levels) nlevels = max_ob_levels
 246       if ( nlevels < 1 ) then
 247          if ( (kx /= 164) .and. (kx /= 174) .and.   &
 248               (kx /= 165) .and. (kx /= 175) .and. (kx /= 74) ) cycle reports
 249       end if
 250
 251       tdiff = abs(platform%info%dhr-0.1)
 252       dlat_earth = platform%info%lat
 253       dlon_earth = platform%info%lon
 254       if (dlon_earth < 0.0) dlon_earth = dlon_earth + 360.0
 255       if (dlon_earth >= 360.0) dlon_earth = dlon_earth - 360.0
 256       dlat_earth = dlat_earth * deg2rad
 257       dlon_earth = dlon_earth * deg2rad
 258
 259       !---------------------------------------------------------------------------
 260       ! This is basically converting  rh to q i
 261       ! Method :
 262       !  if rh, temp and pr all available computes Qs otherwise sets Qs= missing
 263       !  if rh > 100 sets q = qs otherwise q = rh*Qs/100.0
 264       ! Note: Currently da_obs_proc_station is active only for ob_format_ascii
 265       !      call da_obs_proc_station(platform)
 266       !---------------------------------------------------------------------------
 267
 268       ! Loop over duplicating obs for global
 269       ndup = 1
 270       if (global .and. &
 271          (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
 272       if (test_transforms) ndup = 1
 273
 274       ! It is possible that logic for counting obs is incorrect for the
 275       ! global case with >1 MPI tasks due to obs duplication, halo, etc.
 276       ! TBH:  20050913
 277       dup_loop: do n = 1, ndup
 278          select case(t29)
 279          case (11, 12, 13, 22, 23, 31)
 280             select case (kx)
 281             case (120, 122, 132, 220, 222, 232) ;         ! Sound
 282                if (.not.use_soundobs) cycle reports
 283                if (n==1) iv%info(sound)%ntotal     = iv%info(sound)%ntotal + 1
 284                if (n==1) iv%info(sonde_sfc)%ntotal = iv%info(sonde_sfc)%ntotal + 1
 285                if (outside) cycle reports
 286                if ( thin_conv ) then
 287                   crit = tdiff
 288                   call map2grids_conv(sound,dlat_earth,dlon_earth,crit,iv%info(sound)%nlocal,itx,1,itt,iout,iuse)
 289                   call map2grids_conv(sonde_sfc,dlat_earth,dlon_earth,crit,iv%info(sonde_sfc)%nlocal,itx,1,itt,iout,iuse)
 290                   if ( .not. iuse ) then
 291                      num_thinned = num_thinned + 1
 292                      cycle reports
 293                   end if
 294                else
 295                   iv%info(sound)%nlocal     = iv%info(sound)%nlocal + 1
 296                   iv%info(sonde_sfc)%nlocal = iv%info(sound)%nlocal
 297                end if
 298                fm = 35
 299             case (221) ;                   ! Pilot
 300                if (.not.use_pilotobs) cycle reports
 301                if (n==1) iv%info(pilot)%ntotal = iv%info(pilot)%ntotal + 1
 302                if (outside) cycle reports
 303                if ( thin_conv ) then
 304                   crit = tdiff
 305                   call map2grids_conv(pilot,dlat_earth,dlon_earth,crit,iv%info(pilot)%nlocal,itx,1,itt,iout,iuse)
 306                   if ( .not. iuse ) then
 307                      num_thinned = num_thinned + 1
 308                      cycle reports
 309                   end if
 310                else
 311                   iv%info(pilot)%nlocal = iv%info(pilot)%nlocal + 1
 312                end if
 313                fm = 32
 314             case default
 315                exit dup_loop
 316             end select
 317
 318          case (41)
 319             ! case (130:131, 133, 230:231, 233) ; ! Airep
 320                if (.not.use_airepobs) cycle reports
 321                if (n==1) iv%info(airep)%ntotal = iv%info(airep)%ntotal + 1
 322                if (outside) cycle reports
 323                if ( thin_conv ) then
 324                   crit = tdiff
 325                   call map2grids_conv(airep,dlat_earth,dlon_earth,crit,iv%info(airep)%nlocal,itx,1,itt,iout,iuse)
 326                   if ( .not. iuse ) then
 327                      num_thinned = num_thinned + 1
 328                      cycle reports
 329                   end if
 330                else
 331                   iv%info(airep)%nlocal = iv%info(airep)%nlocal + 1
 332                end if
 333                fm = 42
 334
 335          case (522, 523);        ! Ships
 336                if (.not.use_shipsobs) cycle reports
 337                if (n==1) iv%info(ships)%ntotal = iv%info(ships)%ntotal + 1
 338                if (outside) cycle reports
 339                if ( thin_conv ) then
 340                   crit = tdiff
 341                   call map2grids_conv(ships,dlat_earth,dlon_earth,crit,iv%info(ships)%nlocal,itx,1,itt,iout,iuse)
 342                   if ( .not. iuse ) then
 343                      num_thinned = num_thinned + 1
 344                      cycle reports
 345                   end if
 346                else
 347                   iv%info(ships)%nlocal = iv%info(ships)%nlocal + 1
 348                end if
 349                fm = 13
 350
 351          case (531, 532, 561, 562) ;          ! Buoy
 352                if (.not.use_buoyobs) cycle reports
 353                if (n==1) iv%info(buoy)%ntotal = iv%info(buoy)%ntotal + 1
 354                if (outside) cycle reports
 355                if ( thin_conv ) then
 356                   crit = tdiff
 357                   call map2grids_conv(buoy,dlat_earth,dlon_earth,crit,iv%info(buoy)%nlocal,itx,1,itt,iout,iuse)
 358                   if ( .not. iuse ) then
 359                      num_thinned = num_thinned + 1
 360                      cycle reports
 361                   end if
 362                else
 363                   iv%info(buoy)%nlocal = iv%info(buoy)%nlocal + 1
 364                end if
 365                fm = 18
 366
 367          case (511, 514)
 368             ! case (181, 281) ;                   ! Synop
 369                if (.not.use_synopobs) cycle reports
 370                if (n==1) iv%info(synop)%ntotal = iv%info(synop)%ntotal + 1
 371                if (outside) cycle reports
 372                if ( thin_conv ) then
 373                   crit = tdiff
 374                   call map2grids_conv(synop,dlat_earth,dlon_earth,crit,iv%info(synop)%nlocal,itx,1,itt,iout,iuse)
 375                   if ( .not. iuse ) then
 376                      num_thinned = num_thinned + 1
 377                      cycle reports
 378                   end if
 379                else
 380                   iv%info(synop)%nlocal = iv%info(synop)%nlocal + 1
 381                end if
 382                fm = 12
 383
 384          case (512)
 385             ! case (187, 287) ;                        ! Metar
 386                if (.not.use_metarobs) cycle reports
 387                if (n==1) iv%info(metar)%ntotal = iv%info(metar)%ntotal + 1
 388                if (outside) cycle reports
 389                if ( thin_conv ) then
 390                   crit = tdiff
 391                   call map2grids_conv(metar,dlat_earth,dlon_earth,crit,iv%info(metar)%nlocal,itx,1,itt,iout,iuse)
 392                   if ( .not. iuse ) then
 393                      num_thinned = num_thinned + 1
 394                      cycle reports
 395                   end if
 396                else
 397                   iv%info(metar)%nlocal = iv%info(metar)%nlocal + 1
 398                end if
 399                fm = 15
 400
 401          case (63)
 402             ! case (242:246, 252:253, 255) ;         ! Geo. CMVs
 403                if (.not.use_geoamvobs) cycle reports
 404                if (n==1) iv%info(geoamv)%ntotal = iv%info(geoamv)%ntotal + 1
 405                if (outside) cycle reports
 406                if ( thin_conv ) then
 407                   crit = tdiff
 408                   call map2grids_conv(geoamv,dlat_earth,dlon_earth,crit,iv%info(geoamv)%nlocal,itx,1,itt,iout,iuse)
 409                   if ( .not. iuse ) then
 410                      num_thinned = num_thinned + 1
 411                      cycle reports
 412                   end if
 413                else
 414                   iv%info(geoamv)%nlocal = iv%info(geoamv)%nlocal + 1
 415                end if
 416                fm = 88
 417
 418          case (582, 583)      ! QuikSCAT 582 and WindSat 583
 419                if (.not.use_qscatobs) cycle reports
 420                if (n==1) iv%info(qscat)%ntotal = iv%info(qscat)%ntotal + 1
 421                if (outside) cycle reports
 422                if ( thin_conv ) then
 423                   crit = tdiff
 424                   call map2grids_conv(qscat,dlat_earth,dlon_earth,crit,iv%info(qscat)%nlocal,itx,1,itt,iout,iuse)
 425                   if ( .not. iuse ) then
 426                      num_thinned = num_thinned + 1
 427                      cycle reports
 428                   end if
 429                else
 430                   iv%info(qscat)%nlocal = iv%info(qscat)%nlocal + 1
 431                end if
 432                fm = 281
 433
 434          case (74)       ! GPS PW
 435                if (.not.use_gpspwobs) cycle reports
 436                if (n==1) iv%info(gpspw)%ntotal = iv%info(gpspw)%ntotal + 1
 437                if (outside) cycle reports
 438                if ( thin_conv ) then
 439                   crit = tdiff
 440                   call map2grids_conv(gpspw,dlat_earth,dlon_earth,crit,iv%info(gpspw)%nlocal,itx,1,itt,iout,iuse)
 441                   if ( .not. iuse ) then
 442                      num_thinned = num_thinned + 1
 443                      cycle reports
 444                   end if
 445                else
 446                   iv%info(gpspw)%nlocal = iv%info(gpspw)%nlocal + 1
 447                end if
 448                fm = 111
 449
 450          case (71, 73, 75, 76, 77)    ! Profiler
 451                if (.not.use_profilerobs) cycle reports
 452                if (n==1) iv%info(profiler)%ntotal = iv%info(profiler)%ntotal + 1
 453                if (outside) cycle reports
 454                if ( thin_conv ) then
 455                   crit = tdiff
 456                   call map2grids_conv(profiler,dlat_earth,dlon_earth,crit,iv%info(profiler)%nlocal,itx,1,itt,iout,iuse)
 457                   if ( .not. iuse ) then
 458                      num_thinned = num_thinned + 1
 459                      cycle reports
 460                   end if
 461                else
 462                   iv%info(profiler)%nlocal = iv%info(profiler)%nlocal + 1
 463                end if
 464                fm = 132
 465          case (571, 65)
 466               if (.not. use_ssmiretrievalobs) cycle reports
 467                if (n==1) iv%info(ssmi_rv)%ntotal = iv%info(ssmi_rv)%ntotal + 1
 468                if (outside) cycle reports
 469                if ( thin_conv ) then
 470                   crit = tdiff
 471                   call map2grids_conv(ssmi_rv,dlat_earth,dlon_earth,crit,iv%info(ssmi_rv)%nlocal,itx,1,itt,iout,iuse)
 472                   if ( .not. iuse ) then
 473                      num_thinned = num_thinned + 1
 474                      cycle reports
 475                   end if
 476                else
 477                   iv%info(ssmi_rv)%nlocal = iv%info(ssmi_rv)%nlocal + 1
 478                end if
 479                fm = 125      ! ssmi wind speed & tpw
 480          case default
 481             select case (kx)
 482             case (111 , 210)    ;         !  Tropical Cyclone Bogus
 483                ! Note Tropical cyclone Bougus is given type 135 in Obs-ascii
 484                if (.not.use_bogusobs) cycle reports
 485                if (n==1) iv%info(bogus)%ntotal = iv%info(bogus)%ntotal + 1
 486                if (outside) cycle reports
 487                if ( thin_conv ) then
 488                   crit = tdiff
 489                   call map2grids_conv(bogus,dlat_earth,dlon_earth,crit,iv%info(bogus)%nlocal,itx,1,itt,iout,iuse)
 490                   if ( .not. iuse ) then
 491                      num_thinned = num_thinned + 1
 492                      cycle reports
 493                   end if
 494                else
 495                   iv%info(bogus)%nlocal = iv%info(bogus)%nlocal + 1
 496                end if
 497                fm = 135
 498
 499             case default
 500                if ( print_detail_obs ) then
 501                   write(unit=message(1), fmt='(a, 2i12)') &
 502                      'unsaved obs found with kx & t29= ',kx,t29
 503                   call da_warning(__FILE__,__LINE__,message(1:1))
 504                end if
 505                exit dup_loop
 506             end select
 507          end select
 508          obs_index=fm_index(fm)
 509          iv%info(obs_index)%max_lev = max(iv%info(obs_index)%max_lev, nlevels)
 510       end do dup_loop
 511    end do reports
 512
 513    iv%info(synop)%max_lev     = 1
 514    iv%info(metar)%max_lev     = 1
 515    iv%info(ships)%max_lev     = 1
 516    iv%info(buoy)%max_lev      = 1
 517    iv%info(sonde_sfc)%max_lev = 1
 518
 519    write(unit=message(1),fmt='(A,4(1x,i7))') &
 520       'da_scan_obs_bufr: num_report, num_outside_all, num_outside_time, num_thinned: ', &
 521       num_report, num_outside_all, num_outside_time, num_thinned
 522    call da_message(message(1:1))
 523
 524    if ( thin_conv ) then
 525       do n = 1, num_ob_indexes
 526          call cleangrids_conv(n)
 527       end do
 528    end if
 529
 530    call closbf(iunit)
 531    close(iunit)
 532    call da_free_unit(iunit)
 533
 534    if (trace_use) call da_trace_exit("da_scan_obs_bufr")
 535 #else
 536    call da_error(__FILE__,__LINE__,(/"must compile with BUFR library"/))
 537 #endif
 538
 539 end subroutine da_scan_obs_bufr
 540