var/da/da_obs_io/da_read_obs_bufrgpsro.inc

   1 subroutine da_read_obs_bufrgpsro (iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Read NCEP GPSRO BUFR observation file for input to wrfvar
   5    !---------------------------------------------------------------------------
   6    !   METHOD: use F90 sequantial data structure to avoid read file twice
   7    !            so that da_scan_obs_bufr is not necessary any more.
   8    !            1. read gpsro 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/12/17  - F90 sequantial structure  Peng Xiu
  14    !
  15    !----------------------------------------------------------------------------
  16
  17    use da_control
  18
  19    implicit none
  20
  21    type (iv_type),             intent(inout) :: iv
  22
  23 #ifdef BUFR
  24
  25    real,    parameter   :: r8bfms = 9.0D08  ! BUFR missing value threshold
  26    integer, parameter   :: maxlevs = 500
  27    integer              :: iunit, iost, idate, iret, nlev1, nlev2,k,i,ii
  28    integer              :: num_report, num_outside_all, num_outside_time
  29    integer              :: iyear,imonth,iday,ihour,imin
  30    integer              :: ntotal, nlocal, nlev, ref_qc
  31    real*8               :: obs_time
  32    real*8               :: hdr(10)
  33    real                 :: ntotal_ifgat(0:num_fgat_time)
  34    real*8               :: rdata1(25,maxlevs), rdata2(25,maxlevs)
  35    real                 :: height, ref_data, ref_error
  36    character(len=8)     :: subset
  37    character(len=80)    :: hdstr
  38    character(len=12)    :: filename
  39    logical              :: outside, outside_all
  40    type(info_type)      :: info
  41    type(model_loc_type) :: loc
  42    character(len=5)     :: id
  43    character(len=19)    :: date_char
  44    integer              :: ifgat, kk, num_p
  45    integer              :: err_opt   ! 1: WRF-Var/obsproc, 2: GSI
  46    real                 :: erh90, erh0, err90, err0
  47    integer(i_kind)      :: ireadns
  48    integer              :: num_bufr(7), numbufr, ibufr
  49    type datalink_gpsro              !for gpsro data reading
  50        type (info_type)         :: info
  51        type(model_loc_type)     :: loc
  52        type (gpsref_type)       :: gpsref
  53        integer                  :: ifgat
  54        type (field_type)        :: slp
  55        type (field_type)        :: pw
  56        integer                  :: obs_global_index
  57        type(datalink_gpsro), pointer :: next
  58    end type datalink_gpsro
  59    type(datalink_gpsro),pointer  :: head=>null(), plink=>null()
  60    integer(i_kind)      :: iprof
  61
  62    if (trace_use) call da_trace_entry("da_read_obs_bufrgosro")
  63
  64    iprof  = 0
  65    nlocal = 0
  66    ntotal = 0
  67    num_report       = 0
  68    num_outside_all  = 0
  69    num_outside_time = 0
  70    num_p=0
  71    ntotal_ifgat(0:num_fgat_time)=0
  72
  73    ! open file
  74    !  ---------
  75    num_bufr(:)=0
  76    numbufr=0
  77    if (num_fgat_time>1) then
  78       do i=1,7
  79         call da_get_unit(iunit)
  80         write(filename,fmt='(A,I1,A)') 'gpsro0',i,'.bufr'
  81         open(unit   = iunit, FILE   = trim(filename),iostat =  iost, form = 'unformatted', STATUS = 'OLD')
  82         if (iost == 0) then
  83            numbufr=numbufr+1
  84            num_bufr(numbufr)=i
  85         else
  86            close (iunit)
  87         end if
  88         call da_free_unit(iunit)
  89       end do
  90    else
  91       numbufr=1
  92    end if
  93
  94    if (numbufr==0) numbufr=1
  95
  96 bufrfile:  do ibufr=1,numbufr
  97    if (num_fgat_time==1) then
  98        filename='gpsro.bufr'
  99    else
 100        if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then
 101           filename='gpsro.bufr'
 102       else
 103           write(filename,fmt='(A,I1,A)') 'gpsro0',num_bufr(ibufr),'.bufr'
 104       end if
 105    end if
 106
 107    ! Use specified unit number, so we can control its endian format in environment.
 108    iunit = 96
 109    call closbf(iunit)
 110    open(unit   = iunit, FILE   = trim(filename), &
 111       iostat =  iost, form = 'unformatted', STATUS = 'OLD')
 112    if (iost /= 0) then
 113       write(unit=message(1),fmt='(A,I5,A)') &
 114          "Error",iost," opening PREPBUFR obs file "//trim(filename)
 115       call da_warning(__FILE__,__LINE__,message(1:1))
 116       if (trace_use) call da_trace_exit("da_read_obs_bufrgpsro")
 117       return
 118    end if
 119
 120    !--------------------------------
 121    ! open bufr file then check date
 122    !--------------------------------
 123    call openbf(iunit,'IN',iunit)
 124    call datelen(10)
 125    call readmg(iunit,subset,idate,iret)
 126    if ( iret /= 0 ) then
 127       write(unit=message(1),fmt='(A,I5,A)') &
 128          "Error",iret," reading GPSRO BUFR obs file "//trim(filename)
 129       call da_warning(__FILE__,__LINE__,message(1:1))
 130       call closbf(iunit)
 131       close(iunit)
 132       if (trace_use) call da_trace_exit("da_read_obs_bufrgpsro")
 133       return
 134    end if
 135    write(unit=message(1),fmt='(a,i10)') 'GPSRO BUFR file date is: ', idate
 136    call da_message(message(1:1))
 137    rewind(iunit)
 138
 139    hdstr = 'YEAR MNTH DAYS HOUR MINU PCCF ELRC SAID PTID GEODU'
 140
 141    reports: do while ( ireadns(iunit,subset,idate) == 0 )
 142
 143       num_report = num_report + 1
 144
 145       call ufbint(iunit,hdr,10,1,iret,hdstr)
 146
 147       iyear  = int(hdr(1))
 148       imonth = int(hdr(2))
 149       iday   = int(hdr(3))
 150       ihour  = int(hdr(4))
 151       imin   = int(hdr(5))
 152
 153       write(id, '(i3.3,i2.2)') int(hdr(8)), int(hdr(9)) ! construct id using SAID and PTID
 154       write(date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 155          iyear, '-', imonth, '-', iday, '_', ihour, ':', imin, ':', 0
 156
 157       ! check date
 158       call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 159       if (obs_time < time_slots(0) .or.  &
 160           obs_time >= time_slots(num_fgat_time)) then
 161          num_outside_time = num_outside_time + 1
 162          if ( print_detail_obs ) then
 163             write(unit=stderr,fmt='(a,1x,i4.4,4i2.2,a)')  &
 164                info%id(1:5),iyear,imonth,iday,ihour,imin, '  -> outside_time'
 165          end if
 166          cycle reports
 167       end if
 168
 169       if ( hdr(6) < 100.0 ) then   ! check percentage of confidence PCCF
 170          cycle reports
 171       end if
 172
 173       call ufbseq(iunit,rdata1,25,maxlevs,nlev1,'ROSEQ1')  ! RAOC PROFILE LOCATIONS SEQUENCE
 174       call ufbseq(iunit,rdata2,25,maxlevs,nlev2,'ROSEQ3')  ! RAOC HEIGHT/REFRACTIVITY SEQUENCE
 175
 176       if ( nlev1 /= nlev2 ) then
 177          cycle reports
 178       end if
 179
 180  !--------  determine FGAT index ifgat
 181
 182          do ifgat=1,num_fgat_time
 183             if (obs_time >= time_slots(ifgat-1) .and.  &
 184                 obs_time  < time_slots(ifgat)) exit
 185          end do
 186
 187       iprof = iprof + 1
 188
 189       lev_loop: do k = 1, nlev1
 190
 191          info%lat  = rdata1(1,k)
 192          info%lon  = rdata1(2,k)
 193
 194          ! gpsro.bufr contains missing longitude, occasionally
 195          if ( info%lat > r8bfms .or. info%lon > r8bfms ) then
 196             cycle lev_loop
 197          end if
 198
 199          height    = rdata2(1,k)
 200          ref_data  = rdata2(2,k)
 201
 202          ! check for missing data
 203          if ( height > r8bfms .or. ref_data > r8bfms ) then
 204             cycle lev_loop
 205          end if
 206
 207          ref_qc    = 0                ! initialized to be good
 208          ref_error = ref_data * 0.01
 209
 210          ! check loc
 211          info%lat = max(info%lat, -89.95)
 212          info%lat = min(info%lat,  89.95)
 213          call da_llxy(info, loc, outside, outside_all)
 214          if ( outside_all ) then
 215             num_outside_all = num_outside_all + 1
 216             if ( print_detail_obs ) then
 217                write(unit=stderr,fmt='(a,2(1x,f8.3),a)')  &
 218                   id(1:5), info%lat, info%lon, '  -> outside_domain'
 219             end if
 220             cycle lev_loop
 221          end if
 222          ntotal = ntotal + 1
 223          ntotal_ifgat(ifgat)=ntotal_ifgat(ifgat)+1
 224          if ( outside ) then
 225             cycle lev_loop
 226          end if
 227          nlocal = nlocal + 1
 228
 229          ! check height, only keep data below certain height (default is 30km)
 230          if ( height > top_km_gpsro*1000.0 .or. &
 231               height < bot_km_gpsro*1000.0 ) then
 232             ref_qc = -77
 233          end if
 234
 235          err_opt = 1
 236          ! observation errors  WRF-Var/obsproc
 237          if ( err_opt == 1 ) then
 238             if ( height >= 12000.0 ) then
 239                ref_error = ref_error * 0.3
 240             else
 241                erh90 = (0.3-1.5)*(height-12000.0)/(12000.0-0.0) + 0.3
 242                if ( height >= 5500.0 ) then
 243                   erh0 = (0.3-1.3)*(height-12000.0)/(12000.0-5500.0) + 0.3
 244                else if ( height >= 2500.0) then
 245                   erh0 = (1.3-2.5)*(height-5500.0)/(5500.0-2500.0) + 1.3
 246                else
 247                   erh0 = 2.5
 248                end if
 249                err90 = ref_error * erh90
 250                err0  = ref_error * erh0
 251                ref_error = err90 - (1.0-abs(info%lat)/90.0)*(err90-err0)
 252             end if
 253          end if
 254
 255          ! observation errors  GSI_Q1FY09,  Kuo et al. 2003
 256          if ( err_opt == 2 ) then
 257             if ( (info%lat >= -30.0) .and. (info%lat <= 30.0) ) then   ! tropics
 258                if ( (height >= 7000.0) .and. (height <= 31000.0) ) then
 259                   ref_error = ref_error*(0.1125+(1.25e-5*height))
 260                else if ( height > 31000.0 ) then
 261                   ref_error = ref_error*0.5
 262                else if ( height < 7000.0  ) then
 263                   ref_error = ref_error*(3.0-(4.0e-4*height))
 264                else
 265                   write(unit=message(1),fmt='(a,f8.1,a,f8.2)') 'unable to process with height = ', &
 266                      height, ' at lat = ', info%lat
 267                   call da_error(__FILE__,__LINE__,message(1:1))
 268                end if
 269             else   ! mid-latitudes
 270                if ( (height >= 5000.0) .and. (height <= 25000.0) ) then
 271                   ref_error = ref_error*0.3
 272                else if ( (height >= 25000.0) .and. (height <= 31000.0) ) then
 273                   ref_error = ref_error*(-3.45+(1.5e-4*height))
 274                else if ( height > 31000.0 ) then
 275                   ref_error = ref_error*1.2
 276                else if ( height < 5000.0 ) then
 277                   ref_error = ref_error*(0.75-(9.0e-5*height))
 278                else
 279                   write(unit=message(1),fmt='(a,f8.1,a,f8.2)') 'unable to process with height = ', &
 280                      height, ' at lat = ', info%lat
 281                   call da_error(__FILE__,__LINE__,message(1:1))
 282                end if
 283             end if
 284          end if
 285
 286          !write(info%name, '(a,i6.6,a,a)') 'NCEP_GPSRO_', nlocal, '_', date_char
 287          write(info%name, '(a,i6.6,a,a)') 'NCEP_GPSRO_', iprof, '_', date_char
 288
 289          if ( print_detail_obs ) then
 290             write(unit=stdout,fmt='(a,1x,a,1x,i4.4,4i2.2,2f8.2,f8.1,f8.2,i3,f9.5)')  &
 291                info%name,id(1:5),iyear,imonth,iday,ihour,imin, &
 292                info%lat,info%lon,height,ref_data,ref_qc,ref_error
 293          end if
 294
 295          if (.not. associated(head)) then
 296              nullify ( head )
 297              allocate ( head )
 298              nullify ( head%next )
 299              plink => head
 300          else
 301              allocate ( plink%next )
 302              plink => plink%next
 303              nullify ( plink%next )
 304          end if
 305
 306          plink%info%name      = info%name
 307          plink%info%platform  = 'FM-116 GPSRF'
 308          plink%info%id        = id
 309          plink%info%date_char = date_char
 310          plink%info%levels    = 1              ! each level is treated as separate obs
 311          plink%info%elv       = 0.0            ! not used
 312          plink%info%lat     = info%lat
 313          plink%info%lon     = info%lon
 314
 315          plink%loc%x       = loc%x
 316          plink%loc%y       = loc%y
 317          plink%loc%i       = loc%i
 318          plink%loc%j       = loc%j
 319          plink%loc%dx      = loc%dx
 320          plink%loc%dxm     = loc%dxm
 321          plink%loc%dy      = loc%dy
 322          plink%loc%dym     = loc%dym
 323
 324          plink%slp%inv   = missing_r
 325          plink%slp%qc    = missing_data
 326          plink%slp%error = xmiss
 327          plink%pw%inv    = missing_r
 328          plink%pw%qc     = missing_data
 329          plink%pw%error  = xmiss
 330
 331          plink%obs_global_index = ntotal
 332
 333          nlev = 1
 334          allocate (plink%gpsref%h  (1:nlev))
 335          allocate (plink%gpsref%ref(1:nlev))
 336          allocate (plink%gpsref%p  (1:nlev))
 337          allocate (plink%gpsref%t  (1:nlev))
 338          allocate (plink%gpsref%q  (1:nlev))
 339          do i = 1, nlev
 340             plink%gpsref%h(i)         = height
 341             plink%gpsref%ref(i)%inv   = ref_data
 342             plink%gpsref%ref(i)%qc    = ref_qc
 343             plink%gpsref%ref(i)%error = ref_error
 344             plink%gpsref%p(i)%inv     = missing_r
 345             plink%gpsref%p(i)%qc      = missing_data
 346             plink%gpsref%p(i)%error   = xmiss
 347             plink%gpsref%t(i)%inv     = missing_r
 348             plink%gpsref%t(i)%qc      = missing_data
 349             plink%gpsref%t(i)%error   = xmiss
 350             plink%gpsref%q(i)%inv     = missing_r
 351             plink%gpsref%q(i)%qc      = missing_data
 352             plink%gpsref%q(i)%error   = xmiss
 353          end do
 354          plink%ifgat=ifgat
 355          num_p=num_p+1
 356
 357          end do lev_loop
 358         end do reports
 359
 360         call closbf(iunit)
 361         close(iunit)
 362
 363 end do bufrfile
 364
 365          iv%info(gpsref)%ntotal=ntotal
 366          iv%info(gpsref)%nlocal=nlocal
 367          iv%info(gpsref)%max_lev=1
 368
 369 !allocate iv
 370
 371        if (iv%info(gpsref)%nlocal > 0) allocate(iv%gpsref(1:iv%info(gpsref)%nlocal))
 372
 373        if (iv%info(gpsref)%nlocal > 0) then
 374          allocate (iv%info(gpsref)%name(iv%info(gpsref)%nlocal))
 375          allocate (iv%info(gpsref)%platform(iv%info(gpsref)%nlocal))
 376          allocate (iv%info(gpsref)%id(iv%info(gpsref)%nlocal))
 377          allocate (iv%info(gpsref)%date_char(iv%info(gpsref)%nlocal))
 378          allocate (iv%info(gpsref)%levels(iv%info(gpsref)%nlocal))
 379          allocate (iv%info(gpsref)%lat(iv%info(gpsref)%max_lev,iv%info(gpsref)%nlocal))
 380          allocate (iv%info(gpsref)%lon(iv%info(gpsref)%max_lev,iv%info(gpsref)%nlocal))
 381          allocate (iv%info(gpsref)%elv(iv%info(gpsref)%nlocal))
 382          allocate (iv%info(gpsref)%pstar(iv%info(gpsref)%nlocal))
 383          allocate (iv%info(gpsref)%slp(iv%info(gpsref)%nlocal))
 384          allocate (iv%info(gpsref)%pw(iv%info(gpsref)%nlocal))
 385          allocate (iv%info(gpsref)%x  (kms:kme,iv%info(gpsref)%nlocal))
 386          allocate (iv%info(gpsref)%y  (kms:kme,iv%info(gpsref)%nlocal))
 387          allocate (iv%info(gpsref)%i  (kms:kme,iv%info(gpsref)%nlocal))
 388          allocate (iv%info(gpsref)%j  (kms:kme,iv%info(gpsref)%nlocal))
 389          allocate (iv%info(gpsref)%dx (kms:kme,iv%info(gpsref)%nlocal))
 390          allocate (iv%info(gpsref)%dxm(kms:kme,iv%info(gpsref)%nlocal))
 391          allocate (iv%info(gpsref)%dy (kms:kme,iv%info(gpsref)%nlocal))
 392          allocate (iv%info(gpsref)%dym(kms:kme,iv%info(gpsref)%nlocal))
 393          allocate (iv%info(gpsref)%k  (iv%info(gpsref)%max_lev,iv%info(gpsref)%nlocal))
 394          allocate (iv%info(gpsref)%dz (iv%info(gpsref)%max_lev,iv%info(gpsref)%nlocal))
 395          allocate (iv%info(gpsref)%dzm(iv%info(gpsref)%max_lev,iv%info(gpsref)%nlocal))
 396          allocate (iv%info(gpsref)%zk (iv%info(gpsref)%max_lev,iv%info(gpsref)%nlocal))
 397          allocate (iv%info(gpsref)%proc_domain(iv%info(gpsref)%max_lev,iv%info(gpsref)%nlocal))
 398          allocate (iv%info(gpsref)%thinned(iv%info(gpsref)%max_lev,iv%info(gpsref)%nlocal))
 399          allocate (iv%info(gpsref)%obs_global_index(iv%info(gpsref)%nlocal))
 400          iv%info(gpsref)%proc_domain(:,:)  = .false.
 401          iv%info(gpsref)%thinned(:,:)      = .false.
 402          iv%info(gpsref)%zk(:,:)           = missing_r
 403         end if
 404 !
 405 !sort for 4d var
 406          nlocal=0
 407          do kk=1,num_fgat_time
 408          plink => head
 409          iv%info(gpsref)%ptotal(kk)=0
 410
 411          reports2: do ii=1,num_p
 412
 413          if (plink%ifgat /= kk) then  !sort iv
 414             plink => plink%next
 415             cycle reports2
 416          else
 417          nlocal=nlocal+1
 418          iv%info(gpsref)%name(nlocal)      = plink%info%name
 419          iv%info(gpsref)%platform(nlocal)  = plink%info%platform
 420          iv%info(gpsref)%id(nlocal)        = plink%info%id
 421          iv%info(gpsref)%date_char(nlocal) = plink%info%date_char
 422          iv%info(gpsref)%levels(nlocal)    = plink%info%levels              ! each level is treated as separate obs
 423          iv%info(gpsref)%elv(nlocal)       = plink%info%elv            ! not used
 424          iv%info(gpsref)%lat(:,nlocal)     = plink%info%lat
 425          iv%info(gpsref)%lon(:,nlocal)     = plink%info%lon
 426
 427          iv%info(gpsref)%x(:,nlocal)       = plink%loc%x
 428          iv%info(gpsref)%y(:,nlocal)       = plink%loc%y
 429          iv%info(gpsref)%i(:,nlocal)       = plink%loc%i
 430          iv%info(gpsref)%j(:,nlocal)       = plink%loc%j
 431          iv%info(gpsref)%dx(:,nlocal)      = plink%loc%dx
 432          iv%info(gpsref)%dxm(:,nlocal)     = plink%loc%dxm
 433          iv%info(gpsref)%dy(:,nlocal)      = plink%loc%dy
 434          iv%info(gpsref)%dym(:,nlocal)     = plink%loc%dym
 435
 436          iv%info(gpsref)%slp(nlocal)%inv   = plink%slp%inv
 437          iv%info(gpsref)%slp(nlocal)%qc    = plink%slp%qc
 438          iv%info(gpsref)%slp(nlocal)%error = plink%slp%error
 439          iv%info(gpsref)%pw(nlocal)%inv    = plink%pw%inv
 440          iv%info(gpsref)%pw(nlocal)%qc     = plink%pw%qc
 441          iv%info(gpsref)%pw(nlocal)%error  = plink%pw%error
 442
 443          iv%info(gpsref)%obs_global_index(nlocal) = plink%obs_global_index
 444
 445          nlev = 1
 446          allocate (iv%gpsref(nlocal)%h  (1:nlev))
 447          allocate (iv%gpsref(nlocal)%ref(1:nlev))
 448          allocate (iv%gpsref(nlocal)%p  (1:nlev))
 449          allocate (iv%gpsref(nlocal)%t  (1:nlev))
 450          allocate (iv%gpsref(nlocal)%q  (1:nlev))
 451
 452          do i = 1, nlev
 453             iv%gpsref(nlocal)%h(i)         = plink%gpsref%h(i)
 454             iv%gpsref(nlocal)%ref(i)%inv   = plink%gpsref%ref(i)%inv
 455             iv%gpsref(nlocal)%ref(i)%qc    = plink%gpsref%ref(i)%qc
 456             iv%gpsref(nlocal)%ref(i)%error = plink%gpsref%ref(i)%error
 457
 458             iv%gpsref(nlocal)%p(i)%inv     = plink%gpsref%p(i)%inv
 459             iv%gpsref(nlocal)%p(i)%qc      = plink%gpsref%p(i)%qc
 460             iv%gpsref(nlocal)%p(i)%error   = plink%gpsref%p(i)%error
 461             iv%gpsref(nlocal)%t(i)%inv     = plink%gpsref%t(i)%inv
 462             iv%gpsref(nlocal)%t(i)%qc      = plink%gpsref%t(i)%qc
 463             iv%gpsref(nlocal)%t(i)%error   = plink%gpsref%t(i)%error
 464             iv%gpsref(nlocal)%q(i)%inv     = plink%gpsref%q(i)%inv
 465             iv%gpsref(nlocal)%q(i)%qc      = plink%gpsref%q(i)%qc
 466             iv%gpsref(nlocal)%q(i)%error   = plink%gpsref%q(i)%error
 467          end do
 468        end if
 469        plink => plink%next
 470       end do reports2
 471
 472       ntotal_ifgat(kk)=ntotal_ifgat(kk)+ntotal_ifgat(kk-1)
 473       iv%info(gpsref)%ptotal(kk)=ntotal_ifgat(kk)
 474       iv%info(gpsref)%plocal(kk)=nlocal
 475       ! thinning is not applied to GPSRO BUFR, simply make a copy from total number
 476       iv%info(gpsref)%thin_ptotal(kk)=ntotal_ifgat(kk)
 477       iv%info(gpsref)%thin_plocal(kk)=nlocal
 478    end do
 479
 480    write(unit=message(1),fmt='(A,3(1x,i7))') &
 481       'da_read_obs_bufrgpsro: num_report, num_outside_all, num_outside_time: ', &
 482       num_report, num_outside_all, num_outside_time
 483    call da_message(message(1:1))
 484
 485    if ( nlocal /= iv%info(gpsref)%nlocal ) then
 486       call da_error(__FILE__,__LINE__,(/"numbers mismatch between scanning and reading NCEP GSPRO BUFR file"/))
 487    end if
 488
 489    ! Release the linked list
 490    plink => head
 491
 492    DO WHILE ( ASSOCIATED (plink) )
 493       head => plink%next
 494       deallocate (plink%gpsref%h)
 495       deallocate (plink%gpsref%ref)
 496       deallocate (plink%gpsref%p)
 497       deallocate (plink%gpsref%t)
 498       deallocate (plink%gpsref%q)
 499       deallocate ( plink )
 500       plink => head
 501    ENDDO
 502
 503    NULLIFY (head)
 504
 505    if (trace_use) call da_trace_exit("da_read_obs_bufrgpsro")
 506 #else
 507    call da_error(__FILE__,__LINE__,(/"must compile with BUFR library"/))
 508 #endif
 509
 510 end subroutine da_read_obs_bufrgpsro