var/da/da_radiance/da_read_obs_bufrseviri.inc

   1 subroutine da_read_obs_bufrseviri (obstype,iv,infile)
   2
   3 ! subprogram:    read_seviri                  read bufr format seviri data
   4    !--------------------------------------------------------
   5    !  Purpose: read in NCEP bufr eos AIRS/AMSUA/HSB 1b data
   6    !            to innovation structure
   7    !
   8    !   METHOD: use F90 sequantial data structure to avoid read file twice
   9    !            so that da_scan_bufrairs is not necessary any more.
  10    !            1. read file radiance data in sequential data structure
  11    !            2. do gross QC check
  12    !            3. assign sequential data structure to innovation structure
  13    !                  and deallocate sequential data structure
  14    !
  15    !  HISTORY: 2013/03/26 - Creation            Hongli Wang
  16    !
  17    !------------------------------------------------------------------------------
  18
  19   implicit none
  20
  21   real(r_kind)     :: POinT001 =   0.001_r_kind
  22   real(r_kind)     :: POinT01  =   0.01_r_kind
  23   real(r_kind)     :: TEN      =  10.0_r_kind
  24   real(r_kind)     :: R45      =  45.0_r_kind
  25   real(r_kind)     :: R60      =  60.0_r_kind
  26   real(r_kind)     :: R90      =  90.0_r_kind
  27   real(r_kind)     :: R180     = 180.0_r_kind
  28   real(r_kind)     :: R360     = 360.0_r_kind
  29
  30   character(9)      ,  intent (in)  :: obstype
  31   type (iv_type)    ,intent (inout) :: iv
  32 #ifdef BUFR
  33   real(kind=8)    ::  obs_time
  34   type (datalink_type), pointer  :: head, p, current, prev
  35   type(info_type)                :: info
  36   type(model_loc_type)           :: loc
  37   type(model_loc_type)           :: loc_fov
  38
  39 !!! for seviri
  40   character(80),  intent (in) :: infile
  41
  42   character(8)       :: subset,subcsr,subasr
  43   character(80)      :: hdrsevi                             ! seviri header
  44
  45   integer(i_kind)    :: nchanl,ilath,ilonh,ilzah,iszah,irec
  46   integer(i_kind)    :: nhdr,nchn,ncld,nbrst !,idate,lnbufr
  47   integer(i_kind)    :: ireadmg,ireadsb,iret
  48
  49   real(r_double),allocatable,dimension(:)   :: hdr
  50   real(r_double),allocatable,dimension(:,:) :: datasev1,datasev2
  51
  52   logical            :: clrsky,allsky,allchnmiss
  53   real               :: rclrsky
  54   integer :: kidsat
  55   integer(i_kind)    :: idate5(6)
  56   integer (i_kind), allocatable :: ptotal(:), nread(:)
  57   integer(i_kind)    :: idate, im, iy, idd, ihh
  58 !!! end for seviri
  59
  60   ! Number of channels for sensors in BUFR
  61   integer(i_kind),parameter :: nchan = 8
  62   integer(i_kind),parameter :: n_totchan  = 12
  63   integer(i_kind),parameter :: maxinfo    =  33
  64   integer(i_kind)   :: inst,platform_id,satellite_id,sensor_id
  65   real(r_kind)      :: crit
  66   integer(i_kind)   :: ifgat, iout, iobs
  67   logical           :: outside, outside_all, iuse,outside_fov
  68
  69   integer           :: numbufr,ibufr,jj
  70   logical :: found, head_found
  71   real(r_kind)      :: step, start,step_adjust
  72   character(len=4)  :: senname
  73   integer(i_kind)   :: size,size_tmp
  74   character(10)     :: date
  75   real(r_kind)      :: sstime, tdiff, t4dv
  76   integer(i_kind)   :: nmind
  77
  78   ! Other work variables
  79   real(r_kind)     :: clr_amt,piece
  80   real(r_kind)     :: dlon, dlat
  81   real(r_kind)     :: dlon_earth,dlat_earth,dlon_earth_deg,dlat_earth_deg
  82   real(r_kind)     :: lza, lzaest,sat_height_ratio
  83   real(r_kind)     :: pred, crit1, dist1
  84   real(r_kind)     :: sat_zenang
  85   real(r_kind)     :: radi
  86   real(r_kind)     :: tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr
  87   real(r_kind),dimension(0:4) :: rlndsea
  88   real(r_kind),dimension(0:3) :: sfcpct
  89   real(r_kind),dimension(0:3) :: ts
  90   real(r_kind),dimension(10) :: sscale
  91   real(r_kind),dimension(n_totchan) :: temperature
  92   real(r_kind),allocatable,dimension(:):: data_all
  93   real(r_kind) disterr,disterrmax,rlon00,rlat00
  94
  95   logical          :: assim,valid
  96   logical          :: seviri
  97   logical          :: data_on_edges,luse
  98   integer(i_kind)  :: nreal, ityp,ityp2, isflg
  99   integer(i_kind)  :: ifov, instr, ioff, ilat, ilon, sensorindex
 100   integer(i_kind)  :: num_seviri_file
 101   integer(i_kind)  :: num_seviri_local, num_seviri_global, num_seviri_used, num_seviri_thinned
 102   integer(i_kind)  :: num_seviri_used_tmp
 103   integer(i_kind)  :: i, j, l, iskip, ifovn, bad_line
 104
 105   integer(i_kind)  :: itx, k, nele, itt, n
 106   integer(i_kind)  :: iexponent
 107   integer(i_kind)  :: idomsfc
 108   integer(i_kind)  :: ntest
 109   integer(i_kind)  :: error_status
 110   integer(i_kind)  :: num_bufr(7)
 111
 112   integer          :: iost, lnbufr
 113   character(20)    ::filename
 114   real, allocatable :: in(:), out(:)
 115
 116   ! Set standard parameters
 117   integer(i_kind),parameter:: ichan=-999  ! fov-based surface code is not channel specific for seviri
 118   real(r_kind),parameter:: expansion=one  ! exansion factor for fov-based location code.
 119   real(r_kind),parameter:: tbmin  = 50._r_kind
 120   real(r_kind),parameter:: tbmax  = 550._r_kind
 121   real(r_kind),parameter:: earth_radius = 6371000._r_kind
 122
 123   ilath=8                      ! the position of latitude in the header
 124   ilonh=9                      ! the position of longitude in the header
 125   ilzah=10                     ! satellite zenith angle
 126   iszah=11                     ! solar zenith angle
 127   subcsr='NC021043'            ! sub message
 128   subasr='NC021042'            ! sub message
 129
 130    if(trace_use) call da_trace_entry("da_read_obs_bufrseviri")
 131
 132   !  0.0  Initialize variables
 133   !-----------------------------------
 134   sensor_id = 21 !seviri
 135   disterrmax=zero
 136   ntest=0
 137   nreal  = maxinfo
 138   seviri= obstype == 'seviri'
 139
 140   bad_line=-1
 141   step  = 1.0
 142   start = 0.0
 143   step_adjust = 1.0_r_kind
 144   senname = 'SEVIRI'
 145   num_bufr(:)=0
 146   numbufr=0
 147   allocate(nread(1:rtminit_nsensor))
 148   allocate(ptotal(0:num_fgat_time))
 149   nread(1:rtminit_nsensor) = 0
 150   ptotal(0:num_fgat_time) = 0
 151   iobs = 0                 ! for thinning, argument is inout
 152   num_seviri_file    = 0
 153   num_seviri_local   = 0
 154   num_seviri_global  = 0
 155   num_seviri_used    = 0
 156   num_seviri_thinned = 0
 157
 158   ! 1.0 Assign file names and prepare to read bufr files
 159   !-------------------------------------------------------------------------
 160
 161    if (num_fgat_time>1) then
 162       do i=1,7
 163          call da_get_unit(lnbufr)
 164          write(filename,fmt='(A,2I1,A)') trim(infile),0,i,'.bufr'
 165          open(unit   = lnbufr, FILE   = trim(filename),iostat =  iost, form = 'unformatted', STATUS = 'OLD')
 166          if (iost == 0) then
 167             numbufr=numbufr+1
 168             num_bufr(numbufr)=i
 169          else
 170             close (lnbufr)
 171          end if
 172          call da_free_unit(lnbufr)
 173       end do
 174    else
 175      numbufr=1
 176    end if
 177
 178
 179    if (numbufr==0) numbufr=1
 180
 181   bufrfile:  do ibufr=1,numbufr
 182    if (num_fgat_time==1) then
 183       filename=trim(infile)//'.bufr'
 184    else
 185       if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then
 186          filename=trim(infile)//'.bufr'
 187       else
 188          write(filename,fmt='(A,2I1,A)') trim(infile),0,num_bufr(ibufr),'.bufr'
 189       end if
 190    end if
 191
 192    ! dont change, WRFDA uses specified units to read radiance data
 193    lnbufr = 92
 194    open(unit=lnbufr,file=trim(filename),form='unformatted', &
 195       iostat = iost, status = 'old' ) !,convert='little_endian')
 196    if (iost /= 0) then
 197       call da_warning(__FILE__,__LINE__, &
 198          (/"Cannot open file "//infile/))
 199       if(trace_use) call da_trace_exit("da_read_obs_bufrsevri")
 200       return
 201    end if
 202
 203    ! Open BUFR table
 204    call openbf(lnbufr,'IN',lnbufr)
 205    call datelen(10)
 206    call readmg(lnbufr,subset,idate,iret)
 207
 208    ! Check the data set
 209   if( iret/=0) then
 210      write(message(1),fmt='(A)')'SKIP PROCESSING OF SEVIRI FILE'
 211      write(message(2),fmt='(A,I2,A)')'infile = ', lnbufr, infile
 212      call da_warning(__FILE__,__LINE__,message(1:2))
 213      if(trace_use) call da_trace_exit("da_read_obs_bufrseviri")
 214      return
 215   endif
 216
 217   clrsky=.false.
 218   allsky=.false.
 219   if(subset == subcsr) then
 220      clrsky=.true.
 221      write(message(1),fmt='(A)')'READ SEVIRI FILE'
 222      write(message(2),fmt='(A,L1,2A)')'clrsky= ', clrsky,' subset= ', subset
 223      call da_message(message(1:2))
 224   elseif(subset == subasr) then
 225      allsky=.true.
 226      write(message(1),fmt='(A)')'READ SEVIRI FILE'
 227      write(message(2),fmt='(A,L1,2A)')'allsky= ', allsky,' subset= ', subset
 228      call da_message(message(1:2))
 229   else
 230      write(message(1),fmt='(A)')'SKIP PROCESSING OF UNKNOWN SEVIRI FILE'
 231      write(message(2),fmt='(A,I2,3A)')'infile=', lnbufr, infile,' subset=', subset
 232      call da_warning(__FILE__,__LINE__,message(1:2))
 233      if(trace_use) call da_trace_exit("da_read_obs_bufrseviri")
 234      return
 235   endif
 236
 237   ! Set BUFR string based on seviri data set
 238   if (clrsky) then
 239      hdrsevi='SAID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAZA SOZA'
 240      nhdr=11
 241      nchn=12
 242      ncld=nchn
 243      nbrst=nchn
 244   else if (allsky) then
 245      hdrsevi='SAID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH'
 246      nhdr=9
 247      nchn=11
 248      ncld=2
 249      nbrst=nchn*6                ! channel dependent: all, clear, cloudy, low,
 250                                  !middle and high clouds
 251   endif
 252   allocate(datasev1(1,ncld))     ! not channel dependent
 253   allocate(datasev2(1,nbrst))    ! channel dependent: all, clear, cloudy, low,
 254                                  !middle and high clouds
 255   allocate(hdr(nhdr))
 256
 257
 258
 259    iy=0
 260    im=0
 261    idd=0
 262    ihh=0
 263
 264    sensorindex=1
 265    write(unit=date,fmt='( i10)') idate
 266    read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
 267    write(unit=stdout,fmt='(a,4i4,2x,a)') &
 268       'Bufr file date is ',iy,im,idd,ihh,trim(infile)
 269
 270    ! 2.0 Loop to read bufr file and assign information to a sequential structure
 271    !-------------------------------------------------------------------------
 272
 273    ! Allocate arrays to hold data
 274    nele=nreal+nchan
 275    allocate(data_all(nele))
 276    if ( ibufr == 1 ) then
 277       allocate (head)
 278       nullify  ( head % next )
 279       p => head
 280    end if
 281
 282
 283 ! Big loop to read data file
 284
 285   do while(ireadmg(lnbufr,subset,idate)>=0)
 286
 287      read_loop: do while (ireadsb(lnbufr)==0)
 288          num_seviri_file = num_seviri_file + 1
 289
 290     ! Read SEVIRI information
 291          call ufbint(lnbufr,hdr,nhdr,1,iret,hdrsevi)
 292
 293         kidsat = nint(hdr(1))
 294         ! SAID 55 is meteosat-8  or msg-1
 295         ! SAID 56 is meteosat-9  or msg-2
 296         ! SAID 57 is meteosat-10 or msg-3
 297         ! SAID 70 is meteosat-11 or msg-4
 298         if ( ( kidsat > 70) .or. ( kidsat < 55) ) then
 299            write(unit=message(1),fmt='(A,I6)') 'Unknown platform: ', kidsat
 300            call da_warning(__FILE__,__LINE__,message(1:1))
 301         end if
 302         platform_id  = 12  ! MSG - Meteosat Second Generation
 303         if ( kidsat == 55 ) then
 304             satellite_id = 1
 305         else if ( kidsat == 56 ) then
 306             satellite_id = 2
 307         else if ( kidsat == 57 ) then
 308             satellite_id = 3
 309         else if ( kidsat == 70 ) then
 310             satellite_id = 4
 311         end if
 312
 313         if (clrsky) then     ! asr bufr has no sza
 314         ! remove the obs whose satellite zenith angles larger than 65 degree
 315            if ( hdr(ilzah) > 65.0 ) cycle read_loop
 316         end if
 317
 318         call ufbint(lnbufr,datasev1,1,ncld,iret,'NCLDMNT')
 319
 320         if(iret /= 1) cycle read_loop
 321         do n=1,ncld
 322            if(datasev1(1,n)>= 0.0 .and. datasev1(1,n) <= 100.0 ) then
 323               rclrsky=datasev1(1,n)
 324                ! first QC filter out data with less clear sky fraction
 325                if ( rclrsky < 70.0 ) cycle read_loop
 326                !if ( rclrsky < 90.0 ) cycle read_loop
 327            end if
 328         end do
 329
 330         call ufbrep(lnbufr,datasev2,1,nbrst,iret,'TMBRST')
 331
 332         ! Check if data of channel 1-11 are missing
 333
 334         allchnmiss=.true.
 335         do n=1,11
 336            if(datasev2(1,n)<500.)  then
 337               allchnmiss=.false.
 338            end if
 339         end do
 340         if(allchnmiss) cycle read_loop
 341
 342         ! Check observing position
 343          info%lat  =  hdr(ilath)           ! latitude
 344          info%lon  =  hdr(ilonh)           ! longitude
 345          if( abs(info%lat) > R90  .or. abs(info%lon) > R360 .or. &
 346          (abs(info%lat) == R90 .and. info%lon /= ZERO) )then
 347          write(unit=stdout,fmt=*) &
 348          'READ_SEVIRI:  ### ERROR IN READING ', senname, ' BUFR DATA:', &
 349                ' STRANGE OBS POINT (LAT,LON):', info%lat, info%lon
 350             cycle read_loop
 351          end if
 352
 353          call da_llxy (info, loc, outside, outside_all)
 354          if (outside_all) cycle
 355              inst = 0
 356          do i = 1, rtminit_nsensor
 357             if (platform_id  == rtminit_platform(i) &
 358                .and. satellite_id == rtminit_satid(i)    &
 359                .and. sensor_id    == rtminit_sensor(i)) then
 360                inst = i
 361                exit
 362             end if
 363          end do
 364          if (inst == 0) cycle read_loop
 365
 366         ! Check obs time
 367          idate5(1) = nint(hdr(2)) ! year
 368          idate5(2) = nint(hdr(3)) ! month
 369          idate5(3) = nint(hdr(4)) ! day
 370          idate5(4) = nint(hdr(5)) ! hour
 371          idate5(5) = nint(hdr(6)) ! minute
 372          idate5(6) = nint(hdr(7)) ! second
 373
 374          if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. &
 375              idate5(2) < 1    .or. idate5(2) >   12 .or. &
 376              idate5(3) < 1    .or. idate5(3) >   31 .or. &
 377              idate5(4) < 0    .or. idate5(4) >   24 .or. &
 378              idate5(5) < 0    .or. idate5(5) >   60 ) then
 379
 380             write(message(1),fmt='(A)')'ERROR IN READING SEVIRI BUFR DATA'
 381             write(message(2),fmt='(A,5I8)')'STRANGE OBS TIME (YMDHM):', idate5(1:5)
 382             call da_warning(__FILE__,__LINE__,message(1:2))
 383             cycle read_loop
 384          end if
 385
 386          call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
 387          if ( obs_time < time_slots(0) .or.  &
 388            obs_time >= time_slots(num_fgat_time) ) cycle read_loop
 389          do ifgat=1,num_fgat_time
 390             if ( obs_time >= time_slots(ifgat-1) .and.  &
 391                 obs_time  < time_slots(ifgat) ) exit
 392          end do
 393          num_seviri_global = num_seviri_global + 1
 394          ptotal(ifgat) = ptotal(ifgat) + 1
 395
 396          if (outside) cycle ! No good for this PE
 397          num_seviri_local = num_seviri_local + 1
 398
 399          write(unit=info%date_char, &
 400          fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 401          idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
 402          ':', idate5(5), ':', idate5(6)
 403          info%elv = 0.0  !aquaspot%selv
 404
 405          ! 3.0  Make Thinning
 406          ! Map obs to thinning grid
 407          !-------------------------------------------------------------------
 408          if (thinning) then
 409             dlat_earth = info%lat
 410             dlon_earth = info%lon
 411             if (dlon_earth<zero) dlon_earth = dlon_earth+r360
 412             if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 413             dlat_earth = dlat_earth*deg2rad
 414             dlon_earth = dlon_earth*deg2rad
 415             crit = 1.
 416             call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 417             if (.not. iuse) then
 418                num_seviri_thinned=num_seviri_thinned+1
 419                cycle
 420             end if
 421          end if
 422
 423          ! data SEVIRI channel number(CHNM) and radiance (SCRA)
 424
 425          num_seviri_used = num_seviri_used + 1
 426          nread(inst) = nread(inst) + 1
 427          iskip = 0
 428          do i=1,n_totchan
 429             ! check that tb is within reasonal bound
 430             if ( datasev2(1,i) < tbmin .or. datasev2(1,i) > tbmax ) then
 431                temperature(i) = missing_r
 432             else
 433                temperature(i) = datasev2(1,i)
 434             end if
 435          end do
 436
 437          if(iskip > 0)write(6,*) ' READ_SEVIRI : iskip > 0 ',iskip
 438
 439          do l=1,nchan
 440             data_all(l+nreal) = temperature(l+3)   ! brightness temerature
 441          end do
 442
 443     ! 4.0 assign information to sequential radiance structure
 444     !--------------------------------------------------------------------------
 445     !        iscan = nint(hdr(ilzah))+1.001_r_kind
 446          allocate ( p % tb_inv (1:nchan ))
 447          p%info             = info
 448          p%loc              = loc
 449          p%landsea_mask     = 1
 450          p%scanpos          = nint(hdr(ilzah))+1.001_r_kind
 451          p%satzen           = hdr(ilzah)               ! satellite zenith angle (deg)
 452          p%satazi           = 0.0                      ! satellite azimuth angle (deg)
 453          p%solzen           = 0.0                      ! solar zenith angle (deg)
 454          p%solazi           = 0.0                      ! solar azimuth angle (deg)
 455          p%tb_inv(1:nchan)  = data_all(nreal+1:nreal+nchan)
 456          p%sensor_index     = inst
 457          p%ifgat            = ifgat
 458
 459          allocate (p%next)   ! add next data
 460          p => p%next
 461          nullify (p%next)
 462      end do read_loop
 463   end do
 464   call closbf(lnbufr)
 465
 466   !Deallocate temporary arrays for next bufrfile do loop
 467   deallocate(datasev1)
 468   deallocate(datasev2)
 469   deallocate(hdr)
 470   deallocate(data_all)
 471 end do bufrfile
 472
 473    if (thinning .and. num_seviri_global > 0 ) then
 474
 475 #ifdef DM_PARALLEL
 476
 477       ! Get minimum crit and associated processor index.
 478       j = 0
 479       do ifgat = 1, num_fgat_time
 480          do n = 1, iv%num_inst
 481             j = j + thinning_grid(n,ifgat)%itxmax
 482          end do
 483       end do
 484
 485       allocate ( in  (j) )
 486       allocate ( out (j) )
 487       j = 0
 488       do ifgat = 1, num_fgat_time
 489          do n = 1, iv%num_inst
 490             do i = 1, thinning_grid(n,ifgat)%itxmax
 491                j = j + 1
 492                in(j) = thinning_grid(n,ifgat)%score_crit(i)
 493             end do
 494          end do
 495       end do
 496       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 497
 498       call wrf_dm_bcast_real (out, j)
 499
 500       j = 0
 501       do ifgat = 1, num_fgat_time
 502          do n = 1, iv%num_inst
 503             do i = 1, thinning_grid(n,ifgat)%itxmax
 504                j = j + 1
 505                if ( ABS(out(j)-thinning_grid(n,ifgat)%score_crit(i)) > 1.0E-10 ) thinning_grid(n,ifgat)%ibest_obs(i) = 0
 506             end do
 507          end do
 508       end do
 509
 510       deallocate( in  )
 511       deallocate( out )
 512
 513 #endif
 514
 515       ! Delete the nodes which being thinning out
 516       p => head
 517       prev => head
 518       head_found = .false.
 519       num_seviri_used_tmp = num_seviri_used
 520       do j = 1, num_seviri_used_tmp
 521          n = p%sensor_index
 522          ifgat = p%ifgat
 523          found = .false.
 524
 525          do i = 1, thinning_grid(n,ifgat)%itxmax
 526             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 527                found = .true.
 528                exit
 529             end if
 530          end do
 531
 532          ! free current data
 533          if ( .not. found ) then
 534             current => p
 535             p => p%next
 536             if ( head_found ) then
 537                prev%next => p
 538             else
 539                head => p
 540                prev => p
 541             end if
 542             deallocate ( current % tb_inv )
 543             deallocate ( current )
 544             num_seviri_thinned = num_seviri_thinned + 1
 545             num_seviri_used = num_seviri_used - 1
 546             nread(n) = nread(n) - 1
 547             continue
 548          end if
 549
 550          if ( found .and. head_found ) then
 551             prev => p
 552             p => p%next
 553             continue
 554          end if
 555
 556          if ( found .and. .not. head_found ) then
 557             head_found = .true.
 558             head => p
 559             prev => p
 560             p => p%next
 561          end if
 562
 563       end do
 564
 565    end if  ! End of thinning
 566
 567    iv%total_rad_pixel   = iv%total_rad_pixel + num_seviri_used
 568    iv%total_rad_channel = iv%total_rad_channel + num_seviri_used*nchan
 569
 570    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_seviri_used
 571    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_seviri_global
 572
 573    do i = 1, num_fgat_time
 574       ptotal(i) = ptotal(i) + ptotal(i-1)
 575       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 576    end do
 577    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 578       write(unit=message(1),fmt='(A,I10,A,I10)') &
 579           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 580       call da_warning(__FILE__,__LINE__,message(1:1))
 581    endif
 582
 583    write(unit=stdout,fmt='(a)') '   num_seviri_file num_seviri_global  num_seviri_local   num_seviri_used num_seviri_thinned'
 584    write(stdout,'(5(8x,i10))') num_seviri_file, num_seviri_global, num_seviri_local, num_seviri_used, num_seviri_thinned
 585
 586
 587
 588    !  5.0 allocate innovation radiance structure
 589    !----------------------------------------------------------------
 590
 591
 592    do i = 1, iv%num_inst
 593
 594       if (nread(i) < 1) cycle
 595       iv%instid(i)%num_rad  = nread(i)
 596       iv%instid(i)%info%nlocal = nread(i)
 597       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 598         'Allocating space for radiance innov structure', &
 599          i, iv%instid(i)%rttovid_string, iv%instid(i)%num_rad
 600       call da_allocate_rad_iv (i, nchan, iv)
 601
 602    end do
 603
 604    !  6.0 assign sequential structure to innovation structure
 605    !-------------------------------------------------------------
 606    nread(1:rtminit_nsensor) = 0
 607    p => head
 608
 609    do n = 1, num_seviri_used
 610       i = p%sensor_index
 611       nread(i) = nread(i) + 1
 612
 613       call da_initialize_rad_iv (i, nread(i), iv, p)
 614
 615       current => p
 616       p => p%next
 617       ! free current data
 618       deallocate ( current % tb_inv )
 619       deallocate ( current )
 620    end do
 621    deallocate ( p )
 622    deallocate (nread)
 623    deallocate (ptotal)
 624
 625    call closbf(lnbufr)
 626    close(lnbufr)
 627
 628    call da_free_unit(lnbufr)
 629
 630    if(trace_use) call da_trace_exit("da_read_obs_bufrseviri")
 631 #else
 632    call da_error(__FILE__,__LINE__,(/"Needs to be compiled with a BUFR library"/))
 633 #endif
 634
 635
 636 end subroutine da_read_obs_bufrseviri