var/da/da_radiance/da_read_obs_hdf5ahi.inc

   1 subroutine da_read_obs_hdf5ahi (iv,infile_tb,infile_clp)
   2    !--------------------------------------------------------
   3    !  Purpose: read in CMA AHI Level-1B and Level-2 data in HDF5 format
   4    !           and form innovation structure
   5    !
   6    !   METHOD: use F90 sequantial data structure to avoid read the file twice
   7    !            1. read file radiance data in sequential data structure
   8    !            2. do gross QC check
   9    !            3. assign sequential data structure to innovation structure
  10    !                  and deallocate sequential data structure
  11    !
  12    !  HISTORY: 2016/10/12 - Creation         Yuanbing Wang, NUIST/CAS, NCAR/NESL/MMM/DAS
  13    !
  14    !  To be developed:
  15    !            1. more general variable names;
  16    !            2. get time and dimension information from file
  17    !            3. more readable and efficient programm
  18    !------------------------------------------------------------------------------
  19
  20    implicit none
  21
  22    character(len=*), intent(in)    :: infile_tb, infile_clp
  23    type(iv_type),    intent(inout) :: iv
  24
  25 #if defined(HDF5)
  26 ! fixed parameter values
  27    integer,parameter::nlatitude=600          ! Maximum allowed NumberOfScans
  28    integer,parameter::nlongitude=700          ! low  resolution pixel width
  29    integer,parameter::time_dims=6       ! Time dimension
  30    integer,parameter::nfile_max = 8     ! each hdf file contains ~50min of data
  31                                         ! at most 8 files for a 6-h time window
  32 ! interface variable
  33    integer iret                         ! return status
  34    integer(HID_T) fhnd1       ! file handle
  35    integer(HID_T) dhnd1         ! dataset handle
  36    integer(HSIZE_T) sz1(2)              ! array size 1
  37
  38 ! array data
  39    real(4) :: vlatitude(nlongitude,nlatitude)  !  value for latitude
  40    real(4) :: vlongitude(nlongitude,nlatitude) !  value for longitude
  41
  42    real(4) :: tb07(nlongitude,nlatitude)  ! tb for band 7
  43    real(4) :: tb08(nlongitude,nlatitude)  ! tb for band 8
  44    real(4) :: tb09(nlongitude,nlatitude)  ! tb for band 9
  45    real(4) :: tb10(nlongitude,nlatitude)  ! tb for band 10
  46    real(4) :: tb11(nlongitude,nlatitude)  ! tb for band 11
  47    real(4) :: tb12(nlongitude,nlatitude)  ! tb for band 12
  48    real(4) :: tb13(nlongitude,nlatitude)  ! tb for band 13
  49    real(4) :: tb14(nlongitude,nlatitude)  ! tb for band 14
  50    real(4) :: tb15(nlongitude,nlatitude)  ! tb for band 15
  51    real(4) :: tb16(nlongitude,nlatitude)  ! tb for band 16
  52
  53    real(4) :: sat_zenith(nlongitude,nlatitude)  ! satellite_zenith_angle
  54    integer(4) :: cloud_mask(nlongitude,nlatitude)     !obs cloud mask
  55
  56    real(r_kind),parameter  :: tbmin  = 50._r_kind
  57    real(r_kind),parameter  :: tbmax  = 550._r_kind
  58    real(r_kind),parameter  :: tb_scale  = 100._r_kind
  59
  60    real(kind=8)                   :: obs_time
  61    type (datalink_type),pointer   :: head, p, current, prev
  62    type(info_type)                :: info
  63    type(model_loc_type)           :: loc
  64
  65    integer(i_kind)   :: idate5(6)
  66    integer(i_kind)   :: inst,platform_id,satellite_id,sensor_id
  67    real(r_kind)      :: tb, crit
  68    integer(i_kind)   :: ifgat, iout, iobs
  69    logical           :: outside, outside_all, iuse
  70
  71    integer           :: i,j,k,l,m,n, ifile,landsea_mask
  72    logical           :: found, head_found, head_allocated
  73
  74 ! Other work variables
  75    real(r_kind)     :: dlon_earth,dlat_earth
  76    integer(i_kind)  :: num_ahi_local, num_ahi_global, num_ahi_used, num_ahi_thinned
  77    integer(i_kind)  :: num_ahi_used_tmp, num_ahi_file
  78    integer(i_kind)  :: num_ahi_local_local, num_ahi_global_local, num_ahi_file_local
  79    integer(i_kind)  :: itx, itt
  80    character(80)    :: filename1, filename2
  81    integer          :: nchan,ilongitude,ilatitude,ichannels
  82    integer          :: nfile, ahi_info_unit
  83    character(80)    :: fname_tb(nfile_max)
  84    character(80)    :: fname_clp(nfile_max)
  85    logical          :: fexist, got_clp_file
  86
  87 ! Allocatable arrays
  88    integer(i_kind),allocatable  :: ptotal(:)
  89    real,allocatable             :: in(:), out(:)
  90    real(r_kind),allocatable     :: data_all(:)
  91
  92    if (trace_use) call da_trace_entry("da_read_obs_hdf5ahi")
  93
  94 !  0.0  Initialize variables
  95 !-----------------------------------
  96    head_allocated = .false.
  97    platform_id  = 31  ! Table-2 Col 1 corresponding to 'himawari'
  98    satellite_id = 8   ! Table-2 Col 3
  99    sensor_id    = 56  ! Table-3 Col 2 corresponding to 'ahi'
 100
 101    allocate(ptotal(0:num_fgat_time))
 102    ptotal(0:num_fgat_time) = 0
 103    iobs = 0                 ! for thinning, argument is inout
 104    num_ahi_file    = 0
 105    num_ahi_local   = 0
 106    num_ahi_global  = 0
 107    num_ahi_used    = 0
 108    num_ahi_thinned = 0
 109
 110    sz1(1)=nlongitude
 111    sz1(2)=nlatitude
 112
 113    do i = 1, rtminit_nsensor
 114       if (platform_id  == rtminit_platform(i) &
 115           .and. satellite_id == rtminit_satid(i)    &
 116           .and. sensor_id    == rtminit_sensor(i)) then
 117          inst = i
 118          exit
 119       end if
 120    end do
 121    if (inst == 0) then
 122       call da_warning(__FILE__,__LINE__, &
 123           (/"The combination of Satellite_Id and Sensor_Id for AHI is not found"/))
 124       if (trace_use) call da_trace_exit("da_read_obs_hdf5ahi")
 125       return
 126    end if
 127
 128 ! Initialize HDF5 library and Fortran90 interface
 129    call H5open_f(iret)
 130    if(iret.lt.0)then
 131       call da_warning(__FILE__,__LINE__,(/"Problems initializing HDF5 Lib, can't read AHI data."/))
 132       if (trace_use) call da_trace_exit("da_read_obs_hdf5ahi")
 133       return
 134    endif
 135
 136    nchan = iv%instid(inst)%nchan
 137    write(unit=stdout,fmt=*)'AHI nchan: ',nchan
 138    allocate(data_all(1:nchan))
 139
 140 ! 1.0 Assign file names and prepare to read ahi files
 141 !-------------------------------------------------------------------------
 142    nfile       = 0  !initialize
 143    fname_tb(:) = '' !initialize
 144    ! first check if hdf file is available
 145    filename1 = trim(infile_tb)
 146    filename2 = trim(infile_clp)
 147    inquire (file=filename1, exist=fexist)
 148    if ( fexist ) then
 149       nfile = 1
 150       fname_tb(nfile)  = filename1
 151       fname_clp(nfile) = filename2
 152    else
 153       ! check if L1SGRTBR-0x.h5 is available for multiple input files
 154       ! here 0x is the input file sequence number
 155       ! do not confuse it with fgat time slot index
 156       do i = 1, nfile_max
 157          write(filename1,fmt='(A,A,I2.2,A)') trim(infile_tb),'-',i
 158          write(filename2,fmt='(A,A,I2.2,A)') trim(infile_clp),'-',i
 159          inquire (file=filename1, exist=fexist)
 160          if ( fexist ) then
 161             nfile = nfile + 1
 162             fname_tb(nfile)  = filename1
 163             fname_clp(nfile) = filename2
 164          else
 165             exit
 166          end if
 167       end do
 168    end if
 169
 170    if ( nfile == 0 ) then
 171       call da_warning(__FILE__,__LINE__,(/"No valid AHI file found."/))
 172       if (trace_use) call da_trace_exit("da_read_obs_hdf5ahi")
 173       return
 174    end if
 175
 176   !open the data info file
 177    call da_get_unit(ahi_info_unit)
 178    open(unit=ahi_info_unit,file='ahi_info',status='old',iostat=iret)
 179    if(iret /= 0)then
 180       call da_warning(__FILE__,__LINE__,(/"data_info file read error"/))
 181    endif
 182    read(ahi_info_unit,*)
 183    read(ahi_info_unit,*)
 184    read(ahi_info_unit,*)
 185    read(ahi_info_unit,*)
 186    read(ahi_info_unit,*)
 187    read(ahi_info_unit,*)
 188    read(ahi_info_unit,*)
 189
 190    infile_loop:  do ifile = 1, nfile
 191       num_ahi_file_local    = 0
 192       num_ahi_local_local   = 0
 193       num_ahi_global_local  = 0
 194
 195    ! open infile_tb HDF5 file for read
 196       call H5Fopen_f(fname_tb(ifile),H5F_ACC_RDONLY_F,fhnd1,iret,H5P_DEFAULT_F)
 197       if(iret.lt.0)then
 198          call da_warning(__FILE__,__LINE__,(/"Cannot open HDF5 file"//trim(fname_tb(ifile))/))
 199          cycle infile_loop
 200       endif
 201
 202    ! read lat
 203       call H5Dopen_f(fhnd1,'pixel_latitude',dhnd1,iret)
 204       call H5Dread_f(dhnd1,H5T_IEEE_F32LE,vlatitude,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 205       if(iret.lt.0)then
 206          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Latitude"/))
 207       endif
 208       call H5Dclose_f(dhnd1,iret)
 209
 210    ! read lon
 211       call H5Dopen_f(fhnd1,'pixel_longitude',dhnd1,iret)
 212       call H5Dread_f(dhnd1,H5T_IEEE_F32LE,vlongitude,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 213       if(iret.lt.0)then
 214           call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Longitude"/))
 215           call da_trace_exit("da_read_obs_hdf5ahi")
 216       endif
 217       call H5Dclose_f(dhnd1,iret)
 218    ! sample display
 219       write(unit=stdout,fmt=*)'latitude,longitude(pixel=1,scan=1): ',vlatitude(1,1),vlongitude(1,1)
 220
 221    ! read tb for band 7
 222       call H5Dopen_f(fhnd1,'NOMChannelIRX0390_2000',dhnd1,iret)
 223       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb07,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 224       if(iret.lt.0)then
 225          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Band 7"/))
 226       endif
 227       call H5Dclose_f(dhnd1,iret)
 228    ! sample display
 229       write(unit=stdout,fmt=*) 'tb07(pixel=1,scan=1): ',tb07(1,1)
 230
 231    ! read tb for band 8
 232       call H5Dopen_f(fhnd1,'NOMChannelIRX0620_2000',dhnd1,iret)
 233       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb08,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 234       if(iret.lt.0)then
 235          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Band 8"/))
 236       endif
 237       call H5Dclose_f(dhnd1,iret)
 238    ! sample display
 239       write(unit=stdout,fmt=*) 'tb08(pixel=1,scan=1): ',tb08(1,1)
 240
 241    ! read tb for band 9
 242       call H5Dopen_f(fhnd1,'NOMChannelIRX0700_2000',dhnd1,iret)
 243       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb09,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 244       if(iret.lt.0)then
 245          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Band 9"/))
 246       endif
 247       call H5Dclose_f(dhnd1,iret)
 248    ! sample display
 249       write(unit=stdout,fmt=*) 'tb09(pixel=1,scan=1): ',tb09(1,1)
 250
 251    ! read tb for band 10
 252       call H5Dopen_f(fhnd1,'NOMChannelIRX0730_2000',dhnd1,iret)
 253       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb10,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 254       if(iret.lt.0)then
 255          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Band 10"/))
 256       endif
 257       call H5Dclose_f(dhnd1,iret)
 258    ! sample display
 259       write(unit=stdout,fmt=*) 'tb10(pixel=1,scan=1): ',tb10(1,1)
 260
 261    ! read tb for band 11
 262       call H5Dopen_f(fhnd1,'NOMChannelIRX0860_2000',dhnd1,iret)
 263       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb11,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 264       if(iret.lt.0)then
 265          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Band 11"/))
 266       endif
 267       call H5Dclose_f(dhnd1,iret)
 268    ! sample display
 269       write(unit=stdout,fmt=*) 'tb11(pixel=1,scan=1): ',tb11(1,1)
 270
 271    ! read tb for band 12
 272       call H5Dopen_f(fhnd1,'NOMChannelIRX0960_2000',dhnd1,iret)
 273       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb12,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 274       if(iret.lt.0)then
 275          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Band 12"/))
 276       endif
 277       call H5Dclose_f(dhnd1,iret)
 278    ! sample display
 279       write(unit=stdout,fmt=*) 'tb12(pixel=1,scan=1): ',tb12(1,1)
 280
 281    ! read tb for band 13
 282       call H5Dopen_f(fhnd1,'NOMChannelIRX1040_2000',dhnd1,iret)
 283       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb13,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 284       if(iret.lt.0)then
 285          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: 13"/))
 286       endif
 287       call H5Dclose_f(dhnd1,iret)
 288    ! sample display
 289       write(unit=stdout,fmt=*) 'tb13(pixel=1,scan=1): ',tb13(1,1)
 290
 291    ! read tb for band 14
 292       call H5Dopen_f(fhnd1, 'NOMChannelIRX1120_2000',dhnd1,iret)
 293       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb14,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 294       if(iret.lt.0)then
 295          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Band 14"/))
 296       endif
 297       call H5Dclose_f(dhnd1,iret)
 298    ! sample display
 299       write(unit=stdout,fmt=*) 'tb14(pixel=1,scan=1): ',tb14(1,1)
 300
 301    ! read tb for band 15
 302       call H5Dopen_f(fhnd1,'NOMChannelIRX1230_2000',dhnd1,iret)
 303       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb15,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 304       if(iret.lt.0)then
 305          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: 15"/))
 306       endif
 307       call H5Dclose_f(dhnd1,iret)
 308    ! sample display
 309       write(unit=stdout,fmt=*) 'tb15(pixel=1,scan=1): ',tb15(1,1)
 310
 311    ! read tb for band 16
 312       call H5Dopen_f(fhnd1,'NOMChannelIRX1330_2000',dhnd1,iret)
 313       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,tb16,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 314       if(iret.lt.0)then
 315          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: 16"/))
 316       endif
 317       call H5Dclose_f(dhnd1,iret)
 318    ! sample display
 319       write(unit=stdout,fmt=*) 'tb16(pixel=1,scan=1): ',tb16(1,1)
 320
 321    ! read array: satellite_zenith_angle
 322    ! read
 323       call H5Dopen_f(fhnd1,'pixel_satellite_zenith_angle',dhnd1,iret)
 324       call H5Dread_f(dhnd1,H5T_IEEE_F32LE,sat_zenith,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 325       if(iret.lt.0)then
 326          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: satellite_zenith_angle"/))
 327       endif
 328       call H5Dclose_f(dhnd1,iret)
 329    ! sample display
 330       write(unit=stdout,fmt=*)'sat_zenith(pixel=1,scan=1): ',sat_zenith(1,1)
 331
 332    ! close infile_tb and HDF5
 333       call H5Fclose_f(fhnd1,iret)
 334
 335      !open infile_clw file and HDF5
 336       got_clp_file = .false.
 337       call H5Fopen_f(fname_clp(ifile),H5F_ACC_RDONLY_F,fhnd1,iret,H5P_DEFAULT_F)
 338       if ( iret == 0 ) then
 339          got_clp_file = .true.
 340       endif
 341       ! to do: when got_clp_file=.true., need to check GranuleID for consistency
 342       ! betweee tb and clw files
 343
 344       if ( got_clp_file ) then
 345
 346       ! read CLOUD_MASK from infile_clw:
 347          call H5Dopen_f(fhnd1,'cloud_mask',dhnd1,iret)
 348          call H5Dread_f(dhnd1,H5T_NATIVE_INTEGER,cloud_mask,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 349          if(iret.lt.0)then
 350             call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: CLOUD_MASK data"/))
 351          endif
 352          call H5Dclose_f(dhnd1,iret)
 353       ! sample display
 354          write(unit=stdout,fmt=*)'cloud_mask(pixel=1,scan=1): ',cloud_mask(1,1)
 355
 356       ! close infile_clw file and HDF5
 357                 call H5Fclose_f(fhnd1,iret)
 358       end if
 359
 360          !read date information
 361           read(ahi_info_unit,*) idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),idate5(6)
 362
 363 ! 2.0 Loop to read hdf file and assign information to a sequential structure
 364 !-------------------------------------------------------------------------
 365
 366    ! Allocate arrays to hold data
 367       if ( .not. head_allocated ) then
 368          allocate (head)
 369          nullify  ( head % next )
 370          p => head
 371          head_allocated = .true.
 372       end if
 373
 374    ! start scan_loop
 375       scan_loop:     do ilatitude=1, nlatitude
 376
 377          call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
 378
 379          if ( obs_time < time_slots(0) .or. obs_time >= time_slots(num_fgat_time) ) cycle scan_loop
 380          do ifgat=1,num_fgat_time
 381             if ( obs_time >= time_slots(ifgat-1) .and. obs_time  < time_slots(ifgat) ) exit
 382          end do
 383
 384       ! start fov_loop: longitude
 385          fov_loop:      do ilongitude=1, nlongitude
 386
 387                         if ( sat_zenith(ilongitude,ilatitude) > 65.0 ) cycle fov_loop
 388
 389             num_ahi_file       = num_ahi_file + 1
 390             num_ahi_file_local = num_ahi_file_local + 1
 391             info%lat  =  vlatitude(ilongitude,ilatitude)
 392             info%lon  =  vlongitude(ilongitude,ilatitude)
 393
 394             call da_llxy (info, loc, outside, outside_all)
 395             if (outside_all) cycle fov_loop
 396
 397             num_ahi_global       = num_ahi_global + 1
 398             num_ahi_global_local = num_ahi_global_local + 1
 399             ptotal(ifgat) = ptotal(ifgat) + 1
 400             if (outside) cycle fov_loop   ! No good for this PE
 401
 402             num_ahi_local       = num_ahi_local + 1
 403             num_ahi_local_local = num_ahi_local_local + 1
 404             write(unit=info%date_char, &
 405             fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 406                idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
 407                ':', idate5(5), ':', idate5(6)
 408             info%elv = 0.0
 409
 410 ! 3.0  Make Thinning
 411 ! Map obs to thinning grid
 412 !-------------------------------------------------------------------
 413             if (thinning) then
 414                dlat_earth = info%lat !degree
 415                dlon_earth = info%lon
 416                if (dlon_earth<zero)  dlon_earth = dlon_earth+r360
 417                if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 418                dlat_earth = dlat_earth*deg2rad !radian
 419                dlon_earth = dlon_earth*deg2rad
 420                crit = 1.
 421                call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 422                if (.not. iuse) then
 423                   num_ahi_thinned = num_ahi_thinned+1
 424                   cycle fov_loop
 425                end if
 426             end if
 427
 428             num_ahi_used = num_ahi_used + 1
 429             data_all = missing_r
 430
 431             tb = tb07(ilongitude,ilatitude) / tb_scale
 432             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 433             data_all(1)= tb
 434
 435             tb = tb08(ilongitude,ilatitude) / tb_scale
 436             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 437             data_all(2)= tb
 438
 439                         tb = tb09(ilongitude,ilatitude) / tb_scale
 440             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 441             data_all(3)= tb
 442
 443             tb = tb10(ilongitude,ilatitude) / tb_scale
 444             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 445             data_all(4)= tb
 446
 447             tb = tb11(ilongitude,ilatitude) / tb_scale
 448             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 449             data_all(5)= tb
 450
 451                         tb = tb12(ilongitude,ilatitude) / tb_scale
 452             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 453             data_all(6)= tb
 454
 455             tb = tb13(ilongitude,ilatitude) / tb_scale
 456             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 457             data_all(7)= tb
 458
 459             tb = tb14(ilongitude,ilatitude) / tb_scale
 460             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 461             data_all(8)= tb
 462
 463                         tb = tb15(ilongitude,ilatitude) / tb_scale
 464             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 465             data_all(9)= tb
 466
 467                         tb = tb16(ilongitude,ilatitude) / tb_scale
 468             if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 469             data_all(10)= tb
 470
 471 ! 4.0 assign information to sequential radiance structure
 472 !--------------------------------------------------------------------------
 473             allocate ( p % tb_inv (1:nchan ))
 474             p%info             = info
 475             p%loc              = loc
 476             p%landsea_mask     = 1
 477             p%scanpos          = ilongitude !nint(sat_zenith(ilongitude,ilatitude))+1.001_r_kind !
 478             p%satzen           = sat_zenith(ilongitude,ilatitude)
 479             p%satazi           = 0
 480             p%solzen           = 0
 481             p%solazi           = 0
 482                         p%cloudflag        = cloud_mask(ilongitude,ilatitude)
 483             p%tb_inv(1:nchan)  = data_all(1:nchan)
 484             p%sensor_index     = inst
 485             p%ifgat            = ifgat
 486
 487             allocate (p%next)   ! add next data
 488             p => p%next
 489             nullify (p%next)
 490          end do fov_loop
 491       end do scan_loop
 492
 493       write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_file    : ',num_ahi_file_local
 494       write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_global  : ',num_ahi_global_local
 495       write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_local   : ',num_ahi_local_local
 496    end do infile_loop
 497
 498    close(ahi_info_unit) !close date information file
 499    call da_free_unit(ahi_info_unit)
 500    call H5close_f(iret)
 501    deallocate(data_all) ! Deallocate data arrays
 502
 503    if (thinning .and. num_ahi_global > 0 ) then
 504 #ifdef DM_PARALLEL
 505    ! Get minimum crit and associated processor index.
 506       j = 0
 507       do ifgat = 1, num_fgat_time
 508          j = j + thinning_grid(inst,ifgat)%itxmax
 509       end do
 510
 511       allocate ( in  (j) )
 512       allocate ( out (j) )
 513       j = 0
 514       do ifgat = 1, num_fgat_time
 515          do i = 1, thinning_grid(inst,ifgat)%itxmax
 516             j = j + 1
 517             in(j) = thinning_grid(inst,ifgat)%score_crit(i)
 518          end do
 519       end do
 520       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 521
 522       call wrf_dm_bcast_real (out, j)
 523
 524       j = 0
 525       do ifgat = 1, num_fgat_time
 526          do i = 1, thinning_grid(inst,ifgat)%itxmax
 527             j = j + 1
 528             if ( ABS(out(j)-thinning_grid(inst,ifgat)%score_crit(i)) > 1.0E-10 ) &
 529             thinning_grid(inst,ifgat)%ibest_obs(i) = 0
 530          end do
 531       end do
 532
 533       deallocate( in  )
 534       deallocate( out )
 535
 536 #endif
 537
 538    ! Delete the nodes which being thinning out
 539       p => head
 540       prev => head
 541       head_found = .false.
 542       num_ahi_used_tmp = num_ahi_used
 543       do j = 1, num_ahi_used_tmp
 544          n = p%sensor_index
 545          ifgat = p%ifgat
 546          found = .false.
 547
 548          do i = 1, thinning_grid(n,ifgat)%itxmax
 549             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 550                found = .true.
 551                exit
 552             end if
 553          end do
 554
 555       ! free current data
 556          if ( .not. found ) then
 557             current => p
 558             p => p%next
 559             if ( head_found ) then
 560                prev%next => p
 561             else
 562                head => p
 563                prev => p
 564             end if
 565             deallocate ( current % tb_inv )
 566             deallocate ( current )
 567             num_ahi_thinned = num_ahi_thinned + 1
 568             num_ahi_used = num_ahi_used - 1
 569             continue
 570          end if
 571
 572          if ( found .and. head_found ) then
 573             prev => p
 574             p => p%next
 575             continue
 576          end if
 577
 578          if ( found .and. .not. head_found ) then
 579             head_found = .true.
 580             head => p
 581             prev => p
 582             p => p%next
 583          end if
 584
 585       end do
 586
 587    end if  ! End of thinning
 588
 589    iv%total_rad_pixel   = iv%total_rad_pixel + num_ahi_used
 590    iv%total_rad_channel = iv%total_rad_channel + num_ahi_used*nchan
 591
 592    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_ahi_used
 593    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_ahi_global
 594
 595    do i = 1, num_fgat_time
 596       ptotal(i) = ptotal(i) + ptotal(i-1)
 597       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 598    end do
 599    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 600       write(unit=message(1),fmt='(A,I10,A,I10)') &
 601           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 602       call da_warning(__FILE__,__LINE__,message(1:1))
 603    endif
 604
 605    write(unit=stdout,fmt='(a)') 'AHI data counts: '
 606    write(stdout,fmt='(a,i7)') ' In file: ',num_ahi_file
 607    write(stdout,fmt='(a,i7)') ' Global : ',num_ahi_global
 608    write(stdout,fmt='(a,i7)') ' Local  : ',num_ahi_local
 609    write(stdout,fmt='(a,i7)') ' Used   : ',num_ahi_used
 610    write(stdout,fmt='(a,i7)') ' Thinned: ',num_ahi_thinned
 611
 612
 613 !  5.0 allocate innovation radiance structure
 614 !----------------------------------------------------------------
 615
 616    if (num_ahi_used > 0) then
 617       iv%instid(inst)%num_rad  = num_ahi_used
 618       iv%instid(inst)%info%nlocal = num_ahi_used
 619       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 620          'Allocating space for radiance innov structure', &
 621          inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
 622       call da_allocate_rad_iv (inst, nchan, iv)
 623    end if
 624
 625 !  6.0 assign sequential structure to innovation structure
 626 !-------------------------------------------------------------
 627    p => head
 628
 629    do n = 1, num_ahi_used
 630       i = p%sensor_index
 631       call da_initialize_rad_iv (i, n, iv, p)
 632       current => p
 633       p => p%next
 634    ! free current data
 635       deallocate ( current % tb_inv )
 636       deallocate ( current )
 637    end do
 638    deallocate ( p )
 639    deallocate (ptotal)
 640
 641    if (trace_use) call da_trace_exit("da_read_obs_hdf5ahi")
 642 #else
 643    call da_error(__FILE__,__LINE__,(/"Needs to be compiled with HDF5 library"/))
 644 #endif
 645 end subroutine da_read_obs_hdf5ahi