var/da/da_radiance/da_read_obs_bufriasi.inc

   1 subroutine da_read_obs_bufriasi (obstype,iv,infile)
   2 !$$$  subprogram documentation block
   3 !                .      .    .                                       .
   4 ! subprogram:    read_iasi                  read bufr format iasi data
   5
   6   implicit none
   7
   8
   9   character(9)      ,  intent (in)  :: obstype
  10   character(100)    ,  intent (in)  :: infile
  11   type (iv_type)    ,intent (inout) :: iv
  12 #ifdef BUFR
  13   real(kind=8)    ::  obs_time
  14   type (datalink_type), pointer  :: head, p, current, prev
  15   type(info_type)                :: info
  16   type(model_loc_type)           :: loc
  17   type(model_loc_type)           :: loc_fov
  18
  19 ! Subroutine constants
  20   real(r_kind)                   :: ten      =  10.0_r_kind
  21   real(r_kind)                   :: r90      =  90.0_r_kind
  22   real(r_kind)                   :: r360     = 360.0_r_kind
  23
  24 ! Number of channels for sensors in BUFR
  25   integer(i_kind),parameter      :: nchan = 616        !--- 616 subset ch out of 8078 ch for AIRS
  26   integer(i_kind),parameter      :: n_totchan  = 616
  27   integer(i_kind),parameter      :: maxinfo    =  33
  28   integer(i_kind)                :: inst,platform_id,satellite_id,sensor_id
  29   integer(i_kind), allocatable   :: ptotal(:), nread(:)
  30   real(r_kind)                   :: crit
  31   integer(i_kind)                :: ifgat, iout, iobs
  32   logical                        :: outside, outside_all, iuse
  33
  34 ! BUFR functions
  35   integer(i_kind)   :: ireadsb,ireadmg
  36
  37 ! Variables for BUFR IO
  38   real(r_double),dimension(5)           :: linele
  39   real(r_double),dimension(14)          :: allspot
  40   real(r_double),dimension(2,n_totchan) :: allchan
  41   real(r_double),dimension(3,10)        :: cscale
  42   real(r_double),dimension(6)           :: cloud_frac
  43   integer           :: numbufr,ibufr
  44   logical           :: found, head_found
  45   real(r_kind)      :: step, start,step_adjust
  46   character(len=8)  :: subset
  47   character(len=4)  :: senname
  48   character(len=80) :: allspotlist
  49   integer(i_kind)   :: jstart
  50   integer(i_kind)   :: iret
  51   character(10)     :: date
  52
  53 ! Work variables for time
  54   integer(i_kind)   :: idate, im, iy, idd, ihh
  55   integer(i_kind)   :: idate5(6)
  56
  57 ! Other work variables
  58   real(r_kind)                          :: piece
  59   real(r_kind)                          :: dlon_earth,dlat_earth
  60   real(r_kind)                          :: lza, lzaest,sat_height_ratio
  61   real(r_kind)                          :: sat_zenang
  62   real(r_kind)                          :: radi
  63   real(r_kind),dimension(10)            :: sscale
  64   real(r_kind),dimension(n_totchan)     :: temperature
  65   real(r_kind),allocatable,dimension(:) :: data_all
  66   real(r_kind)                          :: Effective_Temperature
  67
  68   logical          :: iasi
  69   integer(i_kind)  :: nreal, ksatid
  70   integer(i_kind)  :: ifov, iscn
  71   integer(i_kind)  :: num_iasi_file
  72   integer(i_kind)  :: num_iasi_local, num_iasi_global, num_iasi_used, num_iasi_thinned
  73   integer(i_kind)  :: num_iasi_used_tmp
  74   integer(i_kind)  :: i, j, l, iskip, ifovn, bad_line
  75   integer(i_kind)  :: itx, k, nele, itt, n
  76   integer(i_kind)  :: iexponent
  77   integer(i_kind)  :: num_bufr(7)
  78   integer          :: iost, lnbufr
  79   character(20)    :: filename
  80   real,allocatable :: in(:), out(:)
  81
  82 ! Set standard parameters
  83   integer(i_kind),parameter :: ichan=-999  ! fov-based surface code is not channel specific for iasi
  84   real(r_kind),parameter    :: expansion=one         ! exansion factor for fov-based location code.                                                 ! use one for ir sensors.
  85   real(r_kind),parameter    :: tbmin  = 50._r_kind
  86   real(r_kind),parameter    :: tbmax  = 550._r_kind
  87   real(r_kind),parameter    :: earth_radius = 6371000._r_kind
  88
  89
  90    if (trace_use) call da_trace_entry("da_read_obs_bufriasi")
  91 !  0.0  Initialize variables
  92 !-----------------------------------
  93   platform_id  = 10   ! metop series
  94   sensor_id = 16 !iasi
  95   nreal  = maxinfo
  96   iasi=      obstype == 'iasi'
  97   bad_line=-1
  98   step  = 3.334
  99   start = -48.33
 100   step_adjust = 0.625_r_kind
 101   senname = 'IASI'
 102   allspotlist= &
 103    'SIID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAZA BEARAZ SOZA SOLAZI SAID'
 104   num_bufr(:)=0
 105   numbufr=0
 106   allocate(nread(1:rtminit_nsensor))
 107   allocate(ptotal(0:num_fgat_time))
 108   nread(1:rtminit_nsensor) = 0
 109   ptotal(0:num_fgat_time) = 0
 110   iobs = 0                 ! for thinning, argument is inout
 111   num_iasi_file    = 0
 112   num_iasi_local   = 0
 113   num_iasi_global  = 0
 114   num_iasi_used    = 0
 115   num_iasi_thinned = 0
 116    if (num_fgat_time>1) then
 117       do i=1,7
 118          call da_get_unit(lnbufr)
 119          write(filename,fmt='(A,2I1,A)') trim(infile),0,i,'.bufr'
 120          open(unit   = lnbufr, FILE   = trim(filename),iostat =  iost, form = 'unformatted', STATUS = 'OLD')
 121          if (iost == 0) then
 122             numbufr=numbufr+1
 123             num_bufr(numbufr)=i
 124          else
 125             close (lnbufr)
 126          end if
 127          call da_free_unit(lnbufr)
 128       end do
 129    else
 130      numbufr=1
 131    end if
 132
 133    if (numbufr==0) numbufr=1
 134
 135   bufrfile:  do ibufr=1,numbufr
 136    if (num_fgat_time==1) then
 137       filename=trim(infile)//'.bufr'
 138    else
 139       if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then
 140          filename=trim(infile)//'.bufr'
 141       else
 142          write(filename,fmt='(A,2I1,A)') trim(infile),0,num_bufr(ibufr),'.bufr'
 143       end if
 144    end if
 145    lnbufr = 95
 146    open(unit=lnbufr,file=trim(filename),form='unformatted', &
 147       iostat = iost, status = 'old')
 148    if (iost /= 0) then
 149       call da_warning(__FILE__,__LINE__, &
 150          (/"Cannot open file "//infile/))
 151       if (trace_use) call da_trace_exit("da_read_obs_bufriasi")
 152       return
 153    end if
 154
 155 ! Open BUFR table
 156    call openbf(lnbufr,'IN',lnbufr)
 157    call datelen(10)
 158    call readmg(lnbufr,subset,idate,iret)
 159    iy=0
 160    im=0
 161    idd=0
 162    ihh=0
 163    write(unit=date,fmt='( i10)') idate
 164    read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
 165    write(unit=stdout,fmt='(a,4i4,2x,a)') &
 166       'Bufr file date is ',iy,im,idd,ihh,trim(infile)
 167
 168    ! Loop to read bufr file and assign information to a sequential structure
 169    !-------------------------------------------------------------------------
 170 ! Allocate arrays to hold data
 171   nele=nreal+nchan
 172   allocate(data_all(nele))
 173    if ( ibufr == 1 ) then
 174       allocate (head)
 175       nullify  ( head % next )
 176       p => head
 177    end if
 178
 179
 180 ! Big loop to read data file
 181
 182   do while(ireadmg(lnbufr,subset,idate)>=0)
 183
 184      read_loop: do while (ireadsb(lnbufr)==0)
 185          num_iasi_file = num_iasi_file + 1
 186
 187 !    Read IASI FOV information
 188          call ufbint(lnbufr,linele,5,1,iret,'FOVN SLNM QGFQ MJFC SELV')
 189          if ( linele(3) /= zero) then
 190                  cycle read_loop
 191              end if
 192          if ( bad_line == nint(linele(2))) then
 193 !        zenith angle/scan spot mismatch, reject entire line
 194            cycle read_loop
 195          else
 196            bad_line = -1
 197          end if
 198           call ufbint(lnbufr,allspot,14,1,iret,allspotlist)
 199          if(iret /= 1) cycle read_loop
 200
 201          ksatid = nint(allspot(14))
 202          ! SAID 3 is metop-b (metop-1)
 203          ! SAID 4 is metop-a (metop-2)
 204          ! SAID 5 is metop-c (metop-3)
 205          if ( ( ksatid > 5) .or. ( ksatid < 3) ) then
 206             write(unit=message(1),fmt='(A,I6)') 'Unknown platform: ', ksatid
 207             call da_warning(__FILE__,__LINE__,message(1:1))
 208          end if
 209          if ( ksatid == 3 ) then
 210              satellite_id = 1
 211          else if ( ksatid == 4 ) then
 212              satellite_id = 2
 213          else if ( ksatid == 5 ) then
 214              satellite_id = 3
 215          end if
 216
 217 !    Check observing position
 218          info%lat  =  allspot(8)  ! latitude
 219          info%lon  =  allspot(9)  ! longitude)
 220          if( abs(info%lat) > r90  .or. abs(info%lon) > r360 .or. &
 221          (abs(info%lat) == r90 .and. info%lon /= zero) )then
 222          write(unit=stdout,fmt=*) &
 223          'READ_IASI:  ### ERROR IN READING ', senname, ' BUFR DATA:', &
 224                ' STRANGE OBS POINT (LAT,LON):', info%lat, info%lon
 225             cycle read_loop
 226          end if
 227
 228          call da_llxy (info, loc, outside, outside_all)
 229          if (outside_all) cycle
 230              inst = 0
 231          do i = 1, rtminit_nsensor
 232             if (platform_id  == rtminit_platform(i) &
 233                .and. satellite_id == rtminit_satid(i)    &
 234                .and. sensor_id    == rtminit_sensor(i)) then
 235                inst = i
 236                exit
 237             end if
 238          end do
 239          if (inst == 0) cycle read_loop
 240
 241 !    Check obs time
 242          idate5(1) = nint(allspot(2)) ! year
 243          idate5(2) = nint(allspot(3)) ! month
 244          idate5(3) = nint(allspot(4)) ! day
 245          idate5(4) = nint(allspot(5)) ! hour
 246          idate5(5) = nint(allspot(6)) ! minute
 247          idate5(6) = nint(allspot(7)) ! second
 248
 249          if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. &
 250              idate5(2) < 1    .or. idate5(2) >   12 .or. &
 251              idate5(3) < 1    .or. idate5(3) >   31 .or. &
 252              idate5(4) <0     .or. idate5(4) >   24 .or. &
 253              idate5(5) <0     .or. idate5(5) >   60 )then
 254
 255             write(6,*)'READ_IASI:  ### ERROR IN READING ', 'IASI', ' BUFR DATA:', &
 256                 ' STRANGE OBS TIME (YMDHM):', idate5(1:5)
 257              cycle read_loop
 258          end if
 259          call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
 260          if ( obs_time < time_slots(0) .or.  &
 261            obs_time >= time_slots(num_fgat_time) ) cycle read_loop
 262          do ifgat=1,num_fgat_time
 263             if ( obs_time >= time_slots(ifgat-1) .and.  &
 264                 obs_time  < time_slots(ifgat) ) exit
 265          end do
 266          num_iasi_global = num_iasi_global + 1
 267          ptotal(ifgat) = ptotal(ifgat) + 1
 268
 269 !    Observational info
 270          sat_zenang  = allspot(10)            ! satellite zenith angle
 271          ifov = nint(linele(1))               ! field of view
 272
 273 !    IASI fov ranges from 1 to 120.   Current angle dependent bias
 274 !    correction has a maximum of 90 scan positions.   Geometry
 275 !    of IASI scan allows us to remap 1-120 to 1-60.   Variable
 276 !    ifovn below contains the remapped IASI fov.  This value is
 277 !    passed on to and used in setuprad
 278          ifovn = (ifov-1)/2 + 1
 279          iscn = nint(linele(2))               ! scan line
 280
 281 !    Check field of view (FOVN) and satellite zenith angle (SAZA)
 282          if( ifov <= 0 .or. ifov > 120 .or. sat_zenang > 90._r_kind ) then
 283          write(unit=stdout,fmt=*) &
 284          'READ_IASI:  ### ERROR IN READING ', senname, ' BUFR DATA:', &
 285                 ' STRANGE OBS INFO(FOVN,SLNM,SAZA):', ifov, iscn, allspot(10)
 286            cycle read_loop
 287          end if
 288          if ( ifov <= 60 ) sat_zenang = -sat_zenang
 289
 290 !    Compare IASI satellite scan angle and zenith angle
 291          piece = -step_adjust
 292          if ( mod(ifovn,2) == 1) piece = step_adjust
 293          lza = ((start + float((ifov-1)/4)*step) + piece)*deg2rad
 294          sat_height_ratio = (earth_radius + linele(5))/earth_radius
 295          lzaest = asin(sat_height_ratio*sin(lza))*rad2deg
 296          if (abs(sat_zenang - lzaest) > one) then
 297                  bad_line = iscn
 298           write(unit=stdout,fmt=*) 'abs(sat_zenang - lzaest) > one',abs(sat_zenang - lzaest)
 299            cycle read_loop
 300          end if
 301 !   Clear Amount  (percent clear)
 302          call ufbrep(lnbufr,cloud_frac,1,6,iret,'FCPH')
 303          if (outside) cycle ! No good for this PE
 304          num_iasi_local = num_iasi_local + 1
 305
 306          write(unit=info%date_char, &
 307          fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 308          idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
 309          ':', idate5(5), ':', idate5(6)
 310          info%elv = 0.0  !aquaspot%selv
 311
 312          !  Make Thinning
 313          !  Map obs to thinning grid
 314          !-------------------------------------------------------------------
 315          if (thinning) then
 316             dlat_earth = info%lat
 317             dlon_earth = info%lon
 318             if (dlon_earth<zero) dlon_earth = dlon_earth+r360
 319             if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 320             dlat_earth = dlat_earth*deg2rad
 321             dlon_earth = dlon_earth*deg2rad
 322             crit = 1.
 323             call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 324             if (.not. iuse) then
 325                num_iasi_thinned=num_iasi_thinned+1
 326                cycle
 327             end if
 328          end if
 329          call ufbrep(lnbufr,cscale,3,10,iret,'STCH ENCH CHSF')
 330          if(iret /= 10) then
 331            write(unit=stdout,fmt=*)  'READ_IASI  read scale error ',iret
 332            cycle read_loop
 333          end if
 334
 335 ! The scaling factors are as follows, cscale(1) is the start channel number,
 336 !                                     cscale(2) is the end channel number,
 337 !                                     cscale(3) is the exponent scaling factor
 338 ! In our case (616 channels) there are 10 groups of cscale (dimension :: cscale(3,10))
 339 !  The units are W/m2..... you need to convert to mW/m2.... (subtract 5 from cscale(3)
 340          do i=1,10  ! convert exponent scale factor to int and change units
 341              iexponent = -(nint(cscale(3,i)) - 5)
 342              sscale(i)=ten**iexponent
 343          end do
 344
 345 !    Read IASI channel number(CHNM) and radiance (SCRA)
 346
 347          call ufbint(lnbufr,allchan,2,n_totchan,iret,'SCRA CHNM')
 348          if( iret /= n_totchan)then
 349               write(unit=stdout,fmt=*) &
 350                           'READ_IASI:  ### ERROR IN READING ', senname, ' BUFR DATA:', &
 351                 iret, ' CH DATA IS READ INSTEAD OF ',n_totchan
 352            cycle read_loop
 353          end if
 354          num_iasi_used = num_iasi_used + 1
 355          nread(inst) = nread(inst) + 1
 356          iskip = 0
 357          jstart=1
 358          do i=1,n_totchan
 359 !     check that channel number is within reason
 360             if (( allchan(1,i) > zero .and. allchan(1,i) < 99999._r_kind) .and. &  ! radiance bounds
 361                (allchan(2,i) < 8462._r_kind .and. allchan(2,i) > zero )) then     ! chan # bounds
 362 !         radiance to BT calculation
 363               radi = allchan(1,i)
 364               scaleloop: do j=jstart,10
 365                  if(allchan(2,i) >= cscale(1,j) .and. allchan(2,i) <= cscale(2,j))then
 366                     radi = allchan(1,i)*sscale(j)
 367                     jstart=j
 368                     exit scaleloop
 369                  end if
 370               end do scaleloop
 371               if (rtm_option == rtm_option_crtm) then
 372 #ifdef CRTM
 373                  call CRTM_Planck_Temperature(inst,i,radi,temperature(i))
 374 #endif
 375               else if (rtm_option == rtm_option_rttov) then
 376 #ifdef RTTOV
 377                  Effective_Temperature = LOG( ( coefs(inst)%coef%planck1(i) / radi ) + 1 )
 378                  Effective_Temperature = coefs(inst)%coef%planck2(i)/Effective_Temperature
 379                  temperature(i) = ( Effective_Temperature - coefs(inst)%coef%ff_bco(i) ) / &
 380                                   coefs(inst)%coef%ff_bcs(i)
 381 #endif
 382               end if
 383               if(temperature(i) < tbmin .or. temperature(i) > tbmax ) then
 384                  temperature(i) = min(tbmax,max(zero,temperature(i)))
 385
 386               end if
 387             else           ! error with channel number or radiance
 388 !             write(6,*)'READ_IASI:  iasi chan error',i,allchan(1,i), allchan(2,i)
 389               temperature(i) = min(tbmax,max(zero,temperature(i)))
 390
 391             end if
 392          end do
 393          if(iskip > 0)write(6,*) ' READ_IASI : iskip > 0 ',iskip
 394 !         if( iskip >= 10 )cycle read_loop
 395          data_all(5) = sat_zenang             ! satellite zenith angle (deg)
 396          data_all(6) = allspot(11)            ! satellite azimuth angle (deg)
 397          data_all(9) = allspot(12)            ! solar zenith angle (deg)
 398          data_all(10)= allspot(13)            ! solar azimuth angle (deg)
 399          do l=1,nchan
 400             data_all(l+nreal) = temperature(l)   ! brightness temerature
 401          end do
 402
 403 !  4.0   assign information to sequential radiance structure
 404 !--------------------------------------------------------------------------
 405          allocate ( p % tb_inv (1:nchan) )
 406          p%info             = info
 407          p%loc              = loc
 408          p%landsea_mask     = 1
 409          p%scanpos          = ifovn
 410          p%satzen           = data_all(5)
 411          p%satazi           = data_all(6)  ! look angle (deg) ! airsspot%bearaz
 412          p%solzen           = data_all(9)
 413          p%solazi           = data_all(10)
 414          p%tb_inv(1:nchan)  = data_all(nreal+1:nreal+nchan)
 415          p%sensor_index     = inst
 416          p%ifgat            = ifgat
 417
 418          allocate (p%next)   ! add next data
 419          p => p%next
 420          nullify (p%next)
 421      end do read_loop
 422   end do
 423   call closbf(lnbufr)
 424
 425   !Deallocate temporary array for next bufrfile do loop
 426   deallocate(data_all)
 427 end do bufrfile
 428
 429    if (thinning .and. num_iasi_global > 0 ) then
 430
 431 #ifdef DM_PARALLEL
 432
 433       ! Get minimum crit and associated processor index.
 434       j = 0
 435       do ifgat = 1, num_fgat_time
 436          do n = 1, iv%num_inst
 437             j = j + thinning_grid(n,ifgat)%itxmax
 438          end do
 439       end do
 440
 441       allocate ( in  (j) )
 442       allocate ( out (j) )
 443       j = 0
 444       do ifgat = 1, num_fgat_time
 445          do n = 1, iv%num_inst
 446             do i = 1, thinning_grid(n,ifgat)%itxmax
 447                j = j + 1
 448                in(j) = thinning_grid(n,ifgat)%score_crit(i)
 449             end do
 450          end do
 451       end do
 452       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 453
 454       call wrf_dm_bcast_real (out, j)
 455
 456       j = 0
 457       do ifgat = 1, num_fgat_time
 458          do n = 1, iv%num_inst
 459             do i = 1, thinning_grid(n,ifgat)%itxmax
 460                j = j + 1
 461                if ( ABS(out(j)-thinning_grid(n,ifgat)%score_crit(i)) > 1.0E-10 ) thinning_grid(n,ifgat)%ibest_obs(i) = 0
 462             end do
 463          end do
 464       end do
 465
 466       deallocate( in  )
 467       deallocate( out )
 468
 469 #endif
 470
 471       ! Delete the nodes which being thinning out
 472       p => head
 473       prev => head
 474       head_found = .false.
 475       num_iasi_used_tmp = num_iasi_used
 476       do j = 1, num_iasi_used_tmp
 477          n = p%sensor_index
 478          ifgat = p%ifgat
 479          found = .false.
 480
 481          do i = 1, thinning_grid(n,ifgat)%itxmax
 482             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 483                found = .true.
 484                exit
 485             end if
 486          end do
 487
 488          ! free current data
 489          if ( .not. found ) then
 490             current => p
 491             p => p%next
 492             if ( head_found ) then
 493                prev%next => p
 494             else
 495                head => p
 496                prev => p
 497             end if
 498             deallocate ( current % tb_inv )
 499             deallocate ( current )
 500             num_iasi_thinned = num_iasi_thinned + 1
 501             num_iasi_used = num_iasi_used - 1
 502             nread(n) = nread(n) - 1
 503             continue
 504          end if
 505
 506          if ( found .and. head_found ) then
 507             prev => p
 508             p => p%next
 509             continue
 510          end if
 511
 512          if ( found .and. .not. head_found ) then
 513             head_found = .true.
 514             head => p
 515             prev => p
 516             p => p%next
 517          end if
 518
 519       end do
 520
 521    end if  ! End of thinning
 522
 523    iv%total_rad_pixel   = iv%total_rad_pixel + num_iasi_used
 524    iv%total_rad_channel = iv%total_rad_channel + num_iasi_used*nchan
 525
 526    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_iasi_used
 527    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_iasi_global
 528
 529    do i = 1, num_fgat_time
 530       ptotal(i) = ptotal(i) + ptotal(i-1)
 531       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 532    end do
 533    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 534       write(unit=message(1),fmt='(A,I10,A,I10)') &
 535           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 536       call da_warning(__FILE__,__LINE__,message(1:1))
 537    endif
 538
 539    write(unit=message(1),fmt='(a)') '   num_iasi_file num_iasi_global  num_iasi_local   num_iasi_used num_iasi_thinned'
 540    write(unit=message(2),fmt='(5(6x,i10))') num_iasi_file, num_iasi_global, num_iasi_local, num_iasi_used, num_iasi_thinned
 541    call da_message(message(1:2))
 542
 543    !  5.0 allocate innovation radiance structure
 544    !----------------------------------------------------------------
 545
 546
 547    do i = 1, iv%num_inst
 548
 549       if (nread(i) < 1) cycle
 550       iv%instid(i)%num_rad  = nread(i)
 551       iv%instid(i)%info%nlocal = nread(i)
 552       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 553         'Allocating space for radiance innov structure', &
 554          i, iv%instid(i)%rttovid_string, iv%instid(i)%num_rad
 555       call da_allocate_rad_iv (i, nchan, iv)
 556
 557    end do
 558
 559    !  6.0 assign sequential structure to innovation structure
 560    !-------------------------------------------------------------
 561    nread(1:rtminit_nsensor) = 0
 562    p => head
 563
 564    do n = 1, num_iasi_used
 565       i = p%sensor_index
 566       nread(i) = nread(i) + 1
 567
 568       call da_initialize_rad_iv (i, nread(i), iv, p)
 569
 570       current => p
 571       p => p%next
 572       ! free current data
 573       deallocate ( current % tb_inv )
 574       deallocate ( current )
 575    end do
 576    deallocate ( p )
 577    deallocate (nread)
 578    deallocate (ptotal)
 579
 580    call closbf(lnbufr)
 581    close(lnbufr)
 582
 583 !   call da_free_unit(lnbufr)
 584
 585    if (trace_use) call da_trace_exit("da_read_obs_bufriasi")
 586 #else
 587    call da_error(__FILE__,__LINE__,(/"Needs to be compiled with a BUFR library"/))
 588 #endif
 589
 590
 591 end subroutine da_read_obs_bufriasi