var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc

   1 subroutine da_read_obs_netcdf4ahi_jaxa (iv, infile_tb, infile_tbp)
   2    !--------------------------------------------------------
   3    !  Purpose: read in JAXA AHI Level-1 data in NETCDF4 format for different times
   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    !  HISTORY: 2020/03/01 - Add clear sky cloud detection procedures   Dongmei Xu, NUIST, NCAR/MMM
  14    !  To be devoloped: 1.time information; 2.dimension sequence
  15    !------------------------------------------------------------------------------
  16
  17    use netcdf
  18
  19    implicit none
  20
  21    character(len=*), intent(in)    :: infile_tb, infile_tbp
  22    type(iv_type),    intent(inout) :: iv
  23
  24 ! fixed parameter values
  25    integer,parameter::time_dims=6       ! Time dimension
  26    integer,parameter::nfile_max = 8     ! each netcdf file contains
  27    real,parameter::scale_factor_tb=0.01          ! Maximum allowed NumberOfScans
  28    real,parameter::add_offset_tb=273.15          ! low  resolution pixel width
  29
  30 ! interface variable
  31    integer iret, rcode, ncid                      ! return status
  32
  33 ! array data
  34    real(4), allocatable :: vlatitude(:)  !  value for latitude
  35    real(4), allocatable :: vlongitude(:) !  value for longitude
  36
  37    real(4), allocatable    :: tbb(:,:,:)  ! tb for band 7-16
  38    real(4), allocatable    :: sat_zenith(:,:),sat_azi(:,:),sol_azi(:,:),sol_zenith(:,:)
  39
  40    real(r_kind),parameter  :: tbmin  = 50._r_kind
  41    real(r_kind),parameter  :: tbmax  = 550._r_kind
  42
  43    real(kind=8)                   :: obs_time
  44    type (datalink_type),pointer   :: head, p, current, prev
  45    type(info_type)                :: info
  46    type(model_loc_type)           :: loc
  47
  48    integer(i_kind)    :: idate5(6)
  49    character(len=80) :: filename,str1,str2
  50
  51    integer(i_kind)   :: inst,platform_id,satellite_id,sensor_id
  52    real(r_kind)      :: tb, crit
  53    integer(i_kind)   :: ifgat, iout, iobs
  54    logical           :: outside, outside_all, iuse
  55
  56    integer           :: i,j,k,l,m,n, ifile, landsea_mask
  57
  58    logical           :: found, head_found, head_allocated
  59
  60 ! Other work variables
  61    real(r_kind)     :: dlon_earth,dlat_earth
  62    integer(i_kind)  :: num_ahi_local, num_ahi_global, num_ahi_used, num_ahi_thinned
  63    integer(i_kind)  :: num_ahi_used_tmp, num_ahi_file
  64    integer(i_kind)  :: num_ahi_local_local, num_ahi_global_local, num_ahi_file_local
  65    integer(i_kind)  :: itx, itt
  66    character(80)    :: filename1,filename2
  67    integer          :: nchan,nlongitude,nlatitude,ilongitude,ilatitude,ichannels
  68    integer          :: lonstart,latstart,ahi_info_unit
  69    integer          :: LatDimID,LonDimID
  70    integer          :: latid,lonid,tbb_id,sazid,cltyid,ter_id
  71    integer          :: nfile
  72    character(80)    :: fname_tb(nfile_max),fname_tbp(nfile_max)
  73    logical          :: fexist, got_clp_file
  74
  75 ! for qc stuff
  76    integer, parameter :: ch7  = 1
  77    integer, parameter :: ch10 = 4
  78    integer, parameter :: ch13 = 7
  79    integer, parameter :: ch14 = 8
  80    integer, parameter :: ch15 = 9
  81    integer            :: superob_width ! Must be ≥ 0
  82    integer            :: first_boxcenter, last_boxcenter_x, last_boxcenter_y, box_bottom, box_upper, box_left, box_right
  83    integer            :: isup, jsup, ixsup, iysup, ibox, jbox, nkeep
  84    real(4), allocatable    :: tbbp(:,:,:)  ! tb for band 7-16
  85    real(r_kind), allocatable       :: tbb_temp1(:,:), tbb_temp2(:,:)
  86    real(r_kind), allocatable       :: rtct(:,:), rtct2(:,:),rtmt(:,:)
  87    real(r_kind), allocatable       :: tempir(:,:),cwvt1(:,:)
  88    real(r_kind), allocatable    :: ter_sat(:,:),SDob_10(:,:),SDob_13(:,:),SDob_14(:,:)
  89    real(r_kind), allocatable    :: cloud_mask(:,:) ! For reading L2 cloud information
  90    character(200)    :: filenameStatic
  91 ! Allocatable arrays
  92    integer(i_kind),allocatable  :: ptotal(:)
  93    real,allocatable             :: in(:), out(:)
  94    real(r_kind)                 :: temp1 = 0.0
  95    character(len=80) tbb_name(10)
  96    data tbb_name/'tbb_07','tbb_08','tbb_09','tbb_10','tbb_11', &
  97                  'tbb_12','tbb_13','tbb_14','tbb_15','tbb_16'/
  98
  99    if (trace_use) call da_trace_entry("da_read_obs_netcdf4ahi_jaxa")
 100
 101 !  0.0  Initialize variables
 102 !-----------------------------------
 103    head_allocated = .false.
 104    platform_id  = 31  ! Table-2 Col 1 corresponding to 'himawari'
 105    satellite_id = 8   ! Table-2 Col 3
 106    sensor_id    = 56  ! Table-3 Col 2 corresponding to 'ahi'
 107
 108    allocate(ptotal(0:num_fgat_time))
 109    ptotal(0:num_fgat_time) = 0
 110    iobs = 0                 ! for thinning, argument is inout
 111    num_ahi_file    = 0
 112    num_ahi_local   = 0
 113    num_ahi_global  = 0
 114    num_ahi_used    = 0
 115    num_ahi_thinned = 0
 116
 117    inst = 0
 118    do i = 1, rtminit_nsensor
 119       if (platform_id  == rtminit_platform(i) &
 120           .and. satellite_id == rtminit_satid(i)    &
 121           .and. sensor_id    == rtminit_sensor(i)) then
 122          inst = i
 123          exit
 124       end if
 125    end do
 126    if (inst == 0) then
 127       call da_warning(__FILE__,__LINE__, &
 128           (/"The combination of Satellite_Id and Sensor_Id for AHI is not found"/))
 129       if (trace_use) call da_trace_exit("da_read_obs_netcdf4ahi_jaxa")
 130       return
 131    end if
 132
 133    nchan = iv%instid(inst)%nchan
 134    write(unit=stdout,fmt=*)'AHI nchan: ',nchan
 135    superob_width = 2*ahi_superob_halfwidth+1
 136
 137 ! 1.0 Assign file names and prepare to read ahi files
 138 !-------------------------------------------------------------------------
 139    nfile       = 0  !initialize
 140    fname_tb(:) = '' !initialize
 141    fname_tbp(:) = '' ! initialize
 142
 143    ! first check if ahi nc file is available
 144    filename1 = trim(infile_tb)
 145    filename2 = trim(infile_tbp)
 146    inquire (file=filename1, exist=fexist)
 147    if ( fexist ) then
 148       nfile = 1
 149       fname_tb(nfile)  = filename1
 150           fname_tbp(nfile) = filename2
 151    else
 152       ! check if netcdf4 files are available for multiple input files
 153       ! here 0x is the input file sequence number
 154       ! do not confuse it with fgat time slot index
 155       do i = 1, nfile_max
 156          write(filename1,fmt='(A,A,I2.2,A)') trim(infile_tb),'-',i
 157                  write(filename2,fmt='(A,A,I2.2,A)') trim(infile_tbp),'-',i
 158          inquire (file=filename1, exist=fexist)
 159          if ( fexist ) then
 160             nfile = nfile + 1
 161             fname_tb(nfile)  = filename1
 162                         fname_tbp(nfile) = filename2
 163          else
 164             exit
 165          end if
 166       end do
 167    end if
 168
 169    if ( nfile == 0 ) then
 170       call da_warning(__FILE__,__LINE__, &
 171          (/"No valid AHI file found."/))
 172       if (trace_use) call da_trace_exit("da_read_obs_netcdf4ahi_jaxa")
 173       return
 174    end if
 175
 176   ! Added to get area information from the file, in case you don't
 177   !  want to read the entire area, which is large and memory intensive
 178   !open the data area info file
 179    call da_get_unit(ahi_info_unit)
 180    open(unit=ahi_info_unit,file='ahi_info',status='old',iostat=iret)
 181    if(iret /= 0)then
 182        call da_warning(__FILE__,__LINE__,(/"area_info file read error"/))
 183    endif
 184   !read date information
 185    read(ahi_info_unit,*)
 186    read(ahi_info_unit,*)
 187    read(ahi_info_unit,*)
 188    read(ahi_info_unit,*)
 189    read(ahi_info_unit,*)
 190    read(ahi_info_unit,*) lonstart,latstart,nlongitude,nlatitude
 191    close(ahi_info_unit)
 192    call da_free_unit(ahi_info_unit)
 193
 194    infile_loop:  do ifile = 1, nfile
 195       num_ahi_file_local    = 0
 196       num_ahi_local_local   = 0
 197       num_ahi_global_local  = 0
 198
 199    ! open NETCDF4 L1 file for read for the current time
 200       iret = nf90_open(fname_tb(ifile), NF90_NOWRITE, ncid)
 201       if(iret /= 0)then
 202          call da_warning(__FILE__,__LINE__, &
 203               (/"Cannot open NETCDF4 file "//trim(fname_tb(ifile))/))
 204          cycle infile_loop
 205       endif
 206
 207       ! read dimensions:        latitude and longitude
 208       ! If any value in the file is < 0, that means to ignore them
 209       !     and get the number of lats/lons from file, and use everything
 210       if ( nlongitude < 0 .or. nlatitude < 0 .or. lonstart < 0 .or. latstart < 0 ) then
 211          iret = nf90_inq_dimid(ncid, "latitude", LatDimID)
 212          iret = nf90_inquire_dimension(ncid, LatDimID, len=nlatitude)
 213
 214          iret = nf90_inq_dimid(ncid, "longitude", LonDimID)
 215          iret = nf90_inquire_dimension(ncid, LonDimID, len=nlongitude)
 216
 217          lonstart = 1
 218          latstart = 1
 219       endif
 220
 221       write(unit=stdout,fmt=*)'nlongitude,nlatitude: ',nlongitude,nlatitude
 222       write(unit=stdout,fmt=*)'lonstart,latstart: ',lonstart,latstart
 223
 224  ! read array: time
 225       iret = nf90_get_att(ncid, nf90_global, "id", filename)
 226       if(iret /= 0)then
 227          call da_warning(__FILE__,__LINE__, &
 228              (/"NETCDF4 read error for: observation date"/))
 229       endif
 230           read(filename,"(A7,I4,I2,I2,A1,I2,I2)") str1,idate5(1),idate5(2),idate5(3), &
 231                                                                       str2,idate5(4),idate5(5)
 232           idate5(6)=00
 233           write(unit=stdout,fmt=*)'observation date: ', idate5
 234
 235    ! read array: latlon
 236    ! read lat
 237           iret = nf90_inq_varid(ncid, 'latitude', latid)
 238           allocate( vlatitude(nlatitude))
 239       iret = nf90_get_var(ncid, latid, vlatitude,start=(/latstart/),count=(/nlatitude/))
 240       if(iret /= 0)then
 241          call da_warning(__FILE__,__LINE__, &
 242              (/"NETCDF4 read error for: Latitude of Observation Point"/))
 243       endif
 244    ! read lon
 245           iret = nf90_inq_varid(ncid, 'longitude', lonid)
 246           allocate( vlongitude(nlongitude))
 247       iret = nf90_get_var(ncid, lonid, vlongitude,start=(/lonstart/),count=(/nlongitude/))
 248       if(iret /= 0)then
 249           call da_warning(__FILE__,__LINE__, &
 250               (/"NETCDF4 read error for: Longitude of Observation Point"/))
 251           call da_trace_exit("da_read_obs_netcdf4ahi_jaxa")
 252       endif
 253    ! sample display
 254       write(unit=stdout,fmt=*)'latitude,longitude(pixel=1,scan=1): ',vlatitude(1),vlongitude(1)
 255
 256    ! read array: tb for band 7-16
 257    ! read
 258       allocate( tbb(nlongitude,nlatitude,nchan))
 259       do k=1,nchan
 260                 iret = nf90_inq_varid(ncid, tbb_name(k), tbb_id)
 261                 iret = nf90_get_var(ncid, tbb_id, tbb(:,:,k),start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/))
 262                 if(iret /= 0)then
 263                   call da_warning(__FILE__,__LINE__, &
 264                            (/"NETCDF4 read error for: Brightness Temperature"/))
 265                 endif
 266                 do j=1,nlatitude
 267                   do i=1,nlongitude
 268                           tbb(i,j,k)=tbb(i,j,k) * scale_factor_tb + add_offset_tb
 269               if( tbb(i,j,k) < tbmin .or. tbb(i,j,k) > tbmax ) tbb(i,j,k) = missing_r
 270                   end do
 271                 end do
 272            ! sample display
 273                    write(unit=stdout,fmt=*)&
 274                         'tbb(pixel=1,scan=1,chan=',k,'): ', tbb(1,1,k)
 275       end do
 276
 277    ! read array: satellite zenith angle
 278    ! read
 279           iret = nf90_inq_varid(ncid, 'SAZ', sazid)
 280           allocate( sat_zenith(nlongitude,nlatitude))
 281           iret = nf90_get_var(ncid, sazid, sat_zenith,start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/))
 282       if(iret /= 0)then
 283          call da_warning(__FILE__,__LINE__, &
 284             (/"NETCDF4 read error for: satellite zenith angle"/))
 285       endif
 286                 do j=1,nlatitude
 287                         do i=1,nlongitude
 288                                 sat_zenith(i,j)=sat_zenith(i,j) * scale_factor_tb
 289                         end do
 290                 end do
 291    ! sample display
 292       write(unit=stdout,fmt=*)&
 293          'satellite zenith angle(pixel=1,scan=1): ',sat_zenith(1,1)
 294
 295
 296    ! read array: satellite azimuth angle
 297    ! read
 298           iret = nf90_inq_varid(ncid, 'SAA', sazid)
 299           allocate( sat_azi(nlongitude,nlatitude))
 300           iret = nf90_get_var(ncid, sazid, sat_azi,start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/))
 301       if(iret /= 0)then
 302          call da_warning(__FILE__,__LINE__, &
 303             (/"NETCDF4 read error for: satellite azimuth angle"/))
 304       endif
 305                 do j=1,nlatitude
 306                         do i=1,nlongitude
 307                                 sat_azi(i,j)=sat_azi(i,j) * scale_factor_tb
 308                         end do
 309                 end do
 310    ! sample display
 311       write(unit=stdout,fmt=*)&
 312          'satellite azimuth angle(pixel=1,scan=1): ',sat_azi(1,1)
 313
 314    ! read array: solar zenith angle
 315    ! read
 316           iret = nf90_inq_varid(ncid, 'SOZ', sazid)
 317           allocate( sol_zenith(nlongitude,nlatitude))
 318           iret = nf90_get_var(ncid, sazid, sol_zenith,start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/))
 319       if(iret /= 0)then
 320          call da_warning(__FILE__,__LINE__, &
 321             (/"NETCDF4 read error for: satellite zenith angle"/))
 322       endif
 323                 do j=1,nlatitude
 324                         do i=1,nlongitude
 325                                 sol_zenith(i,j)=sol_zenith(i,j) * scale_factor_tb
 326                         end do
 327                 end do
 328    ! sample display
 329       write(unit=stdout,fmt=*)&
 330          'solar zenith angle(pixel=1,scan=1): ',sol_zenith(1,1)
 331
 332
 333    ! read array: solar azimuth angle
 334    ! read
 335           iret = nf90_inq_varid(ncid, 'SOA', sazid)
 336           allocate( sol_azi(nlongitude,nlatitude))
 337           iret = nf90_get_var(ncid, sazid, sol_azi,start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/))
 338       if(iret /= 0)then
 339          call da_warning(__FILE__,__LINE__, &
 340             (/"NETCDF4 read error for: satellite azimuth angle"/))
 341       endif
 342                 do j=1,nlatitude
 343                         do i=1,nlongitude
 344                                 sol_azi(i,j)=sol_azi(i,j) * scale_factor_tb
 345                         end do
 346                 end do
 347    ! sample display
 348       write(unit=stdout,fmt=*)&
 349          'solar azimuth angle(pixel=1,scan=1): ',sol_azi(1,1)
 350    ! close infile_tb file
 351       iret = nf90_close(ncid)
 352
 353       !  Read the level-2 cloud information file from JAXA
 354       !  This must be on the same grid as everything else
 355       allocate(cloud_mask(nlongitude,nlatitude))
 356       cloud_mask = -9999.0
 357       inquire(file='L2AHICLP', exist=fexist)
 358       if ( fexist ) then
 359          iret = nf90_open('L2AHICLP', nf90_NOWRITE, ncid)
 360          if(iret /= 0)then
 361             call da_warning(__FILE__,__LINE__,(/"NETCDF4 read error for: CLTYPE data"/))
 362          endif
 363
 364          iret = nf90_inq_varid(ncid, "CLTYPE", cltyid)
 365          iret = nf90_get_var(ncid,cltyid,cloud_mask,start=(/lonstart,latstart/), count=(/nlongitude,nlatitude/))
 366          if(iret /= 0)then
 367             call da_warning(__FILE__,__LINE__,(/"NETCDF4 read error for: CLTYPE data"/))
 368          endif
 369          iret = nf90_close(ncid)
 370       end if
 371
 372       ! Note, these are only deallocated if use_clddet_zz = true
 373       !  so their allocation should probably be moved inside the conditional...
 374       allocate( tbbp(nlongitude,nlatitude,nchan))
 375       allocate( ter_sat(nlongitude,nlatitude))
 376       allocate( SDob_10(nlongitude,nlatitude))
 377       allocate( SDob_13(nlongitude,nlatitude))
 378       allocate( SDob_14(nlongitude,nlatitude))
 379       allocate( rtct(nlongitude,nlatitude))
 380       allocate( rtct2(nlongitude,nlatitude))
 381       allocate( rtmt(nlongitude,nlatitude))
 382       allocate( tempir(nlongitude,nlatitude))
 383       allocate( cwvt1(nlongitude,nlatitude))
 384       allocate( tbb_temp1(nlongitude,nlatitude))
 385       allocate( tbb_temp2(nlongitude,nlatitude))
 386
 387    if (use_clddet_zz) then
 388 ! --------------------------
 389 ! 0-call read static data: terrain and landmask   ?????
 390 ! --------------------------
 391    ! open NETCDF4 L1 file for read for the previous time
 392       iret = nf90_open(fname_tbp(ifile), NF90_NOWRITE, ncid)
 393       if(iret /= 0)then
 394          call da_warning(__FILE__,__LINE__, &
 395               (/"Cannot open NETCDF4 file "//trim(fname_tbp(ifile))/))
 396          cycle infile_loop
 397       endif
 398
 399    ! read array: tb for band 7-16
 400    ! read
 401
 402       do k=1,nchan
 403                 iret = nf90_inq_varid(ncid, tbb_name(k), tbb_id)
 404                 iret = nf90_get_var(ncid, tbb_id, tbbp(:,:,k),start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/))
 405                 if(iret /= 0)then
 406                   call da_warning(__FILE__,__LINE__, &
 407                            (/"NETCDF4 read error for: Brightness Temperature"/))
 408                 endif
 409                 do j=1,nlatitude
 410                   do i=1,nlongitude
 411                           tbbp(i,j,k)=tbbp(i,j,k) * scale_factor_tb + add_offset_tb
 412               if( tbbp(i,j,k) < tbmin .or. tbbp(i,j,k) > tbmax ) tbbp(i,j,k) = missing_r
 413                   end do
 414                 end do
 415            ! sample display
 416                    write(unit=stdout,fmt=*)&
 417                         'tbbp(pixel=1,scan=1,chan=',k,'): ', tbbp(1,1,k)
 418       end do
 419
 420    ! close infile_tbp file
 421       iret = nf90_close(ncid)
 422
 423 ! --------------------------
 424 ! 1-call read static data: terrain and landmask   ?????
 425 ! --------------------------
 426
 427
 428       filenameStatic="static_ahi.nc"
 429       iret = nf90_open(filenameStatic, nf90_NOWRITE, ncid)
 430       if(iret /= 0)then
 431          print*,"Cannot open NETCDF4 file "//trim(filenameStatic)
 432       endif
 433       iret = nf90_inq_varid(ncid, "AhiTer", ter_id)
 434       iret = nf90_get_var(ncid, ter_id, ter_sat(:,:),start=(/lonstart,latstart/),count=(/nlongitude,nlatitude/))
 435       if(iret /= 0) print*,"NETCDF4 read error for: SatTer"
 436       where (ter_sat == -999.) ter_sat = missing_r
 437       iret = nf90_close(ncid)
 438     !  print *, 'ter_sat',maxval(ter_sat),      minval(ter_sat)
 439
 440 !!2- ZZ test0-SDob
 441
 442           tbb_temp1 = missing_r
 443       tbb_temp1 = tbb(:,:,ch13)   ! 10.4um
 444           SDob_13=missing_r
 445       call qc_SDob(nlongitude,nlatitude,tbb_temp1,SDob_13)
 446       !print *, 'SDob_13',maxval(SDob_13),      minval(SDob_13)
 447
 448           !!!!!!!!!! for ZZ test6
 449           tbb_temp1 = missing_r
 450       tbb_temp1 = tbb(:,:,ch10)   ! 7.3um
 451           SDob_10=missing_r
 452       call qc_SDob(nlongitude,nlatitude,tbb_temp1,SDob_10)
 453 !print *, 'SDob_10',maxval(SDob_10),    minval(SDob_10)
 454           tbb_temp1 = missing_r
 455       tbb_temp1 = tbb(:,:,ch14)   ! 11.2umum
 456           SDob_14=missing_r
 457       call qc_SDob(nlongitude,nlatitude,tbb_temp1,SDob_14)
 458 !print *, 'SDob_14',maxval(SDob_14),    minval(SDob_14)
 459
 460 !! ZZ test1- RTCT
 461           tbb_temp1 = missing_r
 462       tbb_temp1 = tbb(:,:,ch14)   ! 11.2um
 463       rtct=0.1
 464       call qc_RTCT(nlongitude,nlatitude,tbb_temp1,ter_sat,rtct,rtct2)
 465
 466 !! ZZ test5-RFMFT
 467
 468           tbb_temp1 = missing_r
 469           tbb_temp2 = missing_r
 470       tbb_temp1 = tbb(:,:,ch14);  tbb_temp2 = tbb(:,:,ch15)   ! 11.2um and 12.4um
 471           rtmt=missing_r
 472       call qc_rtmt(nlongitude,nlatitude,tbb_temp1,tbb_temp2,rtmt)
 473 !print *, 'rtmt',maxval(rtmt),  minval(rtmt)
 474
 475 !! ZZ test6-CIRH2O
 476           tbb_temp1 = missing_r
 477           tbb_temp2 = missing_r
 478           tbb_temp1 = tbb(:,:,ch10)
 479           tbb_temp2 = tbb(:,:,ch14)   ! 7.3um and 11.2um
 480           cwvt1=missing_r
 481       call qc_cwvt(nlongitude,nlatitude,tbb_temp1,tbb_temp2,cwvt1)
 482 !print *, 'cwvt1',maxval(cwvt1),        minval(cwvt1)
 483
 484
 485 !! ZZ test10- TEMPIR
 486           tbb_temp1 = missing_r
 487           tbb_temp2 = missing_r
 488       tbb_temp1 = tbbp(:,:,ch14);  tbb_temp2 = tbb(:,:,ch14)   ! 11.2um
 489       where(tbb_temp1 /=missing_r .and. tbb_temp2 /= missing_r) tempir = tbb_temp1-tbb_temp2
 490       where(tbb_temp1 ==missing_r .or.  tbb_temp2 == missing_r) tempir = missing_r
 491 !print *, 'tempir',maxval(tempir),      minval(tempir),nlongitude,nlatitude
 492 end if
 493
 494
 495
 496
 497 ! 2.0 Loop to read netcdf and assign information to a sequential structure
 498 !-------------------------------------------------------------------------
 499    ! Allocate arrays to hold data
 500       if ( .not. head_allocated ) then
 501          allocate (head)
 502          nullify  ( head % next )
 503          p => head
 504          head_allocated = .true.
 505       end if
 506    ! start scan_loop
 507
 508       first_boxcenter = ahi_superob_halfwidth + 1
 509       last_boxcenter_y = superob_width * (nlatitude / superob_width) - ahi_superob_halfwidth
 510       last_boxcenter_x = superob_width * (nlongitude / superob_width) - ahi_superob_halfwidth
 511
 512       scan_loop:     do ilatitude=first_boxcenter, nlatitude, superob_width
 513          call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
 514          if ( obs_time < time_slots(0) .or.  &
 515             obs_time >= time_slots(num_fgat_time) ) cycle scan_loop
 516          do ifgat=1,num_fgat_time
 517             if ( obs_time >= time_slots(ifgat-1) .and.  &
 518                obs_time  < time_slots(ifgat) ) exit
 519          end do
 520       ! start fov_loop
 521
 522          jbox = ilatitude/superob_width + 1
 523          if ( ahi_superob_halfwidth .gt. 0 ) then
 524             box_bottom = superob_width * (jbox-1) +1
 525             box_upper  = superob_width * jbox
 526          else
 527             box_bottom = superob_width * (jbox-1)
 528             box_upper  = superob_width * (jbox-1)
 529          end if
 530
 531          fov_loop:      do ilongitude=first_boxcenter, nlongitude, superob_width
 532
 533             if ( sat_zenith(ilongitude,ilatitude) > 65.0 ) cycle fov_loop
 534
 535             num_ahi_file       = num_ahi_file + 1
 536             num_ahi_file_local = num_ahi_file_local + 1
 537             info%lat  =  vlatitude(ilatitude)
 538             info%lon  =  vlongitude(ilongitude)
 539             call da_llxy (info, loc, outside, outside_all)
 540             if (outside_all) cycle fov_loop
 541             num_ahi_global       = num_ahi_global + 1
 542             num_ahi_global_local = num_ahi_global_local + 1
 543             ptotal(ifgat) = ptotal(ifgat) + 1
 544             if (outside) cycle fov_loop   ! No good for this PE
 545             num_ahi_local       = num_ahi_local + 1
 546             num_ahi_local_local = num_ahi_local_local + 1
 547             write(unit=info%date_char, &
 548             fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 549                idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
 550                ':', idate5(5), ':', idate5(6)
 551             info%elv = 0.0
 552
 553             ! 3.0  Make Thinning
 554             ! Map obs to thinning grid
 555             !-------------------------------------------------------------------
 556             if (thinning) then
 557                dlat_earth = info%lat !degree
 558                dlon_earth = info%lon
 559                if (dlon_earth<zero)  dlon_earth = dlon_earth+r360
 560                if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 561                dlat_earth = dlat_earth*deg2rad !radian
 562                dlon_earth = dlon_earth*deg2rad
 563                crit = 1.
 564                call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 565                if (.not. iuse) then
 566                   num_ahi_thinned = num_ahi_thinned+1
 567                   cycle fov_loop
 568                end if
 569             end if
 570             num_ahi_used = num_ahi_used + 1
 571
 572             ibox = ilongitude/superob_width + 1
 573             if ( ahi_superob_halfwidth .gt. 0 ) then
 574                box_left  = superob_width * (ibox-1) +1 ! will exceed nlatitude/nlongitude if ahi_superob_halfwidth = 0
 575                box_right = superob_width * ibox
 576             else
 577                box_left  = superob_width * (ibox-1)
 578                box_right = superob_width * (ibox-1)
 579             end if
 580
 581             ! Super-ob BT for this channel
 582
 583             allocate ( p % tb_inv (1:nchan) )
 584             p % tb_inv = missing_r
 585
 586             do k = 1, nchan
 587                nkeep = count(tbb(box_left:box_right,box_bottom:box_upper,k) > 0.0 )
 588                temp1 = sum ( tbb(box_left:box_right,box_bottom:box_upper,k), &
 589                              tbb(box_left:box_right,box_bottom:box_upper,k) > 0.0 )
 590
 591                if (ahi_superob_halfwidth .gt.0 .and. nkeep .gt. 0) then
 592                   tb = temp1 / real(nkeep,8)
 593                else
 594                   ! Extract single pixel BT and radiance value for this channel
 595                   tb = tbb(ilongitude,ilatitude,k)
 596                end if
 597
 598                ! Apply clear-sky bias correction if requested
 599                if ( ahi_apply_clrsky_bias ) then
 600                   tb = tb - satinfo(inst)%clearSkyBias(k)
 601                  !write(*,*)'channel, clearSkyBias = ',satinfo(inst)%ichan(k),satinfo(inst)%clearSkyBias(k)
 602                endif
 603
 604                if (tb > 0.0)   p % tb_inv(k) = tb
 605
 606              end do
 607
 608             ! Preprocessing for cloud detection including
 609             ! extracting Tb values from cloud QC buffer
 610
 611             if (use_clddet_zz) then
 612
 613                if (.not. allocated(p % superob)) then
 614                   allocate( p % superob(superob_width,superob_width) )
 615                end if
 616
 617                ! Loops over superob pixels
 618                do jsup = 1, superob_width
 619                do isup = 1, superob_width
 620                   iysup = ilatitude + jsup-1-ahi_superob_halfwidth
 621                   ixsup = ilongitude + isup-1-ahi_superob_halfwidth
 622                   if (ixsup < 1 .or. ixsup > nlongitude .or. iysup < 1 .or. iysup > nlatitude ) then
 623                      cycle
 624                   end if
 625                   if ( .not. allocated(p % superob(isup,jsup) % tb_obs ) ) then
 626                      allocate ( p % superob(isup,jsup) % tb_obs  (1:nchan,1) )
 627                      allocate ( p % superob(isup,jsup) % cld_qc(1) )
 628                   end if
 629                   p % superob(isup,jsup) % tb_obs(1:nchan,1)      = tbb( ixsup, iysup, 1:nchan )
 630                   p % superob(isup,jsup) % cld_qc(1) % TEMPIR        = tempir(ixsup, iysup)
 631                   p % superob(isup,jsup) % cld_qc(1) % RTCT          = rtct(ixsup, iysup)
 632                   p % superob(isup,jsup) % cld_qc(1) % RFMFT         = rtmt(ixsup, iysup)
 633                   p % superob(isup,jsup) % cld_qc(1) % CIRH2O        = cwvt1(ixsup, iysup)
 634                   p % superob(isup,jsup) % cld_qc(1) % tb_stddev_10  = SDob_10(ixsup, iysup)
 635                   p % superob(isup,jsup) % cld_qc(1) % tb_stddev_13  = SDob_13(ixsup, iysup)
 636                   p % superob(isup,jsup) % cld_qc(1) % tb_stddev_14  = SDob_14(ixsup, iysup)
 637                   p % superob(isup,jsup) % cld_qc(1) % terr_hgt      = ter_sat(ixsup, iysup)
 638
 639                end do !isup
 640                end do !jsup
 641
 642             end if
 643             ! 4.0 assign information to sequential radiance structure
 644             !--------------------------------------------------------------------------
 645             p%info             = info
 646             p%loc              = loc
 647             p%landsea_mask     = 1
 648             p%scanpos          = ilongitude !nint(sat_zenith(ilongitude,ilatitude))+1.001_r_kind !
 649             p%satzen           = sat_zenith(ilongitude,ilatitude)
 650             p%satazi           = sat_azi(ilongitude,ilatitude)
 651             p%solzen           = sol_zenith(ilongitude,ilatitude)
 652             p%solazi           = sol_azi(ilongitude,ilatitude)
 653             p%sensor_index     = inst
 654             p%ifgat            = ifgat
 655             p%cloudflag        = int(cloud_mask(ilongitude,ilatitude)) ! AHI L-2 cloud flag
 656             allocate (p%next)   ! add next data
 657             p => p%next
 658             nullify (p%next)
 659          end do fov_loop
 660       end do scan_loop
 661       write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_file    : ',num_ahi_file_local
 662       write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_global  : ',num_ahi_global_local
 663       write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_ahi_local   : ',num_ahi_local_local
 664
 665
 666    deallocate(vlatitude)
 667    deallocate(vlongitude)
 668    deallocate(tbb)
 669    deallocate(sat_zenith)
 670    deallocate(sol_zenith)
 671    deallocate(sat_azi)
 672    deallocate(sol_azi)
 673    deallocate(cloud_mask)
 674
 675    if( use_clddet_zz ) deallocate(tbbp)
 676    if( use_clddet_zz ) deallocate(ter_sat)
 677    if( use_clddet_zz ) deallocate(SDob_10)
 678    if( use_clddet_zz ) deallocate(SDob_13)
 679    if( use_clddet_zz ) deallocate(SDob_14)
 680    if( use_clddet_zz ) deallocate(rtct)
 681    if( use_clddet_zz ) deallocate(rtct2)
 682    if( use_clddet_zz ) deallocate(rtmt)
 683    if( use_clddet_zz ) deallocate(tempir)
 684    if( use_clddet_zz ) deallocate(cwvt1)
 685    if( use_clddet_zz ) deallocate(tbb_temp1)
 686    if( use_clddet_zz ) deallocate(tbb_temp2)
 687
 688    end do infile_loop
 689
 690    if (thinning .and. num_ahi_global > 0 ) then
 691 #ifdef DM_PARALLEL
 692    ! Get minimum crit and associated processor index.
 693       j = 0
 694       do ifgat = 1, num_fgat_time
 695          j = j + thinning_grid(inst,ifgat)%itxmax
 696       end do
 697
 698       allocate ( in  (j) )
 699       allocate ( out (j) )
 700       j = 0
 701       do ifgat = 1, num_fgat_time
 702          do i = 1, thinning_grid(inst,ifgat)%itxmax
 703             j = j + 1
 704             in(j) = thinning_grid(inst,ifgat)%score_crit(i)
 705          end do
 706       end do
 707       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 708
 709       call wrf_dm_bcast_real (out, j)
 710
 711       j = 0
 712       do ifgat = 1, num_fgat_time
 713          do i = 1, thinning_grid(inst,ifgat)%itxmax
 714             j = j + 1
 715             if ( ABS(out(j)-thinning_grid(inst,ifgat)%score_crit(i)) > 1.0E-10 ) &
 716             thinning_grid(inst,ifgat)%ibest_obs(i) = 0
 717          end do
 718       end do
 719
 720       deallocate( in  )
 721       deallocate( out )
 722
 723 #endif
 724
 725    ! Delete the nodes which being thinning out
 726       p => head
 727       prev => head
 728       head_found = .false.
 729       num_ahi_used_tmp = num_ahi_used
 730       do j = 1, num_ahi_used_tmp
 731          n = p%sensor_index
 732          ifgat = p%ifgat
 733          found = .false.
 734
 735          do i = 1, thinning_grid(n,ifgat)%itxmax
 736             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 737                found = .true.
 738                exit
 739             end if
 740          end do
 741
 742       ! free current data
 743          if ( .not. found ) then
 744             current => p
 745             p => p%next
 746             if ( head_found ) then
 747                prev%next => p
 748             else
 749                head => p
 750                prev => p
 751             end if
 752             deallocate ( current % tb_inv )
 753             if ( allocated ( current % superob ) ) then
 754                do jsup = 1, superob_width
 755                do isup = 1, superob_width
 756                   deallocate ( current % superob(isup,jsup) % tb_obs )
 757                   deallocate ( current % superob(isup,jsup) % cld_qc )
 758                end do
 759                end do
 760                deallocate ( current % superob )
 761             end if
 762             deallocate ( current )
 763             num_ahi_thinned = num_ahi_thinned + 1
 764             num_ahi_used = num_ahi_used - 1
 765             continue
 766          end if
 767
 768          if ( found .and. head_found ) then
 769             prev => p
 770             p => p%next
 771             continue
 772          end if
 773
 774          if ( found .and. .not. head_found ) then
 775             head_found = .true.
 776             head => p
 777             prev => p
 778             p => p%next
 779          end if
 780
 781       end do
 782
 783    end if  ! End of thinning
 784    iv%total_rad_pixel   = iv%total_rad_pixel + num_ahi_used
 785    iv%total_rad_channel = iv%total_rad_channel + num_ahi_used*nchan
 786
 787    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_ahi_used
 788    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_ahi_global
 789
 790    do i = 1, num_fgat_time
 791       ptotal(i) = ptotal(i) + ptotal(i-1)
 792       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 793    end do
 794    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 795       write(unit=message(1),fmt='(A,I10,A,I10)') &
 796           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 797       call da_warning(__FILE__,__LINE__,message(1:1))
 798    endif
 799
 800    write(unit=stdout,fmt='(a)') 'AHI data counts: '
 801    write(stdout,fmt='(a,i7)') ' In file: ',num_ahi_file
 802    write(stdout,fmt='(a,i7)') ' Global : ',num_ahi_global
 803    write(stdout,fmt='(a,i7)') ' Local  : ',num_ahi_local
 804    write(stdout,fmt='(a,i7)') ' Used   : ',num_ahi_used
 805    write(stdout,fmt='(a,i7)') ' Thinned: ',num_ahi_thinned
 806
 807 !  5.0 allocate innovation radiance structure
 808 !----------------------------------------------------------------
 809
 810    if (num_ahi_used > 0) then
 811       iv%instid(inst)%num_rad  = num_ahi_used
 812       iv%instid(inst)%info%nlocal = num_ahi_used
 813       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 814          'Allocating space for radiance innov structure', &
 815          inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
 816       call da_allocate_rad_iv (inst, nchan, iv)
 817    end if
 818
 819 !  6.0 assign sequential structure to innovation structure
 820 !-------------------------------------------------------------
 821    p => head
 822
 823    do n = 1, num_ahi_used
 824       i = p%sensor_index
 825       call da_initialize_rad_iv (i, n, iv, p)
 826       current => p
 827       p => p%next
 828    ! free current data
 829       deallocate ( current % tb_inv )
 830       if ( allocated ( current % superob ) ) then
 831          do jsup = 1, superob_width
 832          do isup = 1, superob_width
 833             deallocate ( current % superob(isup,jsup) % tb_obs )
 834             deallocate ( current % superob(isup,jsup) % cld_qc )
 835          end do
 836          end do
 837          deallocate ( current % superob )
 838       end if
 839       deallocate ( current )
 840    end do
 841    deallocate ( p )
 842    deallocate (ptotal)
 843
 844    if (trace_use) call da_trace_exit("da_read_obs_netcdf4ahi_jaxa")
 845
 846 end subroutine da_read_obs_netcdf4ahi_jaxa
 847