var/da/da_radiance/da_read_obs_bufratms.inc

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