var/da/da_radiance/da_read_obs_bufrtovs.inc

   1 subroutine da_read_obs_bufrtovs (obstype,iv, infile)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: read in NCEP bufr tovs 1b data to innovation structure
   5    !
   6    !   METHOD: use F90 sequential data structure to avoid reading file twice
   7    !            so that da_scan_bufrtovs is not necessary any more.
   8    !            1. read file radiance 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    !   File History:
  13    !      Peiming Dong Added NPP atms, 2012/4/17
  14    !---------------------------------------------------------------------------
  15
  16    implicit none
  17
  18    character(5)      ,  intent (in)    :: obstype
  19    character(20)     ,  intent (in)    :: infile
  20    type (iv_type)    ,  intent (inout) :: iv
  21
  22 #ifdef BUFR
  23
  24    integer          :: iost
  25    integer(i_kind), allocatable :: nread(:)
  26
  27    integer(i_kind),parameter:: n1bhdr=15
  28    integer(i_kind),parameter:: n2bhdr=2
  29    integer(i_kind),parameter:: maxinfo=12
  30    integer(i_kind),parameter:: maxchanl=100
  31
  32    logical hirs2,hirs3,hirs4,msu,amsua,amsub,mhs,atms
  33    logical outside, outside_all, iuse
  34    integer :: inst
  35
  36    character(10) date
  37    character(8) subset,subfgn
  38    character(80) hdr1b
  39    character(80) hdr2b
  40    integer(i_kind) ihh,i,j,k,ifov,idd,ireadmg,ireadsb
  41    integer(i_kind) iret,idate,im,iy,nchan
  42    integer :: num_bufr(7), numbufr, ibufr
  43    character(20) :: filename
  44
  45    ! thinning variables
  46    integer(i_kind) itt,itx,iobs,iout
  47    real(r_kind) terrain,timedif,crit,dist
  48    real(r_kind) dlon_earth,dlat_earth
  49
  50    real(r_kind) tbmin,tbmax, tbbad
  51    real(r_kind) panglr,rato
  52    ! real(r_kind) rmask
  53    real(r_kind) step,start
  54
  55    real(r_double),dimension(maxinfo+maxchanl):: data1b8
  56    real(r_double),dimension(n1bhdr):: bfr1bhdr
  57    real(r_double),dimension(n2bhdr):: bfr2bhdr
  58
  59    ! Instrument triplet, follow the convension of RTTOV
  60    integer   :: platform_id, satellite_id, sensor_id
  61
  62    ! pixel information
  63    integer   ::  year,month,day,hour,minute,second  ! observation time
  64    real*8    ::  obs_time
  65    ! real      ::  rlat, rlon                         !  lat/lon in degrees   for Anfovs
  66    real      ::  satzen, satazi, solzen ,solazi       !  scan angles for Anfovs
  67    integer   ::  landsea_mask
  68    real      ::  srf_height
  69    ! channels' bright temperature
  70    real , allocatable ::   tb_inv(:)                    !  bright temperatures
  71    !  end type bright_temperature
  72
  73    type (datalink_type), pointer    :: head, p, current, prev
  74
  75    integer                        ::  ifgat
  76    type(info_type)                ::  info
  77    type(model_loc_type)           ::  loc
  78
  79    data hdr1b /'SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SAZA SOZA HOLS LSQL SOLAZI'/
  80    data hdr2b /'CLATH CLONH'/
  81    !  data hdr1b /'FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SAZA SOZA HOLS LSQL SLNM BEARAZ'/
  82
  83    data tbmin,tbmax,tbbad / 50.0_r_kind, 550.0_r_kind, -9.99e11_r_kind /
  84    integer :: num_tovs_local, num_tovs_file, num_tovs_global, num_tovs_selected
  85    integer :: num_tovs_thinned, num_tovs_used, num_tovs_used_tmp
  86    integer :: lnbufr
  87    integer :: n
  88    integer(i_kind), allocatable :: ptotal(:)
  89    real , allocatable :: in(:), out(:)
  90    logical :: found, head_found
  91
  92    call da_trace_entry("da_read_obs_bufrtovs")
  93
  94    ! Initialize variables
  95
  96    nchan = 22
  97    allocate(nread(1:rtminit_nsensor))
  98    allocate(ptotal(0:num_fgat_time))
  99    nread(1:rtminit_nsensor) = 0
 100    ptotal(0:num_fgat_time) = 0
 101
 102    ! Set various variables depending on type of data to be read
 103
 104    ! platform_id  = 1                 !! for NOAA series
 105    ! platform_id  = 10                !! for METOP series
 106
 107    hirs2=     obstype == 'hirs2'
 108    hirs3=     obstype == 'hirs3'
 109    hirs4=     obstype == 'hirs4'
 110    msu=       obstype == 'msu  '
 111    amsua=     obstype == 'amsua'
 112    amsub=     obstype == 'amsub'
 113    mhs=       obstype == 'mhs  '
 114    atms=      obstype == 'atms '
 115
 116    if (hirs2) then
 117       sensor_id    =  0
 118       step   = 1.80_r_kind
 119       start  = -49.5_r_kind
 120       nchan=nchan_hirs2
 121       subfgn='NC021021'
 122       rato=1.1363987_r_kind
 123    else if (hirs3) then
 124       sensor_id    =  0
 125       step   = 1.80_r_kind
 126       start  = -49.5_r_kind
 127       nchan=nchan_hirs3
 128       subfgn='NC021025'
 129    else if (hirs4) then
 130       sensor_id    =  0
 131       step   = 1.80_r_kind
 132       start  = -49.5_r_kind
 133       nchan=nchan_hirs4
 134       subfgn='NC021028'
 135    else if (mhs) then
 136       sensor_id    =  15
 137       step   = 10.0_r_kind/9.0_r_kind
 138       start  = -445.0_r_kind/9.0_r_kind
 139       nchan=nchan_mhs
 140       subfgn='NC021027'
 141    else if (msu) then
 142       sensor_id    =  1
 143       step   = 9.474_r_kind
 144       start  = -47.37_r_kind
 145       nchan=nchan_msu
 146       subfgn='NC021022'
 147       rato=1.1363987_r_kind
 148    else if (amsua) then
 149       sensor_id    =  3
 150       step   = three + one/three
 151       start  = -48.33_r_kind
 152       nchan=nchan_amsua
 153       subfgn='NC021023'
 154    else if (amsub)  then
 155       sensor_id    =  4
 156       step   = 1.1_r_kind
 157       start  = -48.95_r_kind
 158       nchan=nchan_amsub
 159       subfgn='NC021024'
 160    else if (atms) then
 161       sensor_id    =  19
 162       step   = 1.11_r_kind
 163       start  = -52.725_r_kind
 164       nchan=22
 165       subfgn='NC021203'
 166       hdr1b='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SAZA SOZA HMSL LSQL SOLAZI'
 167    end if
 168
 169    allocate (tb_inv(nchan))
 170
 171    num_tovs_file     = 0    ! number of obs in file
 172    num_tovs_global   = 0    ! number of obs within whole domain
 173    num_tovs_local    = 0    ! number of obs within tile
 174    num_tovs_thinned  = 0    ! number of obs rejected by thinning
 175    num_tovs_used     = 0    ! number of obs entered into innovation computation
 176    num_tovs_selected = 0    ! number of obs limited for debuging
 177    iobs = 0                 ! for thinning, argument is inout
 178
 179    ! 0.0  Open unit to satellite bufr file and read file header
 180    !--------------------------------------------------------------
 181
 182    num_bufr(:)=0
 183    numbufr=0
 184    if (num_fgat_time>1) then
 185       do i=1,7
 186          call da_get_unit(lnbufr)
 187          write(filename,fmt='(A,2I1,A)') trim(infile),0,i,'.bufr'
 188          open(unit   = lnbufr, FILE   = trim(filename),iostat =  iost, form = 'unformatted', STATUS = 'OLD')
 189          if (iost == 0) then
 190             numbufr=numbufr+1
 191             num_bufr(numbufr)=i
 192          else
 193             close (lnbufr)
 194          end if
 195          call da_free_unit(lnbufr)
 196       end do
 197    else
 198      numbufr=1
 199    end if
 200
 201    if (numbufr==0) numbufr=1
 202
 203 bufrfile:  do ibufr=1,numbufr
 204    if (num_fgat_time==1) then
 205       filename=trim(infile)//'.bufr'
 206    else
 207       if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then
 208          filename=trim(infile)//'.bufr'
 209       else
 210          write(filename,fmt='(A,2I1,A)') trim(infile),0,num_bufr(ibufr),'.bufr'
 211       end if
 212    end if
 213
 214 !  We want to use specific unit number for bufr data, so we can control the endian format in environment.
 215    lnbufr = 99
 216
 217    open(unit=lnbufr,file=trim(filename),form='unformatted', &
 218       iostat = iost, status = 'old')
 219    if (iost /= 0) then
 220       call da_warning(__FILE__,__LINE__, &
 221          (/"Cannot open file "//infile/))
 222       call da_trace_exit("da_read_obs_bufrtovs")
 223       return
 224    end if
 225
 226    call openbf(lnbufr,'IN',lnbufr)
 227    call datelen(10)
 228    call readmg(lnbufr,subset,idate,iret)
 229    if (subset /= subfgn) then
 230       call closbf(lnbufr)
 231       close(lnbufr)
 232       message(1)='The file title does not match the data subset'
 233       write(unit=message(2),fmt=*) &
 234          'infile=', lnbufr, infile,' subset=', subset, ' subfgn=',subfgn
 235       call da_error(__FILE__,__LINE__,message(1:2))
 236    end if
 237
 238    iy=0
 239    im=0
 240    idd=0
 241    ihh=0
 242    write(unit=date,fmt='( i10)') idate
 243    read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
 244    write(unit=stdout,fmt=*) &
 245       'Bufr file date is ',iy,im,idd,ihh,infile
 246
 247    ! Loop to read bufr file and assign information to a sequential structure
 248    !-------------------------------------------------------------------------
 249
 250    if ( ibufr == 1 ) then
 251       allocate (head)
 252       !  allocate ( head % tb_inv (1:nchan) )
 253       nullify  ( head % next )
 254       p => head
 255    endif
 256
 257    if (tovs_start > 1) then
 258       write (unit=stdout,fmt='(A,I6)') "   Skipping tovs obs before", tovs_start
 259    end if
 260
 261
 262    obs: do while (ireadmg(lnbufr,subset,idate)==0 .and. subset==subfgn)
 263       do while (ireadsb(lnbufr)==0)
 264
 265          ! 1.0     Read header record and data record
 266          call ufbint(lnbufr,bfr1bhdr,n1bhdr,1,iret,hdr1b)
 267          call ufbint(lnbufr,bfr2bhdr,n2bhdr,1,iret,hdr2b)
 268          call ufbrep(lnbufr,data1b8,1,nchan,iret,'TMBR')
 269          ! call ufbrep(lnbufr,data1b8,1,1,iret,'BEARAZ')
 270
 271          ! check if observation outside range
 272
 273          num_tovs_file = num_tovs_file + 1
 274
 275          ! 2.0     Extract observation location and other required information
 276          !     QC1:  judge if data is in the domain, read next record if not
 277          !------------------------------------------------------------------------
 278          ! rlat = bfr1bhdr(bufr_lat)
 279          ! rlon = bfr1bhdr(bufr_lat)
 280          ! if (rlon < 0.0) rlon = rlon+360.0
 281
 282          if(abs(bfr2bhdr(1)) <= 90. .and. abs(bfr2bhdr(2)) <= 360.)then
 283             info%lat = bfr2bhdr(1)
 284             info%lon = bfr2bhdr(2)
 285          elseif(abs(bfr1bhdr(9)) <= 90. .and. abs(bfr1bhdr(10)) <= 360.)then
 286             info%lat = bfr1bhdr(9)
 287             info%lon = bfr1bhdr(10)
 288          endif
 289
 290          call da_llxy (info, loc, outside, outside_all)
 291
 292          if (outside_all) cycle
 293
 294          !  3.0     Extract other information
 295          !------------------------------------------------------
 296          !  3.1     Extract satellite id and scan position.
 297
 298          if ( nint(bfr1bhdr(bufr_satellite_id)) >= 200 .and. nint(bfr1bhdr(bufr_satellite_id)) <= 204 ) then
 299             platform_id = 1
 300             satellite_id = nint(bfr1bhdr(bufr_satellite_id))-192
 301          else if ( nint(bfr1bhdr(bufr_satellite_id)) >= 205 .and. nint(bfr1bhdr(bufr_satellite_id)) <= 209 ) then
 302             platform_id = 1
 303             satellite_id = nint(bfr1bhdr(bufr_satellite_id))-191
 304          else if ( nint(bfr1bhdr(bufr_satellite_id)) == 223 ) then ! noaa-19
 305             platform_id = 1
 306             satellite_id = nint(bfr1bhdr(bufr_satellite_id))-204
 307          else if ( nint(bfr1bhdr(bufr_satellite_id)) >= 3 .and. nint(bfr1bhdr(bufr_satellite_id)) <= 5 ) then
 308             platform_id = 10
 309             satellite_id = nint(bfr1bhdr(bufr_satellite_id))-2
 310          else if ( nint(bfr1bhdr(bufr_satellite_id)) == 708 ) then ! tiros-n
 311             platform_id = 19
 312             satellite_id = 0
 313          else if ( nint(bfr1bhdr(bufr_satellite_id)) == 224 ) then ! npp
 314             platform_id = 17
 315             satellite_id = 0
 316          else if ( nint(bfr1bhdr(bufr_satellite_id)) >= 701 .and. nint(bfr1bhdr(bufr_satellite_id)) <= 707 ) then
 317             platform_id = 1
 318             satellite_id = nint(bfr1bhdr(bufr_satellite_id))-700
 319          end if
 320          ifov = nint(bfr1bhdr(bufr_ifov))
 321
 322
 323          !  QC2:  limb pixel rejected (not implemented)
 324
 325          !  3.2     Extract date information.
 326
 327          year   = bfr1bhdr(bufr_year)
 328          month  = bfr1bhdr(bufr_month)
 329          day    = bfr1bhdr(bufr_day)
 330          hour   = bfr1bhdr(bufr_hour)
 331          minute = bfr1bhdr(bufr_minute)
 332          second = bfr1bhdr(bufr_second)
 333
 334          write(unit=info%date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 335             year, '-', month, '-', day, '_', hour, ':', minute, ':', second
 336
 337          !  QC3: time consistency check with the background date
 338
 339          if (year <= 99) then
 340             if (year < 78) then
 341                year = year + 2000
 342             else
 343                year = year + 1900
 344             end if
 345          end if
 346
 347          call da_get_julian_time(year,month,day,hour,minute,obs_time)
 348
 349          if (obs_time < time_slots(0) .or.  &
 350             obs_time >= time_slots(num_fgat_time)) cycle
 351
 352          ! 3.2.1 determine FGAT index ifgat
 353
 354          do ifgat=1,num_fgat_time
 355             if (obs_time >= time_slots(ifgat-1) .and.  &
 356                 obs_time  < time_slots(ifgat)) exit
 357          end do
 358
 359          ! 3.3 Find wrfvar instrument index from RTTOV instrument triplet
 360          !     go to next data if id is not in the lists
 361
 362          inst = 0
 363          do i = 1, rtminit_nsensor
 364             if (platform_id  == rtminit_platform(i) &
 365                .and. satellite_id == rtminit_satid(i)    &
 366                .and. sensor_id    == rtminit_sensor(i)) then
 367                inst = i
 368                exit
 369             end if
 370          end do
 371          if (inst == 0) cycle
 372
 373          ! 3.4 extract satellite and solar angle
 374
 375          panglr=(start+float(ifov-1)*step)*deg2rad
 376          if (hirs2 .or. msu) then
 377             satzen = asin(rato*sin(panglr))*rad2deg
 378             satzen = abs(satzen)
 379          else
 380             satzen = bfr1bhdr(bufr_satzen) !*deg2rad   ! local zenith angle
 381             satzen = abs(satzen)
 382
 383             ! if (amsua .and. ifov .le. 15) satzen=-satzen
 384             ! if (amsub .and. ifov .le. 45) satzen=-satzen
 385             ! if (hirs3 .and. ifov .le. 28) satzen=-satzen
 386          end if
 387          if ( satzen > 65.0 ) cycle   ! CRTM has a satzen > 65.0 check
 388          satazi = panglr*rad2deg            ! look angle
 389          ! if (satazi<0.0) satazi = satazi+360.0
 390          solzen = bfr1bhdr(bufr_solzen)              ! solar zenith angle
 391          solazi = bfr1bhdr(bufr_solazi)              !RTTOV9_3
 392
 393          num_tovs_global = num_tovs_global + 1
 394          ptotal(ifgat) = ptotal(ifgat) + 1
 395
 396          if (num_tovs_global < tovs_start) then
 397             cycle
 398          end if
 399
 400          if (num_tovs_global > tovs_end) then
 401             write (unit=stdout,fmt='(A,I6)') "   Skipping radiance obs after", tovs_end
 402             exit obs
 403          end if
 404
 405          num_tovs_selected = num_tovs_selected + 1
 406
 407          if (num_tovs_selected > max_tovs_input) then
 408             write(unit=message(1),fmt='(A,I10,A)') &
 409                "Max number of tovs",max_tovs_input," reached"
 410             call da_warning(__FILE__,__LINE__,message(1:1))
 411             num_tovs_selected = num_tovs_selected-1
 412             num_tovs_global   = num_tovs_global-1
 413             ptotal(ifgat) = ptotal(ifgat) - 1
 414             exit obs
 415          end if
 416
 417          if (outside) cycle ! No good for this PE
 418          num_tovs_local = num_tovs_local + 1
 419
 420          !  Make Thinning
 421          !  Map obs to thinning grid
 422          !-------------------------------------------------------------------
 423          if (thinning) then
 424             dlat_earth = info%lat
 425             dlon_earth = info%lon
 426             if (dlon_earth<zero) dlon_earth = dlon_earth+r360
 427             if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 428             dlat_earth = dlat_earth*deg2rad
 429             dlon_earth = dlon_earth*deg2rad
 430             timedif = 0.0 !2.0_r_kind*abs(tdiff)        ! range:  0 to 6
 431             terrain = 0.01_r_kind*abs(bfr1bhdr(13))
 432             crit = 1.0 !0.01_r_kind+terrain + timedif !+ 10.0_r_kind*float(iskip)
 433             call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 434             if (.not. iuse) then
 435                num_tovs_thinned=num_tovs_thinned+1
 436                cycle
 437             end if
 438          end if
 439
 440          num_tovs_used = num_tovs_used + 1
 441          nread(inst) = nread(inst) + 1
 442
 443          ! 3.5 extract surface information
 444          srf_height = bfr1bhdr(bufr_station_height)          ! station height
 445          if (srf_height < 8888.0 .AND. srf_height > -416.0) then
 446          else
 447             srf_height = 0.0
 448          endif
 449
 450          landsea_mask = nint(bfr1bhdr(bufr_landsea_mask))  ! 0:land ; 1:sea (same as RTTOV)
 451          ! rmask=one                          ! land
 452          ! if (nint(bfr1bhdr(bufr_landsea_mask))==1) rmask=0.0_r_kind   ! reverse the land/sea mask in bufr
 453          ! landsea_mask = rmask+.001_r_kind             ! land sea mask
 454
 455          info%elv = srf_height
 456
 457          ! 3.6 extract channel bright temperature
 458
 459          tb_inv(1:nchan) = data1b8(1:nchan)
 460          do k = 1, nchan
 461             if ( tb_inv(k) < tbmin .or. tb_inv(k) > tbmax) &
 462                tb_inv(k) = missing_r
 463          end do
 464          if ( all(tb_inv<0.0) ) then
 465             num_tovs_local = num_tovs_local -1
 466             num_tovs_used = num_tovs_used - 1
 467             nread(inst) = nread(inst) - 1
 468             cycle
 469          end if
 470
 471          !  4.0   assign information to sequential radiance structure
 472          !--------------------------------------------------------------------------
 473          allocate (p % tb_inv (1:nchan))
 474          p%info             = info
 475          p%loc              = loc
 476          p%landsea_mask     = landsea_mask
 477          p%scanpos          = ifov
 478          p%satzen           = satzen
 479          p%satazi           = satazi
 480          p%solzen           = solzen
 481          p%tb_inv(1:nchan)  = tb_inv(1:nchan)
 482          p%sensor_index     = inst
 483          p%ifgat            = ifgat
 484 !RTTOV9_3
 485          p%solazi           = solazi
 486  !end of RTTOV9_3
 487          allocate (p%next)   ! add next data
 488
 489          p => p%next
 490          nullify (p%next)
 491       end do
 492    end do obs
 493
 494    call closbf(lnbufr)
 495    close(lnbufr)
 496
 497 end do bufrfile
 498
 499    if (thinning .and. num_tovs_global > 0 ) then
 500
 501 #ifdef DM_PARALLEL
 502
 503       ! Get minimum crit and associated processor index.
 504       j = 0
 505       do ifgat = 1, num_fgat_time
 506          do n = 1, iv%num_inst
 507             j = j + thinning_grid(n,ifgat)%itxmax
 508          end do
 509       end do
 510
 511       allocate ( in  (j) )
 512       allocate ( out (j) )
 513
 514       j = 0
 515       do ifgat = 1, num_fgat_time
 516          do n = 1, iv%num_inst
 517             do i = 1, thinning_grid(n,ifgat)%itxmax
 518                j = j + 1
 519                in(j) = thinning_grid(n,ifgat)%score_crit(i)
 520             end do
 521          end do
 522       end do
 523       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 524
 525       call wrf_dm_bcast_real (out, j)
 526
 527       j = 0
 528       do ifgat = 1, num_fgat_time
 529          do n = 1, iv%num_inst
 530             do i = 1, thinning_grid(n,ifgat)%itxmax
 531                j = j + 1
 532                if ( ABS(out(j)-thinning_grid(n,ifgat)%score_crit(i)) > 1.0E-10 ) thinning_grid(n,ifgat)%ibest_obs(i)  = 0
 533             end do
 534          end do
 535       end do
 536
 537       deallocate( in  )
 538       deallocate( out )
 539
 540 #endif
 541
 542       ! Delete the nodes which being thinning out
 543       p => head
 544       prev => head
 545       head_found = .false.
 546       num_tovs_used_tmp = num_tovs_used
 547       do j = 1, num_tovs_used_tmp
 548          n = p%sensor_index
 549          ifgat = p%ifgat
 550          found = .false.
 551
 552          do i = 1, thinning_grid(n,ifgat)%itxmax
 553             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 554                found = .true.
 555                exit
 556             endif
 557          end do
 558
 559          ! free current data
 560          if ( .not. found ) then
 561             current => p
 562             p => p%next
 563             if ( head_found ) then
 564                prev%next => p
 565             else
 566                head => p
 567                prev => p
 568             endif
 569             deallocate ( current % tb_inv )
 570             deallocate ( current )
 571             num_tovs_thinned = num_tovs_thinned + 1
 572             num_tovs_used = num_tovs_used - 1
 573             nread(n) = nread(n) - 1
 574             continue
 575          endif
 576
 577          if ( found .and. head_found ) then
 578             prev => p
 579             p => p%next
 580             continue
 581          endif
 582
 583          if ( found .and. .not. head_found ) then
 584             head_found = .true.
 585             head => p
 586             prev => p
 587             p => p%next
 588          endif
 589
 590       end do
 591
 592    endif  ! End of thinning
 593
 594    iv%total_rad_pixel   = iv%total_rad_pixel   + num_tovs_used
 595    iv%total_rad_channel = iv%total_rad_channel + num_tovs_used*nchan
 596
 597    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_tovs_used
 598    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_tovs_global
 599
 600    do i = 1, num_fgat_time
 601       ptotal(i) = ptotal(i) + ptotal(i-1)
 602       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 603    end do
 604    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 605       write(unit=message(1),fmt='(A,I10,A,I10)') &
 606           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 607       call da_warning(__FILE__,__LINE__,message(1:1))
 608    endif
 609
 610    write(unit=stdout,fmt='(a)') 'num_tovs_file num_tovs_global num_tovs_local num_tovs_used num_tovs_thinned'
 611    write(unit=stdout,fmt='(5i10)') num_tovs_file,num_tovs_global, num_tovs_local,num_tovs_used,num_tovs_thinned
 612
 613    deallocate(tb_inv)
 614
 615    !  5.0 allocate innovation radiance structure
 616    !----------------------------------------------------------------
 617
 618    do i = 1, iv%num_inst
 619       if (nread(i) < 1) cycle
 620       iv%instid(i)%num_rad = nread(i)
 621       iv%instid(i)%info%nlocal = nread(i)
 622       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 623         'Allocating space for radiance innov structure', &
 624          i, iv%instid(i)%rttovid_string, iv%instid(i)%num_rad
 625
 626       call da_allocate_rad_iv(i,nchan,iv)
 627
 628    end do
 629
 630    !  6.0 assign sequential structure to innovation structure
 631    !-------------------------------------------------------------
 632    nread(1:rtminit_nsensor) = 0
 633    p => head
 634    ! do while ( associated(p) )
 635
 636    do n = 1, num_tovs_used
 637       i = p%sensor_index
 638       nread(i) = nread(i) + 1
 639
 640       call da_initialize_rad_iv (i, nread(i), iv, p)
 641
 642       current => p
 643       p => p%next
 644
 645       ! free current data
 646       deallocate ( current % tb_inv )
 647       deallocate ( current )
 648    end do
 649
 650    deallocate ( p )
 651
 652    deallocate (nread)
 653    deallocate (ptotal)
 654
 655    ! check if sequential structure has been freed
 656    !
 657    ! p => head
 658    ! do i = 1, num_rad_selected
 659    !    write (unit=stdout,fmt=*)  i, p%tb_inv(1:nchan)
 660    !    p => p%next
 661    ! end do
 662
 663    call da_trace_exit("da_read_obs_bufrtovs")
 664 #else
 665    call da_error(__FILE__,__LINE__,(/"Needs to be compiled with a BUFR library"/))
 666 #endif
 667
 668
 669 end subroutine da_read_obs_bufrtovs
 670
 671