var/da/da_obs_io/da_read_obs_bufrgpsro_eph.inc

   1 subroutine da_read_obs_bufrgpsro_eph (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 reading the file twice
   7    !            1. read gpsro data in sequential data structure
   8    !            2. do gross QC check
   9    !            3. assign sequential data structure to innovation structure
  10    !                  and deallocate sequential data structure
  11    !
  12    !  HISTORY: 2009/12/17  - F90 sequantial structure  Peng Xiu
  13    !           Nov. 2015   - added the option for the assimilation of
  14    !                         gps excess phase which requires the input data
  15    !                         as profiles
  16    !
  17    !----------------------------------------------------------------------------
  18
  19    use da_control
  20
  21    implicit none
  22
  23    type (iv_type),             intent(inout) :: iv
  24
  25 #ifdef BUFR
  26
  27    real,    parameter   :: r8bfms = 9.0D08  ! BUFR missing value threshold
  28    integer, parameter   :: maxlevs = 500
  29    integer              :: iunit, iost, idate, iret, nlev1, nlev2,k,i, nlev_bad
  30    integer              :: num_report, num_outside_all, num_outside_time
  31    integer              :: iyear,imonth,iday,ihour,imin
  32    integer              :: ntotal, nlocal, nlev, ref_qc, iobs
  33    real*8               :: obs_time
  34    real*8               :: hdr(12)
  35    real                 :: ntotal_ifgat(0:num_fgat_time)
  36    real*8               :: rdata1(25,maxlevs), rdata2(25,maxlevs)
  37    real                 :: height, ref_data, ref_error
  38    real                 :: azim
  39    character(len=8)     :: subset
  40    character(len=80)    :: hdstr
  41    logical              :: outside, outside_all
  42    type(info_type)      :: info !for input to llxy
  43    type(model_loc_type) :: loc  !for output from llxy
  44    character(len=5)     :: id
  45    character(len=19)    :: date_char
  46    character(len=12)    :: platform_name
  47    integer              :: ifgat, kk
  48    integer              :: err_opt   ! 1: WRF-Var/obsproc, 2: GSI
  49    real                 :: erh90, erh0, err90, err0
  50    integer(i_kind)      :: ireadns
  51    type datalink_gpsro  !for gpsro data reading
  52        type (info_type)         :: info
  53        type (model_loc_type)    :: loc
  54        type (gpsref_type)       :: gpsref
  55        integer                  :: ifgat
  56        real                     :: rfict   !to be used for gpseph
  57        real, allocatable        :: azim(:) !to be used for gpseph
  58        real, allocatable        :: lat(:)  !to be used for gpseph
  59        real, allocatable        :: lon(:)  !to be used for gpseph
  60        integer                  :: obs_global_index
  61        type(datalink_gpsro), pointer :: next
  62    end type datalink_gpsro
  63    type(datalink_gpsro),pointer  :: head=>null(), plink=>null()
  64    type(ob_in_mean_h)   :: pseudo_ob
  65    integer              :: iprof
  66    integer              :: obs_index
  67    integer              :: lowest_level
  68    integer,parameter    :: nfile_max = 8
  69    integer              :: nfile, ifile, max_lev
  70    character(80)        :: fname(nfile_max)
  71    character(80)        :: infile_prefix, fname_tmp
  72    logical              :: fexist
  73    logical              :: treat_as_profile
  74
  75    if (trace_use) call da_trace_entry("da_read_obs_bufrgpsro_eph")
  76
  77    treat_as_profile = .false.
  78    platform_name    = 'FM-116 GPSRF'
  79    if ( use_gpsephobs ) then
  80       treat_as_profile = .true.
  81       platform_name    = 'FM-118 GPSEP'
  82    end if
  83
  84    max_lev = 1
  85    infile_prefix = 'gpsro'
  86    iprof  = 0
  87    nlocal = 0
  88    ntotal = 0
  89    num_report       = 0
  90    num_outside_all  = 0
  91    num_outside_time = 0
  92    ntotal_ifgat(0:num_fgat_time)=0
  93
  94    ! checking for input files
  95    nfile    = 0  !initialize
  96    fname(:) = '' !initialize
  97    ! first check if gpsro.bufr is available
  98    fname_tmp = trim(infile_prefix)//'.bufr'
  99    inquire (file=fname_tmp, exist=fexist)
 100    if ( fexist ) then
 101       nfile = 1
 102       fname(nfile)  = fname_tmp
 103    else
 104       ! check if gpsro-0x.bufr is available for multiple input files
 105       ! here 0x is the input file sequence number
 106       ! do not confuse it with fgat time slot index
 107       do i = 1, nfile_max
 108          write(fname_tmp,fmt='(A,A,I2.2,A)') trim(infile_prefix),'-',i,'.bufr'
 109          inquire (file=fname_tmp, exist=fexist)
 110          if ( fexist ) then
 111             nfile = nfile + 1
 112             fname(nfile)  = fname_tmp
 113          else
 114             exit
 115          end if
 116       end do
 117    end if
 118
 119    if ( nfile == 0 ) then
 120       call da_warning(__FILE__,__LINE__, &
 121          (/"No valid NCEP BUFR gpsro.bufr or gpsro-01.bufr file found."/))
 122       if (trace_use) call da_trace_exit("da_read_obs_bufrgpsro_eph")
 123       return
 124    end if
 125
 126 bufrfile:  do ifile = 1, nfile
 127
 128    ! Use specified unit number, so we can control its endian format in environment.
 129    iunit = 96
 130    call closbf(iunit)
 131    open(unit=iunit, file=trim(fname(ifile)), &
 132       iostat=iost, form='unformatted', status='OLD')
 133    if (iost /= 0) then
 134       write(unit=message(1),fmt='(A,I5,A)') &
 135          "Error",iost," opening PREPBUFR obs file "//trim(fname(ifile))
 136       call da_warning(__FILE__,__LINE__,message(1:1))
 137       if (trace_use) call da_trace_exit("da_read_obs_bufrgpsro_eph")
 138       return
 139    end if
 140
 141    !--------------------------------
 142    ! open bufr file then check date
 143    !--------------------------------
 144    call openbf(iunit,'IN',iunit)
 145    call datelen(10)
 146    call readmg(iunit,subset,idate,iret)
 147    if ( iret /= 0 ) then
 148       write(unit=message(1),fmt='(A,I5,A)') &
 149          "Error",iret," reading GPSRO BUFR obs file "//trim(fname(ifile))
 150       call da_warning(__FILE__,__LINE__,message(1:1))
 151       call closbf(iunit)
 152       close(iunit)
 153       if (trace_use) call da_trace_exit("da_read_obs_bufrgpsro_eph")
 154       return
 155    end if
 156    write(unit=message(1),fmt='(a,i10)') 'GPSRO BUFR file date is: ', idate
 157    call da_message(message(1:1))
 158    rewind(iunit)
 159
 160    hdstr = 'YEAR MNTH DAYS HOUR MINU PCCF ELRC SAID PTID GEODU CLATH CLONH'
 161
 162    reports: do while ( ireadns(iunit,subset,idate) == 0 )
 163
 164       num_report = num_report + 1
 165
 166       call ufbint(iunit,hdr,12,1,iret,hdstr)
 167
 168       iyear  = int(hdr(1))
 169       imonth = int(hdr(2))
 170       iday   = int(hdr(3))
 171       ihour  = int(hdr(4))
 172       imin   = int(hdr(5))
 173
 174       write(id, '(i3.3,i2.2)') int(hdr(8)), int(hdr(9)) ! construct id using SAID and PTID
 175       write(date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 176          iyear, '-', imonth, '-', iday, '_', ihour, ':', imin, ':', 0
 177
 178       ! check date
 179       call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 180       if (obs_time < time_slots(0) .or.  &
 181           obs_time >= time_slots(num_fgat_time)) then
 182          num_outside_time = num_outside_time + 1
 183          if ( print_detail_obs ) then
 184             write(unit=stderr,fmt='(a,1x,i4.4,4i2.2,a)')  &
 185                info%id(1:5),iyear,imonth,iday,ihour,imin, '  -> outside_time'
 186          end if
 187          cycle reports
 188       end if
 189
 190       ! determine FGAT index ifgat
 191       do ifgat=1,num_fgat_time
 192          if (obs_time >= time_slots(ifgat-1) .and.  &
 193              obs_time  < time_slots(ifgat)) exit
 194       end do
 195
 196       if ( hdr(6) < 100.0 ) then   ! check percentage of confidence PCCF
 197          cycle reports
 198       end if
 199
 200       call ufbseq(iunit,rdata1,25,maxlevs,nlev1,'ROSEQ1')  ! RAOC PROFILE LOCATIONS SEQUENCE
 201       call ufbseq(iunit,rdata2,25,maxlevs,nlev2,'ROSEQ3')  ! RAOC HEIGHT/REFRACTIVITY SEQUENCE
 202
 203       if ( nlev1 /= nlev2 ) then
 204          cycle reports
 205       end if
 206
 207       if ( treat_as_profile ) then
 208          !only one lat/lon (perigee point at occultation point)
 209          !for each report that contains multiple levels of data
 210          info%lat  = hdr(11)
 211          info%lon  = hdr(12)
 212          ! gpsro.bufr contains missing longitude, occasionally
 213          if ( info%lat > r8bfms .or. info%lon > r8bfms ) then
 214             cycle reports
 215          end if
 216          ! check loc
 217          info%lat = max(info%lat, -89.95)
 218          info%lat = min(info%lat,  89.95)
 219          call da_llxy(info, loc, outside, outside_all)
 220          if ( outside_all ) then
 221             num_outside_all = num_outside_all + 1
 222             if ( print_detail_obs ) then
 223                write(unit=stderr,fmt='(a,2(1x,f8.3),a)')  &
 224                   id(1:5), info%lat, info%lon, '  -> outside_domain'
 225             end if
 226             cycle reports
 227          end if
 228
 229          ! check all levels
 230          nlev_bad = 0
 231          do k = 1, nlev1
 232             !info%lat  = rdata1(1,k)
 233             !info%lon  = rdata1(2,k)
 234             !azim      = rdata1(3,k)
 235             !height    = rdata2(1,k)
 236             !ref_data  = rdata2(2,k)
 237             if ( abs(rdata1(1,k)) > 90.0 .or. abs(rdata1(2,k)) > 180.0 .or. &
 238                  abs(rdata1(3,k)) > 180.0 .or. rdata2(1,k) > r8bfms .or. &
 239                  rdata2(2,k) > r8bfms ) then
 240                nlev_bad = nlev_bad + 1
 241             end if
 242          end do
 243          if ( nlev_bad == nlev1 ) cycle reports
 244
 245          ntotal = ntotal + 1
 246          ntotal_ifgat(ifgat)=ntotal_ifgat(ifgat)+1
 247          if ( use_gpsephobs .and. gpseph_loadbalance ) then
 248             if ( myproc /= MOD((ntotal-1), num_procs) ) cycle reports
 249          else
 250             if ( outside ) cycle reports
 251          end if
 252
 253          nlocal = nlocal + 1
 254
 255          if (.not. associated(head)) then
 256              nullify ( head )
 257              allocate ( head )
 258              nullify ( head%next )
 259              plink => head
 260          else
 261              allocate ( plink%next )
 262              plink => plink%next
 263              nullify ( plink%next )
 264          end if
 265
 266          write(info%name, '(a,i6.6,a,a)') 'NCEP_GPSRO_', nlocal, '_', date_char
 267
 268          nlev = nlev1
 269          max_lev = max(max_lev, nlev)
 270
 271          plink%info%name      = info%name
 272          plink%info%platform  = platform_name
 273          plink%info%id        = id
 274          plink%info%date_char = date_char
 275          plink%info%levels    = nlev
 276          plink%info%elv       = 0.0            ! not used
 277          plink%info%lat       = info%lat
 278          plink%info%lon       = info%lon
 279          plink%loc%x          = loc%x
 280          plink%loc%y          = loc%y
 281          plink%loc%i          = loc%i
 282          plink%loc%j          = loc%j
 283          plink%loc%dx         = loc%dx
 284          plink%loc%dxm        = loc%dxm
 285          plink%loc%dy         = loc%dy
 286          plink%loc%dym        = loc%dym
 287
 288          plink%obs_global_index = ntotal
 289
 290          plink%ifgat = ifgat
 291          plink%rfict = hdr(7)*0.001 !ELRC: EARTH LOC RADIUS OF CURVATURE
 292
 293          allocate (plink%gpsref%h  (1:nlev))
 294          allocate (plink%gpsref%ref(1:nlev))
 295          if ( use_gpsephobs ) then
 296             allocate (plink%azim(1:nlev))
 297             allocate (plink%lat(1:nlev))
 298             allocate (plink%lon(1:nlev))
 299          end if
 300       end if !treat_as_profile
 301
 302       iprof = iprof + 1
 303
 304       lev_loop: do k = 1, nlev1
 305
 306          info%lat  = rdata1(1,k)
 307          info%lon  = rdata1(2,k)
 308          azim      = rdata1(3,k)
 309          height    = rdata2(1,k)
 310          ref_data  = rdata2(2,k)
 311          ref_qc    = 0                ! initialized to be good
 312
 313          if ( .not. treat_as_profile ) then
 314
 315             ! gpsro.bufr contains missing longitude, occasionally
 316             if ( info%lat > r8bfms .or. info%lon > r8bfms ) then
 317                cycle lev_loop
 318             end if
 319
 320             ! check for missing data
 321             if ( height > r8bfms .or. ref_data > r8bfms ) then
 322                cycle lev_loop
 323             end if
 324
 325             ! check loc
 326             info%lat = max(info%lat, -89.95)
 327             info%lat = min(info%lat,  89.95)
 328             call da_llxy(info, loc, outside, outside_all)
 329             if ( outside_all ) then
 330                num_outside_all = num_outside_all + 1
 331                if ( print_detail_obs ) then
 332                   write(unit=stderr,fmt='(a,2(1x,f8.3),a)')  &
 333                      id(1:5), info%lat, info%lon, '  -> outside_domain'
 334                end if
 335                cycle lev_loop
 336             end if
 337             ntotal = ntotal + 1
 338             ntotal_ifgat(ifgat)=ntotal_ifgat(ifgat)+1
 339             if ( outside ) then
 340                cycle lev_loop
 341             end if
 342             nlocal = nlocal + 1
 343          end if ! not profile
 344
 345          if ( treat_as_profile ) then
 346             plink%azim(k) = azim
 347             plink%lat(k)  = missing_r
 348             plink%lon(k)  = missing_r
 349             plink%gpsref%h(k)         = missing_r
 350             plink%gpsref%ref(k)%inv   = missing_r
 351             plink%gpsref%ref(k)%qc    = missing_data
 352             plink%gpsref%ref(k)%error = xmiss
 353             if ( use_gpsephobs ) then
 354                if ( info%lat < r8bfms .and. info%lon < r8bfms .and. &
 355                     azim < r8bfms ) then
 356                   plink%azim(k) = azim
 357                   if ( gpsro_drift == 0 ) then
 358                      plink%lat(k)  = plink%info%lat
 359                      plink%lon(k)  = plink%info%lon
 360                   else
 361                      plink%lat(k)  = info%lat
 362                      plink%lon(k)  = info%lon
 363                   end if
 364                end if
 365             end if
 366             if ( height < r8bfms .and. ref_data < r8bfms ) then
 367                plink%gpsref%h(k)         = height
 368                plink%gpsref%ref(k)%inv   = ref_data
 369                plink%gpsref%ref(k)%qc    = ref_qc
 370                plink%gpsref%ref(k)%error = ref_error
 371             end if
 372          end if !treat_as_profile
 373
 374          if ( ref_data < r8bfms ) then ! non-missing data
 375             ref_error = ref_data * 0.01
 376          end if
 377
 378          if ( height < r8bfms .and. info%lat < r8bfms ) then ! non-missing height
 379
 380             ! check height, only keep data below certain height (default is 30km)
 381             if ( height > top_km_gpsro*1000.0 .or. &
 382                  height < bot_km_gpsro*1000.0 ) then
 383                ref_qc = -77
 384             end if
 385
 386             err_opt = 1
 387             ! observation errors  WRF-Var/obsproc
 388             if ( err_opt == 1 ) then
 389                if ( height >= 12000.0 ) then
 390                   ref_error = ref_error * 0.3
 391                else
 392                   erh90 = (0.3-1.5)*(height-12000.0)/(12000.0-0.0) + 0.3
 393                   if ( height >= 5500.0 ) then
 394                      erh0 = (0.3-1.3)*(height-12000.0)/(12000.0-5500.0) + 0.3
 395                   else if ( height >= 2500.0) then
 396                      erh0 = (1.3-2.5)*(height-5500.0)/(5500.0-2500.0) + 1.3
 397                   else
 398                      erh0 = 2.5
 399                   end if
 400                   err90 = ref_error * erh90
 401                   err0  = ref_error * erh0
 402                   ref_error = err90 - (1.0-abs(info%lat)/90.0)*(err90-err0)
 403                end if
 404             end if
 405
 406             ! observation errors  GSI_Q1FY09,  Kuo et al. 2003
 407             if ( err_opt == 2 ) then
 408                if ( (info%lat >= -30.0) .and. (info%lat <= 30.0) ) then   ! tropics
 409                   if ( (height >= 7000.0) .and. (height <= 31000.0) ) then
 410                      ref_error = ref_error*(0.1125+(1.25e-5*height))
 411                   else if ( height > 31000.0 ) then
 412                      ref_error = ref_error*0.5
 413                   else if ( height < 7000.0  ) then
 414                      ref_error = ref_error*(3.0-(4.0e-4*height))
 415                   else
 416                      write(unit=message(1),fmt='(a,f8.1,a,f8.2)') 'unable to process with height = ', &
 417                         height, ' at lat = ', info%lat
 418                      call da_error(__FILE__,__LINE__,message(1:1))
 419                   end if
 420                else   ! mid-latitudes
 421                   if ( (height >= 5000.0) .and. (height <= 25000.0) ) then
 422                      ref_error = ref_error*0.3
 423                   else if ( (height >= 25000.0) .and. (height <= 31000.0) ) then
 424                      ref_error = ref_error*(-3.45+(1.5e-4*height))
 425                   else if ( height > 31000.0 ) then
 426                      ref_error = ref_error*1.2
 427                   else if ( height < 5000.0 ) then
 428                      ref_error = ref_error*(0.75-(9.0e-5*height))
 429                   else
 430                      write(unit=message(1),fmt='(a,f8.1,a,f8.2)') 'unable to process with height = ', &
 431                         height, ' at lat = ', info%lat
 432                      call da_error(__FILE__,__LINE__,message(1:1))
 433                   end if
 434                end if
 435             end if
 436
 437          end if ! non-missing height
 438
 439          !write(info%name, '(a,i6.6,a,a)') 'NCEP_GPSRO_', nlocal, '_', date_char
 440          write(info%name, '(a,i6.6,a,a)') 'NCEP_GPSRO_', iprof, '_', date_char
 441
 442          if ( print_detail_obs ) then
 443             write(unit=stdout,fmt='(a,1x,a,1x,i4.4,4i2.2,2f8.2,f8.1,f8.2,i3,f9.5)')  &
 444                info%name,id(1:5),iyear,imonth,iday,ihour,imin, &
 445                info%lat,info%lon,height,ref_data,ref_qc,ref_error
 446          end if
 447
 448          if ( .not. treat_as_profile ) then
 449             if (.not. associated(head)) then
 450                 nullify ( head )
 451                 allocate ( head )
 452                 nullify ( head%next )
 453                 plink => head
 454             else
 455                 allocate ( plink%next )
 456                 plink => plink%next
 457                 nullify ( plink%next )
 458             end if
 459
 460             plink%info%name      = info%name
 461             plink%info%platform  = platform_name
 462             plink%info%id        = id
 463             plink%info%date_char = date_char
 464             plink%info%levels    = 1              ! each level is treated as separate obs
 465             plink%info%elv       = 0.0            ! not used
 466             plink%info%lat       = info%lat
 467             plink%info%lon       = info%lon
 468             plink%loc%x          = loc%x
 469             plink%loc%y          = loc%y
 470             plink%loc%i          = loc%i
 471             plink%loc%j          = loc%j
 472             plink%loc%dx         = loc%dx
 473             plink%loc%dxm        = loc%dxm
 474             plink%loc%dy         = loc%dy
 475             plink%loc%dym        = loc%dym
 476
 477             plink%obs_global_index = ntotal
 478
 479             ! 1 level
 480             allocate (plink%gpsref%h  (1:1))
 481             allocate (plink%gpsref%ref(1:1))
 482             plink%gpsref%h(1)         = height
 483             plink%gpsref%ref(1)%inv   = ref_data
 484             plink%gpsref%ref(1)%qc    = ref_qc
 485             plink%gpsref%ref(1)%error = ref_error
 486             plink%ifgat=ifgat
 487          end if !not treat_as_profile
 488
 489       end do lev_loop
 490    end do reports
 491
 492    call closbf(iunit)
 493    close(iunit)
 494
 495 end do bufrfile
 496
 497    ! assign data in either iv%gpsref or iv%gpseph
 498    if ( use_gpsephobs ) then
 499       obs_index = gpseph
 500       !re-set mex_lev for pseudo_ob on model mean height levels
 501       max_lev = kde-kds+1
 502       allocate (pseudo_ob%ref(kds:kde))
 503       allocate (pseudo_ob%lat(kds:kde))
 504       allocate (pseudo_ob%lon(kds:kde))
 505       allocate (pseudo_ob%azim(kds:kde))
 506    else
 507       obs_index = gpsref
 508    end if
 509
 510    iv%info(obs_index)%max_lev = max_lev
 511    iv%info(obs_index)%ntotal  = ntotal
 512    iv%info(obs_index)%nlocal  = nlocal
 513
 514    !allocate iv%info
 515
 516    if (iv%info(obs_index)%nlocal > 0) then
 517       allocate (iv%info(obs_index)%name(iv%info(obs_index)%nlocal))
 518       allocate (iv%info(obs_index)%platform(iv%info(obs_index)%nlocal))
 519       allocate (iv%info(obs_index)%id(iv%info(obs_index)%nlocal))
 520       allocate (iv%info(obs_index)%date_char(iv%info(obs_index)%nlocal))
 521       allocate (iv%info(obs_index)%levels(iv%info(obs_index)%nlocal))
 522       allocate (iv%info(obs_index)%lat(iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
 523       allocate (iv%info(obs_index)%lon(iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
 524       allocate (iv%info(obs_index)%elv(iv%info(obs_index)%nlocal))
 525       allocate (iv%info(obs_index)%pstar(iv%info(obs_index)%nlocal))
 526       allocate (iv%info(obs_index)%slp(iv%info(obs_index)%nlocal))
 527       allocate (iv%info(obs_index)%pw(iv%info(obs_index)%nlocal))
 528       allocate (iv%info(obs_index)%x  (kms:kme,iv%info(obs_index)%nlocal))
 529       allocate (iv%info(obs_index)%y  (kms:kme,iv%info(obs_index)%nlocal))
 530       allocate (iv%info(obs_index)%i  (kms:kme,iv%info(obs_index)%nlocal))
 531       allocate (iv%info(obs_index)%j  (kms:kme,iv%info(obs_index)%nlocal))
 532       allocate (iv%info(obs_index)%dx (kms:kme,iv%info(obs_index)%nlocal))
 533       allocate (iv%info(obs_index)%dxm(kms:kme,iv%info(obs_index)%nlocal))
 534       allocate (iv%info(obs_index)%dy (kms:kme,iv%info(obs_index)%nlocal))
 535       allocate (iv%info(obs_index)%dym(kms:kme,iv%info(obs_index)%nlocal))
 536       allocate (iv%info(obs_index)%k  (iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
 537       allocate (iv%info(obs_index)%dz (iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
 538       allocate (iv%info(obs_index)%dzm(iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
 539       allocate (iv%info(obs_index)%zk (iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
 540       allocate (iv%info(obs_index)%proc_domain(iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
 541       allocate (iv%info(obs_index)%thinned(iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
 542       allocate (iv%info(obs_index)%obs_global_index(iv%info(obs_index)%nlocal))
 543       iv%info(obs_index)%proc_domain(:,:)  = .false.
 544       iv%info(obs_index)%thinned(:,:)      = .false.
 545       iv%info(obs_index)%zk(:,:)           = missing_r
 546    end if
 547
 548    if ( use_gpsephobs ) then
 549       if (iv%info(gpseph)%nlocal > 0) allocate(iv%gpseph(1:iv%info(gpseph)%nlocal))
 550    else
 551       if (iv%info(gpsref)%nlocal > 0) allocate(iv%gpsref(1:iv%info(gpsref)%nlocal))
 552    end if
 553
 554    ! sort obs and assign to iv structure
 555
 556    iobs=0
 557    do kk = 1, num_fgat_time
 558       plink => head
 559       iv%info(obs_index)%ptotal(kk)=0
 560
 561       reports2: do i = 1, iv%info(obs_index)%nlocal
 562
 563          if ( plink%ifgat /= kk ) then  ! not in current time slot
 564             plink => plink%next
 565             cycle reports2
 566          else
 567
 568             iobs=iobs+1
 569
 570             iv%info(obs_index)%name(iobs)      = plink%info%name
 571             iv%info(obs_index)%platform(iobs)  = plink%info%platform
 572             iv%info(obs_index)%id(iobs)        = plink%info%id
 573             iv%info(obs_index)%date_char(iobs) = plink%info%date_char
 574             iv%info(obs_index)%levels(iobs)    = plink%info%levels
 575             iv%info(obs_index)%elv(iobs)       = plink%info%elv            ! not used
 576
 577             iv%info(obs_index)%lat(:,iobs)     = plink%info%lat
 578             iv%info(obs_index)%lon(:,iobs)     = plink%info%lon
 579             iv%info(obs_index)%x(:,iobs)       = plink%loc%x
 580             iv%info(obs_index)%y(:,iobs)       = plink%loc%y
 581             iv%info(obs_index)%i(:,iobs)       = plink%loc%i
 582             iv%info(obs_index)%j(:,iobs)       = plink%loc%j
 583             iv%info(obs_index)%dx(:,iobs)      = plink%loc%dx
 584             iv%info(obs_index)%dxm(:,iobs)     = plink%loc%dxm
 585             iv%info(obs_index)%dy(:,iobs)      = plink%loc%dy
 586             iv%info(obs_index)%dym(:,iobs)     = plink%loc%dym
 587
 588             iv%info(obs_index)%slp(iobs)%inv   = missing_r
 589             iv%info(obs_index)%slp(iobs)%qc    = missing_data
 590             iv%info(obs_index)%slp(iobs)%error = xmiss
 591             iv%info(obs_index)%pw(iobs)%inv    = missing_r
 592             iv%info(obs_index)%pw(iobs)%qc     = missing_data
 593             iv%info(obs_index)%pw(iobs)%error  = xmiss
 594
 595             iv%info(obs_index)%obs_global_index(iobs) = plink%obs_global_index
 596
 597             if ( use_gpsephobs ) then
 598                nlev = plink%info%levels
 599                !create pseudo_ob on grid mean altitude for gpseph
 600                call da_gpseph_create_ob(nlev, plink%gpsref%h(1:nlev),       &
 601                                               plink%gpsref%ref(1:nlev)%inv, &
 602                                               plink%lat(1:nlev),            &
 603                                               plink%lon(1:nlev),            &
 604                                               plink%azim(1:nlev),           &
 605                                               pseudo_ob, lowest_level)
 606                if ( lowest_level < 0 ) then
 607                   ! no valid levels found for creating ob
 608                   iv%info(gpseph)%levels(iobs) = 1
 609                   iv%gpseph(iobs)%level1 = 0
 610                   iv%gpseph(iobs)%level2 = 0
 611                   iv%gpseph(iobs)%rfict  = plink%rfict
 612                   !hcl tmp hack: still allocate for 1 level
 613                   allocate (iv%gpseph(iobs)%h  (1))
 614                   allocate (iv%gpseph(iobs)%ref(1))
 615                   allocate (iv%gpseph(iobs)%eph(1))
 616                   allocate (iv%gpseph(iobs)%lat(1))
 617                   allocate (iv%gpseph(iobs)%lon(1))
 618                   allocate (iv%gpseph(iobs)%azim(1))
 619                   iv%gpseph(iobs)%h(1)         = global_h_mean(1)*1000.0 !km to m
 620                   iv%gpseph(iobs)%lat(1)       = missing_r
 621                   iv%gpseph(iobs)%lon(1)       = missing_r
 622                   iv%gpseph(iobs)%azim(1)      = missing_r
 623                   iv%gpseph(iobs)%ref(1)%inv   = missing_r
 624                   iv%gpseph(iobs)%ref(1)%qc    = missing_data
 625                   iv%gpseph(iobs)%ref(1)%error = xmiss
 626                   iv%gpseph(iobs)%eph(1)%inv   = missing_r
 627                   iv%gpseph(iobs)%eph(1)%qc    = missing_data
 628                   iv%gpseph(iobs)%eph(1)%error = xmiss
 629                else
 630                   iv%gpseph(iobs)%level1 = lowest_level+1
 631                   iv%gpseph(iobs)%level2 = kde-1
 632                   iv%info(gpseph)%levels(iobs) = max_lev !it is kde-kds+1 for gpseph
 633                   iv%gpseph(iobs)%rfict  = plink%rfict
 634
 635                   allocate (iv%gpseph(iobs)%h  (kds:kde))
 636                   allocate (iv%gpseph(iobs)%ref(kds:kde))
 637                   allocate (iv%gpseph(iobs)%eph(kds:kde))
 638                   ! iv%gpseph includes additional azim/lat/lon info of each level
 639                   allocate (iv%gpseph(iobs)%lat(kds:kde))
 640                   allocate (iv%gpseph(iobs)%lon(kds:kde))
 641                   allocate (iv%gpseph(iobs)%azim(kds:kde))
 642                   do k = kts, kte
 643                      iv%gpseph(iobs)%h(k)         = global_h_mean(k)*1000.0 !km to m
 644                      iv%gpseph(iobs)%lat(k)       = pseudo_ob%lat(k)
 645                      iv%gpseph(iobs)%lon(k)       = pseudo_ob%lon(k)
 646                      iv%gpseph(iobs)%azim(k)      = pseudo_ob%azim(k)
 647                      iv%gpseph(iobs)%ref(k)%inv   = pseudo_ob%ref(k)
 648                      iv%gpseph(iobs)%ref(k)%qc    = missing_data !plink%gpsref%ref(i)%qc
 649                      iv%gpseph(iobs)%ref(k)%error = xmiss        !plink%gpsref%ref(i)%error
 650                      iv%gpseph(iobs)%eph(k)%inv   = missing_r
 651                      iv%gpseph(iobs)%eph(k)%qc    = missing_data
 652                      iv%gpseph(iobs)%eph(k)%error = xmiss
 653                   end do
 654                end if
 655
 656             else !not use_gpsephobs
 657
 658                nlev = plink%info%levels
 659                allocate (iv%gpsref(iobs)%h  (1:nlev))
 660                allocate (iv%gpsref(iobs)%ref(1:nlev))
 661                allocate (iv%gpsref(iobs)%p  (1:nlev))
 662                allocate (iv%gpsref(iobs)%t  (1:nlev))
 663                allocate (iv%gpsref(iobs)%q  (1:nlev))
 664                do k = 1, nlev
 665                   iv%gpsref(iobs)%h(k)         = plink%gpsref%h(k)
 666                   iv%gpsref(iobs)%ref(k)%inv   = plink%gpsref%ref(k)%inv
 667                   iv%gpsref(iobs)%ref(k)%qc    = plink%gpsref%ref(k)%qc
 668                   iv%gpsref(iobs)%ref(k)%error = plink%gpsref%ref(k)%error
 669
 670                   iv%gpsref(iobs)%p(k)%inv     = missing_r
 671                   iv%gpsref(iobs)%p(k)%qc      = missing_data
 672                   iv%gpsref(iobs)%p(k)%error   = xmiss
 673                   iv%gpsref(iobs)%t(k)%inv     = missing_r
 674                   iv%gpsref(iobs)%t(k)%qc      = missing_data
 675                   iv%gpsref(iobs)%t(k)%error   = xmiss
 676                   iv%gpsref(iobs)%q(k)%inv     = missing_r
 677                   iv%gpsref(iobs)%q(k)%qc      = missing_data
 678                   iv%gpsref(iobs)%q(k)%error   = xmiss
 679                end do
 680
 681             end if !use_gpsephobs
 682          end if ! if current time slot
 683          plink => plink%next
 684       end do reports2
 685
 686       ntotal_ifgat(kk)=ntotal_ifgat(kk)+ntotal_ifgat(kk-1)
 687       iv%info(obs_index)%ptotal(kk)=ntotal_ifgat(kk)
 688       iv%info(obs_index)%plocal(kk)=iobs
 689       ! thinning is not applied to GPSRO BUFR, simply make a copy from total number
 690       iv%info(obs_index)%thin_ptotal(kk)=ntotal_ifgat(kk)
 691       iv%info(obs_index)%thin_plocal(kk)=iobs
 692    end do
 693
 694    write(unit=message(1),fmt='(A,3(1x,i7))') &
 695       'da_read_obs_bufrgpsro_eph: num_report, num_outside_all, num_outside_time: ', &
 696       num_report, num_outside_all, num_outside_time
 697    call da_message(message(1:1))
 698
 699    if ( iobs /= iv%info(obs_index)%nlocal ) then
 700       call da_error(__FILE__,__LINE__,(/"numbers mismatch between scanning and reading NCEP GSPRO BUFR file"/))
 701    end if
 702
 703    if ( use_gpsephobs ) then
 704       deallocate (pseudo_ob%ref)
 705       deallocate (pseudo_ob%lat)
 706       deallocate (pseudo_ob%lon)
 707       deallocate (pseudo_ob%azim)
 708    end if
 709
 710    ! Release the linked list
 711    plink => head
 712
 713    DO WHILE ( ASSOCIATED (plink) )
 714       head => plink%next
 715       deallocate (plink%gpsref%h)
 716       deallocate (plink%gpsref%ref)
 717       if ( use_gpsephobs ) then
 718          deallocate (plink%azim)
 719          deallocate (plink%lat)
 720          deallocate (plink%lon)
 721       end if
 722       deallocate ( plink )
 723       plink => head
 724    ENDDO
 725
 726    NULLIFY (head)
 727
 728    if (trace_use) call da_trace_exit("da_read_obs_bufrgpsro_eph")
 729 #else
 730    call da_error(__FILE__,__LINE__,(/"must compile with BUFR library"/))
 731 #endif
 732
 733 end subroutine da_read_obs_bufrgpsro_eph