var/da/da_radiance/da_read_obs_hdf5gmi.inc

   1 subroutine da_read_obs_hdf5gmi (iv, infile_tb, infile_clw)
   2    !--------------------------------------------------------
   3    !  Purpose: read in NASA/JAXA GPM GMI Level-1B 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: 2013/10/22 - Creation         Syed RH Rizvi, NCAR/NESL/MMM/DAS
  13    !           2017/02/17   Modification     Wenying He
  14    !           2021/11/10   Modification     Dongmei Xu
  15    ! METHODS,
  16    !Shen, et al.，2021: Assimilation of GPM Microwave Imager Radiance data with
  17    ! the WRF Hybrid 3DEnVar System for the Prediction of Typhoon Chan-hom (2015),
  18    ! Atmospheric Research. 251, 105422.
  19    !------------------------------------------------------------------------------
  20
  21    implicit none
  22
  23    character(len=*), intent(in)    :: infile_tb, infile_clw
  24    type(iv_type),    intent(inout) :: iv
  25
  26 #if defined(HDF5)
  27 ! fixed parameter values
  28    integer,parameter::max_scan=2963     ! Maximum allowed NumberOfScans
  29    integer,parameter::ovr=20            ! Number of OverlapScans
  30    integer,parameter::hi_rez_fov=221    ! high resolution pixel width
  31    integer,parameter::lo_rez_fov=221    ! low  resolution pixel width
  32    integer,parameter::time_dims=6       ! Time dimension
  33    integer,parameter::nfile_max = 54     ! each hdf file contains ~50min of data
  34                                         ! at most 8 files for a 6-h time window
  35 ! interface variable
  36    integer iret                         ! return status
  37    integer(HID_T) fhnd1, fhnd2, fhnd3   ! file handle
  38    integer(HID_T) ahnd1, ahnd2, ahnd3   ! attribute handle
  39    integer(HID_T) dhnd1, dhnd2, dhnd3   ! dataset handle
  40    integer(HID_T) shnd1, shnd2, shnd3   ! dataspace handle
  41    integer(HSIZE_T) sz1(3)              ! array size 1
  42    integer(HSIZE_T) sz2(3)              ! array size 2
  43
  44    integer(4) :: nscan                       ! NumberOfScans
  45    integer(4) :: year(max_scan),month(max_scan),day(max_scan),hour(max_scan),minute(max_scan),second(max_scan)  !ScanTime
  46
  47    type AM2_COMMON_SCANTIME
  48       real(8)    tai93sec
  49       integer(2) year
  50       integer(2) month
  51       integer(2) day
  52       integer(2) hour
  53       integer(2) minute
  54       integer(2) second
  55       integer(2) ms
  56       integer(2) reserve
  57    endtype
  58
  59 ! array data
  60    type(AM2_COMMON_SCANTIME) st(max_scan)  ! scantime
  61    real(4) :: lathi(hi_rez_fov,max_scan)   ! lat for high-fre
  62    real(4) :: latlo(lo_rez_fov,max_scan)   ! lat for low-fre
  63    real(4) :: lonhi(hi_rez_fov,max_scan)   ! lon for high-fre
  64    real(4) :: lonlo(lo_rez_fov,max_scan)   ! lon for low-fre
  65
  66    real(4)    :: ear_in(lo_rez_fov,max_scan)  ! earth incidence
  67    real(4)    :: ear_az(lo_rez_fov,max_scan)  ! earth azimuth
  68    real(4)    :: sun_az(lo_rez_fov,max_scan)  ! sun_azimuth
  69    real(4)    :: sun_el(lo_rez_fov,max_scan)  ! sun_elevation
  70    real(4)    :: sun_zen(lo_rez_fov,max_scan) ! sun_zenith
  71    real(4)    :: sca                          ! ScaleFactor
  72
  73    real(4)    :: clw(lo_rez_fov,max_scan)     ! obs retrieved cloud liquid water
  74    integer(4) :: surftype(lo_rez_fov,max_scan)! obs retrieved surfacetype
  75    integer(4) :: qcflg(lo_rez_fov,max_scan)   ! obs retrieved qualityflag
  76
  77
  78    real(r_kind)            :: R90    = 90.0_r_kind
  79    real(r_kind),parameter  :: tbmin  = 50._r_kind
  80    real(r_kind),parameter  :: tbmax  = 550._r_kind
  81
  82    real(4)          :: sat_zenith(lo_rez_fov,max_scan)
  83    real(4)          :: sat_azimuth(lo_rez_fov,max_scan)
  84
  85    real(kind=8)                   :: obs_time
  86    type (datalink_type),pointer   :: head, p, current, prev
  87    type(info_type)                :: info
  88    type(model_loc_type)           :: loc
  89
  90    integer(i_kind)    :: idate5(6)
  91
  92    integer(i_kind)   :: inst,platform_id,satellite_id,sensor_id
  93    real(r_kind)      :: tb, crit
  94    integer(i_kind)   :: ifgat, iout, iobs
  95    logical           :: outside, outside_all, iuse
  96
  97    integer           :: i,j,k,l,m,n, ifile, landsea_mask
  98    logical           :: found, head_found, head_allocated
  99
 100 ! Other work variables
 101    real(r_kind)     :: dlon_earth,dlat_earth
 102    integer(i_kind)  :: num_gmi_local, num_gmi_global, num_gmi_used, num_gmi_thinned
 103    integer(i_kind)  :: num_gmi_used_tmp, num_gmi_file
 104    integer(i_kind)  :: num_gmi_local_local, num_gmi_global_local, num_gmi_file_local
 105    integer(i_kind)  :: itx, itt
 106    character(80)    :: filename1, filename2
 107    integer          :: nchan,ifov,iscan,ichannels
 108    integer          :: nfile
 109    character(80)    :: fname_tb(nfile_max)
 110    character(80)    :: fname_clw(nfile_max)
 111    logical          :: fexist, got_clw_file
 112
 113 ! Allocatable arrays
 114    integer(i_kind),allocatable  :: ptotal(:)
 115    real,allocatable             :: in(:), out(:)
 116    real(r_kind),allocatable     :: data_all(:)
 117
 118    real,allocatable             :: obstime(:,:)
 119
 120    real(r_kind)    :: sun_zenith, sun_azimuth
 121
 122    integer,parameter  :: num_low_freq_chan=9
 123    real(4)            :: tb_low(num_low_freq_chan,lo_rez_fov,max_scan)
 124    integer,parameter  :: num_hig_freq_chan=4
 125    real(4)            :: tb_hig(num_hig_freq_chan,hi_rez_fov,max_scan)
 126
 127
 128    if (trace_use) call da_trace_entry("da_read_obs_hdf5gmi")
 129
 130 !  0.0  Initialize variables
 131 !-----------------------------------
 132    head_allocated = .false.
 133    platform_id  = 37  ! Table-2 Col 1 corresponding to 'gpm'
 134    satellite_id = 1   ! Table-2 Col 3
 135    sensor_id    = 71  ! Table-3 Col 2 corresponding to 'gmi'
 136
 137    allocate(ptotal(0:num_fgat_time))
 138    ptotal(0:num_fgat_time) = 0
 139    iobs = 0                 ! for thinning, argument is inout
 140    num_gmi_file    = 0
 141    num_gmi_local   = 0
 142    num_gmi_global  = 0
 143    num_gmi_used    = 0
 144    num_gmi_thinned = 0
 145
 146    nscan=max_scan
 147    sz1(1)=max_scan
 148    sz1(2)=lo_rez_fov
 149
 150    do i = 1, rtminit_nsensor
 151       if (platform_id  == rtminit_platform(i) &
 152           .and. satellite_id == rtminit_satid(i)    &
 153           .and. sensor_id    == rtminit_sensor(i)) then
 154          inst = i
 155          exit
 156       end if
 157    end do
 158    if (inst == 0) then
 159       call da_warning(__FILE__,__LINE__, &
 160           (/"The combination of Satellite_Id and Sensor_Id for GMI is not found"/))
 161       if (trace_use) call da_trace_exit("da_read_obs_hdf5gmi")
 162       return
 163    end if
 164
 165 ! Initialize HDF5 library and Fortran90 interface
 166    call H5open_f(iret)
 167    if(iret.lt.0)then
 168       call da_warning(__FILE__,__LINE__, &
 169            (/"Problems in Initializing HDF5 library. Can not read GMI HDF5 data. "/))
 170       if (trace_use) call da_trace_exit("da_read_obs_hdf5gmi")
 171       return
 172    endif
 173
 174    nchan = iv%instid(inst)%nchan
 175    write(unit=stdout,fmt=*)'GMI nchan: ',nchan
 176    allocate(data_all(1:nchan))
 177
 178 ! 1.0 Assign file names and prepare to read gmi files
 179 !-------------------------------------------------------------------------
 180    nfile       = 0  !initialize
 181    fname_tb(:) = '' !initialize
 182    ! first check if L1SGRTBR.h5 is available
 183    filename1 = trim(infile_tb)//'.HDF5'
 184    filename2 = trim(infile_clw)//'.HDF5'
 185    inquire (file=filename1, exist=fexist)
 186    if ( fexist ) then
 187       nfile = 1
 188       fname_tb(nfile)  = filename1
 189       fname_clw(nfile) = filename2
 190    else
 191       ! check if 1B.GPM.GMI.TB0x.HDF5 is available for multiple input files
 192       ! here 0x is the input file sequence number
 193       ! do not confuse it with fgat time slot index
 194       do i = 1, nfile_max
 195          write(filename1,fmt='(A,A,I2.2,A)') trim(infile_tb),'-',i,'.HDF5'
 196          write(filename2,fmt='(A,A,I2.2,A)') trim(infile_clw),'-',i,'.HDF5'
 197          inquire (file=filename1, exist=fexist)
 198          if ( fexist ) then
 199             nfile = nfile + 1
 200             fname_tb(nfile)  = filename1
 201             fname_clw(nfile) = filename2
 202          write(unit=stdout,fmt=*)'fname_tb=',nfile,filename1
 203          else
 204             exit
 205          end if
 206       end do
 207    end if
 208
 209    if ( nfile == 0 ) then
 210       call da_warning(__FILE__,__LINE__, &
 211          (/"No valid GMI L1B.HDF5 or L1B.HDF5 file found."/))
 212       if (trace_use) call da_trace_exit("da_read_obs_hdf5gmi")
 213       return
 214    end if
 215
 216    infile_loop:  do ifile = 1, nfile
 217       num_gmi_file_local    = 0
 218       num_gmi_local_local   = 0
 219       num_gmi_global_local  = 0
 220
 221    ! open HDF5 file for read
 222       call H5Fopen_f(fname_tb(ifile),H5F_ACC_RDONLY_F,fhnd1,iret,H5P_DEFAULT_F)
 223       if(iret.lt.0)then
 224          call da_warning(__FILE__,__LINE__, &
 225               (/"Cannot open HDF5 file "//trim(fname_tb(ifile))/))
 226          cycle infile_loop
 227       endif
 228       got_clw_file = .false.
 229       call H5Fopen_f(fname_clw(ifile),H5F_ACC_RDONLY_F,fhnd2,iret,H5P_DEFAULT_F)
 230       if ( iret == 0 ) then
 231          got_clw_file = .true.
 232       endif
 233       ! to do: when got_clw_file=.true., need to check GranuleID for consistency
 234       ! betweee tb and clw files
 235
 236     ! read lat
 237       call H5Dopen_f(fhnd1,'/S1/Latitude',dhnd1,iret)
 238       sz1(1)=max_scan
 239       sz1(2)=lo_rez_fov
 240       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,latlo,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 241       if(iret.lt.0)then
 242          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Latitude"/))
 243       endif
 244       call H5Dclose_f(dhnd1,iret)
 245
 246     ! read lon
 247       call H5Dopen_f(fhnd1,'/S1/Longitude',dhnd1,iret)
 248       sz1(1)=max_scan
 249       sz1(2)=lo_rez_fov
 250
 251       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,lonlo,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 252       if(iret.lt.0)then
 253          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: Latitude"/))
 254       endif
 255       call H5Dclose_f(dhnd1,iret)
 256
 257    !   write(unit=stdout,fmt=*)'NumberOfScans(RETRIEVE BY ARRAY SIZE): ',nscan
 258    !   write(unit=stdout,fmt=*)'OverlapScans(FIXED VALUE): ',ovr
 259
 260       ! sample display
 261       write(unit=stdout,fmt=*)'latitude,longitude(pixel=1,scan=1): ',&
 262       latlo(1,1),lonlo(1,1),latlo(2,1),lonlo(2,1),latlo(1,2),lonlo(1,2)
 263
 264
 265     ! read array: satellite_zenith_angle
 266       call H5Dopen_f(fhnd1,'/S1/incidenceAngle',dhnd1,iret)
 267       sz1(1)=max_scan
 268       sz1(2)=lo_rez_fov
 269       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,ear_in,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 270       if(iret.lt.0)then
 271          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: satellite_zenith_angle"/))
 272       endif
 273       call H5Dclose_f(dhnd1,iret)
 274    ! sample display
 275       write(unit=stdout,fmt=*)'sat_zenith(pixel=1,scan=1): ',ear_in(1,1)
 276
 277    ! read array: satellite_azimuth_angle
 278       call H5Dopen_f(fhnd1,'/S1/satAzimuthAngle',dhnd1,iret)
 279       sz1(1)=max_scan
 280       sz1(2)=lo_rez_fov
 281
 282       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,ear_az,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 283       if(iret.lt.0)then
 284          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: satellite_Azimuth_angle"/))
 285       endif
 286       call H5Dclose_f(dhnd1,iret)
 287    ! sample display
 288       write(unit=stdout,fmt=*)'sat_azimuth(pixel=1,scan=1): ',ear_az(1,1),ear_az(1,2),ear_az(2,1)
 289
 290     ! read array: sun_azimuth_angle
 291       call H5Dopen_f(fhnd1,'/S1/solarAzimuthAngle',dhnd1,iret)
 292       sz1(1)=max_scan
 293       sz1(2)=lo_rez_fov
 294
 295       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,sun_az,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 296       if(iret.lt.0)then
 297          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: sun_Azimuth_angle"/))
 298       endif
 299       call H5Dclose_f(dhnd1,iret)
 300    ! sample display
 301       write(unit=stdout,fmt=*)'sun_azimuth(pixel=1,scan=1): ',sun_az(1,1)
 302
 303     ! read array: solar_ZenAngle
 304       call H5Dopen_f(fhnd1,'/S1/solarZenAngle',dhnd1,iret)
 305       sz1(1)=max_scan
 306       sz1(2)=lo_rez_fov
 307
 308        call H5Dread_f(dhnd1,H5T_NATIVE_REAL,sun_el,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 309       if(iret.lt.0)then
 310          call da_warning(__FILE__,__LINE__,(/"HDF5 read error for: sun_elevation_angle"/))
 311       endif
 312       call H5Dclose_f(dhnd1,iret)
 313    ! sample display
 314       write(unit=stdout,fmt=*)'sun_Zenangle(pixel=1,scan=1): ',sun_el(1,1)
 315
 316       sun_zen(:,:)    =sun_el(:,:)
 317       sat_zenith(:,:) =ear_in(:,:)
 318       sat_azimuth(:,:)=ear_az(:,:)
 319
 320    ! read array: scantime
 321    ! read Year
 322       call H5Dopen_f(fhnd1,'/S1/ScanTime/Year',dhnd1,iret)
 323       sz1(1)=max_scan
 324       call H5Dread_f(dhnd1,H5T_NATIVE_INTEGER,year,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 325       if(iret.lt.0)then
 326          call da_warning(__FILE__,__LINE__, &
 327              (/"HDF5 read error for: ScanTime/Year"/))
 328       endif
 329       call H5Dclose_f(dhnd1,iret)
 330
 331    ! read Month
 332       call H5Dopen_f(fhnd1,'/S1/ScanTime/Month',dhnd1,iret)
 333       sz1(1)=max_scan
 334       call H5Dread_f(dhnd1,H5T_NATIVE_INTEGER,month,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 335       if(iret.lt.0)then
 336          call da_warning(__FILE__,__LINE__, &
 337              (/"HDF5 read error for: ScanTime/Month"/))
 338       endif
 339       call H5Dclose_f(dhnd1,iret)
 340
 341    ! read Day
 342       call H5Dopen_f(fhnd1,'/S1/ScanTime/DayOfMonth',dhnd1,iret)
 343       sz1(1)=max_scan
 344       call H5Dread_f(dhnd1,H5T_NATIVE_INTEGER,day,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 345       if(iret.lt.0)then
 346          call da_warning(__FILE__,__LINE__, &
 347              (/"HDF5 read error for: ScanTime/Day"/))
 348       endif
 349       call H5Dclose_f(dhnd1,iret)
 350
 351    ! read Hour
 352       call H5Dopen_f(fhnd1,'/S1/ScanTime/Hour',dhnd1,iret)
 353       sz1(1)=max_scan
 354       call H5Dread_f(dhnd1,H5T_NATIVE_INTEGER,hour,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 355       if(iret.lt.0)then
 356          call da_warning(__FILE__,__LINE__, &
 357              (/"HDF5 read error for: ScanTime/Hour"/))
 358       endif
 359       call H5Dclose_f(dhnd1,iret)
 360
 361      ! read Minute
 362       call H5Dopen_f(fhnd1,'/S1/ScanTime/Minute',dhnd1,iret)
 363       sz1(1)=max_scan
 364       call H5Dread_f(dhnd1,H5T_NATIVE_INTEGER,minute,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 365       if(iret.lt.0)then
 366          call da_warning(__FILE__,__LINE__, &
 367              (/"HDF5 read error for: ScanTime/Minute"/))
 368       endif
 369       call H5Dclose_f(dhnd1,iret)
 370
 371      ! read Second
 372       call H5Dopen_f(fhnd1,'/S1/ScanTime/Second',dhnd1,iret)
 373       sz1(1)=max_scan
 374       call H5Dread_f(dhnd1,H5T_NATIVE_INTEGER,second,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 375       if(iret.lt.0)then
 376          call da_warning(__FILE__,__LINE__, &
 377              (/"HDF5 read error for: ScanTime/Second"/))
 378       endif
 379       call H5Dclose_f(dhnd1,iret)
 380
 381    ! sample display
 382       allocate  (obstime(1:time_dims,1:nscan))  ! year, month, day, hour, min, sec
 383       do j = 1, nscan
 384          obstime(1,j) = year(j)
 385          obstime(2,j) = month(j)
 386          obstime(3,j) = day(j)
 387          obstime(4,j) = hour(j)
 388          obstime(5,j) = minute(j)
 389          obstime(6,j) = second(j)
 390       end do
 391       write(unit=stdout,fmt=*)'time(scan=1) year: ',year(1),' month:',month(1),' day: ',day(1),&
 392          ' hour: ',hour(1),' minute: ',minute(1),' second: ',second(1)
 393
 394     ! read Tb_low--S1---9 low-frequence channel
 395       call H5Dopen_f(fhnd1,'/S1/Tb',dhnd1,iret)
 396
 397       sz1(1)=max_scan
 398       sz1(2)=lo_rez_fov
 399       sz1(3)=num_low_freq_chan
 400       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,Tb_low,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 401       if(iret.lt.0)then
 402          call da_warning(__FILE__,__LINE__, &
 403              (/"HDF5 read error for: Tb_low"/))
 404       endif
 405       call H5Dclose_f(dhnd1,iret)
 406
 407      ! sample display
 408      ! write(unit=stdout,fmt=*)'tb_low(pixel=1,scan=1,chan=1): ',&
 409      ! tb_low(1,1,1),tb_low(1,1,2),tb_low(2,1,1),tb_low(2,1,2),tb_low(1,2,1),tb_low(1,2,2)
 410
 411      ! read Tb_hig--S2---4 hig-frequence channel
 412       call H5Dopen_f(fhnd1,'/S2/Tb',dhnd1,iret)
 413       sz1(1)=max_scan
 414       sz1(2)=hi_rez_fov
 415       sz1(3)=num_hig_freq_chan
 416       call H5Dread_f(dhnd1,H5T_NATIVE_REAL,Tb_hig,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 417       if(iret.lt.0)then
 418          call da_warning(__FILE__,__LINE__, &
 419              (/"HDF5 read error for: Tb_hig"/))
 420       endif
 421       call H5Dclose_f(dhnd1,iret)
 422
 423      ! sample display
 424      ! write(unit=stdout,fmt=*)'tb_hig(pixel=1,scan=1,chan=1): ',&
 425      ! tb_hig(1,1,1),tb_hig(1,1,2),tb_hig(2,1,1),tb_hig(2,1,2),tb_hig(1,2,1),tb_hig(1,2,2)
 426
 427      ! close infile_tb and HDF5
 428       call H5Fclose_f(fhnd1,iret)
 429
 430       if ( got_clw_file ) then
 431       ! read CLW from infile_clw:
 432          call H5Dopen_f(fhnd2,'/S1/cloudWaterPath',dhnd2,iret)
 433          sz1(1)=max_scan
 434          sz1(2)=lo_rez_fov
 435          call H5Dread_f(dhnd2,H5T_NATIVE_REAL,clw,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 436          if(iret.lt.0)then
 437             call da_warning(__FILE__,__LINE__, &
 438                (/"HDF5 read error for: CLW data"/))
 439          endif
 440          call H5Dclose_f(dhnd2,iret)
 441       ! sample display
 442          write(unit=stdout,fmt=*)'clw(pixel=1,scan=1): ',clw(100,1),clw(2,100),clw(100,2)
 443
 444       ! read SurfaceType from infile_clw:
 445          call H5Dopen_f(fhnd2,'/S1/surfaceTypeIndex',dhnd2,iret)
 446          sz1(1)=max_scan
 447          sz1(2)=lo_rez_fov
 448          call H5Dread_f(dhnd2,H5T_NATIVE_INTEGER,surftype,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 449          if(iret.lt.0)then
 450             call da_warning(__FILE__,__LINE__, &
 451                (/"HDF5 read error for: CLW data"/))
 452          endif
 453          call H5Dclose_f(dhnd2,iret)
 454       ! sample display
 455       !   write(unit=stdout,fmt=*)'surftype(pixel=1,scan=1): ',&
 456       !  surftype(1,1),surftype(2,1),surftype(1,2)
 457
 458       ! read QualityFlag from infile_clw:
 459          call H5Dopen_f(fhnd2,'/S1/qualityFlag',dhnd2,iret)
 460          sz1(1)=max_scan
 461          sz1(2)=lo_rez_fov
 462          call H5Dread_f(dhnd2,H5T_NATIVE_INTEGER,qcflg,sz1,iret,H5S_ALL_F,H5S_ALL_F)
 463          if(iret.lt.0)then
 464             call da_warning(__FILE__,__LINE__, &
 465                (/"HDF5 read error for: CLW data"/))
 466          endif
 467          call H5Dclose_f(dhnd2,iret)
 468       ! sample display
 469       !   write(unit=stdout,fmt=*)'qcflg(pixel=1,scan=1): ',&
 470       !  qcflg(1,1),qcflg(2,1),qcflg(1,2)
 471
 472       ! close infile_clw and HDF5
 473          call H5Fclose_f(fhnd2,iret)
 474       end if
 475 !****display all the hdf5 data with asc format****
 476 !      do i = 1, nscan
 477 !      do j = 1, lo_rez_fov
 478 !      write(65,fmt='(2i4,2f8.3,13f10.3)') j,i,latlo(j,i),lonlo(j,i),sun_zen(j,i),&
 479 !      sat_zenith(j,i),sat_azimuth(j,i),(tb_low(k,j,i),k=1,num_low_freq_chan)
 480 !       write(64,fmt='(2i4,2f8.3,f10.3,2i4)') j,i,latlo(j,i),lonlo(j,i),clw(j,i),surftype(j,i),qcflg(j,k)
 481 !      if(abs(sat_zenith(j,i)).gt.80)write(unit=stdout,fmt='(a,2i4,f12.2)')'bad_sat_zenith',i,j,sat_zenith(j,i)
 482 !      if(abs(sat_azimuth(j,i)).gt.180)write(unit=stdout,fmt='(a,2i4,f12.2)'),'bad_sat_azimuth',i,j,sat_azimuth(j,i)
 483 !      if(abs(sun_zen(j,i)).gt.180)write(unit=stdout,fmt='(a,2i4,f12.2)')'bad_sun_zenith',i,j,sun_zen(j,i)
 484 !      if(abs(sun_az(j,i)).gt.180)write(unit=stdout,fmt='(a,2i4,f12.2)')'bad_sun_azi',i,j,sun_az(j,i)
 485 !      end do
 486 !      end do
 487
 488 ! 2.0 Loop to read hdf file and assign information to a sequential structure
 489 !-------------------------------------------------------------------------
 490
 491    ! Allocate arrays to hold data
 492       if ( .not. head_allocated ) then
 493          allocate (head)
 494          nullify  ( head % next )
 495          p => head
 496          head_allocated = .true.
 497       end if
 498    ! start scan_loop
 499       scan_loop:     do iscan=1, nscan
 500          do i = 1, 6
 501             idate5(i)=obstime(i, iscan)
 502          end do
 503          call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
 504          if ( obs_time < time_slots(0) .or.  &
 505             obs_time >= time_slots(num_fgat_time) ) cycle scan_loop
 506          do ifgat=1,num_fgat_time
 507             if ( obs_time >= time_slots(ifgat-1) .and.  &
 508                obs_time  < time_slots(ifgat) ) exit
 509          end do
 510
 511       ! start fov_loop
 512          fov_loop:      do ifov=1, lo_rez_fov
 513             num_gmi_file       = num_gmi_file + 1
 514             num_gmi_file_local = num_gmi_file_local + 1
 515             info%lat  =  latlo(ifov,iscan)
 516             info%lon  =  lonlo(ifov,iscan)
 517
 518             call da_llxy (info, loc, outside, outside_all)
 519             if (outside_all) cycle fov_loop
 520
 521             num_gmi_global       = num_gmi_global + 1
 522             num_gmi_global_local = num_gmi_global_local + 1
 523             ptotal(ifgat) = ptotal(ifgat) + 1
 524             if (outside) cycle fov_loop   ! No good for this PE
 525
 526             landsea_mask = 0    !0:land; 1:sea !juded from amsr2
 527             if(surftype(ifov,iscan) < 3.0 .or. surftype(ifov,iscan) > 12.0) &  !surftype 3-12:land
 528             landsea_mask = 1
 529
 530             info%elv = 0.0
 531
 532             num_gmi_local       = num_gmi_local + 1
 533             num_gmi_local_local = num_gmi_local_local + 1
 534             write(unit=info%date_char, &
 535             fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 536                idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
 537                ':', idate5(5), ':', idate5(6)
 538
 539 ! 3.0  Make Thinning
 540 ! Map obs to thinning grid
 541 !-------------------------------------------------------------------
 542             if (thinning) then
 543                dlat_earth = info%lat !degree
 544                dlon_earth = info%lon
 545                if (dlon_earth<zero)  dlon_earth = dlon_earth+r360
 546                if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 547                dlat_earth = dlat_earth*deg2rad !radian
 548                dlon_earth = dlon_earth*deg2rad
 549                crit = 1.
 550                call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 551                if (.not. iuse) then
 552                   num_gmi_thinned = num_gmi_thinned+1
 553                   cycle fov_loop
 554                end if
 555             end if
 556
 557             num_gmi_used = num_gmi_used + 1
 558             data_all = missing_r
 559
 560             do k=1,num_low_freq_chan
 561                tb = tb_low(k,ifov,iscan)
 562                if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 563                data_all(k)= tb
 564             enddo
 565
 566             do k=1,num_hig_freq_chan
 567                tb = tb_hig(k,ifov,iscan)
 568                if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 569                data_all(9+k)= tb
 570             enddo
 571 !           write(unit=stdout,fmt=*)'data_all: '
 572
 573 ! 4.0 assign information to sequential radiance structure
 574 !--------------------------------------------------------------------------
 575             allocate ( p % tb_inv (1:nchan ))
 576             p%info             = info
 577             p%loc              = loc
 578             p%landsea_mask     = landsea_mask
 579             p%scanpos          = ifov
 580             p%satzen           = sat_zenith(ifov,iscan)
 581             p%satazi           = sat_azimuth(ifov,iscan)
 582             p%solzen           = sun_zen(ifov,iscan)
 583             p%solazi           = sun_az(ifov,iscan)
 584             p%clw              = clw(ifov,iscan)
 585             p%tb_inv(1:nchan)  = data_all(1:nchan)
 586             p%sensor_index     = inst
 587             p%ifgat            = ifgat
 588
 589             allocate (p%next)   ! add next data
 590             p => p%next
 591             nullify (p%next)
 592          end do fov_loop
 593       end do scan_loop
 594
 595    ! Dellocate arrays
 596       deallocate  (obstime)
 597
 598       write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_gmi_file    : ',num_gmi_file_local
 599       write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_gmi_global  : ',num_gmi_global_local
 600       write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_gmi_local   : ',num_gmi_local_local
 601    end do infile_loop
 602
 603    call H5close_f(iret)
 604
 605    deallocate(data_all) ! Deallocate data arrays
 606
 607    if (thinning .and. num_gmi_global > 0 ) then
 608 #ifdef DM_PARALLEL
 609    ! Get minimum crit and associated processor index.
 610       j = 0
 611       do ifgat = 1, num_fgat_time
 612          j = j + thinning_grid(inst,ifgat)%itxmax
 613       end do
 614
 615       allocate ( in  (j) )
 616       allocate ( out (j) )
 617       j = 0
 618       do ifgat = 1, num_fgat_time
 619          do i = 1, thinning_grid(inst,ifgat)%itxmax
 620             j = j + 1
 621             in(j) = thinning_grid(inst,ifgat)%score_crit(i)
 622          end do
 623       end do
 624       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 625
 626       call wrf_dm_bcast_real (out, j)
 627
 628       j = 0
 629       do ifgat = 1, num_fgat_time
 630          do i = 1, thinning_grid(inst,ifgat)%itxmax
 631             j = j + 1
 632             if ( ABS(out(j)-thinning_grid(inst,ifgat)%score_crit(i)) > 1.0E-10 ) &
 633             thinning_grid(inst,ifgat)%ibest_obs(i) = 0
 634          end do
 635       end do
 636
 637       deallocate( in  )
 638       deallocate( out )
 639
 640 #endif
 641
 642    ! Delete the nodes which being thinning out
 643       p => head
 644       prev => head
 645       head_found = .false.
 646       num_gmi_used_tmp = num_gmi_used
 647       do j = 1, num_gmi_used_tmp
 648          n = p%sensor_index
 649          ifgat = p%ifgat
 650          found = .false.
 651
 652          do i = 1, thinning_grid(n,ifgat)%itxmax
 653             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 654                found = .true.
 655                exit
 656             end if
 657          end do
 658
 659       ! free current data
 660          if ( .not. found ) then
 661             current => p
 662             p => p%next
 663             if ( head_found ) then
 664                prev%next => p
 665             else
 666                head => p
 667                prev => p
 668             end if
 669             deallocate ( current % tb_inv )
 670             deallocate ( current )
 671             num_gmi_thinned = num_gmi_thinned + 1
 672             num_gmi_used = num_gmi_used - 1
 673             continue
 674          end if
 675
 676          if ( found .and. head_found ) then
 677             prev => p
 678             p => p%next
 679             continue
 680          end if
 681
 682          if ( found .and. .not. head_found ) then
 683             head_found = .true.
 684             head => p
 685             prev => p
 686             p => p%next
 687          end if
 688
 689       end do
 690
 691    end if  ! End of thinning
 692
 693    iv%total_rad_pixel   = iv%total_rad_pixel + num_gmi_used
 694    iv%total_rad_channel = iv%total_rad_channel + num_gmi_used*nchan
 695
 696    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_gmi_used
 697    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_gmi_global
 698
 699    do i = 1, num_fgat_time
 700       ptotal(i) = ptotal(i) + ptotal(i-1)
 701       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 702    end do
 703    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 704       write(unit=message(1),fmt='(A,I10,A,I10)') &
 705           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 706       call da_warning(__FILE__,__LINE__,message(1:1))
 707    endif
 708
 709    write(unit=stdout,fmt='(a)') 'GMI data counts: '
 710    write(stdout,fmt='(a,i7)') ' In file: ',num_gmi_file
 711    write(stdout,fmt='(a,i7)') ' Global : ',num_gmi_global
 712    write(stdout,fmt='(a,i7)') ' Local  : ',num_gmi_local
 713    write(stdout,fmt='(a,i7)') ' Used   : ',num_gmi_used
 714    write(stdout,fmt='(a,i7)') ' Thinned: ',num_gmi_thinned
 715
 716 !  5.0 allocate innovation radiance structure
 717 !----------------------------------------------------------------
 718
 719    if (num_gmi_used > 0) then
 720       iv%instid(inst)%num_rad  = num_gmi_used
 721       iv%instid(inst)%info%nlocal = num_gmi_used
 722       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 723          'Allocating space for radiance innov structure', &
 724          inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
 725
 726       call da_allocate_rad_iv (inst, nchan, iv)
 727
 728    end if
 729
 730    write(unit=stdout,fmt=*)'for test: ',inst, nchan
 731
 732 !  6.0 assign sequential structure to innovation structure
 733 !-------------------------------------------------------------
 734    p => head
 735
 736    do n = 1, num_gmi_used
 737       i = p%sensor_index
 738       call da_initialize_rad_iv (i, n, iv, p)
 739       current => p
 740       p => p%next
 741    ! free current data
 742       deallocate ( current % tb_inv )
 743       deallocate ( current )
 744    end do
 745    deallocate ( p )
 746    deallocate (ptotal)
 747
 748    if (trace_use) call da_trace_exit("da_read_obs_hdf5gmi")
 749 #else
 750    call da_error(__FILE__,__LINE__,(/"Needs to be compiled with HDF5 library"/))
 751 #endif
 752 end subroutine da_read_obs_hdf5gmi