var/da/da_radiance/da_read_obs_bufrairs.inc

   1 subroutine da_read_obs_bufrairs(obstype,iv,infile)
   2    !--------------------------------------------------------
   3    !  Purpose: read in NCEP bufr eos AIRS/AMSUA/HSB 1b data
   4    !            to innovation structure
   5    !
   6    !   METHOD: use F90 sequantial data structure to avoid read file twice
   7    !            so that da_scan_bufrairs 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    !
  13    !  HISTORY: 2006/01/03 - Creation            Zhiquan Liu
  14    !           2008/03/01 - VISNIR Cloud Cover  Tom Auligne
  15    !           2008/04/01 - Warmest FoV         Tom Auligne
  16    !           2008/07/15 - BUFR format change  Hui-Chuan Lin
  17    !            NCEP msg type NC021250 (center FOV data) discontinued after Aug 2007
  18    !            replaced by NC021249 (every FOV data)
  19    !
  20    !------------------------------------------------------------------------------
  21
  22
  23
  24
  25  implicit none
  26
  27   character(9)      ,  intent (in)  :: obstype
  28   character(100)    ,  intent (in)  :: infile
  29   type (iv_type)    ,intent (inout) :: iv
  30
  31 #ifdef BUFR
  32
  33 ! Number of channels for sensors in BUFR
  34   integer(i_kind),parameter :: N_AIRSCHAN = 281  !- 281 subset ch out of 2378 ch for AIRS
  35   integer(i_kind),parameter :: N_AMSUCHAN =  15
  36   integer(i_kind),parameter :: N_HSBCHAN  =   4
  37   integer(i_kind),parameter :: N_MAXCHAN  = 350
  38   integer(i_kind),parameter :: maxinfo    =  12
  39
  40 ! BUFR format for AQUASPOT (SPITSEQN)
  41   integer(i_kind),parameter :: N_AQUASPOT_LIST = 25
  42   type aquaspot_list
  43      sequence
  44      real(r_double) :: said   ! Satellite identifier
  45      real(r_double) :: orbn   ! Orbit number
  46      real(r_double) :: slnm   ! Scan line number
  47      real(r_double) :: mjfc   ! Major frame count
  48      real(r_double) :: selv   ! Height of station
  49      real(r_double) :: soza   ! Solar zenith angle
  50      real(r_double) :: solazi ! Solar azimuth angle
  51      real(r_double) :: intms(2,9) ! SATELLITE inSTRUMENT TEMPERATURES
  52   end type aquaspot_list
  53   real(r_double), dimension(1:N_AQUASPOT_LIST) :: aquaspot_list_array
  54
  55
  56 ! BUFR format for AIRSSPOT (SITPSEQN)
  57   integer(i_kind),parameter :: N_AIRSSPOT_LIST = 12
  58   type airsspot_list
  59      sequence
  60      real(r_double) :: siid  ! Satellite instruments
  61      real(r_double) :: year
  62      real(r_double) :: mnth
  63      real(r_double) :: days
  64      real(r_double) :: hour
  65      real(r_double) :: minu
  66      real(r_double) :: seco
  67      real(r_double) :: clath ! Latitude (high accuracy)
  68      real(r_double) :: clonh ! Longitude (high accuracy)
  69      real(r_double) :: saza  ! Satellite zenith angle
  70      real(r_double) :: bearaz ! Bearing or azimuth
  71      real(r_double) :: fovn  ! Field of view number
  72   end type airsspot_list
  73   real(r_double), dimension(1:N_AIRSSPOT_LIST) :: airsspot_list_array
  74
  75
  76 ! BUFR format for AIRSCHAN (SCBTSEQN)
  77   integer(i_kind),parameter :: N_AIRSCHAN_LIST = 4
  78   type airschan_list
  79      sequence
  80      real(r_double) :: chnm    ! Channel number
  81      real(r_double) :: logrcw  ! Log-10 of (Temperature-radiance central wavenumber
  82      real(r_double) :: acqf    ! Channel quality flags for ATOVS
  83      real(r_double) :: tmbrst  ! Brightness temperature
  84   end type airschan_list
  85   real(r_double), dimension(1:N_AIRSCHAN_LIST,1:N_MAXCHAN) :: airschan_list_array
  86
  87 ! BUFR talble file sequencial number
  88   character(len=512)  :: table_file
  89
  90
  91 ! Variables for BUFR IO
  92   type(aquaspot_list) :: aquaspot
  93   type(airsspot_list) :: airsspot
  94   type(airschan_list) :: airschan(N_MAXCHAN)
  95
  96   real(r_kind)      :: step, start
  97   real(r_kind)      :: airdata(N_AIRSCHAN+maxinfo)
  98   character(len=8)  :: subset
  99   character(len=4)  :: senname
 100   character(len=8)  :: spotname
 101   character(len=8)  :: channame
 102   integer(i_kind)   :: nchan,nchanr
 103   integer(i_kind)   :: iret
 104
 105
 106 ! Work variables for time
 107   integer(i_kind)   :: idate, ifgat
 108   integer(i_kind)   :: idate5(6)
 109   character(len=10) :: date
 110   integer(i_kind)   :: nmind
 111   integer(i_kind)   :: iy, im, idd, ihh
 112   real*8            :: obs_time
 113
 114
 115 ! Other work variables
 116   integer(i_kind)  :: nreal
 117   integer(i_kind)  :: k, iobsout
 118   integer(i_kind),dimension(19)::icount
 119   real(r_kind)     :: rlat, rlon, dx, dy, dx1, dy1, sstx, dlon, dlat
 120   real(r_kind)     :: sza, timedif, pred, crit1
 121   integer(i_kind)  :: klat1, klon1, klatp1, klonp1
 122   integer(i_kind)  :: ifov,size
 123   integer(i_kind)  :: num_eos_file, num_eos_global, num_eos_local, num_eos_kept, num_eos_thinned
 124   integer(i_kind)  :: inst,platform_id,satellite_id,sensor_id
 125   logical          :: iflag,outside,outside_all
 126   integer(i_kind)  :: i, l, n, error, airs_table_unit
 127   integer          :: iost, lnbufr
 128   real(r_kind),allocatable,dimension(:,:):: airdata_all
 129   real(kind=8)     :: tocc(1,1)
 130   integer          ::num_bufr(7),numbufr,ibufr
 131   character(20)    ::filename
 132   integer(i_kind), allocatable :: ptotal(:)
 133
 134 ! Set standard parameters
 135   real(r_kind)     :: POinT001 =   0.001_r_kind
 136   real(r_kind)     :: POinT01  =   0.01_r_kind
 137   real(r_kind)     :: TEN      =  10.0_r_kind
 138   real(r_kind)     :: R45      =  45.0_r_kind
 139   real(r_kind)     :: R60      =  60.0_r_kind
 140   real(r_kind)     :: R90      =  90.0_r_kind
 141   real(r_kind)     :: R180     = 180.0_r_kind
 142   real(r_kind)     :: R360     = 360.0_r_kind
 143
 144 ! Thinning variables
 145   integer(i_kind) itt,itx,iobs,iout,size_tmp, j
 146   real(r_kind) crit,dist
 147   real(r_kind) dlon_earth,dlat_earth
 148   logical luse
 149   real , allocatable :: in(:), out(:)
 150   logical :: found, head_found
 151
 152   logical           :: airs, eos_amsua, hsb, airstab
 153   type(info_type)            :: info
 154   type(model_loc_type)       :: loc
 155   type (datalink_type), pointer   :: head, p, current, prev
 156
 157   if (trace_use) call da_trace_entry("da_read_obs_bufrairs")
 158
 159 !  0.0  Initialize variables
 160 !-----------------------------------
 161   platform_id  = 9   ! eos series
 162   satellite_id = 2   ! eos-2
 163   nreal  = maxinfo
 164   allocate(ptotal(0:num_fgat_time))
 165   ptotal(0:num_fgat_time) = 0
 166   airs=      obstype == 'airs     '
 167   eos_amsua= obstype == 'eos_amsua'
 168   hsb=       obstype == 'hsb      '
 169
 170   icount=0
 171   if(airs)then
 172      sensor_id = 11
 173      step   = 1.1_r_kind
 174      start = -48.9_r_kind
 175      senname = 'AIRS'
 176      nchan  = N_AIRSCHAN
 177      nchanr = N_AIRSCHAN
 178   else if(eos_amsua)then
 179      sensor_id = 3
 180      step   = three + one/three
 181      start  = -48.33_r_kind
 182      senname = 'AMSU'
 183      nchan  = N_AMSUCHAN
 184      nchanr = N_AMSUCHAN
 185   else if(hsb)then
 186      sensor_id = 12
 187      step   = 1.1_r_kind
 188      start  = -48.95_r_kind
 189      senname = 'HSB'
 190      nchan  = N_HSBCHAN
 191      nchanr = N_HSBCHAN+1
 192   end if
 193   spotname = trim(senname)//'SPOT'
 194   channame = trim(senname)//'CHAN'
 195
 196       do inst = 1, rtminit_nsensor
 197         if (    platform_id  == rtminit_platform(inst) &
 198           .and. satellite_id == rtminit_satid(inst)    &
 199           .and. sensor_id    == rtminit_sensor(inst)    ) then
 200             exit
 201         end if
 202       end do
 203
 204       if ( inst == rtminit_nsensor .and.           &
 205            platform_id  /= rtminit_platform(inst)  &
 206           .or. satellite_id /= rtminit_satid(inst) &
 207           .or. sensor_id /= rtminit_sensor(inst)  ) then
 208          if (trace_use) call da_trace_exit("da_read_obs_bufrairs")
 209          return
 210       end if
 211
 212   num_eos_file    = 0  ! number of obs in file
 213   num_eos_global  = 0  ! number of obs within whole domain
 214   num_eos_local   = 0  ! number of obs within tile
 215   num_eos_kept    = 0  ! number of obs kept for innovation computation
 216   num_eos_thinned = 0  ! number of obs rejected by thinning
 217   size = 0     !
 218   iobs = 0     ! for thinning, argument is inout
 219
 220 !    1.0  Open BUFR table and BUFR file
 221 !--------------------------------------------------------------
 222   table_file = 'gmao_airs_bufr.tbl'      ! make table file name
 223   inquire(file=table_file,exist=airstab)
 224   if (airstab) then
 225       if (print_detail_rad) then
 226          write(unit=message(1),fmt=*) &
 227             'Reading BUFR Table A file: ',trim(table_file)
 228          call da_message(message(1:1))
 229       end if
 230       call da_get_unit(airs_table_unit)
 231       open(unit=airs_table_unit,file=table_file,iostat = iost)
 232       if (iost /= 0) then
 233          call da_error(__FILE__,__LINE__, &
 234             (/"Cannot open file "//table_file/))
 235       end if
 236   end if
 237
 238 ! Open BUFR file
 239 !  call da_get_unit(lnbufr)
 240
 241    num_bufr(:)=0
 242    numbufr=0
 243    if (num_fgat_time>1) then
 244       do i=1,7
 245          call da_get_unit(lnbufr)
 246          write(filename,fmt='(A,2I1,A)') trim(infile),0,i,'.bufr'
 247          open(unit   = lnbufr, FILE   = trim(filename),iostat =  iost, form = 'unformatted', STATUS = 'OLD')
 248          if (iost == 0) then
 249             numbufr=numbufr+1
 250             num_bufr(numbufr)=i
 251          else
 252             close (lnbufr)
 253          end if
 254          call da_free_unit(lnbufr)
 255       end do
 256    else
 257      numbufr=1
 258    end if
 259
 260    if (numbufr==0) numbufr=1
 261
 262 bufrfile:  do ibufr=1,numbufr
 263    if (num_fgat_time==1) then
 264       filename=trim(infile)//'.bufr'
 265    else
 266       if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then
 267          filename=trim(infile)//'.bufr'
 268       else
 269          write(filename,fmt='(A,2I1,A)') trim(infile),0,num_bufr(ibufr),'.bufr'
 270       end if
 271    end if
 272
 273   lnbufr=97
 274
 275   open(unit=lnbufr,file=trim(filename),form='unformatted',iostat = iost, status = 'old')
 276   if (iost /= 0) then
 277      call da_warning(__FILE__,__LINE__, &
 278         (/"Cannot open file "//infile/))
 279      close(lnbufr)
 280      if (airstab) then
 281         close(airs_table_unit)
 282         call da_free_unit(airs_table_unit)
 283      end if
 284      if (trace_use) call da_trace_exit("da_read_obs_bufrairs")
 285      return
 286   end if
 287   if ( airstab ) then
 288      call openbf(lnbufr,'IN',airs_table_unit)
 289   else
 290      call openbf(lnbufr,'IN',lnbufr)
 291   end if
 292   call datelen(10)
 293
 294
 295 !   2.0  Read header
 296 !---------------------------
 297   call readmg(lnbufr,subset,idate,iret)
 298
 299   iy = 0
 300   im = 0
 301   idd = 0
 302   ihh = 0
 303   if( iret /= 0 ) goto 1000     ! no data?
 304
 305   write(unit=date,fmt='( i10)') idate
 306   read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
 307   write(unit=stdout,fmt='(a,4i4,2x,a)') &
 308       'Bufr file date is ',iy,im,idd,ihh,trim(infile)
 309
 310 !   3.0 Loop over observations
 311 !----------------------------
 312      if ( ibufr == 1 ) then
 313         allocate ( head )
 314       ! allocate ( head % tb (1:nchan) )
 315         nullify  ( head % next )
 316         p => head
 317      endif
 318
 319   loop_obspoints: do
 320
 321 !   3.1 Read headder
 322 !-------------------------------
 323      call readsb(lnbufr,iret)
 324
 325      if( iret /=0 )then
 326         call readmg(lnbufr,subset,idate,iret)
 327         if( iret /= 0 ) exit loop_obspoints     ! end of file
 328         cycle loop_obspoints
 329      end if
 330
 331      num_eos_file = num_eos_file + 1
 332
 333 !   3.2 Read AQUASPOT (SPITSEQN)
 334 !------------------------
 335      call ufbseq(lnbufr,aquaspot_list_array,N_AQUASPOT_LIST,1,iret,'SPITSEQN')
 336      aquaspot = aquaspot_list( aquaspot_list_array(1), &
 337                                aquaspot_list_array(2), &
 338                                aquaspot_list_array(3), &
 339                                aquaspot_list_array(4), &
 340                                aquaspot_list_array(5), &
 341                                aquaspot_list_array(6), &
 342                                aquaspot_list_array(7), &
 343                                RESHAPE(aquaspot_list_array(8:25), (/2,9/)) )
 344
 345 !   3.3 Read AIRSSPOT or AMSUSPOT or HSBSPOT
 346 !-------------------------------------------------
 347      if ( trim(senname) == 'AIRS' ) then
 348         call ufbseq(lnbufr,airsspot_list_array,N_AIRSSPOT_LIST,1,iret,'SITPSEQN')
 349      else
 350         call ufbseq(lnbufr,airsspot_list_array,N_AIRSSPOT_LIST,1,iret,spotname)
 351      end if
 352      airsspot = airsspot_list( airsspot_list_array(1), &
 353                                airsspot_list_array(2), &
 354                                airsspot_list_array(3), &
 355                                airsspot_list_array(4), &
 356                                airsspot_list_array(5), &
 357                                airsspot_list_array(6), &
 358                                airsspot_list_array(7), &
 359                                airsspot_list_array(8), &
 360                                airsspot_list_array(9), &
 361                                airsspot_list_array(10), &
 362                                airsspot_list_array(11), &
 363                                airsspot_list_array(12) )
 364
 365 !   3.4 Read AIRSCHAN or AMSUCHAN or HSBCHAN
 366 !-------------------------------------------
 367      if ( trim(senname) == 'AIRS' ) then
 368         call ufbseq(lnbufr,airschan_list_array,N_AIRSCHAN_LIST,N_MAXCHAN,iret,'SCBTSEQN')
 369      else
 370         call ufbseq(lnbufr,airschan_list_array,N_AIRSCHAN_LIST,N_MAXCHAN,iret,channame)
 371      end if
 372      do l = 1 , N_MAXCHAN
 373         airschan(l) = airschan_list( airschan_list_array(1,l), &
 374                                      airschan_list_array(2,l), &
 375                                      airschan_list_array(3,l), &
 376                                      airschan_list_array(4,l)  )
 377      end do
 378
 379      if (iret /= nchanr) then
 380         write(unit=message(1),fmt=*) &
 381             'Cannot read ', channame, &
 382             ' bufr data:', &
 383             iret, ' ch data is read instead of', nchanr
 384         call da_warning(__FILE__,__LINE__,message(1:1))
 385         cycle loop_obspoints
 386      end if
 387
 388 !   3.5 Read Cloud Cover from AIRS/VISNIR
 389 !-------------------------------------------
 390      call ufbint(lnbufr,tocc,1,1,iret,'TOCC')
 391
 392 !   4.0  Check observing position (lat/lon)
 393 !      QC1:  juge if data is in the domain,
 394 !            read next record if not
 395 !------------------------------------------
 396      if( abs(airsspot%clath) > R90  .or. &
 397           abs(airsspot%clonh) > R360 .or. &
 398           (abs(airsspot%clath) == R90 .and. airsspot%clonh /= ZERO) )then
 399         cycle loop_obspoints
 400      end if
 401
 402 !    Retrieve observing position
 403      if(airsspot%clonh >= R360) then
 404         airsspot%clonh = airsspot%clonh - R360
 405 !     else if(airsspot%clonh < ZERO) then
 406 !        airsspot%clonh = airsspot%clonh + R360
 407      end if
 408
 409         info%lat  = airsspot%clath
 410         info%lon  = airsspot%clonh
 411         call da_llxy (info, loc, outside, outside_all )
 412
 413         if ( outside_all ) cycle loop_obspoints
 414
 415 !  4.1  Check obs time
 416 !-------------------------------------
 417      idate5(1) = airsspot%year ! year
 418      idate5(2) = airsspot%mnth ! month
 419      idate5(3) = airsspot%days ! day
 420      idate5(4) = airsspot%hour ! hour
 421      idate5(5) = airsspot%minu ! minute
 422      idate5(6) = airsspot%seco ! second
 423
 424      if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. &
 425           idate5(2) <    1 .or. idate5(2) >   12 .or. &
 426           idate5(3) <    1 .or. idate5(3) >   31 .or. &
 427           idate5(4) <    0 .or. idate5(4) >   24 .or. &
 428           idate5(5) <    0 .or. idate5(5) >   60 .or. &
 429           idate5(6) <    0 .or. idate5(6) >   60 ) then
 430         cycle loop_obspoints
 431      end if
 432
 433 !  QC3: time consistency check with the background date
 434
 435       if (idate5(1) .LE. 99) then
 436         if (idate5(1) .LT. 78) then
 437           idate5(1) = idate5(1) + 2000
 438         else
 439           idate5(1) = idate5(1) + 1900
 440         end if
 441       end if
 442
 443       call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
 444
 445       if ( obs_time < time_slots(0) .or.  &
 446            obs_time >= time_slots(num_fgat_time) ) cycle
 447
 448 !  3.2.1   determine FGAT index ifgat
 449 !
 450        do ifgat=1,num_fgat_time
 451            if ( obs_time >= time_slots(ifgat-1) .and.  &
 452                 obs_time  < time_slots(ifgat) ) exit
 453        end do
 454
 455        num_eos_global = num_eos_global + 1
 456        ptotal(ifgat) = ptotal(ifgat) + 1
 457
 458        if (outside) cycle loop_obspoints
 459
 460        num_eos_local = num_eos_local + 1
 461
 462        write(unit=info%date_char, &
 463          fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 464          idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
 465          ':', idate5(5), ':', idate5(6)
 466
 467        info%elv = 0.0  !aquaspot%selv
 468
 469 !  4.2  Check observational info
 470 !-------------------------------------------------------
 471      if( airsspot%fovn <    0.0_r_kind .or. airsspot%fovn > 100.0_r_kind .or. &
 472           airsspot%saza < -360.0_r_kind .or. airsspot%saza > 360.0_r_kind .or. &
 473           aquaspot%soza < -180.0_r_kind .or. aquaspot%soza > 180.0_r_kind )then
 474          write(unit=message(1),fmt=*) &
 475             'Cannot read ', channame, ' bufr data:', &
 476             ' strange obs info(fov,saza,soza):', &
 477             airsspot%fovn, airsspot%saza, aquaspot%soza
 478         call da_warning(__FILE__,__LINE__,message(1:1))
 479         cycle loop_obspoints
 480      end if
 481
 482 !    Retrieve observing info
 483      ifov = int( airsspot%fovn + POinT001 )
 484      sza  = abs(airsspot%saza)
 485 !     if( ((airs .or. hsb) .and. ifov <= 45) .or. &
 486 !          ( eos_amsua     .and. ifov <= 15) )then
 487 !        sza = - sza
 488 !     end if
 489
 490 !     airdata(6) = (start + float(ifov-1)*step)  ! look angle (deg)
 491 !     airdata(9) = ZERO                          ! surface height
 492 !     airdata(10)= POinT001                      ! land sea mask
 493
 494 !  4.1 Make Thinning
 495 !  Map obs to thinning grid
 496 !-------------------------------------------------------------------
 497        if (thinning) then
 498           dlat_earth = info%lat
 499           dlon_earth = info%lon
 500           if (dlon_earth<zero) dlon_earth = dlon_earth+r360
 501           if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 502           dlat_earth = dlat_earth*deg2rad
 503           dlon_earth = dlon_earth*deg2rad
 504           crit = 1.0 ! 0.01_r_kind+terrain + timedif !+ 10.0_r_kind*float(iskip)
 505           if (airs_warmest_fov) &
 506              crit = 1E10 * exp(-(airschan(129)%tmbrst-220.0)/2)   ! warmest bt for window channel (10.36 micron)
 507           call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,luse)
 508           if (.not. luse) then
 509              num_eos_thinned = num_eos_thinned + 1
 510              cycle loop_obspoints
 511           end if
 512        end if
 513
 514
 515 !   4.3 Retrieve Tb
 516 !-----------------------
 517      iflag = .false.
 518
 519      do l=1,nchan
 520         airdata(l+nreal) = airschan(l)%tmbrst            ! brightness temperature
 521         if( airdata(l+nreal) > 0.0_r_kind .and. airdata(l+nreal) < 500.0_r_kind )then
 522            iflag = .true.
 523         else
 524            airdata(l+nreal) = missing_r
 525         end if
 526      end do
 527
 528      if ( .not. iflag )then
 529         write(unit=message(1),fmt=*) &
 530           'Error in reading ', channame, ' bufr data:', &
 531           ' all tb data is missing'
 532         call da_warning(__FILE__,__LINE__,message(1:1))
 533         num_eos_local = num_eos_local - 1
 534         cycle loop_obspoints
 535      end if
 536
 537      num_eos_kept = num_eos_kept + 1
 538
 539 !  4.0   assign information to sequential radiance structure
 540 !--------------------------------------------------------------------------
 541    allocate ( p % tb_inv (1:nchan) )
 542    p%info             = info
 543    p%loc              = loc
 544    p%landsea_mask     = POinT001
 545    p%scanline         = int(aquaspot%slnm + POinT001)
 546    p%scanpos          = ifov
 547    p%satzen           = sza
 548    p%satazi           = (start + float(ifov-1)*step)  ! look angle (deg) ! airsspot%bearaz
 549    p%solzen           = aquaspot%soza
 550    p%solazi           = aquaspot%solazi
 551    p%tb_inv(1:nchan)  = airdata(nreal+1:nreal+nchan)
 552    p%sensor_index     = inst
 553    p%ifgat            = ifgat
 554
 555    size = size + 1
 556    allocate ( p%next, stat=error)   ! add next data
 557    if (error /= 0 ) then
 558       call da_error(__FILE__,__LINE__, &
 559           (/"Cannot allocate radiance structure"/))
 560    end if
 561
 562    p => p%next
 563    nullify (p%next)
 564
 565   end do loop_obspoints
 566
 567   call closbf(lnbufr)
 568   close(lnbufr)
 569
 570 end do bufrfile
 571
 572    if (thinning .and. num_eos_global > 0) then
 573
 574 #ifdef DM_PARALLEL
 575
 576       ! Get minimum crit and associated processor index.
 577       j = 0
 578       do ifgat = 1, num_fgat_time
 579          do n = 1, iv%num_inst
 580             j = j + thinning_grid(n,ifgat)%itxmax
 581          end do
 582       end do
 583
 584       allocate ( in  (j) )
 585       allocate ( out (j) )
 586
 587       j = 0
 588       do ifgat = 1, num_fgat_time
 589          do n = 1, iv%num_inst
 590             do i = 1, thinning_grid(n,ifgat)%itxmax
 591                j = j + 1
 592                in(j) = thinning_grid(n,ifgat)%score_crit(i)
 593             end do
 594          end do
 595       end do
 596       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 597
 598       call wrf_dm_bcast_real (out, j)
 599
 600       j = 0
 601       do ifgat = 1, num_fgat_time
 602          do n = 1, iv%num_inst
 603             do i = 1, thinning_grid(n,ifgat)%itxmax
 604                j = j + 1
 605                if ( ABS(out(j)-thinning_grid(n,ifgat)%score_crit(i)) > 1.0E-10 ) thinning_grid(n,ifgat)%ibest_obs(i)  = 0
 606             end do
 607          end do
 608       end do
 609
 610       deallocate( in  )
 611       deallocate( out )
 612
 613 #endif
 614
 615       ! Delete the nodes which being thinning out
 616       if ( size > 0 ) then
 617          p => head
 618          prev => head
 619          head_found = .false.
 620          size_tmp = size
 621          do j = 1, size_tmp
 622             n = p%sensor_index
 623             ifgat = p%ifgat
 624             found = .false.
 625
 626             do i = 1, thinning_grid(n,ifgat)%itxmax
 627                if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 628                   found = .true.
 629                   exit
 630                endif
 631             end do
 632
 633             ! free current data
 634             if ( .not. found ) then
 635                current => p
 636                p => p%next
 637                if ( head_found ) then
 638                   prev%next => p
 639                else
 640                   head => p
 641                   prev => p
 642                endif
 643                deallocate ( current % tb_inv )
 644                deallocate ( current )
 645                num_eos_thinned = num_eos_thinned + 1
 646                num_eos_kept = num_eos_kept - 1
 647                size = size - 1
 648                continue
 649             endif
 650
 651             if ( found .and. head_found ) then
 652                prev => p
 653                p => p%next
 654                continue
 655             endif
 656
 657             if ( found .and. .not. head_found ) then
 658                head_found = .true.
 659                head => p
 660                prev => p
 661                p => p%next
 662             endif
 663          end do
 664       endif
 665
 666    endif  ! End of thinning
 667
 668    iv%total_rad_pixel   = iv%total_rad_pixel + size
 669    iv%total_rad_channel = iv%total_rad_channel + size*nchan
 670
 671    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + size
 672    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_eos_global
 673    iv%instid(inst)%info%nlocal = size
 674
 675    do ifgat = 1, num_fgat_time
 676       ptotal(ifgat) = ptotal(ifgat) + ptotal(ifgat-1)
 677       iv%info(radiance)%ptotal(ifgat) = iv%info(radiance)%ptotal(ifgat) + ptotal(ifgat)
 678    end do
 679
 680    write(unit=stdout,fmt='(a)') '   num_eos_file num_eos_global  num_eos_local   num_eos_kept num_eos_thinned'
 681    write(unit=stdout,fmt='(5(5x,i10))') num_eos_file,num_eos_global, num_eos_local,num_eos_kept,num_eos_thinned
 682
 683 !  5.0 allocate innovation radiance structure
 684 !----------------------------------------------------------------
 685 !  do i = 1, iv%num_inst
 686    if ( size > 0 ) then
 687       iv%instid(inst)%num_rad = size
 688       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 689         'Allocating space for radiance innov structure', &
 690          inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
 691    end if
 692 !  end do
 693
 694 !   6.0 assign sequential structure to innovation structure
 695 !-------------------------------------------------------------
 696
 697   n = 0
 698   p => head
 699
 700   call da_allocate_rad_iv (inst, nchan, iv)
 701
 702   do i = 1, size
 703 !   inst = p%sensor_index
 704    n = n + 1
 705
 706    call da_initialize_rad_iv (inst, n, iv, p)
 707    iv%instid(inst)%rain_flag(n) = tocc(1,1) ! Temporary dumping of AIRS/VISNIR cloud cover
 708
 709    current => p
 710    p => p%next
 711
 712 ! free current data
 713    deallocate ( current % tb_inv )
 714    deallocate ( current )
 715
 716  end do
 717
 718  deallocate (p)
 719  deallocate (ptotal)
 720
 721 1000 continue
 722    call closbf(lnbufr)
 723    close(lnbufr)
 724 !   call da_free_unit(lnbufr)
 725    if (airstab) then
 726       close(airs_table_unit)
 727       call da_free_unit(airs_table_unit)
 728    end if
 729
 730    if (trace_use) call da_trace_exit("da_read_obs_bufrairs")
 731 #else
 732    call da_error(__FILE__,__LINE__,(/"Needs to be compiled with a BUFR library"/))
 733 #endif
 734
 735 end subroutine da_read_obs_bufrairs