var/da/da_radiance/da_read_obs_netcdf4ahi_geocat.inc

   1 subroutine da_read_obs_netcdf4ahi_geocat (iv, infile_tb, infile_clp)
   2    !--------------------------------------------------------
   3    !  Purpose: read in GEOCAT AHI Level-1 and Level-2 data in NETCDF4 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/22 - Creation         Yuanbing Wang, NUIST/CAS, NCAR/NESL/MMM/DAS
  13    !  To be devoloped: 1.time information; 2.dimension sequence
  14    !------------------------------------------------------------------------------
  15
  16    use netcdf
  17    implicit none
  18
  19    character(len=*), intent(in)    :: infile_tb, infile_clp
  20    type(iv_type),    intent(inout) :: iv
  21
  22 ! fixed parameter values
  23    integer,parameter::time_dims=6       ! Time dimension
  24    integer,parameter::nfile_max = 8     ! each netcdf file contains
  25    real,parameter::add_offset_tb1=285.0
  26    real,parameter::add_offset_tb2=235.0
  27    real,parameter::add_offset_tb3=260.0
  28    real,parameter::add_offset_tb4=240.0
  29    real,parameter::add_offset_saz=90.0
  30    real,parameter::scale_factor_tb1=0.00350962858973968
  31    real,parameter::scale_factor_tb2=0.00198370311593982
  32    real,parameter::scale_factor_tb3=0.00274666585283975
  33    real,parameter::scale_factor_tb4=0.00213629566331980
  34    real,parameter::scale_factor_lat=0.00274666585283975
  35    real,parameter::scale_factor_lon=0.00549333170567949
  36    real,parameter::scale_factor_saz=0.00274666585283975
  37
  38 ! interface variable
  39    integer iret, rcode, ncid                      ! return status
  40
  41 ! array data
  42    real(4), allocatable :: vlatitude(:,:)  !  value for latitude
  43    real(4), allocatable :: vlongitude(:,:) !  value for longitude
  44
  45    real(4), allocatable    :: tbb(:,:,:)  ! tb for band 7-16
  46    real(4), allocatable    :: sat_zenith(:,:)
  47
  48    byte, allocatable    ::cloud_mask(:,:)
  49
  50    real(r_kind),parameter  :: tbmin  = 50._r_kind
  51    real(r_kind),parameter  :: tbmax  = 550._r_kind
  52
  53    real(kind=8)                   :: obs_time
  54    type (datalink_type),pointer   :: head, p, current, prev
  55    type(info_type)                :: info
  56    type(model_loc_type)           :: loc
  57
  58    integer(i_kind)    :: idate5(6)
  59    character(len=80)  :: filename,str_tmp
  60
  61    integer(i_kind)   :: inst,platform_id,satellite_id,sensor_id
  62    real(r_kind)      :: tb, crit
  63    integer(i_kind)   :: ifgat, iout, iobs
  64    logical           :: outside, outside_all, iuse
  65
  66    integer           :: i,j,k,l,m,n, ifile, landsea_mask
  67    logical           :: found, head_found, head_allocated
  68
  69 ! Other work variables
  70    real(r_kind)     :: dlon_earth,dlat_earth
  71    integer(i_kind)  :: num_ahi_local, num_ahi_global, num_ahi_used, num_ahi_thinned
  72    integer(i_kind)  :: num_ahi_used_tmp, num_ahi_file
  73    integer(i_kind)  :: num_ahi_local_local, num_ahi_global_local, num_ahi_file_local
  74    integer(i_kind)  :: itx, itt
  75    character(80)    :: filename1,filename2
  76    integer          :: nchan,nlongitude,nlatitude,ilongitude,ilatitude,ichannels
  77    integer          :: lonstart,latstart
  78    integer          :: LatDimID,LonDimID
  79    integer          :: latid,lonid,tbb_id,sazid,cltyid
  80    integer          :: nfile
  81    character(80)    :: fname_tb(nfile_max),fname_clp(nfile_max)
  82    integer          :: vtype
  83    character(80)    :: vname
  84    logical          :: fexist,got_clp_file
  85    integer          :: ahi_info_unit
  86
  87 ! Allocatable arrays
  88    integer(i_kind),allocatable  :: ptotal(:)
  89    real,allocatable             :: in(:), out(:)
  90    real(r_kind),allocatable     :: data_all(:)
  91
  92    character(len=80) tbb_name(10)
  93    data tbb_name/'himawari_8_ahi_channel_7_brightness_temperature',&
  94                  'himawari_8_ahi_channel_8_brightness_temperature',&
  95                  'himawari_8_ahi_channel_9_brightness_temperature',&
  96                  'himawari_8_ahi_channel_10_brightness_temperature',&
  97                  'himawari_8_ahi_channel_11_brightness_temperature',&
  98                  'himawari_8_ahi_channel_12_brightness_temperature',&
  99                  'himawari_8_ahi_channel_13_brightness_temperature',&
 100                  'himawari_8_ahi_channel_14_brightness_temperature',&
 101                  'himawari_8_ahi_channel_15_brightness_temperature',&
 102                  'himawari_8_ahi_channel_16_brightness_temperature'/
 103
 104    if (trace_use) call da_trace_entry("da_read_obs_netcdf4ahi_geocat")
 105
 106 !  0.0  Initialize variables
 107 !-----------------------------------
 108    head_allocated = .false.
 109    platform_id  = 31  ! Table-2 Col 1 corresponding to 'himawari'
 110    satellite_id = 8   ! Table-2 Col 3
 111    sensor_id    = 56  ! Table-3 Col 2 corresponding to 'ahi'
 112
 113    allocate(ptotal(0:num_fgat_time))
 114    ptotal(0:num_fgat_time) = 0
 115    iobs = 0                 ! for thinning, argument is inout
 116    num_ahi_file    = 0
 117    num_ahi_local   = 0
 118    num_ahi_global  = 0
 119    num_ahi_used    = 0
 120    num_ahi_thinned = 0
 121
 122    do i = 1, rtminit_nsensor
 123       if (platform_id  == rtminit_platform(i) &
 124           .and. satellite_id == rtminit_satid(i)    &
 125           .and. sensor_id    == rtminit_sensor(i)) then
 126          inst = i
 127          exit
 128       end if
 129    end do
 130    if (inst == 0) then
 131       call da_warning(__FILE__,__LINE__, &
 132           (/"The combination of Satellite_Id and Sensor_Id for AHI is not found"/))
 133       if (trace_use) call da_trace_exit("da_read_obs_netcdf4ahi_geocat")
 134       return
 135    end if
 136
 137    nchan = iv%instid(inst)%nchan
 138    write(unit=stdout,fmt=*)'AHI nchan: ',nchan
 139    allocate(data_all(1:nchan))
 140
 141 ! 1.0 Assign file names and prepare to read ahi files
 142 !-------------------------------------------------------------------------
 143    nfile       = 0  !initialize
 144    fname_tb(:) = '' !initialize
 145
 146    ! first check if ahi nc file is available
 147    filename1 = trim(infile_tb)
 148    filename2 = trim(infile_clp)
 149    inquire (file=filename1, exist=fexist)
 150    if ( fexist ) then
 151       nfile = 1
 152       fname_tb(nfile)  = filename1
 153       fname_clp(nfile) = filename2
 154    else
 155       ! check if netcdf4 files are available for multiple input files
 156       ! here 0x is the input file sequence number
 157       ! do not confuse it with fgat time slot index
 158       do i = 1, nfile_max
 159          write(filename1,fmt='(A,A,I2.2,A)') trim(infile_tb),'-',i
 160          write(filename2,fmt='(A,A,I2.2,A)') trim(infile_clp),'-',i
 161          inquire (file=filename1, exist=fexist)
 162          if ( fexist ) then
 163             nfile = nfile + 1
 164             fname_tb(nfile)  = filename1
 165             fname_clp(nfile) = filename2
 166          else
 167             exit
 168          end if
 169       end do
 170    end if
 171
 172    if ( nfile == 0 ) then
 173       call da_warning(__FILE__,__LINE__, &
 174          (/"No valid AHI file found."/))
 175       if (trace_use) call da_trace_exit("da_read_obs_netcdf4ahi_geocat")
 176       return
 177    end if
 178
 179
 180   !open the data area info file
 181    call da_get_unit(ahi_info_unit)
 182    open(unit=ahi_info_unit,file='ahi_info',status='old',iostat=iret)
 183    if(iret /= 0)then
 184        call da_warning(__FILE__,__LINE__,(/"area_info file read error"/))
 185    endif
 186   !read date information
 187    read(ahi_info_unit,*)
 188    read(ahi_info_unit,*)
 189    read(ahi_info_unit,*)
 190    read(ahi_info_unit,*)
 191    read(ahi_info_unit,*)
 192    read(ahi_info_unit,*) lonstart,latstart,nlongitude,nlatitude
 193    close(ahi_info_unit)
 194    call da_free_unit(ahi_info_unit)
 195
 196    write(*,*) lonstart,latstart,nlongitude,nlatitude
 197
 198    infile_loop:  do ifile = 1, nfile
 199       num_ahi_file_local    = 0
 200       num_ahi_local_local   = 0
 201       num_ahi_global_local  = 0
 202
 203    ! open NETCDF4 L1 file for read
 204       iret = nf90_open(fname_tb(ifile), nf90_NOWRITE, ncid)
 205       if(iret /= 0)then
 206          call da_warning(__FILE__,__LINE__, &
 207               (/"Cannot open NETCDF4 file "//trim(fname_tb(ifile))/))
 208          cycle infile_loop
 209       endif
 210
 211    ! read dimensions:   latitude and longitude
 212      ! iret = nf90_inq_dimid(ncid, "lines", LatDimID)
 213      ! iret = nf90_inquire_dimension(ncid, LatDimID, len=nlatitude)
 214
 215      ! iret = nf90_inq_dimid(ncid, "elements", LonDimID)
 216      ! iret = nf90_inquire_dimension(ncid, LonDimID, len=nlongitude)
 217
 218      ! write(unit=stdout,fmt=*) nlongitude,nlatitude
 219
 220    ! read array: time
 221       iret = nf90_get_att(ncid, nf90_global, "Image_Date_Time", filename)
 222       if(iret /= 0)then
 223          call da_warning(__FILE__,__LINE__, &
 224              (/"NETCDF4 read error for: observation date"/))
 225       end if
 226          read(filename,"(I4,A1,I2,A1,I2,A1,I2,A1,I2,A1,I2,A1)") idate5(1),str_tmp,idate5(2),str_tmp,&
 227            idate5(3),str_tmp,idate5(4),str_tmp,idate5(5),str_tmp,idate5(6),str_tmp
 228
 229          write(unit=stdout,fmt=*)'observation date: ', idate5
 230
 231    ! read array: lat
 232       iret = nf90_inq_varid(ncid, 'pixel_latitude', latid)
 233       allocate(vlatitude(nlongitude,nlatitude))
 234       iret = nf90_get_var(ncid,latid,vlatitude,start=(/lonstart,latstart/), &
 235                           count=(/nlongitude,nlatitude/)) !
 236       if(iret /= 0)then
 237          call da_warning(__FILE__,__LINE__, &
 238              (/"NETCDF4 read error for: Latitude of Observation Point"/))
 239       endif
 240       do j=1,nlatitude
 241       do i=1,nlongitude
 242           vlatitude(i,j)=vlatitude(i,j) * scale_factor_lat
 243       end do
 244       end do
 245    ! sample display
 246       write(unit=stdout,fmt=*)'vlatitude(pixel=1,scan=1): ',vlatitude(1,1)
 247
 248    ! read lon
 249       iret = nf90_inq_varid(ncid, 'pixel_longitude', lonid)
 250       allocate(vlongitude(nlongitude,nlatitude))
 251       iret = nf90_get_var(ncid,lonid,vlongitude,start=(/lonstart,latstart/), &
 252                           count=(/nlongitude,nlatitude/))
 253       if(iret /= 0)then
 254           call da_warning(__FILE__,__LINE__, &
 255               (/"NETCDF4 read error for: Longitude of Observation Point"/))
 256           call da_trace_exit("da_read_obs_netcdf4ahi_geocat")
 257       endif
 258       do j=1,nlatitude
 259       do i=1,nlongitude
 260          vlongitude(i,j)=vlongitude(i,j) * scale_factor_lon
 261       end do
 262       end do
 263    ! sample display
 264       write(unit=stdout,fmt=*)'vlongitude(pixel=1,scan=1): ',vlongitude(1,1)
 265
 266    ! read array: tb for band 7-16
 267       allocate(tbb(nlongitude,nlatitude,nchan))
 268       do k=1,nchan
 269          iret = nf90_inq_varid(ncid, tbb_name(k), tbb_id)
 270          iret = nf90_get_var(ncid,tbb_id,tbb(:,:,k),start=(/lonstart,latstart/), &
 271                              count=(/nlongitude,nlatitude/))
 272          if(iret /= 0)  then
 273             call da_warning(__FILE__,__LINE__, &
 274                    (/"NETCDF4 read error for: Brightness Temperature"/))
 275          endif
 276          do j=1,nlatitude
 277          do i=1,nlongitude
 278             if(k==1) then
 279               tbb(i,j,k)=tbb(i,j,k) * scale_factor_tb1 + add_offset_tb1
 280             end if
 281             if(k>=2 .and. k<=4) then
 282               tbb(i,j,k)=tbb(i,j,k) * scale_factor_tb2 + add_offset_tb2
 283             end if
 284             if(k>=5 .and. k<=9) then
 285               tbb(i,j,k)=tbb(i,j,k) * scale_factor_tb3 + add_offset_tb3
 286             end if
 287             if(k==10) then
 288               tbb(i,j,k)=tbb(i,j,k) * scale_factor_tb4 + add_offset_tb4
 289             end if
 290          end do
 291          end do
 292          ! sample display
 293          write(unit=stdout,fmt=*) 'tbb(pixel=1,scan=1,chan=',k,'): ', tbb(1,1,k)
 294       end do
 295
 296    ! read array: satellite zenith angle
 297       iret = nf90_inq_varid(ncid, 'pixel_satellite_zenith_angle', sazid)
 298       allocate(sat_zenith(nlongitude,nlatitude))
 299       iret = nf90_get_var(ncid,sazid,sat_zenith,start=(/lonstart,latstart/), &
 300                           count=(/nlongitude,nlatitude/))
 301       if(iret /= 0)then
 302          call da_warning(__FILE__,__LINE__, &
 303             (/"NETCDF4 read error for: satellite zenith angle"/))
 304       endif
 305       do j=1,nlatitude
 306       do i=1,nlongitude
 307          sat_zenith(i,j)=sat_zenith(i,j) * scale_factor_saz + add_offset_saz
 308       end do
 309       end do
 310    ! sample display
 311       write(unit=stdout,fmt=*) 'satellite zenith angle(pixel=1,scan=1): ',sat_zenith(1,1)
 312
 313    ! close infile_tb file
 314       iret = nf90_close(ncid)
 315
 316      !open infile_clp file
 317       got_clp_file = .false.
 318       iret = nf90_open(fname_clp(ifile), nf90_NOWRITE, ncid)
 319       if ( iret == 0 ) then
 320          got_clp_file = .true.
 321       endif
 322
 323       if ( got_clp_file ) then
 324         ! read array: eps_cmask_ahi_cloud_mask
 325           rcode = nf90_inq_varid(ncid, "eps_cmask_ahi_cloud_mask", cltyid)
 326           allocate(cloud_mask(nlongitude,nlatitude))
 327           iret = nf90_get_var(ncid,cltyid,cloud_mask,start=(/lonstart,latstart/), &
 328                               count=(/nlongitude,nlatitude/))
 329           if(iret /= 0)then
 330             call da_warning(__FILE__,__LINE__,(/"NETCDF4 read error for: CLTYPE data"/))
 331           endif
 332         ! sample display
 333           write(unit=stdout,fmt=*)'cloud_mask(pixel=1,scan=1): ',cloud_mask(1,1)
 334         ! close infile_clp file
 335           iret = nf90_close(ncid)
 336       end if
 337
 338 ! 2.0 Loop to read netcdf and assign information to a sequential structure
 339 !-------------------------------------------------------------------------
 340
 341    ! Allocate arrays to hold data
 342       if ( .not. head_allocated ) then
 343          allocate (head)
 344          nullify  ( head % next )
 345          p => head
 346          head_allocated = .true.
 347       end if
 348
 349    ! start scan_loop
 350       scan_loop:     do ilatitude=1, nlatitude
 351
 352          call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
 353          if ( obs_time < time_slots(0) .or.  &
 354             obs_time >= time_slots(num_fgat_time) ) cycle scan_loop
 355          do ifgat=1,num_fgat_time
 356             if ( obs_time >= time_slots(ifgat-1) .and.  &
 357                obs_time  < time_slots(ifgat) ) exit
 358          end do
 359
 360       ! start fov_loop
 361          fov_loop:      do ilongitude=1, nlongitude
 362
 363             if ( sat_zenith(ilongitude,ilatitude) > 65.0 ) cycle fov_loop
 364
 365             num_ahi_file       = num_ahi_file + 1
 366             num_ahi_file_local = num_ahi_file_local + 1
 367             info%lat  =  vlatitude(ilongitude,ilatitude)
 368             info%lon  =  vlongitude(ilongitude,ilatitude)
 369
 370             call da_llxy (info, loc, outside, outside_all)
 371             if (outside_all) cycle fov_loop
 372
 373             num_ahi_global       = num_ahi_global + 1
 374             num_ahi_global_local = num_ahi_global_local + 1
 375             ptotal(ifgat) = ptotal(ifgat) + 1
 376             if (outside) cycle fov_loop   ! No good for this PE
 377
 378             num_ahi_local       = num_ahi_local + 1
 379             num_ahi_local_local = num_ahi_local_local + 1
 380             write(unit=info%date_char, &
 381             fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 382                idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
 383                ':', idate5(5), ':', idate5(6)
 384             info%elv = 0.0
 385
 386 ! 3.0  Make Thinning
 387 ! Map obs to thinning grid
 388 !-------------------------------------------------------------------
 389             if (thinning) then
 390                dlat_earth = info%lat !degree
 391                dlon_earth = info%lon
 392                if (dlon_earth<zero)  dlon_earth = dlon_earth+r360
 393                if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 394                dlat_earth = dlat_earth*deg2rad !radian
 395                dlon_earth = dlon_earth*deg2rad
 396                crit = 1.
 397                call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 398                if (.not. iuse) then
 399                   num_ahi_thinned = num_ahi_thinned+1
 400                   cycle fov_loop
 401                end if
 402             end if
 403
 404             num_ahi_used = num_ahi_used + 1
 405             data_all = missing_r
 406
 407             do k=1,nchan
 408                tb = tbb(ilongitude,ilatitude,k)
 409                if( tb < tbmin .or. tb > tbmax ) tb = missing_r
 410                data_all(k)= tb
 411             enddo
 412
 413 ! 4.0 assign information to sequential radiance structure
 414 !--------------------------------------------------------------------------
 415             allocate ( p % tb_inv (1:nchan ))
 416             p%info             = info
 417             p%loc              = loc
 418             p%landsea_mask     = 1
 419             p%scanpos          = ilongitude !nint(sat_zenith(ilongitude,ilatitude))+1.001_r_kind !
 420             p%satzen           = sat_zenith(ilongitude,ilatitude)
 421             p%satazi           = 0
 422             p%solzen           = 0
 423             p%solazi           = 0
 424             p%tb_inv(1:nchan)  = data_all(1:nchan)
 425             p%sensor_index     = inst
 426             p%ifgat            = ifgat
 427                         p%cloudflag        = cloud_mask(ilongitude,ilatitude)
 428
 429             allocate (p%next)   ! add next data
 430             p => p%next
 431             nullify (p%next)
 432          end do fov_loop
 433       end do scan_loop
 434
 435       write(stdout,fmt='(3a,i10)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_file    : ',num_ahi_file_local
 436       write(stdout,fmt='(3a,i10)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_global  : ',num_ahi_global_local
 437       write(stdout,fmt='(3a,i10)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_local   : ',num_ahi_local_local
 438    end do infile_loop
 439
 440    deallocate(data_all) ! Deallocate data arrays
 441    !deallocate(cloudflag)
 442    deallocate(vlatitude)
 443    deallocate(vlongitude)
 444    deallocate(tbb)
 445    deallocate(sat_zenith)
 446    if( got_clp_file ) deallocate(cloud_mask)
 447
 448    if (thinning .and. num_ahi_global > 0 ) then
 449 #ifdef DM_PARALLEL
 450    ! Get minimum crit and associated processor index.
 451       j = 0
 452       do ifgat = 1, num_fgat_time
 453          j = j + thinning_grid(inst,ifgat)%itxmax
 454       end do
 455
 456       allocate ( in  (j) )
 457       allocate ( out (j) )
 458       j = 0
 459       do ifgat = 1, num_fgat_time
 460          do i = 1, thinning_grid(inst,ifgat)%itxmax
 461             j = j + 1
 462             in(j) = thinning_grid(inst,ifgat)%score_crit(i)
 463          end do
 464       end do
 465       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 466
 467       call wrf_dm_bcast_real (out, j)
 468
 469       j = 0
 470       do ifgat = 1, num_fgat_time
 471          do i = 1, thinning_grid(inst,ifgat)%itxmax
 472             j = j + 1
 473             if ( ABS(out(j)-thinning_grid(inst,ifgat)%score_crit(i)) > 1.0E-10 ) &
 474             thinning_grid(inst,ifgat)%ibest_obs(i) = 0
 475          end do
 476       end do
 477
 478       deallocate( in  )
 479       deallocate( out )
 480
 481 #endif
 482
 483    ! Delete the nodes which being thinning out
 484       p => head
 485       prev => head
 486       head_found = .false.
 487       num_ahi_used_tmp = num_ahi_used
 488       do j = 1, num_ahi_used_tmp
 489          n = p%sensor_index
 490          ifgat = p%ifgat
 491          found = .false.
 492
 493          do i = 1, thinning_grid(n,ifgat)%itxmax
 494             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 495                found = .true.
 496                exit
 497             end if
 498          end do
 499
 500       ! free current data
 501          if ( .not. found ) then
 502
 503                         current => p
 504             p => p%next
 505
 506             if ( head_found ) then
 507                prev%next => p
 508             else
 509                head => p
 510                prev => p
 511             end if
 512
 513             deallocate ( current % tb_inv )
 514             deallocate ( current )
 515
 516             num_ahi_thinned = num_ahi_thinned + 1
 517             num_ahi_used = num_ahi_used - 1
 518             continue
 519          end if
 520
 521          if ( found .and. head_found ) then
 522             prev => p
 523             p => p%next
 524             continue
 525          end if
 526
 527          if ( found .and. .not. head_found ) then
 528             head_found = .true.
 529             head => p
 530             prev => p
 531             p => p%next
 532          end if
 533
 534       end do
 535    end if  ! End of thinning
 536
 537    iv%total_rad_pixel   = iv%total_rad_pixel + num_ahi_used
 538    iv%total_rad_channel = iv%total_rad_channel + num_ahi_used*nchan
 539
 540    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_ahi_used
 541    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_ahi_global
 542
 543    do i = 1, num_fgat_time
 544       ptotal(i) = ptotal(i) + ptotal(i-1)
 545       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 546    end do
 547    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 548       write(unit=message(1),fmt='(A,I10,A,I10)') &
 549           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 550       call da_warning(__FILE__,__LINE__,message(1:1))
 551    endif
 552
 553    write(unit=stdout,fmt='(a)') 'AHI data counts: '
 554    write(stdout,fmt='(a,i10)') ' In file: ',num_ahi_file
 555    write(stdout,fmt='(a,i10)') ' Global : ',num_ahi_global
 556    write(stdout,fmt='(a,i10)') ' Local  : ',num_ahi_local
 557    write(stdout,fmt='(a,i10)') ' Used   : ',num_ahi_used
 558    write(stdout,fmt='(a,i10)') ' Thinned: ',num_ahi_thinned
 559
 560 !  5.0 allocate innovation radiance structure
 561 !----------------------------------------------------------------
 562
 563    if (num_ahi_used > 0) then
 564       iv%instid(inst)%num_rad  = num_ahi_used
 565       iv%instid(inst)%info%nlocal = num_ahi_used
 566       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 567          'Allocating space for radiance innov structure', &
 568          inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
 569       call da_allocate_rad_iv (inst, nchan, iv)
 570    end if
 571
 572 !  6.0 assign sequential structure to innovation structure
 573 !-------------------------------------------------------------
 574    p => head
 575
 576    do n = 1, num_ahi_used
 577       i = p%sensor_index
 578       call da_initialize_rad_iv (i, n, iv, p)
 579       current => p
 580       p => p%next
 581    ! free current data
 582       deallocate ( current % tb_inv )
 583       deallocate ( current )
 584    end do
 585    deallocate ( p )
 586    deallocate (ptotal)
 587
 588    if (trace_use) call da_trace_exit("da_read_obs_netcdf4ahi_geocat")
 589
 590 end subroutine da_read_obs_netcdf4ahi_geocat