var/da/da_obs_io/da_read_obs_bufr.inc

   1 subroutine da_read_obs_bufr (iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Read BUFR observation file for input to wrfvar
   5    !---------------------------------------------------------------------------
   6    !   METHOD: use F90 sequential data structure to avoid read file twice
   7    !            so that da_scan_obs_bufr is not necessary any more.
   8    !            1. read prepbufr data in sequential data structure
   9    !            2. do gross QC check
  10    !            3. assign sequential data structure to innovation structure
  11    !                  and deallocate sequential data structure
  12    !
  13    !  HISTORY: 2009/10/13  - F90 sequential structure  Peng Xiu
  14    !
  15    !----------------------------------------------------------------------------
  16
  17    use da_define_structures, only: da_allocate_observations
  18
  19    implicit none
  20
  21    type (iv_type),             intent(inout) :: iv
  22
  23
  24 #ifdef BUFR
  25    character(len=10)    :: filename
  26    real, parameter    :: r8bfms = 9.0D08  ! BUFR missing value threshold
  27    logical                      :: match, end_of_file
  28    character(len=8)             :: subst2, csid, csid2
  29    integer              :: idate2, nlevels2, lv1, lv2
  30    real               :: hdr_save(7)
  31    real*8             :: hdr2(7), r8sid, r8sid2
  32    real               :: pob1, pob2
  33    real               :: temp(8,255)
  34    real               :: obs_save(8,255)
  35    real*8             :: obs2(8,255), qms2(8,255)
  36    real*8             :: oes2(8,255), pco2(8,255)
  37    real               :: pmo_save(2,1)
  38    real*8             :: pmo2(2,1)
  39    equivalence                     (r8sid, csid), (r8sid2, csid2)
  40    ! for thinning
  41    real               :: tdiff                       ! DHR
  42    real               :: dlat_earth,dlon_earth,crit
  43    integer              :: itt,itx,iout
  44    logical                      :: iuse
  45
  46    type (multi_level_type)      :: platform
  47    logical                      :: outside, outside_all, outside_time
  48    integer              :: ilocal(num_ob_indexes)
  49    integer              :: ntotal(num_ob_indexes)
  50    integer              :: nlocal(num_ob_indexes)
  51    integer              :: tp(num_ob_indexes)
  52    character(len=40)            :: obstr,hdstr,qmstr,oestr, pcstr
  53    character(len=8)             :: subset
  54    character(len=14)            :: cdate, dmn, obs_date,platform_name
  55    real*8             :: hdr(7)
  56    real*8             :: pmo(2,1)
  57    real*8             :: obs(8,255),qms(8,255),oes(8,255),pco(8,255)
  58    real               :: woe,toe,qoe,poe,pob,pwe
  59    real*8             :: obs_time
  60    integer              :: iyear, imonth, iday, ihour, imin
  61
  62    integer              :: iost, ndup, n, i, j, k, kk,surface_level, num_report, i1, i2
  63    integer              :: iret, idate, kx, old_nlevels,nlevels, t29,ifgat,ii
  64    integer              :: cat,zqm,pqm,qqm,tqm,wqm,pwq,pmq
  65    integer              :: tpc, wpc,iret2
  66    integer              :: iunit, fm , obs_index
  67
  68    integer              :: qflag           ! Flag for retaining data
  69
  70    real, allocatable  :: in(:), out(:)
  71    integer :: num_bufr(7)
  72    logical                      :: found
  73
  74    logical                      :: use_errtable
  75    integer              :: junit, itype, ivar
  76    real                 :: err_uv, err_t, err_p, err_q, err_pw, coef
  77    integer              :: ibufr
  78    integer              :: num_outside_all, num_outside_time, num_thinned,num_p,numbufr
  79
  80 !added by yw
  81    logical :: fexist, ywflag
  82
  83    type datalink_BUFR              !for PREPBUFR data reading
  84        type (multi_level_type_BUFR)          :: platform_BUFR
  85        integer                               :: fm_BUFR
  86        integer                               :: t29_BUFR
  87        integer                               :: ifgat_BUFR
  88        integer                               :: nlevels_BUFR
  89        integer                               :: kx_BUFR
  90        real                                  :: pco_BUFR(8,255)
  91        type(datalink_BUFR), pointer :: next
  92    end type datalink_BUFR
  93
  94    type(datalink_BUFR),pointer  :: head=>null(), plink=>null()
  95
  96    if (trace_use) call da_trace_entry("da_read_obs_bufr")
  97
  98 !  0.0 Initialize variables
  99 !--------------------------------------------------------------
 100    ilocal(:) = 0
 101    ntotal(:) = 0
 102    nlocal(:) = 0
 103
 104    err_uv = 10.0 ! m/s
 105    err_t  = 5.0  ! degree
 106    err_p  = 200  ! Pa
 107    err_q  = 10   ! RH percent
 108    err_pw = 0.2  ! cm
 109
 110    ! quality marker 0: Keep (always assimilate)
 111    !                1: Good
 112    !                2: Neutral or not checked
 113    !                3: Suspect
 114    if ( anal_type_verify ) then
 115       qflag = min(qmarker_retain,2)
 116    else
 117       qflag = qmarker_retain
 118    end if
 119    write(unit=message(1),fmt='(a,i1,a)') &
 120          "PREPBUFR ob with quality marker <= ", qflag, " will be retained."
 121    call da_message(message(1:1))
 122
 123    num_report       = 0
 124    num_outside_all  = 0
 125    num_outside_time = 0
 126    num_thinned      = 0
 127    num_p            = 0
 128    tp(:)            = 0
 129
 130
 131 ! 1.0  Open file
 132 !----------------------------------------------------------------
 133 !
 134 !check if input file exists
 135 num_bufr(:)=0
 136 numbufr=0
 137 if (num_fgat_time>1) then
 138   do i=1,7
 139      call da_get_unit(iunit)
 140      write(filename,fmt='(A,I1,A)') 'ob0',i,'.bufr'
 141      open(unit   = iunit, FILE   = trim(filename),iostat =  iost, form = 'unformatted', STATUS = 'OLD')
 142      if (iost == 0) then
 143         numbufr=numbufr+1
 144         num_bufr(numbufr)=i
 145      else
 146         close (iunit)
 147      end if
 148      call da_free_unit(iunit)
 149    end do
 150  else
 151    numbufr=1
 152  end if
 153
 154   if (numbufr==0) numbufr=1
 155
 156 !yw added by Wei Yu
 157   inquire (file="ob1.bufr", exist=fexist)
 158   ywflag = .false.
 159   if (fexist ) then
 160      numbufr = 2
 161      ywflag = .true.
 162   endif
 163 !yw end added
 164
 165
 166 bufrfile:  do ibufr=1,numbufr
 167      if (num_fgat_time==1) then
 168          filename='ob.bufr'
 169      else
 170          if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then
 171             filename='ob.bufr'
 172         else
 173             write(filename,fmt='(A,I1,A)') 'ob0',num_bufr(ibufr),'.bufr'
 174         end if
 175      end if
 176
 177 !yw added by Wei Yu
 178      if( (ibufr .eq. 2) .and. ywflag)  then
 179          filename='ob1.bufr'
 180      endif
 181 !yw end added
 182
 183 !
 184 !     We want to use specific unit number to read prepbufr data, which enables us to control its endianness
 185       iunit = 96
 186
 187       open(unit   = iunit, FILE   = trim(filename), &
 188          iostat =  iost, form = 'unformatted', STATUS = 'OLD')
 189       if (iost /= 0) then
 190          write(unit=message(1),fmt='(A,I5,A)') &
 191             "Error",iost," opening PREPBUFR obs file "//trim(filename)
 192          call da_warning(__FILE__,__LINE__,message(1:1))
 193          if ( num_fgat_time == 1 ) then
 194             call da_free_unit(iunit)
 195             if (trace_use) call da_trace_exit("da_read_obs_bufr")
 196             return
 197          else
 198             cycle bufrfile
 199          end if
 200       end if
 201
 202    hdstr='SID XOB YOB DHR TYP ELV T29'
 203    obstr='POB QOB TOB ZOB UOB VOB PWO CAT' ! observation
 204    qmstr='PQM QQM TQM ZQM WQM NUL PWQ NUL' ! quality marker
 205    oestr='POE QOE TOE NUL WOE NUL PWE NUL' ! observation error
 206    pcstr='PPC QPC TPC ZPC WPC NUL PWP NUL' ! program code
 207
 208    call openbf(iunit,'IN',iunit)
 209    call datelen(10)
 210
 211    call readns(iunit,subset,idate,iret)  ! read in the next subset
 212    if ( iret /= 0 ) then
 213       write(unit=message(1),fmt='(A,I5,A)') &
 214          "Error",iret," reading PREPBUFR obs file "//trim(filename)
 215       call da_warning(__FILE__,__LINE__,message(1:1))
 216       call closbf(iunit)
 217       if (trace_use) call da_trace_exit("da_read_obs_bufr")
 218       return
 219    end if
 220    !rewind(iunit)
 221
 222    write(unit=message(1),fmt='(a,i10)') 'BUFR file date is: ', idate
 223    call da_message(message(1:1))
 224
 225    ! oetab(300,33,6) processed in da_setup_obs_structures_bufr.inc
 226    use_errtable = allocated(oetab)
 227
 228 !  2.0 read data
 229 !  scan reports first
 230 !--------------------------------------------------------------
 231
 232    match        = .false.
 233    end_of_file  = .false.
 234    outside_all  = .false.
 235    outside_time = .false.
 236    reports: do while ( .not. end_of_file )
 237
 238       if ( match .or. outside_all .or. outside_time ) then
 239          call readns(iunit,subset,idate,iret)  ! read in the next subset
 240          if ( iret /= 0 ) then
 241             write(unit=message(1),fmt='(A,I3,A,I3)') &
 242                "return code from readns",iret,       &
 243                "reach the end of PREPBUFR obs unit",iunit
 244             !call da_warning(__FILE__,__LINE__,message(1:1))
 245             exit reports
 246          end if
 247       end if
 248
 249       num_report = num_report+1
 250
 251       call ufbint(iunit,hdr,7,1,iret2,hdstr)
 252       call ufbint(iunit,pmo,2,1,nlevels,'PMO PMQ')
 253       call ufbint(iunit,qms,8,255,nlevels,qmstr)
 254       call ufbint(iunit,oes,8,255,nlevels,oestr)
 255       call ufbint(iunit,pco,8,255,nlevels,pcstr)
 256       call ufbint(iunit,obs,8,255,nlevels,obstr)
 257
 258       r8sid = hdr(1)
 259       platform % info % name(1:8)  = subset
 260       platform % info % name(9:40) = '                                '
 261       platform % info % id(1:5)    = csid(1:5)
 262       platform % info % id(6:40)   = '                                   '
 263       platform % info % dhr        = hdr(4)    ! difference in hour
 264       platform % info % elv        = hdr(6)
 265       platform % info % lon        = hdr(2)
 266       platform % info % lat        = hdr(3)
 267
 268       ! blacklisted stations should be handled through an external table.
 269       ! For now, temporary fix is implemented here for known incorrect
 270       ! station info in NCEP PREPBUFR file
 271       if ( trim(platform%info%id) == 'BGQQ' ) then
 272          platform%info%elv = 19
 273          platform%info%lon = -69.21
 274          platform%info%lat = 77.46
 275       end if
 276       if ( trim(platform%info%id) == 'UWKE' ) then
 277          platform%info%elv = 194
 278          platform%info%lon = 52.09
 279          platform%info%lat = 55.56
 280       end if
 281
 282       ! Put a check on Lon and Lat
 283       if ( platform%info%lon >= 180.0 ) platform%info%lon = platform%info%lon - 360.0
 284       ! Fix funny wind direction at Poles
 285       !if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
 286       !   platform%info%lon = 0.0
 287       !end if
 288       platform%info%lat = max(platform%info%lat, -89.95)
 289       platform%info%lat = min(platform%info%lat,  89.95)
 290
 291       ! Restrict to a range of reports, useful for debugging
 292
 293       if (num_report < report_start) cycle reports
 294       if (num_report > report_end)  exit reports
 295
 296       call da_llxy (platform%info, platform%loc,outside, outside_all)
 297
 298       if (outside_all) then
 299          num_outside_all = num_outside_all + 1
 300          if ( print_detail_obs ) then
 301             write(unit=stderr,fmt='(a,1x,a,2(1x,f8.3),a)')  &
 302                platform%info%name(1:8),platform%info%id(1:5), &
 303                platform%info%lat, platform%info%lon, '  -> outside_domain'
 304          end if
 305          cycle reports
 306       end if
 307
 308       ! check date
 309       write(cdate,'(i10)') idate
 310       write(dmn,'(i4,a1)') int(platform%info%dhr*60.0), 'm'
 311       call da_advance_time (cdate(1:10), trim(dmn), obs_date)
 312       if ( obs_date(13:14) /= '00' ) then
 313          write(0,*) 'wrong date: ', trim(cdate), trim(dmn), trim(obs_date)
 314          call da_error(__FILE__,__LINE__,(/"Wrong date"/))
 315       else
 316          read (obs_date(1:12),'(i4,4i2)') iyear, imonth, iday, ihour, imin
 317       end if
 318       call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 319       if (obs_time < time_slots(0) .or.  &
 320          obs_time >= time_slots(num_fgat_time)) then
 321          outside_time = .true.
 322          num_outside_time = num_outside_time + 1
 323          if ( print_detail_obs ) then
 324             write(unit=stderr,fmt='(a,1x,a,1x,a,a)')  &
 325                platform%info%name(1:8),platform%info%id(1:5), &
 326                trim(obs_date), '  -> outside_time'
 327          end if
 328          cycle reports
 329       else
 330          outside_time = .false.
 331       end if
 332
 333 !--------  determine FGAT index ifgat
 334
 335          do ifgat=1,num_fgat_time
 336             if (obs_time >= time_slots(ifgat-1) .and.  &
 337                 obs_time  < time_slots(ifgat)) exit
 338          end do
 339
 340       write(unit=platform%info%date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 341          iyear, '-', imonth, '-', iday, '_', ihour, ':', imin, ':', 0
 342
 343       if (print_detail_obs) then
 344          ! Simplistic approach, could be improved to get it all done on PE 0
 345          if (.NOT. outside) then
 346             write(unit=stdout,fmt='(A,1X,I8,1X,A,F8.2,A,F8.2,A,1X,A,I3,1X,A,2F8.3)') &
 347                "Report",num_report,"at",platform%info%lon,"E",platform%info%lat,"N", &
 348                "on processor", myproc,"position", platform%loc%x,platform%loc%y
 349          end if
 350       end if
 351
 352       t29 = int(0.1 + hdr(7))
 353       kx=int(0.1+hdr(5))
 354
 355       if ( use_errtable ) then
 356          do k = 1, nlevels
 357             pob = obs(1,k)
 358             do lv1 = 1, 32
 359                if ( pob >= oetab(kx,lv1+1,1) .and. pob <= oetab(kx,lv1,1) ) then
 360                   coef = (pob-oetab(kx,lv1,1))/(oetab(kx,lv1,1)-oetab(kx,lv1+1,1))
 361                   oes(1,k) = (1.0+coef)*oetab(kx,lv1,5)-coef*oetab(kx,lv1+1,5) !p
 362                   oes(2,k) = (1.0+coef)*oetab(kx,lv1,3)-coef*oetab(kx,lv1+1,3) !q
 363                   oes(3,k) = (1.0+coef)*oetab(kx,lv1,2)-coef*oetab(kx,lv1+1,2) !t
 364                   oes(5,k) = (1.0+coef)*oetab(kx,lv1,4)-coef*oetab(kx,lv1+1,4) !uv
 365                   oes(7,k) = (1.0+coef)*oetab(kx,lv1,6)-coef*oetab(kx,lv1+1,6) !pw
 366                   exit
 367                end if
 368             end do
 369          end do
 370       end if
 371
 372       call readns(iunit,subst2,idate2,iret)
 373
 374       if ( iret /= 0 ) then
 375          end_of_file = .true.
 376       else
 377          match_check: do
 378             call ufbint(iunit,hdr2,7,1,iret2,hdstr)
 379             ! check if this subset and the previous one are matching mass and wind
 380             match = .true.
 381             if ( subset /= subst2 ) then
 382                match = .false.
 383                exit match_check
 384             end if
 385             r8sid2 = hdr2(1)
 386             if ( csid /= csid2 ) then  ! check SID
 387                match = .false.
 388                exit match_check
 389             end if
 390             do i = 2, 4   ! check XOB, YOB, DHR
 391                if ( hdr(i) /= hdr2(i) ) then
 392                   match = .false.
 393                   exit match_check
 394                end if
 395             end do
 396             if ( hdr(6) /= hdr2(6) ) then   ! check ELV
 397                match = .false.
 398                exit match_check
 399             end if
 400             !The two headers match, now read data from the second subset
 401             call ufbint(iunit,pmo2,2,1,nlevels2,'PMO PMQ')
 402             call ufbint(iunit,qms2,8,255,nlevels2,qmstr)
 403             call ufbint(iunit,oes2,8,255,nlevels2,oestr)
 404             call ufbint(iunit,pco2,8,255,nlevels2,pcstr)
 405             call ufbint(iunit,obs2,8,255,nlevels2,obstr)
 406
 407             if ( use_errtable ) then
 408                kx = nint(hdr2(5))
 409                do k = 1, nlevels2
 410                   pob = obs2(1,k)
 411                   do lv1 = 1, 32
 412                      if ( pob >= oetab(kx,lv1+1,1) .and. pob <= oetab(kx,lv1,1) ) then
 413                         coef = (pob-oetab(kx,lv1,1))/(oetab(kx,lv1,1)-oetab(kx,lv1+1,1))
 414                         oes2(1,k) = (1.0+coef)*oetab(kx,lv1,5)-coef*oetab(kx,lv1+1,5) !p
 415                         oes2(2,k) = (1.0+coef)*oetab(kx,lv1,3)-coef*oetab(kx,lv1+1,3) !q
 416                         oes2(3,k) = (1.0+coef)*oetab(kx,lv1,2)-coef*oetab(kx,lv1+1,2) !t
 417                         oes2(5,k) = (1.0+coef)*oetab(kx,lv1,4)-coef*oetab(kx,lv1+1,4) !uv
 418                         oes2(7,k) = (1.0+coef)*oetab(kx,lv1,6)-coef*oetab(kx,lv1+1,6) !pw
 419                         exit
 420                      end if
 421                   end do
 422                end do
 423             end if
 424
 425             ! If this is a surface report, the wind subset precedes the
 426             ! mass subset - switch the subsets around in order to combine
 427             ! the surface pressure properly
 428             kx = nint(hdr(5))
 429             if ( kx == 280 .or. kx == 281 .or. kx == 284  .or.  &
 430                  kx == 287 .or. kx == 288  ) then
 431                pmo_save = pmo2
 432                pmo2 = pmo
 433                pmo = pmo_save
 434                temp = obs2
 435                obs2 = obs
 436                obs = temp
 437                hdr_save = hdr2
 438                hdr2 = hdr
 439                hdr = hdr_save
 440                temp = qms2
 441                qms2 = qms
 442                qms = temp
 443                temp = oes2
 444                oes2 = oes
 445                oes = temp
 446                temp = pco2
 447                pco2 = pco
 448                pco = temp
 449             end if
 450
 451             ! combine the two matching subsets
 452             do i = 1, 2
 453                if ( pmo(i,1) > r8bfms ) then
 454                   pmo(i,1) = pmo2(i,1)
 455                end if
 456             end do
 457             lev_loop: do lv2 = 1, nlevels2
 458                do lv1 = 1, nlevels
 459                   pob1 = obs(1,lv1)
 460                   pob2 = obs2(1,lv2)
 461                   if ( pob1 == pob2 ) then
 462                      do i = 1, 7   ! skip the CAT
 463                         if ( obs(i,lv1) > r8bfms ) then
 464                            obs(i,lv1) = obs2(i,lv2)
 465                            if ( obs2(i,lv2) <= r8bfms ) then
 466                               obs(8,lv1) = obs2(8,lv2)  ! rewrite CAT
 467                            end if
 468                         end if
 469                         if ( oes(i,lv1) > r8bfms ) then
 470                            oes(i,lv1) = oes2(i,lv2)
 471                         end if
 472                         if ( pco(i,lv1) > r8bfms ) then
 473                            pco(i,lv1) = pco2(i,lv2)
 474                         end if
 475                      end do
 476                      do i = 1, 8
 477                         if ( qms(i,lv1) > r8bfms ) then
 478                            qms(i,lv1) = qms2(i,lv2)
 479                         end if
 480                      end do
 481                      cycle lev_loop
 482                   else if ( (pob2 > pob1) .or. (lv1 .eq. nlevels) ) then
 483                      nlevels = nlevels + 1
 484                      obs(:,nlevels) = obs2(:,lv2)
 485                      qms(:,nlevels) = qms2(:,lv2)
 486                      oes(:,nlevels) = oes2(:,lv2)
 487                      pco(:,nlevels) = pco2(:,lv2)
 488                      cycle lev_loop
 489                   end if
 490                end do
 491             end do lev_loop
 492             ! sort combined report in descending pressure order
 493             do i1 = 1, nlevels-1
 494                do i2 = i1+1, nlevels
 495                   if ( obs(1,i2) .gt. obs(1,i1) ) then
 496                      temp(:,1) = obs(:,i1)
 497                      obs(:,i1) = obs(:,i2)
 498                      obs(:,i2) = temp(:,1)
 499                      temp(:,1) = qms(:,i1)
 500                      qms(:,i1) = qms(:,i2)
 501                      qms(:,i2) = temp(:,1)
 502                      temp(:,1) = oes(:,i1)
 503                      oes(:,i1) = oes(:,i2)
 504                      oes(:,i2) = temp(:,1)
 505                      temp(:,1) = pco(:,i1)
 506                      pco(:,i1) = pco(:,i2)
 507                      pco(:,i2) = temp(:,1)
 508                   end if
 509                end do
 510             end do
 511             exit match_check
 512          end do match_check
 513
 514          if ( .not. match ) then
 515             subset = subst2
 516             idate = idate2
 517          end if
 518
 519       end if
 520
 521       ! skip some types
 522       !  61: Satellite soundings/retrievals/radiances
 523       !  66: SSM/I rain rate product
 524       !  72: NEXTRAD VAD winds
 525       if ( t29 == 61 .or. t29 == 66 .or. t29 == 72 ) cycle reports
 526
 527       platform % info % levels    = nlevels
 528
 529       platform % loc % slp %inv   = missing_r
 530       platform % loc % slp %qc    = missing_data
 531       platform % loc % slp %error = err_p
 532       platform % loc % pw %inv    = missing_r
 533       platform % loc % pw %qc     = missing_data
 534       platform % loc % pw %error  = err_pw
 535       if ( pmo(1,1) < r8bfms ) then
 536          pmq=nint(pmo(2,1))
 537          if ( pmq <= qflag .and. pmq >= 0 ) then
 538             platform % loc % slp % inv =pmo(1,1)*100.0
 539             platform % loc % slp % qc  =pmq
 540             platform % loc % slp % error = err_p          ! hardwired
 541          end if
 542       end if
 543       if ( obs(7,1) < r8bfms ) then
 544          pwq=nint(qms(7,1))
 545 #if ( RWORDSIZE == 4 )
 546          pwe = min(DBLE(err_pw), oes(7,1))
 547 #else
 548          pwe = min(err_pw, oes(7,1))
 549 #endif
 550          if ( pwq <= qflag .and. pwq >= 0 ) then
 551             platform % loc % pw % inv = obs(7,1) * 0.1    ! convert to cm
 552             platform % loc % pw % qc  =  pwq
 553             platform % loc % pw % error = pwe          ! hardwired
 554          end if
 555       end if
 556       !$OMP PARALLEL DO &
 557       !$OMP PRIVATE (i)
 558       do i=1,max_ob_levels
 559          platform % each (i) % height  = missing_r
 560          platform % each (i) % height_qc = missing_data
 561
 562          platform % each (i) % zk = missing_r
 563
 564          platform % each (i) % u % inv = missing_r
 565          platform % each (i) % u % qc  = missing_data
 566          platform % each (i) % u % error = err_uv
 567
 568          platform % each (i) % v = platform % each (i) % u
 569
 570          platform % each (i) % t % inv = missing_r
 571          platform % each (i) % t % qc  = missing_data
 572          platform % each (i) % t % error = err_t
 573
 574          platform % each (i) % p % inv = missing_r
 575          platform % each (i) % p % qc  = missing_data
 576          platform % each (i) % p % error = err_p
 577
 578          platform % each (i) % q % inv = missing_r
 579          platform % each (i) % q % qc  = missing_data
 580          platform % each (i) % q % error = err_q
 581       end do
 582       !$OMP END PARALLEL DO
 583
 584       !!$OMP PARALLEL DO &
 585       !!$OMP PRIVATE (k, tpc, wpc, pqm, qqm, tqm, wqm, zqm, cat, toe, poe, qoe, woe)
 586       do k = 1, platform % info % levels
 587          ! set t units to Kelvin
 588          if (obs(3,k) > -200.0 .and. obs(3,k) < 300.0) then
 589             obs(3,k) = obs(3,k) + t_kelvin
 590          end if
 591
 592          ! scale q and compute t from tv, if they aren't missing
 593          if (obs(2,k) > 0.0 .and. obs(2,k) < r8bfms) then
 594             obs(2,k) = obs(2,k)*1e-6
 595             tpc = nint(pco(3,k))
 596             if (obs(3,k) > -200.0 .and. obs(3,k) < 350.0) then
 597                if ( tpc >= 8 ) then   ! program code 008 VIRTMP
 598                   ! 0.61 is used in NCEP prepdata.f to convert T to Tv
 599                   obs(3,k) = obs(3,k) / (1.0 + 0.61 * obs(2,k))
 600                end if
 601             end if
 602          end if
 603
 604          !Not currently used
 605          !cat=nint(obs(8,k))
 606
 607 #if ( RWORDSIZE == 4 )
 608          toe = min(DBLE(err_t), oes(3,k))
 609          woe = min(DBLE(err_uv), oes(5,k))
 610          qoe = min(DBLE(err_q), oes(2,k)*10.0)  ! convert to % from PREPBUFR percent divided by 10
 611          poe = min(DBLE(err_p), oes(1,k)*100.0) ! convert to Pa
 612 #else
 613          toe = min(err_t, oes(3,k))
 614          woe = min(err_uv, oes(5,k))
 615          qoe = min(err_q, oes(2,k)*10.0)  ! convert to % from PREPBUFR percent divided by 10
 616          poe = min(err_p, oes(1,k)*100.0) ! convert to Pa
 617 #endif
 618
 619          if ( obs(3,k) < r8bfms ) then
 620             tqm=nint(qms(3,k))
 621             if ( tqm <= qflag .and. tqm >= 0 ) then
 622                platform % each (k) % t % inv =obs(3,k)
 623                platform % each (k) % t % qc  =tqm
 624                platform % each (k) % t % error =toe
 625             end if
 626          end if
 627
 628          if (obs(5,k) < r8bfms .and. obs(6,k) < r8bfms ) then
 629             wqm=nint(qms(5,k))
 630             wpc = nint(pco(5,k))
 631             if ( wqm <= qflag .and. wqm >= 0 ) then
 632                platform % each (k) % u % inv =obs(5,k)
 633                platform % each (k) % v % inv =obs(6,k)
 634                platform % each (k) % u % qc  =wqm
 635                platform % each (k) % u % error =woe
 636                platform % each (k) % v % qc  =wqm
 637                platform % each (k) % v % error =woe
 638
 639                ! Convert earth wind to model wind.
 640                ! note on SPSSMI wind: only wspd available (stored in VOB)
 641                ! and direction is initially set to to zero and stored in UOB
 642                ! in wpc = 1 stage. u and v are computed in program wpc = 10 (OIQC).
 643                if ( kx /= 283 .or. ( kx == 283 .and. wpc == 10 ) ) then
 644                   call da_earth_2_model_wind(obs(5,k), obs(6,k), &
 645                      platform % each (k) % u % inv, &
 646                      platform % each (k) % v % inv, &
 647                      platform%info%lon)
 648                end if
 649                if (platform % each (k) % u % inv == 0.0 .and. platform % each (k) % v % inv == 0.0) then
 650                   platform % each (k) % u % inv = missing_r
 651                   platform % each (k) % v % inv = missing_r
 652                   platform % each (k) % u % qc  = missing_data
 653                   platform % each (k) % v % qc  = missing_data
 654                end if
 655             end if
 656          end if
 657          if ( obs(2,k)>0.0 .and. obs(2,k)<r8bfms ) then
 658             qqm=nint(qms(2,k))
 659             if ( qqm<=qflag .and. qqm>=0 ) then
 660                platform % each (k) % q % inv =obs(2,k)
 661                if ( obs(1,k) >= 300.0 .and. obs(1,k) < r8bfms ) then   ! do not use moisture above 300 hPa
 662                   platform % each (k) % q % qc  =qqm
 663                end if
 664                platform % each (k) % q % error = qoe
 665             end if
 666          end if
 667
 668          if ( obs(4,k) < r8bfms )then
 669             zqm=nint(qms(4,k))
 670             if ( zqm <= qflag .and. zqm >= 0 ) then
 671                platform % each (k) % height  = obs(4,k)
 672                platform % each (k) % height_qc =zqm
 673             end if
 674          end if
 675
 676          if ( obs(1,k) > 0.0 .and. obs(1,k) < r8bfms ) then
 677             pqm=nint(qms(1,k))
 678             if ( pqm <= qflag .and. pqm >= 0 ) then
 679                platform % each (k) % p % inv =obs(1,k)*100.0  ! convert to Pa
 680                platform % each (k) % p % qc  =pqm
 681                platform % each (k) % p % error =poe
 682             end if
 683          end if
 684       end do
 685       !!$OMP END PARALLEL DO
 686
 687       ! assign u,v,t,q obs errors for synop and metar
 688       if ( t29 == 512 .or. t29 == 511 .or. t29 == 514 ) then
 689 #if ( RWORDSIZE == 4 )
 690          toe = min(DBLE(err_t), oes(3,1))
 691          woe = min(DBLE(err_uv), oes(5,1))
 692          qoe = min(DBLE(err_q), oes(2,1)*10.0)  ! convert to % from PREPBUFR percent divided by 10
 693 #else
 694          toe = min(err_t, oes(3,1))
 695          woe = min(err_uv, oes(5,1))
 696          qoe = min(err_q, oes(2,1)*10.0)  ! convert to % from PREPBUFR percent divided by 10
 697 #endif
 698          if (obs(5,1) < r8bfms .and. obs(6,1) < r8bfms ) then
 699             wqm=nint(qms(5,1))
 700             if ( wqm == 8 .or. wqm  == 9 .or. wqm  == 15) then
 701                platform%each(1)%u%qc  = 88
 702                platform%each(1)%v%qc  = 88
 703                if ( use_errtable ) then
 704                   platform%each(1)%u%error = woe
 705                   platform%each(1)%v%error = woe
 706                else
 707                   platform%each(1)%u%error = 1.1
 708                   platform%each(1)%v%error = 1.1
 709                end if
 710                ! Convert earth wind to model wind.
 711                call da_earth_2_model_wind(obs(5,1), obs(6,1),     &
 712                   platform%each(1)%u%inv, platform%each(1)%v%inv, &
 713                   platform%info%lon)
 714                if ( platform%each(1)%u%inv == 0.0 .and. platform%each(1)%v%inv == 0.0 ) then
 715                   platform%each(1)%u%inv = missing_r
 716                   platform%each(1)%v%inv = missing_r
 717                   platform%each(1)%u%qc  = missing_data
 718                   platform%each(1)%v%qc  = missing_data
 719                end if
 720             end if
 721          end if
 722          if ( obs(3,1) < r8bfms ) then
 723             tqm=nint(qms(3,1))
 724             if ( tqm == 8 .or. tqm == 9 .or. tqm == 15 ) then
 725                platform%each(1)%t%inv = obs(3,1)
 726                platform%each(1)%t%qc  = 88
 727                if ( use_errtable ) then
 728                   platform%each(1)%t%error = toe
 729                else
 730                   platform%each(1)%t%error = 2.0
 731                end if
 732             end if
 733          end if
 734          if ( obs(2,1) > 0.0 .and. obs(2,1) < r8bfms ) then
 735             qqm=nint(qms(2,1))
 736             if ( qqm == 8 .or. qqm == 9 .or. qqm == 15 ) then
 737                platform%each(1)%q%inv = obs(2,1)
 738                platform%each(1)%q%qc  = 88
 739                if ( use_errtable ) then
 740                   platform%each(1)%q%error = qoe  ! RH percent
 741                else
 742                   platform%each(1)%q%error = 10   ! RH percent
 743                end if
 744             end if
 745          end if
 746       end if
 747       ! assign tpw obs errors for gpspw (t29=74) and SPSSMI retrieved TPW (t29=65)
 748       if ( obs(7,1) > 0.0 .and. obs(7,1) < r8bfms ) then
 749          if ( t29 == 74 .or. t29 == 65 ) then
 750             if ( pwq == 8 .or. pwq  == 9 .or. pwq  == 15) then
 751                platform%loc%pw%inv   = obs(7,1) * 0.1    ! convert to cm
 752                platform%loc%pw%qc    = 88
 753                if ( use_errtable ) then
 754                   platform%loc%pw%error = pwe
 755                else
 756                   platform%loc%pw%error = 0.2     ! hardwired to 0.2 cm
 757                end if
 758             end if
 759          end if
 760       end if
 761
 762       nlevels    = platform%info%levels
 763
 764       if (nlevels > max_ob_levels) then
 765          nlevels = max_ob_levels
 766
 767          write(unit=stderr, fmt='(/a/)') 'Warning: Too many levels.'
 768
 769          write(unit=stderr, fmt='(/2a/2a/2x,a,2f8.2,2(2x,a,f9.2)/2(2x,a,i4)/)') &
 770             'Subset:   ', platform%info%name(1:8), &
 771             'Platfrom: ', trim(platform%info%platform), &
 772             'Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
 773             'elv:   ', platform%info%elv, &
 774             'pstar: ', platform%info%pstar, &
 775             'level: ', platform%info%levels, &
 776             'kx:    ', kx
 777       else if ( nlevels < 1 ) then
 778          if ( (kx /= 164) .and. (kx /= 174) .and.                   &
 779               (kx /= 165) .and. (kx /= 175) .and. (kx /= 74) ) then
 780             write(unit=stderr, fmt='(/a/)') &
 781                'Warning: Too few levels.'
 782
 783             write(unit=stderr, fmt='(/2a/2a/2x,a,2f8.2,2(2x,a,f9.2)/3(2x,a,i4)/)') &
 784                'Subset:   ', platform%info%name(1:8), &
 785                'Platfrom: ', trim(platform%info%platform), &
 786                'Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
 787                'elv:   ', platform%info%elv, &
 788                'pstar: ', platform%info%pstar, &
 789                'level: ', platform%info%levels, &
 790                'kx:    ', kx, &
 791                't29:   ', t29
 792
 793             cycle reports
 794          end if
 795       end if
 796
 797      if (num_fgat_time > 1 ) then !for 4dvar
 798
 799       if (.not. associated(head)) then
 800          nullify ( head )
 801          allocate ( head )
 802          allocate ( head%platform_BUFR%each(1:nlevels) )
 803          nullify ( head%next )
 804          plink => head
 805       else
 806          allocate ( plink%next )
 807          plink => plink%next
 808          allocate ( plink%platform_BUFR%each(1:nlevels) )
 809          nullify ( plink%next )
 810       end if
 811
 812       ! Recalculate the tdiff  based on the base time at each slot
 813       platform%info%dhr = ( obs_time - &
 814              time_slots(0) - (ifgat-1)*(time_slots(num_fgat_time)-time_slots(0))/float(num_fgat_time-1) &
 815                           ) / 60.0
 816       plink%platform_BUFR%info = platform%info
 817       plink%platform_BUFR%loc = platform%loc
 818       plink%platform_BUFR%each(1:nlevels) = platform%each(1:nlevels)
 819       plink%ifgat_BUFR=ifgat
 820       plink%fm_BUFR=0
 821       plink%nlevels_BUFR=nlevels
 822       plink%kx_BUFR=kx
 823       plink%t29_BUFR=t29
 824       plink%pco_BUFR=pco
 825
 826       num_p=num_p+1
 827
 828    else  !3dvar
 829
 830       tdiff = abs(platform%info%dhr-0.02)
 831       dlat_earth = platform%info%lat
 832       dlon_earth = platform%info%lon
 833       if (dlon_earth < 0.0) dlon_earth = dlon_earth + 360.0
 834       if (dlon_earth >= 360.0) dlon_earth = dlon_earth - 360.0
 835       dlat_earth = dlat_earth * deg2rad
 836       dlon_earth = dlon_earth * deg2rad
 837       ndup = 1
 838       if (global .and. &
 839          (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
 840       if (test_transforms) ndup = 1
 841
 842       ! It is possible that logic for counting obs is incorrect for the
 843       ! global case with >1 MPI tasks due to obs duplication, halo, etc.
 844       ! TBH:  20050913
 845       dup_loop: do n = 1, ndup
 846          select case(t29)
 847          case (11, 12, 13, 22, 23, 31)
 848             select case (kx)
 849             case (120, 122, 132, 220, 222, 232) ;         ! Sound
 850                if (.not.use_soundobs) cycle reports
 851                if (n==1) iv%info(sound)%ntotal     = iv%info(sound)%ntotal + 1
 852                if (n==1) iv%info(sonde_sfc)%ntotal = iv%info(sonde_sfc)%ntotal + 1
 853
 854                if (outside) cycle reports
 855                if ( thin_conv_opt(sound) > no_thin .or. &
 856                     thin_conv_opt(sonde_sfc) > no_thin ) then
 857                   crit = tdiff
 858                   call map2grids_conv(sound,dlat_earth,dlon_earth,crit,iv%info(sound)%nlocal,itx,1,itt,iout,iuse)
 859                   call map2grids_conv(sonde_sfc,dlat_earth,dlon_earth,crit,iv%info(sonde_sfc)%nlocal,itx,1,itt,iout,iuse)
 860
 861                   if ( .not. iuse ) then
 862                      num_thinned = num_thinned + 1
 863                      cycle reports
 864                   end if
 865                else
 866                   iv%info(sound)%nlocal     = iv%info(sound)%nlocal + 1
 867                   iv%info(sonde_sfc)%nlocal = iv%info(sound)%nlocal
 868                end if
 869                fm = 35
 870
 871             case (221) ;                   ! Pilot
 872                if (.not.use_pilotobs) cycle reports
 873                if (n==1) iv%info(pilot)%ntotal = iv%info(pilot)%ntotal + 1
 874                if (outside) cycle reports
 875                if ( thin_conv_opt(pilot) > no_thin ) then
 876                   crit = tdiff
 877                   call map2grids_conv(pilot,dlat_earth,dlon_earth,crit,iv%info(pilot)%nlocal,itx,1,itt,iout,iuse)
 878                   if ( .not. iuse ) then
 879                      num_thinned = num_thinned + 1
 880                      cycle reports
 881                   end if
 882                else
 883                   iv%info(pilot)%nlocal = iv%info(pilot)%nlocal + 1
 884                end if
 885                fm = 32
 886
 887             case default
 888                exit dup_loop
 889             end select
 890
 891          case (41)
 892             ! case (130:131, 133, 230:231, 233) ; ! Airep
 893                if (.not.use_airepobs) cycle reports
 894                if (n==1) iv%info(airep)%ntotal = iv%info(airep)%ntotal + 1
 895                if (outside)  cycle reports
 896                if ( thin_conv_opt(airep) > no_thin ) then
 897                   crit = tdiff
 898                   call map2grids_conv(airep,dlat_earth,dlon_earth,crit,iv%info(airep)%nlocal,itx,1,itt,iout,iuse)
 899                   if ( .not. iuse ) then
 900                      num_thinned = num_thinned + 1
 901                      cycle reports
 902                   end if
 903                else
 904                   iv%info(airep)%nlocal = iv%info(airep)%nlocal + 1
 905                end if
 906                fm = 42
 907
 908          case (522, 523);        ! Ships
 909                if (.not.use_shipsobs) cycle reports
 910                if (n==1) iv%info(ships)%ntotal = iv%info(ships)%ntotal + 1
 911                if (outside) cycle reports
 912                if ( thin_conv_opt(ships) > no_thin ) then
 913                   crit = tdiff
 914                   call map2grids_conv(ships,dlat_earth,dlon_earth,crit,iv%info(ships)%nlocal,itx,1,itt,iout,iuse)
 915                   if ( .not. iuse ) then
 916                      num_thinned = num_thinned + 1
 917                      cycle reports
 918                   end if
 919                else
 920                   iv%info(ships)%nlocal = iv%info(ships)%nlocal + 1
 921                end if
 922                fm = 13
 923
 924          case (531, 532, 561, 562) ;          ! Buoy
 925                if (.not.use_buoyobs) cycle reports
 926                if (n==1) iv%info(buoy)%ntotal = iv%info(buoy)%ntotal + 1
 927                if (outside) cycle reports
 928
 929                if ( thin_conv_opt(buoy) > no_thin ) then
 930                   crit = tdiff
 931                   call map2grids_conv(buoy,dlat_earth,dlon_earth,crit,iv%info(buoy)%nlocal,itx,1,itt,iout,iuse)
 932                   if ( .not. iuse ) then
 933                      num_thinned = num_thinned + 1
 934                      cycle reports
 935                   end if
 936                else
 937                   iv%info(buoy)%nlocal = iv%info(buoy)%nlocal + 1
 938                end if
 939                fm = 18
 940
 941          case (511, 514)
 942             ! case (181, 281) ;                   ! Synop
 943                if (.not.use_synopobs) cycle reports
 944
 945                if (n==1) iv%info(synop)%ntotal = iv%info(synop)%ntotal + 1
 946                if (outside) cycle reports
 947                if ( thin_conv_opt(synop) > no_thin ) then
 948                   crit = tdiff
 949                   call map2grids_conv(synop,dlat_earth,dlon_earth,crit,iv%info(synop)%nlocal,itx,1,itt,iout,iuse)
 950                   if ( .not. iuse ) then
 951                      num_thinned = num_thinned + 1
 952                      cycle reports
 953                   end if
 954                else
 955                   iv%info(synop)%nlocal = iv%info(synop)%nlocal + 1
 956                end if
 957                fm = 12
 958
 959          case (512)
 960             ! case (187, 287) ;                        ! Metar
 961                if (.not.use_metarobs) cycle reports
 962
 963                if (n==1) iv%info(metar)%ntotal = iv%info(metar)%ntotal + 1
 964                if (outside) cycle reports
 965
 966                if ( thin_conv_opt(metar) > no_thin ) then
 967                   crit = tdiff
 968                   call map2grids_conv(metar,dlat_earth,dlon_earth,crit,iv%info(metar)%nlocal,itx,1,itt,iout,iuse)
 969                   if ( .not. iuse ) then
 970                      num_thinned = num_thinned + 1
 971                      cycle reports
 972                   end if
 973                else
 974                   iv%info(metar)%nlocal = iv%info(metar)%nlocal + 1
 975                end if
 976                fm = 15
 977
 978          case (63)
 979             ! case (242:246, 252:253, 255) ;         ! Geo. CMVs
 980                if (.not.use_geoamvobs) cycle reports
 981
 982                if (n==1) iv%info(geoamv)%ntotal = iv%info(geoamv)%ntotal + 1
 983                if (outside) cycle reports
 984
 985                if ( thin_conv_opt(geoamv) > no_thin ) then
 986                   crit = tdiff
 987                   call map2grids_conv(geoamv,dlat_earth,dlon_earth,crit,iv%info(geoamv)%nlocal,itx,1,itt,iout,iuse)
 988                   if ( .not. iuse ) then
 989                      num_thinned = num_thinned + 1
 990                      cycle reports
 991                   end if
 992                else
 993                   iv%info(geoamv)%nlocal = iv%info(geoamv)%nlocal + 1
 994                end if
 995                fm = 88
 996
 997          case (581, 582, 583, 584)      ! ERS 581, QuikSCAT 582, WindSat 583, ASCAT 584
 998                if (.not.use_qscatobs) cycle reports
 999
1000                if (n==1) iv%info(qscat)%ntotal = iv%info(qscat)%ntotal + 1
1001                if (outside) cycle reports
1002
1003                if ( thin_conv_opt(qscat) > no_thin ) then
1004                   crit = tdiff
1005                   call map2grids_conv(qscat,dlat_earth,dlon_earth,crit,iv%info(qscat)%nlocal,itx,1,itt,iout,iuse)
1006                   if ( .not. iuse ) then
1007                      num_thinned = num_thinned + 1
1008                      cycle reports
1009                   end if
1010                else
1011                   iv%info(qscat)%nlocal = iv%info(qscat)%nlocal + 1
1012                end if
1013                fm = 281
1014
1015          case (74)       ! GPS PW
1016                if (.not.use_gpspwobs) cycle reports
1017
1018                if (n==1) iv%info(gpspw)%ntotal = iv%info(gpspw)%ntotal + 1
1019                if (outside) cycle reports
1020
1021                if ( thin_conv_opt(gpspw) > no_thin ) then
1022                   crit = tdiff
1023                   call map2grids_conv(gpspw,dlat_earth,dlon_earth,crit,iv%info(gpspw)%nlocal,itx,1,itt,iout,iuse)
1024                   if ( .not. iuse ) then
1025                      num_thinned = num_thinned + 1
1026                      cycle reports
1027                   end if
1028                else
1029                   iv%info(gpspw)%nlocal = iv%info(gpspw)%nlocal + 1
1030                end if
1031                fm = 111
1032
1033          case (71, 73, 75, 76, 77)    ! Profiler
1034                if (.not.use_profilerobs) cycle reports
1035
1036                if (n==1) iv%info(profiler)%ntotal = iv%info(profiler)%ntotal + 1
1037                if (outside) cycle reports
1038
1039                if ( thin_conv_opt(profiler) > no_thin ) then
1040                   crit = tdiff
1041                   call map2grids_conv(profiler,dlat_earth,dlon_earth,crit,iv%info(profiler)%nlocal,itx,1,itt,iout,iuse)
1042                   if ( .not. iuse ) then
1043                      num_thinned = num_thinned + 1
1044                      cycle reports
1045                   end if
1046                else
1047                   iv%info(profiler)%nlocal = iv%info(profiler)%nlocal + 1
1048                end if
1049                fm = 132
1050
1051          case (571, 65)
1052               if (.not. use_ssmiretrievalobs) cycle reports
1053
1054                if (n==1) iv%info(ssmi_rv)%ntotal = iv%info(ssmi_rv)%ntotal + 1
1055                if (outside) cycle reports
1056
1057                if ( thin_conv_opt(ssmi_rv) > no_thin ) then
1058                   crit = tdiff
1059                   call map2grids_conv(ssmi_rv,dlat_earth,dlon_earth,crit,iv%info(ssmi_rv)%nlocal,itx,1,itt,iout,iuse)
1060                   if ( .not. iuse ) then
1061                      num_thinned = num_thinned + 1
1062                      cycle reports
1063                   end if
1064                else
1065                   iv%info(ssmi_rv)%nlocal = iv%info(ssmi_rv)%nlocal + 1
1066                end if
1067                fm = 125      ! ssmi wind speed & tpw
1068
1069          case default
1070             select case (kx)
1071             case (111 , 210)    ;         !  Tropical Cyclone Bogus
1072                ! Note Tropical cyclone Bougus is given type 135 in Obs-ascii
1073                if (.not.use_bogusobs) cycle reports
1074
1075                if (n==1) iv%info(bogus)%ntotal = iv%info(bogus)%ntotal + 1
1076                if (outside) cycle reports
1077
1078                if ( thin_conv_opt(bogus) > no_thin ) then
1079                   crit = tdiff
1080                   call map2grids_conv(bogus,dlat_earth,dlon_earth,crit,iv%info(bogus)%nlocal,itx,1,itt,iout,iuse)
1081                   if ( .not. iuse ) then
1082                      num_thinned = num_thinned + 1
1083                      cycle reports
1084                   end if
1085                else
1086                   iv%info(bogus)%nlocal = iv%info(bogus)%nlocal + 1
1087                end if
1088                fm = 135
1089
1090             case default
1091                if ( print_detail_obs ) then
1092                   write(unit=message(1), fmt='(a, 2i12)') &
1093                      'unsaved obs found with kx & t29= ',kx,t29
1094                   call da_warning(__FILE__,__LINE__,message(1:1))
1095                end if
1096                exit dup_loop
1097             end select
1098          end select
1099
1100       obs_index=fm_index(fm)
1101       iv%info(obs_index)%max_lev = max(iv%info(obs_index)%max_lev, nlevels)
1102
1103       if (.not. associated(head)) then
1104          nullify ( head )
1105          allocate ( head )
1106          allocate ( head%platform_BUFR%each(1:nlevels) )
1107          nullify ( head%next )
1108          plink => head
1109       else
1110          allocate ( plink%next )
1111          plink => plink%next
1112          allocate ( plink%platform_BUFR%each(1:nlevels) )
1113          nullify ( plink%next )
1114       end if
1115
1116       plink%platform_BUFR%info = platform%info
1117       plink%platform_BUFR%loc = platform%loc
1118       plink%platform_BUFR%each(1:nlevels) = platform%each(1:nlevels)
1119       plink%fm_BUFR=fm
1120       plink%ifgat_BUFR=ifgat
1121       plink%nlevels_BUFR=nlevels
1122       plink%kx_BUFR=kx
1123       plink%t29_BUFR=t29
1124       plink%pco_BUFR=pco
1125
1126       num_p=num_p+1
1127          end do dup_loop
1128        end if !3dvar and 4dvar
1129     end do reports
1130
1131 call closbf(iunit)
1132 close(iunit)
1133
1134 end do bufrfile
1135
1136
1137 !  3.0 Thinning based on FGAT
1138 !      Only for 4dvar
1139 !--------------------------------------------------------------
1140
1141  if (num_fgat_time > 1 ) then
1142
1143     do kk=1,num_fgat_time
1144       if ( thin_conv ) then
1145          do n = 1, num_ob_indexes
1146             if ( .not. allocated(thinning_grid_conv(n)%icount) ) cycle
1147             call cleangrids_conv(n)
1148          end do
1149       end if
1150
1151        plink => head
1152        reports2: do ii=1,num_p
1153
1154        if (plink%ifgat_BUFR /= kk) then  !sort iv
1155             if ( .not. associated(plink%next) ) exit reports2
1156             plink => plink%next
1157             cycle reports2
1158       else
1159                  ! for thinning
1160       tdiff = abs(plink%platform_BUFR%info%dhr-0.02)
1161       dlat_earth = plink%platform_BUFR%info%lat
1162       dlon_earth = plink%platform_BUFR%info%lon
1163       if (dlon_earth < 0.0) dlon_earth = dlon_earth + 360.0
1164       if (dlon_earth >= 360.0) dlon_earth = dlon_earth - 360.0
1165       dlat_earth = dlat_earth * deg2rad
1166       dlon_earth = dlon_earth * deg2rad
1167       call da_llxy (plink%platform_BUFR%info, plink%platform_BUFR%loc,outside, outside_all)
1168
1169  ! Loop over duplicating obs for global
1170       ndup = 1
1171       if (global .and. &
1172          (plink%platform_BUFR%loc%i < ids .or. plink%platform_BUFR%loc%i >= ide)) ndup= 2
1173       if (test_transforms) ndup = 1
1174
1175       ! It is possible that logic for counting obs is incorrect for the
1176       ! global case with >1 MPI tasks due to obs duplication, halo, etc.
1177       ! TBH:  20050913
1178       dup_loop2: do n = 1, ndup
1179          select case(plink%t29_BUFR)
1180          case (11, 12, 13, 22, 23, 31)
1181             select case (plink%kx_BUFR)
1182             case (120, 122, 132, 220, 222, 232) ;         ! Sound
1183                if (.not.use_soundobs) then
1184                       plink => plink%next
1185                       cycle reports2
1186                end if
1187                if (n==1) iv%info(sound)%ntotal     = iv%info(sound)%ntotal + 1
1188                if (n==1) iv%info(sonde_sfc)%ntotal = iv%info(sonde_sfc)%ntotal + 1
1189
1190                if (outside) then
1191                     plink => plink%next
1192                     cycle reports2
1193                end if
1194                if ( thin_conv_opt(sound) > no_thin .or. &
1195                     thin_conv_opt(sonde_sfc) > no_thin ) then
1196                   crit = tdiff
1197                   call map2grids_conv(sound,dlat_earth,dlon_earth,crit,iv%info(sound)%nlocal,itx,1,itt,iout,iuse)
1198                   call map2grids_conv(sonde_sfc,dlat_earth,dlon_earth,crit,iv%info(sonde_sfc)%nlocal,itx,1,itt,iout,iuse)
1199
1200                   if ( .not. iuse ) then
1201                      num_thinned = num_thinned + 1
1202                      plink => plink%next
1203                      cycle reports2
1204                   end if
1205                else
1206                   iv%info(sound)%nlocal     = iv%info(sound)%nlocal + 1
1207                   iv%info(sonde_sfc)%nlocal = iv%info(sound)%nlocal
1208                end if
1209                fm = 35
1210
1211             case (221) ;                   ! Pilot
1212                if (.not.use_pilotobs) then
1213                      plink => plink%next
1214                      cycle reports2
1215                end if
1216                if (n==1) iv%info(pilot)%ntotal = iv%info(pilot)%ntotal + 1
1217                if (outside) then
1218                    plink => plink%next
1219                    cycle reports2
1220                end if
1221                if ( thin_conv_opt(pilot) > no_thin ) then
1222                   crit = tdiff
1223                   call map2grids_conv(pilot,dlat_earth,dlon_earth,crit,iv%info(pilot)%nlocal,itx,1,itt,iout,iuse)
1224                   if ( .not. iuse ) then
1225                      num_thinned = num_thinned + 1
1226                      plink => plink%next
1227                      cycle reports2
1228                   end if
1229                else
1230                   iv%info(pilot)%nlocal = iv%info(pilot)%nlocal + 1
1231                end if
1232                fm = 32
1233
1234             case default
1235                exit dup_loop2
1236             end select
1237
1238          case (41)
1239             ! case (130:131, 133, 230:231, 233) ; ! Airep
1240                if (.not.use_airepobs) then
1241                       plink => plink%next
1242                       cycle reports2
1243                end if
1244                if (n==1) iv%info(airep)%ntotal = iv%info(airep)%ntotal + 1
1245                if (outside) then
1246                     plink => plink%next
1247                     cycle reports2
1248                end if
1249                if ( thin_conv_opt(airep) > no_thin ) then
1250                   crit = tdiff
1251                   call map2grids_conv(airep,dlat_earth,dlon_earth,crit,iv%info(airep)%nlocal,itx,1,itt,iout,iuse)
1252                   if ( .not. iuse ) then
1253                      num_thinned = num_thinned + 1
1254                      plink => plink%next
1255                      cycle reports2
1256                   end if
1257                else
1258                   iv%info(airep)%nlocal = iv%info(airep)%nlocal + 1
1259                end if
1260                fm = 42
1261
1262          case (522, 523);        ! Ships
1263                if (.not.use_shipsobs) then
1264                    plink => plink%next
1265                    cycle reports2
1266                end if
1267                if (n==1) iv%info(ships)%ntotal = iv%info(ships)%ntotal + 1
1268                if (outside) then
1269                   plink => plink%next
1270                   cycle reports2
1271                end if
1272                if ( thin_conv_opt(ships) > no_thin ) then
1273                   crit = tdiff
1274                   call map2grids_conv(ships,dlat_earth,dlon_earth,crit,iv%info(ships)%nlocal,itx,1,itt,iout,iuse)
1275                   if ( .not. iuse ) then
1276                      num_thinned = num_thinned + 1
1277                      plink => plink%next
1278                      cycle reports2
1279                   end if
1280                else
1281                   iv%info(ships)%nlocal = iv%info(ships)%nlocal + 1
1282                end if
1283                fm = 13
1284
1285          case (531, 532, 561, 562) ;          ! Buoy
1286                if (.not.use_buoyobs) then
1287                    plink => plink%next
1288                    cycle reports2
1289                end if
1290                if (n==1) iv%info(buoy)%ntotal = iv%info(buoy)%ntotal + 1
1291                if (outside) then
1292                      plink => plink%next
1293                      cycle reports2
1294                end if
1295                if ( thin_conv_opt(buoy) > no_thin ) then
1296                   crit = tdiff
1297                   call map2grids_conv(buoy,dlat_earth,dlon_earth,crit,iv%info(buoy)%nlocal,itx,1,itt,iout,iuse)
1298                   if ( .not. iuse ) then
1299                      num_thinned = num_thinned + 1
1300                      plink => plink%next
1301                      cycle reports2
1302                   end if
1303                else
1304                   iv%info(buoy)%nlocal = iv%info(buoy)%nlocal + 1
1305                end if
1306                fm = 18
1307
1308          case (511, 514)
1309             ! case (181, 281) ;                   ! Synop
1310                if (.not.use_synopobs) then
1311                    plink => plink%next
1312                    cycle reports2
1313                end if
1314                if (n==1) iv%info(synop)%ntotal = iv%info(synop)%ntotal + 1
1315                if (outside) then
1316                      plink => plink%next
1317                      cycle reports2
1318                end if
1319                if ( thin_conv_opt(synop) > no_thin ) then
1320                   crit = tdiff
1321                   call map2grids_conv(synop,dlat_earth,dlon_earth,crit,iv%info(synop)%nlocal,itx,1,itt,iout,iuse)
1322                   if ( .not. iuse ) then
1323                      num_thinned = num_thinned + 1
1324                      plink => plink%next
1325                      cycle reports2
1326                   end if
1327                else
1328                   iv%info(synop)%nlocal = iv%info(synop)%nlocal + 1
1329                end if
1330                fm = 12
1331
1332          case (512)
1333             ! case (187, 287) ;                        ! Metar
1334                if (.not.use_metarobs) then
1335                     plink => plink%next
1336                     cycle reports2
1337                end if
1338                if (n==1) iv%info(metar)%ntotal = iv%info(metar)%ntotal + 1
1339                if (outside) then
1340                     plink => plink%next
1341                     cycle reports2
1342                end if
1343                if ( thin_conv_opt(metar) > no_thin ) then
1344                   crit = tdiff
1345                   call map2grids_conv(metar,dlat_earth,dlon_earth,crit,iv%info(metar)%nlocal,itx,1,itt,iout,iuse)
1346                   if ( .not. iuse ) then
1347                      num_thinned = num_thinned + 1
1348                      plink => plink%next
1349                      cycle reports2
1350                   end if
1351                else
1352                   iv%info(metar)%nlocal = iv%info(metar)%nlocal + 1
1353                end if
1354                fm = 15
1355
1356          case (63)
1357             ! case (242:246, 252:253, 255) ;         ! Geo. CMVs
1358                if (.not.use_geoamvobs) then
1359                        plink => plink%next
1360                        cycle reports2
1361                end if
1362                if (n==1) iv%info(geoamv)%ntotal = iv%info(geoamv)%ntotal + 1
1363                if (outside) then
1364                   plink => plink%next
1365                   cycle reports2
1366                end if
1367                if ( thin_conv_opt(geoamv) > no_thin ) then
1368                   crit = tdiff
1369                   call map2grids_conv(geoamv,dlat_earth,dlon_earth,crit,iv%info(geoamv)%nlocal,itx,1,itt,iout,iuse)
1370                   if ( .not. iuse ) then
1371                      num_thinned = num_thinned + 1
1372                      plink => plink%next
1373                      cycle reports2
1374                   end if
1375                else
1376                   iv%info(geoamv)%nlocal = iv%info(geoamv)%nlocal + 1
1377                end if
1378                fm = 88
1379
1380          case (581, 582, 583, 584)      ! ERS 581, QuikSCAT 582, WindSat 583, ASCAT 584
1381                if (.not.use_qscatobs) then
1382                       plink => plink%next
1383                       cycle reports2
1384                end if
1385                if (n==1) iv%info(qscat)%ntotal = iv%info(qscat)%ntotal + 1
1386                if (outside) then
1387                    plink => plink%next
1388                    cycle reports2
1389                end if
1390                if ( thin_conv_opt(qscat) > no_thin ) then
1391                   crit = tdiff
1392                   call map2grids_conv(qscat,dlat_earth,dlon_earth,crit,iv%info(qscat)%nlocal,itx,1,itt,iout,iuse)
1393                   if ( .not. iuse ) then
1394                      num_thinned = num_thinned + 1
1395                      plink => plink%next
1396                      cycle reports2
1397                   end if
1398                else
1399                   iv%info(qscat)%nlocal = iv%info(qscat)%nlocal + 1
1400                end if
1401                fm = 281
1402
1403          case (74)       ! GPS PW
1404                if (.not.use_gpspwobs) then
1405                     plink => plink%next
1406                     cycle reports2
1407                end if
1408                if (n==1) iv%info(gpspw)%ntotal = iv%info(gpspw)%ntotal + 1
1409                if (outside) then
1410                    plink => plink%next
1411                    cycle reports2
1412                end if
1413                if ( thin_conv_opt(gpspw) > no_thin ) then
1414                   crit = tdiff
1415                   call map2grids_conv(gpspw,dlat_earth,dlon_earth,crit,iv%info(gpspw)%nlocal,itx,1,itt,iout,iuse)
1416                   if ( .not. iuse ) then
1417                      num_thinned = num_thinned + 1
1418                      plink => plink%next
1419                      cycle reports2
1420                   end if
1421                else
1422                   iv%info(gpspw)%nlocal = iv%info(gpspw)%nlocal + 1
1423                end if
1424                fm = 111
1425
1426          case (71, 73, 75, 76, 77)    ! Profiler
1427                if (.not.use_profilerobs) then
1428                      plink => plink%next
1429                      cycle reports2
1430                end if
1431                if (n==1) iv%info(profiler)%ntotal = iv%info(profiler)%ntotal + 1
1432                if (outside) then
1433                   plink => plink%next
1434                   cycle reports2
1435                end if
1436                if ( thin_conv_opt(profiler) > no_thin ) then
1437                   crit = tdiff
1438                   call map2grids_conv(profiler,dlat_earth,dlon_earth,crit,iv%info(profiler)%nlocal,itx,1,itt,iout,iuse)
1439                   if ( .not. iuse ) then
1440                      num_thinned = num_thinned + 1
1441                      plink => plink%next
1442                      cycle reports2
1443                   end if
1444                else
1445                   iv%info(profiler)%nlocal = iv%info(profiler)%nlocal + 1
1446                end if
1447                fm = 132
1448
1449          case (571, 65)
1450               if (.not. use_ssmiretrievalobs) then
1451                          plink => plink%next
1452                          cycle reports2
1453                end if
1454                if (n==1) iv%info(ssmi_rv)%ntotal = iv%info(ssmi_rv)%ntotal + 1
1455                if (outside) then
1456                     plink => plink%next
1457                     cycle reports2
1458                end if
1459                if ( thin_conv_opt(ssmi_rv) > no_thin ) then
1460                   crit = tdiff
1461                   call map2grids_conv(ssmi_rv,dlat_earth,dlon_earth,crit,iv%info(ssmi_rv)%nlocal,itx,1,itt,iout,iuse)
1462                   if ( .not. iuse ) then
1463                      num_thinned = num_thinned + 1
1464                      plink => plink%next
1465                      cycle reports2
1466                   end if
1467                else
1468                   iv%info(ssmi_rv)%nlocal = iv%info(ssmi_rv)%nlocal + 1
1469                end if
1470                fm = 125      ! ssmi wind speed & tpw
1471
1472          case default
1473             select case (plink%kx_BUFR)
1474             case (111 , 210)    ;         !  Tropical Cyclone Bogus
1475                ! Note Tropical cyclone Bougus is given type 135 in Obs-ascii
1476                if (.not.use_bogusobs) then
1477                       plink => plink%next
1478                       cycle reports2
1479                end if
1480                if (n==1) iv%info(bogus)%ntotal = iv%info(bogus)%ntotal + 1
1481                if (outside) then
1482                   plink => plink%next
1483                   cycle reports2
1484                end if
1485                if ( thin_conv_opt(bogus) > no_thin ) then
1486                   crit = tdiff
1487                   call map2grids_conv(bogus,dlat_earth,dlon_earth,crit,iv%info(bogus)%nlocal,itx,1,itt,iout,iuse)
1488                   if ( .not. iuse ) then
1489                      num_thinned = num_thinned + 1
1490                      plink => plink%next
1491                      cycle reports2
1492                   end if
1493                else
1494                   iv%info(bogus)%nlocal = iv%info(bogus)%nlocal + 1
1495                end if
1496                fm = 135
1497
1498             case default
1499                if ( print_detail_obs ) then
1500                   write(unit=message(1), fmt='(a, 2i12)') &
1501                      'unsaved obs found with kx & t29= ',plink%kx_BUFR,plink%t29_BUFR
1502                   call da_warning(__FILE__,__LINE__,message(1:1))
1503                end if
1504                exit dup_loop2
1505             end select
1506          end select
1507
1508       obs_index=fm_index(fm)
1509       iv%info(obs_index)%max_lev = max(iv%info(obs_index)%max_lev, plink%nlevels_BUFR)
1510
1511       plink%fm_BUFR=fm
1512
1513       end do dup_loop2
1514      end if  !sort iv
1515
1516      plink => plink%next
1517      if ( .not. associated(plink) ) exit reports2
1518
1519    end do reports2
1520
1521 end do !kk
1522
1523 end if
1524
1525 !  4.0 Allocate iv
1526 !
1527 !--------------------------------------------------------------
1528
1529    iv%info(synop)%max_lev     = 1
1530    iv%info(metar)%max_lev     = 1
1531    iv%info(ships)%max_lev     = 1
1532    iv%info(buoy)%max_lev      = 1
1533    iv%info(sonde_sfc)%max_lev = 1
1534
1535    call da_allocate_observations (iv)
1536
1537 !  5.0 Transfer p structure into iv structure
1538 !      Also sort iv structure to FGAT time bins
1539 !--------------------------------------------------------------
1540     do kk=1,num_fgat_time
1541
1542        if ( thin_conv ) then
1543           do n = 1, num_ob_indexes
1544              if ( .not. allocated(thinning_grid_conv(n)%icount) ) cycle
1545              call cleangrids_conv(n)
1546           end do
1547        end if
1548
1549        plink => head
1550
1551       reports3: do ii=1,num_p
1552
1553      if (plink%ifgat_BUFR /=  kk) then  !sort iv
1554             if ( .not. associated(plink%next) ) exit reports3
1555             plink => plink%next
1556             cycle reports3
1557       else
1558       ! for thinning
1559       tdiff = abs(plink%platform_BUFR%info%dhr-0.02)
1560       dlat_earth = plink%platform_BUFR%info%lat
1561       dlon_earth = plink%platform_BUFR%info%lon
1562       if (dlon_earth < 0.0) dlon_earth = dlon_earth + 360.0
1563       if (dlon_earth >= 360.0) dlon_earth = dlon_earth - 360.0
1564       dlat_earth = dlat_earth * deg2rad
1565       dlon_earth = dlon_earth * deg2rad
1566       call da_llxy (plink%platform_BUFR%info, plink%platform_BUFR%loc,outside, outside_all)
1567       ndup = 1
1568       if (global .and. &
1569          (plink%platform_BUFR%loc%i < ids .or. plink%platform_BUFR%loc%i >= ide)) ndup= 2
1570       if (test_transforms) ndup = 1
1571       dup_loop3: do n = 1, ndup
1572
1573          select case(plink%t29_BUFR)
1574            case (11, 12, 13, 22, 23, 31)
1575             select case (plink%kx_BUFR)
1576               case (120, 122, 132, 220, 222, 232) ;         ! Sound
1577                if (.not.use_soundobs) then
1578                       plink => plink%next
1579                       cycle reports3
1580                end if
1581                  if (n==1) ntotal(sound) = ntotal(sound) + 1
1582                  if (n==1) ntotal(sonde_sfc) = ntotal(sonde_sfc) + 1
1583                  if (outside) then
1584                     plink => plink%next
1585                     cycle reports3
1586                end if
1587                 if ( thin_conv_opt(sound) > no_thin .or. &
1588                      thin_conv_opt(sonde_sfc) > no_thin ) then
1589                   crit = tdiff
1590                   call map2grids_conv(sound,dlat_earth,dlon_earth,crit,nlocal(sound),itx,1,itt,ilocal(sound),iuse)
1591                   call map2grids_conv(sonde_sfc,dlat_earth,dlon_earth,crit,nlocal(sonde_sfc),itx,1,itt,ilocal(sonde_sfc),iuse)
1592                    if ( .not. iuse ) then
1593                      plink => plink%next
1594                      cycle reports3
1595                   end if
1596                else
1597                   nlocal(sound) = nlocal(sound) + 1
1598                   nlocal(sonde_sfc) = nlocal(sound)
1599                   ilocal(sound) = nlocal(sound)
1600                   ilocal(sonde_sfc) = ilocal(sound)
1601                end if
1602
1603                platform_name ='FM-35 TEMP  '
1604
1605                old_nlevels = plink%nlevels_BUFR
1606
1607                ! Search to see if we have surface obs.
1608
1609                surface_level = 0
1610
1611                do i = 1, plink%nlevels_BUFR
1612                   ! if (elevation and height are the same, it is surface)
1613                   if (abs(plink%platform_BUFR%info%elv - &
1614                      plink%platform_BUFR%each(i)%height) < 0.1) then
1615                      surface_level = i
1616
1617                      ! Save surface pressure.
1618                      iv%sonde_sfc(ilocal(sonde_sfc))%h = plink%platform_BUFR%each(i)%height
1619                      iv%sonde_sfc(ilocal(sonde_sfc))%u = plink%platform_BUFR%each(i)%u
1620                      iv%sonde_sfc(ilocal(sonde_sfc))%v = plink%platform_BUFR%each(i)%v
1621                      iv%sonde_sfc(ilocal(sonde_sfc))%t = plink%platform_BUFR%each(i)%t
1622                      iv%sonde_sfc(ilocal(sonde_sfc))%q = plink%platform_BUFR%each(i)%q
1623                      iv%sonde_sfc(ilocal(sonde_sfc))%p = plink%platform_BUFR%each(i)%p
1624
1625                      exit
1626                   end if
1627                end do
1628
1629                ! processing the sound_sfc data:
1630
1631                if (surface_level > 0) then
1632                   plink%nlevels_BUFR = plink%nlevels_BUFR - 1
1633                else   !missing surface data
1634                   iv%sonde_sfc(ilocal(sonde_sfc))%h = missing_r
1635                   iv%sonde_sfc(ilocal(sonde_sfc))%u%inv   = missing_r
1636                   iv%sonde_sfc(ilocal(sonde_sfc))%u%qc    = missing_data
1637                   iv%sonde_sfc(ilocal(sonde_sfc))%u%error = abs(missing_r)
1638                   iv%sonde_sfc(ilocal(sonde_sfc))%v = iv%sonde_sfc(ilocal(sonde_sfc))%u
1639                   iv%sonde_sfc(ilocal(sonde_sfc))%t = iv%sonde_sfc(ilocal(sonde_sfc))%u
1640                   iv%sonde_sfc(ilocal(sonde_sfc))%p = iv%sonde_sfc(ilocal(sonde_sfc))%u
1641                   iv%sonde_sfc(ilocal(sonde_sfc))%q = iv%sonde_sfc(ilocal(sonde_sfc))%u
1642
1643                end if
1644
1645                if (plink%nlevels_BUFR > 0) then
1646
1647                   if ( ilocal(sound) == nlocal(sound) .or. &
1648                        .not. associated(iv%sound(ilocal(sound))%h) ) then
1649                      allocate (iv%sound(ilocal(sound))%h(1:iv%info(sound)%max_lev))
1650                      allocate (iv%sound(ilocal(sound))%p(1:iv%info(sound)%max_lev))
1651                      allocate (iv%sound(ilocal(sound))%u(1:iv%info(sound)%max_lev))
1652                      allocate (iv%sound(ilocal(sound))%v(1:iv%info(sound)%max_lev))
1653                      allocate (iv%sound(ilocal(sound))%t(1:iv%info(sound)%max_lev))
1654                      allocate (iv%sound(ilocal(sound))%q(1:iv%info(sound)%max_lev))
1655                    endif
1656
1657                   j = 0
1658                   do i = 1, old_nlevels
1659                      if (i == surface_level) cycle
1660                      j=j+1
1661                      iv%sound(ilocal(sound))%h(j) = plink%platform_BUFR%each(i)%height
1662                      iv%sound(ilocal(sound))%p(j) = plink%platform_BUFR%each(i)%p%inv
1663                      iv%sound(ilocal(sound))%u(j) = plink%platform_BUFR%each(i)%u
1664                      iv%sound(ilocal(sound))%v(j) = plink%platform_BUFR%each(i)%v
1665                      iv%sound(ilocal(sound))%t(j) = plink%platform_BUFR%each(i)%t
1666                      iv%sound(ilocal(sound))%q(j) = plink%platform_BUFR%each(i)%q
1667
1668                   end do
1669                end if
1670
1671             case (221) ;           ! Pilot
1672                 if (.not.use_pilotobs) then
1673                      plink => plink%next
1674                      cycle reports3
1675                end if
1676                if (n==1) ntotal(pilot) = ntotal(pilot) + 1
1677                if (outside) then
1678                    plink => plink%next
1679                    cycle reports3
1680                end if
1681                if ( thin_conv_opt(pilot) > no_thin ) then
1682                   crit = tdiff
1683                   call map2grids_conv(pilot,dlat_earth,dlon_earth,crit,nlocal(pilot),itx,1,itt,ilocal(pilot),iuse)
1684                   if ( .not. iuse ) then
1685                      plink => plink%next
1686                      cycle reports3
1687                   end if
1688                else
1689                   nlocal(pilot) = nlocal(pilot) + 1
1690                   ilocal(pilot) = nlocal(pilot)
1691                end if
1692                 platform_name='FM-32 PILOT '
1693
1694                if ( ilocal(pilot) == nlocal(pilot) ) then
1695                   allocate (iv%pilot(ilocal(pilot))%h(1:iv%info(pilot)%max_lev))
1696                   allocate (iv%pilot(ilocal(pilot))%p(1:iv%info(pilot)%max_lev))
1697                   allocate (iv%pilot(ilocal(pilot))%u(1:iv%info(pilot)%max_lev))
1698                   allocate (iv%pilot(ilocal(pilot))%v(1:iv%info(pilot)%max_lev))
1699                endif
1700
1701                do i = 1, plink%nlevels_BUFR
1702                   iv%pilot(ilocal(pilot))%h(i) = plink%platform_BUFR%each(i)%height
1703                   iv%pilot(ilocal(pilot))%p(i) = plink%platform_BUFR%each(i)%p%inv
1704                   iv%pilot(ilocal(pilot))%u(i) = plink%platform_BUFR%each(i)%u
1705                   iv%pilot(ilocal(pilot))%v(i) = plink%platform_BUFR%each(i)%v
1706
1707                end do
1708                 case default
1709                exit dup_loop3
1710             end select
1711          case (41)
1712              if (.not.use_airepobs) then
1713                       plink => plink%next
1714                       cycle reports3
1715                end if
1716
1717             if (n==1) ntotal(airep) = ntotal(airep) + 1
1718              if (outside) then
1719                     plink => plink%next
1720                     cycle reports3
1721                end if
1722             if ( thin_conv_opt(airep) > no_thin ) then
1723                crit = tdiff
1724                call map2grids_conv(airep,dlat_earth,dlon_earth,crit,nlocal(airep),itx,1,itt,ilocal(airep),iuse)
1725                if ( .not. iuse ) then
1726                     plink => plink%next
1727                     cycle reports3
1728                   end if
1729             else
1730                nlocal(airep) = nlocal(airep) + 1
1731                ilocal(airep) = nlocal(airep)
1732             end if
1733             platform_name ='FM-97 AIREP '
1734             if ( ilocal(airep) == nlocal(airep) ) then
1735                allocate (iv%airep(ilocal(airep))%h(1:iv%info(airep)%max_lev))
1736                allocate (iv%airep(ilocal(airep))%p(1:iv%info(airep)%max_lev))
1737                allocate (iv%airep(ilocal(airep))%u(1:iv%info(airep)%max_lev))
1738                allocate (iv%airep(ilocal(airep))%v(1:iv%info(airep)%max_lev))
1739                allocate (iv%airep(ilocal(airep))%t(1:iv%info(airep)%max_lev))
1740                allocate (iv%airep(ilocal(airep))%q(1:iv%info(airep)%max_lev))
1741             end if
1742             do i = 1, plink%nlevels_BUFR
1743                iv % airep (ilocal(airep)) % h(i) = plink%platform_BUFR % each(i) % height
1744                iv % airep (ilocal(airep)) % p(i) = plink%platform_BUFR % each(i) % p % inv
1745                iv % airep (ilocal(airep)) % u(i) = plink%platform_BUFR % each(i) % u
1746                iv % airep (ilocal(airep)) % v(i) = plink%platform_BUFR % each(i) % v
1747                iv % airep (ilocal(airep)) % t(i) = plink%platform_BUFR % each(i) % t
1748                iv % airep (ilocal(airep)) % q(i) = plink%platform_BUFR % each(i) % q
1749             end do
1750
1751          case (522,523);        ! Ships
1752           if (.not.use_shipsobs) then
1753                    plink => plink%next
1754                    cycle reports3
1755            end if
1756
1757             if (n==1) ntotal(ships) = ntotal(ships) + 1
1758              if (outside) then
1759                   plink => plink%next
1760                   cycle reports3
1761                end if
1762             if ( thin_conv_opt(ships) > no_thin ) then
1763                crit = tdiff
1764                call map2grids_conv(ships,dlat_earth,dlon_earth,crit,nlocal(ships),itx,1,itt,ilocal(ships),iuse)
1765                 if ( .not. iuse ) then
1766                    plink => plink%next
1767                    cycle reports3
1768                   end if
1769             else
1770                nlocal(ships) = nlocal(ships) + 1
1771                ilocal(ships) = nlocal(ships)
1772             end if
1773             platform_name ='FM-13 SHIP  '
1774
1775             iv % ships (ilocal(ships)) % h = plink%platform_BUFR % each(1) % height
1776             iv % ships (ilocal(ships)) % u = plink%platform_BUFR % each(1) % u
1777             iv % ships (ilocal(ships)) % v = plink%platform_BUFR % each(1) % v
1778             iv % ships (ilocal(ships)) % t = plink%platform_BUFR % each(1) % t
1779             iv % ships (ilocal(ships)) % p = plink%platform_BUFR % each(1) % p
1780             iv % ships (ilocal(ships)) % q = plink%platform_BUFR % each(1) % q
1781
1782
1783           case (531, 532, 561, 562) ;          ! Buoy
1784            if (.not.use_buoyobs) then
1785                    plink => plink%next
1786                    cycle reports3
1787                end if
1788             if (n==1)  ntotal(buoy) = ntotal(buoy) + 1
1789              if (outside) then
1790                 plink => plink%next
1791                 cycle reports3
1792                end if
1793             if ( thin_conv_opt(buoy) > no_thin ) then
1794                crit = tdiff
1795                call map2grids_conv(buoy,dlat_earth,dlon_earth,crit,nlocal(buoy),itx,1,itt,ilocal(buoy),iuse)
1796                if ( .not. iuse ) then
1797                   plink => plink%next
1798                   cycle reports3
1799                end if
1800             else
1801                nlocal(buoy) = nlocal(buoy) + 1
1802                ilocal(buoy) = nlocal(buoy)
1803             end if
1804             platform_name ='FM-18 BUOY  '
1805
1806             iv%buoy(ilocal(buoy))%h = plink%platform_BUFR%each(1)%height
1807             iv%buoy(ilocal(buoy))%u = plink%platform_BUFR%each(1)%u
1808             iv%buoy(ilocal(buoy))%v = plink%platform_BUFR%each(1)%v
1809             iv%buoy(ilocal(buoy))%t = plink%platform_BUFR%each(1)%t
1810
1811
1812             if ( plink%kx_BUFR == 282 ) then  ! ATLAS BUOY, reported Pstn and Pmslp are both
1813                                    ! missing. Pstn is set to 1013hPa in PREPBUFR
1814                iv%buoy(ilocal(buoy))%p%inv   = missing_r
1815                iv%buoy(ilocal(buoy))%p%qc    = missing_data
1816                iv%buoy(ilocal(buoy))%p%error = err_p
1817             else
1818                iv%buoy(ilocal(buoy))%p = plink%platform_BUFR%each(1)%p
1819             end if
1820             iv%buoy(ilocal(buoy))%q = plink%platform_BUFR%each(1)%q
1821
1822          case (511, 514)
1823               if (.not.use_synopobs) then
1824                  plink => plink%next
1825                  cycle reports3
1826                end if
1827             if (n==1) ntotal(synop) = ntotal(synop) + 1
1828             if (outside) then
1829                plink => plink%next
1830                cycle reports3
1831             end if
1832             if ( thin_conv_opt(synop) > no_thin ) then
1833                crit = tdiff
1834                call map2grids_conv(synop,dlat_earth,dlon_earth,crit,nlocal(synop),itx,1,itt,ilocal(synop),iuse)
1835                 if ( .not. iuse ) then
1836                    plink => plink%next
1837                    cycle reports3
1838                 end if
1839             else
1840                nlocal(synop) = nlocal(synop) + 1
1841                ilocal(synop) = nlocal(synop)
1842             end if
1843             platform_name ='FM-12 SYNOP '
1844
1845             iv % synop (ilocal(synop)) % h = plink%platform_BUFR % each(1) % height
1846             iv % synop (ilocal(synop)) % u = plink%platform_BUFR % each(1) % u
1847             iv % synop (ilocal(synop)) % v = plink%platform_BUFR % each(1) % v
1848             iv % synop (ilocal(synop)) % t = plink%platform_BUFR % each(1) % t
1849             iv % synop (ilocal(synop)) % p = plink%platform_BUFR % each(1) % p
1850             iv % synop (ilocal(synop)) % q = plink%platform_BUFR % each(1) % q
1851
1852             if (iv % synop(ilocal(synop)) % h < plink%platform_BUFR % info % elv) then
1853                iv % synop(ilocal(synop)) % h = plink%platform_BUFR % info % elv
1854             end if
1855
1856          case (512)
1857           if (.not.use_metarobs) then
1858                    plink => plink%next
1859                    cycle reports3
1860              end if
1861             if (n==1) ntotal(metar) = ntotal(metar) + 1
1862              if (outside) then
1863                    plink => plink%next
1864                    cycle reports3
1865                end if
1866             if ( thin_conv_opt(metar) > no_thin ) then
1867                crit = tdiff
1868                call map2grids_conv(metar,dlat_earth,dlon_earth,crit,nlocal(metar),itx,1,itt,ilocal(metar),iuse)
1869                 if ( .not. iuse ) then
1870                    plink => plink%next
1871                    cycle reports3
1872                   end if
1873
1874             else
1875                nlocal(metar) = nlocal(metar) + 1
1876                ilocal(metar) = nlocal(metar)
1877             end if
1878
1879              platform_name='FM-15 METAR '
1880
1881             iv % metar (ilocal(metar)) % h = plink%platform_BUFR % each(1) % height
1882             iv % metar (ilocal(metar)) % u = plink%platform_BUFR % each(1) % u
1883             iv % metar (ilocal(metar)) % v = plink%platform_BUFR % each(1) % v
1884             iv % metar (ilocal(metar)) % t = plink%platform_BUFR % each(1) % t
1885             iv % metar (ilocal(metar)) % p = plink%platform_BUFR % each(1) % p
1886             iv % metar (ilocal(metar)) % q = plink%platform_BUFR % each(1) % q
1887
1888
1889          case (63)
1890               if (.not.use_geoamvobs) then
1891                    plink => plink%next
1892                    cycle reports3
1893                end if
1894                if (n==1) ntotal(geoamv) = ntotal(geoamv) + 1
1895                 if (outside) then
1896                    plink => plink%next
1897                    cycle reports3
1898                end if
1899                if ( thin_conv_opt(geoamv) > no_thin ) then
1900                crit = tdiff
1901                call map2grids_conv(geoamv,dlat_earth,dlon_earth,crit,nlocal(geoamv),itx,1,itt,ilocal(geoamv),iuse)
1902                 if ( .not. iuse ) then
1903                    plink => plink%next
1904                    cycle reports3
1905                   end if
1906             else
1907                nlocal(geoamv) = nlocal(geoamv) + 1
1908                ilocal(geoamv) = nlocal(geoamv)
1909             end if
1910
1911                platform_name ='FM-88 SATOB '
1912             if ( ilocal(geoamv) == nlocal(geoamv) ) then
1913                allocate (iv%geoamv(ilocal(geoamv))%p(1:iv%info(geoamv)%max_lev))
1914                allocate (iv%geoamv(ilocal(geoamv))%u(1:iv%info(geoamv)%max_lev))
1915                allocate (iv%geoamv(ilocal(geoamv))%v(1:iv%info(geoamv)%max_lev))
1916
1917             end if
1918
1919             do i = 1, plink%nlevels_BUFR
1920                iv % geoamv (ilocal(geoamv)) % p(i)  = plink%platform_BUFR % each(i) % p % inv
1921                iv % geoamv (ilocal(geoamv)) % u(i)  = plink%platform_BUFR % each(i) % u
1922                iv % geoamv (ilocal(geoamv)) % v(i)  = plink%platform_BUFR % each(i) % v
1923
1924             end do
1925
1926           case (581, 582, 583, 584)      ! ERS 581, QuikSCAT 582, WindSat 583, ASCAT 584
1927           if (.not.use_qscatobs) then
1928                    plink => plink%next
1929                    cycle reports3
1930                end if
1931             if (n==1) ntotal(qscat) = ntotal(qscat) + 1
1932             if (outside) then
1933                    plink => plink%next
1934                    cycle reports3
1935                end if
1936             if ( thin_conv_opt(qscat) > no_thin ) then
1937                crit = tdiff
1938                call map2grids_conv(qscat,dlat_earth,dlon_earth,crit,nlocal(qscat),itx,1,itt,ilocal(qscat),iuse)
1939                if ( .not. iuse ) then
1940                    plink => plink%next
1941                    cycle reports3
1942                   end if
1943             else
1944                nlocal(qscat) = nlocal(qscat) + 1
1945                ilocal(qscat) = nlocal(qscat)
1946             end if
1947
1948             platform_name ='FM-281 Quiks'
1949
1950             ! prepbufr uses pressure not height, so hardwire height to
1951             ! 0 (sea-level)
1952             iv%qscat(ilocal(qscat))%h = 0.0
1953             iv%qscat(ilocal(qscat))%u = plink%platform_BUFR%each(1)%u
1954             iv%qscat(ilocal(qscat))%v = plink%platform_BUFR%each(1)%v
1955             iv%qscat(ilocal(qscat))%u%error = max(plink%platform_BUFR%each(1)%u%error,1.0)
1956             iv%qscat(ilocal(qscat))%v%error = max(plink%platform_BUFR%each(1)%v%error,1.0)
1957
1958
1959          case (74)       ! GPS PW
1960           if (.not.use_gpspwobs) then
1961                    plink => plink%next
1962                    cycle reports3
1963                end if
1964             if (n==1)  ntotal(gpspw) = ntotal(gpspw) + 1
1965              if (outside) then
1966                    plink => plink%next
1967                    cycle reports3
1968                end if
1969             if ( thin_conv_opt(gpspw) > no_thin ) then
1970                crit = tdiff
1971                call map2grids_conv(gpspw,dlat_earth,dlon_earth,crit,nlocal(gpspw),itx,1,itt,ilocal(gpspw),iuse)
1972                 if ( .not. iuse ) then
1973                    plink => plink%next
1974                    cycle reports3
1975                   end if
1976             else
1977                nlocal(gpspw) = nlocal(gpspw) + 1
1978                ilocal(gpspw) = nlocal(gpspw)
1979             end if
1980             platform_name ='FM-111 GPSPW'
1981
1982             iv%gpspw(ilocal(gpspw))%tpw  = plink%platform_BUFR%loc%pw
1983
1984
1985           case (71, 73, 75, 76, 77)
1986                 if (.not.use_profilerobs) then
1987                    plink => plink%next
1988                    cycle reports3
1989                end if
1990                if (n==1) ntotal(profiler) = ntotal(profiler) + 1
1991                 if (outside) then
1992                    plink => plink%next
1993                    cycle reports3
1994                end if
1995                if ( thin_conv_opt(profiler) > no_thin ) then
1996                crit = tdiff
1997                call map2grids_conv(profiler,dlat_earth,dlon_earth,crit,nlocal(profiler),itx,1,itt,ilocal(profiler),iuse)
1998                 if ( .not. iuse ) then
1999                    plink => plink%next
2000                    cycle reports3
2001                   end if
2002             else
2003                nlocal(profiler) = nlocal(profiler) + 1
2004                ilocal(profiler) = nlocal(profiler)
2005             end if
2006
2007                platform_name ='FM-132 PRFLR'
2008             if ( ilocal(profiler) == nlocal(profiler) ) then
2009                allocate (iv%profiler(ilocal(profiler))%h(1:iv%info(profiler)%max_lev))
2010                allocate (iv%profiler(ilocal(profiler))%p(1:iv%info(profiler)%max_lev))
2011                allocate (iv%profiler(ilocal(profiler))%u(1:iv%info(profiler)%max_lev))
2012                allocate (iv%profiler(ilocal(profiler))%v(1:iv%info(profiler)%max_lev))
2013            end if
2014
2015             do i = 1, plink%nlevels_BUFR
2016                iv%profiler(ilocal(profiler))%h(i) = plink%platform_BUFR%each(i)%height
2017                iv%profiler(ilocal(profiler))%p(i) = plink%platform_BUFR%each(i)%p%inv
2018                iv%profiler(ilocal(profiler))%u(i) = plink%platform_BUFR%each(i)%u
2019                iv%profiler(ilocal(profiler))%v(i) = plink%platform_BUFR%each(i)%v
2020
2021             end do
2022
2023           case (571, 65)   ! SSM/I wind speed & TPW
2024           if (.not. use_ssmiretrievalobs) then
2025                    plink => plink%next
2026                    cycle reports3
2027                end if
2028            if (n==1)  ntotal(ssmi_rv) = ntotal(ssmi_rv) + 1
2029            if (outside) then
2030                    plink => plink%next
2031                    cycle reports3
2032                end if
2033             if ( thin_conv_opt(ssmi_rv) > no_thin ) then
2034                crit = tdiff
2035                call map2grids_conv(ssmi_rv,dlat_earth,dlon_earth,crit,nlocal(ssmi_rv),itx,1,itt,ilocal(ssmi_rv),iuse)
2036                 if ( .not. iuse ) then
2037                    plink => plink%next
2038                    cycle reports3
2039                   end if
2040             else
2041                nlocal(ssmi_rv) = nlocal(ssmi_rv) + 1
2042                ilocal(ssmi_rv) = nlocal(ssmi_rv)
2043             end if
2044
2045             platform_name ='FM-125 SSMI '
2046
2047             select case (plink%kx_BUFR)
2048             case ( 283)  ! wind speed
2049                do i = 1, plink%nlevels_BUFR
2050                   wpc = nint(plink%pco_BUFR(5,i))
2051                   ! if wpc == 1, UOB is set to zero, VOB is speed
2052                   iv%ssmi_rv(ilocal(ssmi_rv))%speed  = plink%platform_BUFR%each(i)%v
2053                   if ( wpc == 10 ) then
2054                      iv%ssmi_rv(ilocal(ssmi_rv))%speed%inv  = sqrt ( &
2055                         plink%platform_BUFR%each(i)%u%inv * plink%platform_BUFR%each(i)%u%inv + &
2056                         plink%platform_BUFR%each(i)%v%inv * plink%platform_BUFR%each(i)%v%inv  )
2057                   end if
2058                   iv%ssmi_rv(ilocal(ssmi_rv))%tpw    = plink%platform_BUFR%loc%pw
2059
2060                end do
2061            case ( 152 )  ! tpw
2062               do i = 1, plink%nlevels_BUFR
2063                  iv%ssmi_rv(ilocal(ssmi_rv))%speed  = plink%platform_BUFR%each(i)%u
2064                  iv%ssmi_rv(ilocal(ssmi_rv))%tpw    = plink%platform_BUFR%loc%pw
2065
2066                end do
2067           case default
2068               exit dup_loop3
2069
2070            end select
2071
2072           case default
2073             select case (plink%kx_BUFR)
2074             case (111 , 210)
2075           if (.not.use_bogusobs) then
2076                    plink => plink%next
2077                    cycle reports3
2078                end if
2079                if (n==1)  ntotal(bogus) = ntotal(bogus) + 1
2080                 if (outside) then
2081                    plink => plink%next
2082                    cycle reports3
2083                end if
2084                 if ( thin_conv_opt(bogus) > no_thin ) then
2085                   crit = tdiff
2086                   call map2grids_conv(bogus,dlat_earth,dlon_earth,crit,nlocal(bogus),itx,1,itt,ilocal(bogus),iuse)
2087                   if ( .not. iuse ) then
2088                      plink => plink%next
2089                      cycle reports3
2090                   end if
2091                else
2092                   nlocal(bogus) = nlocal(bogus) + 1
2093                   ilocal(bogus) = nlocal(bogus)
2094                end if
2095
2096                platform_name ='FM-135 TCBOG'
2097                if ( ilocal(bogus) == nlocal(bogus) ) then
2098                   allocate (iv%bogus(ilocal(bogus))%h(1:iv%info(bogus)%max_lev))
2099                   allocate (iv%bogus(ilocal(bogus))%p(1:iv%info(bogus)%max_lev))
2100                   allocate (iv%bogus(ilocal(bogus))%u(1:iv%info(bogus)%max_lev))
2101                   allocate (iv%bogus(ilocal(bogus))%v(1:iv%info(bogus)%max_lev))
2102                   allocate (iv%bogus(ilocal(bogus))%t(1:iv%info(bogus)%max_lev))
2103                   allocate (iv%bogus(ilocal(bogus))%q(1:iv%info(bogus)%max_lev))
2104                end if
2105
2106                do i = 1, plink%nlevels_BUFR
2107                   iv%bogus(ilocal(bogus))%h(i) = plink%platform_BUFR%each(i)%height
2108                   iv%bogus(ilocal(bogus))%p(i) = plink%platform_BUFR%each(i)%p%inv
2109                   iv%bogus(ilocal(bogus))%u(i) = plink%platform_BUFR%each(i)%u
2110                   iv%bogus(ilocal(bogus))%v(i) = plink%platform_BUFR%each(i)%v
2111                   iv%bogus(ilocal(bogus))%t(i) = plink%platform_BUFR%each(i)%t
2112                   iv%bogus(ilocal(bogus))%q(i) = plink%platform_BUFR%each(i)%q
2113
2114                end do
2115
2116                iv%bogus(ilocal(bogus))%slp    = plink%platform_BUFR%loc%slp
2117               case default
2118               exit dup_loop3
2119              end select
2120          end select
2121
2122          obs_index = fm_index(plink%fm_BUFR)
2123          iv%info(obs_index)%name(ilocal(obs_index))      = plink%platform_BUFR%info%name
2124          iv%info(obs_index)%platform(ilocal(obs_index))  = platform_name
2125          iv%info(obs_index)%id(ilocal(obs_index))        = plink%platform_BUFR%info%id
2126          iv%info(obs_index)%date_char(ilocal(obs_index)) = plink%platform_BUFR%info%date_char
2127          ! nlevels adjusted for some obs types so use that
2128          iv%info(obs_index)%levels(ilocal(obs_index))    = min(iv%info(obs_index)%max_lev, plink%nlevels_BUFR)
2129          iv%info(obs_index)%lat(:,ilocal(obs_index))     = plink%platform_BUFR%info%lat
2130          iv%info(obs_index)%lon(:,ilocal(obs_index))     = plink%platform_BUFR%info%lon
2131          iv%info(obs_index)%elv(ilocal(obs_index))       = plink%platform_BUFR%info%elv
2132          iv%info(obs_index)%pstar(ilocal(obs_index))     = plink%platform_BUFR%info%pstar
2133
2134          iv%info(obs_index)%slp(ilocal(obs_index))           = plink%platform_BUFR%loc%slp
2135          iv%info(obs_index)%pw(ilocal(obs_index))            = plink%platform_BUFR%loc%pw
2136          iv%info(obs_index)%x(:,ilocal(obs_index))           = plink%platform_BUFR%loc%x
2137          iv%info(obs_index)%y(:,ilocal(obs_index))           = plink%platform_BUFR%loc%y
2138          iv%info(obs_index)%i(:,ilocal(obs_index))           = plink%platform_BUFR%loc%i
2139          iv%info(obs_index)%j(:,ilocal(obs_index))           = plink%platform_BUFR%loc%j
2140          iv%info(obs_index)%dx(:,ilocal(obs_index))          = plink%platform_BUFR%loc%dx
2141          iv%info(obs_index)%dxm(:,ilocal(obs_index))         = plink%platform_BUFR%loc%dxm
2142          iv%info(obs_index)%dy(:,ilocal(obs_index))          = plink%platform_BUFR%loc%dy
2143          iv%info(obs_index)%dym(:,ilocal(obs_index))         = plink%platform_BUFR%loc%dym
2144          iv%info(obs_index)%proc_domain(:,ilocal(obs_index)) = plink%platform_BUFR%loc%proc_domain
2145
2146          iv%info(obs_index)%obs_global_index(ilocal(obs_index)) = iv%info(obs_index)%ntotal
2147          ! special case for sonde_sfc, duplicate sound info
2148          if (obs_index == sound) then
2149
2150             iv%info(sonde_sfc)%name(ilocal(sonde_sfc))      = plink%platform_BUFR%info%name
2151             iv%info(sonde_sfc)%platform(ilocal(sonde_sfc))  = platform_name
2152             iv%info(sonde_sfc)%id(ilocal(sonde_sfc))        = plink%platform_BUFR%info%id
2153             iv%info(sonde_sfc)%date_char(ilocal(sonde_sfc)) = plink%platform_BUFR%info%date_char
2154             iv%info(sonde_sfc)%levels(ilocal(sonde_sfc))    = 1
2155             iv%info(sonde_sfc)%lat(:,ilocal(sonde_sfc))     = plink%platform_BUFR%info%lat
2156             iv%info(sonde_sfc)%lon(:,ilocal(sonde_sfc))     = plink%platform_BUFR%info%lon
2157             iv%info(sonde_sfc)%elv(ilocal(sonde_sfc))       = plink%platform_BUFR%info%elv
2158             iv%info(sonde_sfc)%pstar(ilocal(sonde_sfc))     = plink%platform_BUFR%info%pstar
2159
2160             iv%info(sonde_sfc)%slp(ilocal(sonde_sfc))           = plink%platform_BUFR%loc%slp
2161             iv%info(sonde_sfc)%pw(ilocal(sonde_sfc))            = plink%platform_BUFR%loc%pw
2162             iv%info(sonde_sfc)%x(:,ilocal(sonde_sfc))           = plink%platform_BUFR%loc%x
2163             iv%info(sonde_sfc)%y(:,ilocal(sonde_sfc))           = plink%platform_BUFR%loc%y
2164             iv%info(sonde_sfc)%i(:,ilocal(sonde_sfc))           = plink%platform_BUFR%loc%i
2165             iv%info(sonde_sfc)%j(:,ilocal(sonde_sfc))           = plink%platform_BUFR%loc%j
2166             iv%info(sonde_sfc)%dx(:,ilocal(sonde_sfc))          = plink%platform_BUFR%loc%dx
2167             iv%info(sonde_sfc)%dxm(:,ilocal(sonde_sfc))         = plink%platform_BUFR%loc%dxm
2168             iv%info(sonde_sfc)%dy(:,ilocal(sonde_sfc))          = plink%platform_BUFR%loc%dy
2169             iv%info(sonde_sfc)%dym(:,ilocal(sonde_sfc))         = plink%platform_BUFR%loc%dym
2170             iv%info(sonde_sfc)%proc_domain(:,ilocal(sonde_sfc)) = plink%platform_BUFR%loc%proc_domain
2171
2172             iv%info(sonde_sfc)%obs_global_index(ilocal(sonde_sfc)) =iv%info(sonde_sfc)%ntotal
2173
2174          end if
2175
2176        end do dup_loop3
2177      end if !sort iv
2178
2179      plink => plink%next
2180      if ( .not. associated(plink) ) exit reports3
2181
2182    end do reports3
2183
2184    if (num_fgat_time >1 ) then
2185       do n = 1, num_ob_indexes
2186           iv%info(n)%ptotal(kk)=ntotal(n)
2187           iv%info(n)%plocal(kk)=nlocal(n)
2188       end do
2189    else
2190       do n = 1, num_ob_indexes
2191          ntotal(n)=iv%info(n)%ntotal
2192          nlocal(n)=iv%info(n)%nlocal
2193          iv%info(n)%ptotal(1)=iv%info(n)%ntotal
2194          iv%info(n)%plocal(1)=iv%info(n)%nlocal
2195        end do
2196     end if
2197
2198    ! thinning check
2199    if ( thin_conv ) then
2200       do n = 1, num_ob_indexes
2201
2202         if ( thin_conv_opt(n) <= no_thin ) cycle
2203
2204         if ( ntotal(n)>0 ) then
2205
2206             allocate ( in  (thinning_grid_conv(n)%itxmax) )
2207             allocate ( out (thinning_grid_conv(n)%itxmax) )
2208             do i = 1, thinning_grid_conv(n)%itxmax
2209                in(i) = thinning_grid_conv(n)%score_crit(i)
2210             end do
2211 #ifdef DM_PARALLEL
2212             ! Get minimum crit and associated processor index.
2213
2214             call mpi_reduce(in, out, thinning_grid_conv(n)%itxmax, true_mpi_real, mpi_min, root, comm, ierr)
2215             call wrf_dm_bcast_real (out, thinning_grid_conv(n)%itxmax)
2216 #else
2217             out = in
2218 #endif
2219             do i = 1, thinning_grid_conv(n)%itxmax
2220                if ( abs(out(i)-thinning_grid_conv(n)%score_crit(i)) > 1.0E-10 ) then
2221                   thinning_grid_conv(n)%ibest_obs(i) = 0
2222                else
2223                   if(thinning_grid_conv(n)%score_crit(i) < 9.99e6 ) &
2224                   iv%info(n)%thin_ntotal=  iv%info(n)%thin_ntotal + 1
2225                end if
2226             end do
2227
2228             deallocate( in  )
2229             deallocate( out )
2230
2231             do j = (1+tp(n)), nlocal(n)
2232                found = .false.
2233                do i = 1, thinning_grid_conv(n)%itxmax
2234                   if ( thinning_grid_conv(n)%ibest_obs(i) == j .and.         &
2235                        thinning_grid_conv(n)%score_crit(i) < 9.99e6 ) then
2236                      found = .true.
2237                      iv%info(n)%thin_nlocal =  iv%info(n)%thin_nlocal + 1
2238                      exit
2239                   end if
2240                end do
2241                if ( .not. found ) then
2242                   iv%info(n)%thinned(:,j) = .true.
2243                end if
2244             end do
2245
2246          tp(n)=nlocal(n)
2247          end if
2248       end do
2249        if (num_fgat_time >1 ) then
2250           do n = 1, num_ob_indexes
2251               iv%info(n)%thin_plocal(kk)=iv%info(n)%thin_nlocal
2252               iv%info(n)%thin_ptotal(kk)=iv%info(n)%thin_ntotal
2253           end do
2254        else
2255           do n = 1, num_ob_indexes
2256              iv%info(n)%thin_plocal(1)=iv%info(n)%thin_nlocal
2257              iv%info(n)%thin_ptotal(1)=iv%info(n)%thin_ntotal
2258            end do
2259         end if
2260
2261    end if  ! thin_conv
2262
2263  end do  !kk cycle
2264
2265 if ( thin_conv ) then
2266     do n = 1, num_ob_indexes
2267         if ( thin_conv_opt(n) <= no_thin ) cycle
2268         if ( ntotal(n)>0 ) then
2269             if ( nlocal(n) > 0 ) then
2270                if ( ANY(iv%info(n)%thinned(:,:)) ) then
2271                   call da_set_obs_missing(iv,n)  ! assign missing values to those thinned=true data
2272                end if
2273             end if
2274         end if
2275     end do
2276 end if
2277
2278 write(unit=message(1),fmt='(A,4(1x,i7))') &
2279       'da_read_obs_bufr: num_report, num_outside_all, num_outside_time, num_thinned: ', &
2280       num_report, num_outside_all, num_outside_time, num_thinned
2281    call da_message(message(1:1))
2282
2283
2284 ! Release the linked list
2285 plink => head
2286
2287 DO WHILE ( ASSOCIATED (plink) )
2288    head => plink%next
2289    DEALLOCATE ( plink%platform_BUFR%each )
2290    DEALLOCATE ( plink )
2291    plink => head
2292 ENDDO
2293
2294 NULLIFY (head)
2295
2296    if (trace_use) call da_trace_exit("da_read_obs_bufr")
2297 #else
2298    call da_error(__FILE__,__LINE__,(/"must compile with BUFR library"/))
2299 #endif
2300
2301 end subroutine da_read_obs_bufr
2302
2303
2304