var/da/da_radiance/da_read_obs_bufrssmis.inc

   1 subroutine da_read_obs_bufrssmis (obstype,iv,infile)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: read in NCEP bufr SSM/IS data to innovation structure
   5    !
   6    !   METHOD: use F90 sequential data structure to avoid reading 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
  13
  14
  15    use da_control
  16
  17    implicit none
  18
  19    character(5) ,     intent (in)    :: obstype    ! ssmis
  20    character(20),     intent (in)    :: infile     ! ssmis.bufr
  21    type (iv_type),    intent (inout) :: iv
  22
  23 #ifdef BUFR
  24
  25    integer(i_kind), parameter :: bufsat_dmsp16 = 249  ! DMSP16 BUFR identifier
  26    integer(i_kind), parameter :: bufsat_dmsp17 = 285  ! DMSP17 BUFR identifier
  27    integer(i_kind), parameter :: bufsat_dmsp18 = 286  ! DMSP18 BUFR identifier
  28    integer(i_kind), parameter :: n1bhdr = 15
  29    integer(i_kind), parameter :: maxchanl = 24
  30    real(r_kind),    parameter :: tbmin = 70.0_r_kind
  31    real(r_kind),    parameter :: tbmax = 320.0_r_kind
  32
  33    character(80) :: hdr1b
  34    !data hdr1b /'SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SLNM ORBN SELV SURF RAINF'/
  35    data hdr1b /'SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SLNM ORBN SELV SFLG RFLAG'/
  36    character(10) :: date
  37    character(8)  :: subset, subfgn
  38    character(20) :: filename
  39
  40    logical :: outside, outside_all
  41
  42    integer(i_kind) :: iost, inst, lnbufr, ifgat
  43    integer(i_kind) :: num_bufr(7),numbufr,ibufr
  44    integer(i_kind) :: ihh, i, j, n, k, slnm, ifov, idd, ireadmg, ireadsb
  45    integer(i_kind) :: iret, idate, im, iy
  46    integer(i_kind) :: jc, incangl, bch, landsea_mask, rain_flag
  47    integer(i_kind) :: platform_id, satellite_id, sensor_id, nchan, num_ssmis_file
  48    integer(i_kind) :: num_ssmis_local, num_ssmis_global, num_ssmis_used, num_ssmis_thinned
  49    integer(i_kind) :: num_ssmis_used_tmp
  50
  51    real(r_double), dimension(2,maxchanl) :: bufrtbb
  52    real(r_double), dimension(n1bhdr)     :: bfr1bhdr
  53
  54    ! pixel information
  55    integer(i_kind)   ::  year,month,day,hour,minute,second  ! observation time
  56    real(kind=8)    ::  obs_time
  57    real(r_double), allocatable ::  tb_inv(:)          !  bright temperatures
  58
  59    type (datalink_type), pointer  :: head, p, current, prev
  60    type(info_type)                :: info
  61    type(model_loc_type)           :: loc
  62
  63    ! thinning variables
  64    integer(i_kind) :: itt,itx,iobs,iout
  65    real(r_kind)    :: terrain,timedif,crit,dist
  66    real(r_kind)    :: dlon_earth,dlat_earth
  67    logical         :: iuse
  68    real, allocatable :: in(:), out(:)
  69    logical           :: found, head_found
  70
  71    integer(i_kind), allocatable :: ptotal(:), nread(:)
  72
  73    call da_trace_entry("da_read_obs_bufrssmis")
  74
  75    allocate(nread(1:rtminit_nsensor))
  76    allocate(ptotal(0:num_fgat_time))
  77    nread(1:rtminit_nsensor) = 0
  78    ptotal(0:num_fgat_time) = 0
  79
  80    platform_id  = 2                 ! for DMSP series
  81    sensor_id    = 10                ! for SSMIS
  82    nchan        = nchan_ssmis
  83
  84    allocate (tb_inv(nchan))
  85    num_ssmis_file    = 0
  86    num_ssmis_local   = 0
  87    num_ssmis_global  = 0
  88    num_ssmis_used    = 0
  89    num_ssmis_thinned = 0
  90    iobs = 0                 ! for thinning, argument is inout
  91
  92    ! 0.0  Open unit to satellite bufr file and read file header
  93    !--------------------------------------------------------------
  94
  95 !   call da_get_unit(lnbufr)
  96
  97    num_bufr(:)=0
  98    numbufr=0
  99    if (num_fgat_time>1) then
 100       do i=1,7
 101          call da_get_unit(lnbufr)
 102          write(filename,fmt='(A,2I1,A)') trim(infile),0,i,'.bufr'
 103          open(unit   = lnbufr, FILE   = trim(filename),iostat =  iost, form = 'unformatted', STATUS = 'OLD')
 104          if (iost == 0) then
 105             numbufr=numbufr+1
 106             num_bufr(numbufr)=i
 107          else
 108             close (lnbufr)
 109          end if
 110          call da_free_unit(lnbufr)
 111       end do
 112    else
 113      numbufr=1
 114    end if
 115
 116    if (numbufr==0) numbufr=1
 117
 118 bufrfile:  do ibufr=1,numbufr
 119    if (num_fgat_time==1) then
 120       filename=trim(infile)//'.bufr'
 121    else
 122       if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then
 123          filename=trim(infile)//'.bufr'
 124       else
 125          write(filename,fmt='(A,2I1,A)') trim(infile),0,num_bufr(ibufr),'.bufr'
 126       end if
 127    end if
 128
 129    lnbufr=98
 130    open(unit=lnbufr,file=trim(filename),form='unformatted', &
 131       iostat = iost, status = 'old')
 132    if (iost /= 0) then
 133       call da_warning(__FILE__,__LINE__, &
 134          (/"Cannot open file "//infile/))
 135       call da_trace_exit("da_read_obs_bufrssmis")
 136       return
 137    end if
 138
 139    call openbf(lnbufr,'IN',lnbufr)
 140    call datelen(10)
 141    call readmg(lnbufr,subset,idate,iret)
 142
 143    iy=0
 144    im=0
 145    idd=0
 146    ihh=0
 147    write(unit=date,fmt='( i10)') idate
 148    read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
 149    write(unit=stdout,fmt=*) &
 150       'Bufr file date is ',iy,im,idd,ihh,infile
 151
 152    ! Loop to read bufr file and assign information to a sequential structure
 153    !-------------------------------------------------------------------------
 154    if ( ibufr == 1 ) then
 155       allocate (head)
 156       nullify  ( head % next )
 157       p => head
 158    end if
 159
 160 ! Set various variables depending on type of data to be read
 161
 162    !subfgn = 'NC003003'
 163    subfgn = 'NC021201'
 164    incangl = 53.2_r_kind
 165
 166    subset_loop: do while (ireadmg(lnbufr,subset,idate)==0)
 167
 168       read_loop: do while (ireadsb(lnbufr)==0 .and. subset==subfgn)
 169
 170          num_ssmis_file = num_ssmis_file + 1
 171          ! 1.0     Read header record and data record
 172
 173          call ufbint(lnbufr,bfr1bhdr,n1bhdr,1,iret,hdr1b)
 174          call ufbrep(lnbufr,bufrtbb,2,maxchanl,iret,"CHNM TMBR" )
 175
 176          ! check if observation outside range
 177
 178          ! 2.0     Extract observation location and other required information
 179          !     QC1:  judge if data is in the domain, read next record if not
 180          !------------------------------------------------------------------------
 181
 182          info%lat  =  bfr1bhdr(bufr_lat)
 183          info%lon  =  bfr1bhdr(bufr_lon)
 184          call da_llxy (info, loc, outside, outside_all)
 185
 186
 187          if (outside_all) cycle
 188
 189          !  3.0     Extract other information
 190
 191          info%elv  = 0.0
 192          landsea_mask = nint(bfr1bhdr(bufr_landsea_mask))   ! ssmis surface flag
 193                                                             ! 0:land, 2:near coast, 3:ice,
 194                                                             ! 4:possible ice, 5:ocean, 6:coast
 195          ! RTTOV surftype: 0:land, 1:sea, 2:sea ice
 196          if ( landsea_mask == 5 ) then
 197             landsea_mask = 1
 198          else if ( landsea_mask == 2 .or. landsea_mask == 6 ) then
 199             landsea_mask = 0
 200          else if ( landsea_mask == 3 .or. landsea_mask == 4 ) then
 201             landsea_mask = 2
 202          end if
 203          rain_flag = nint(bfr1bhdr(15))    ! 0:no rain, 1:rain
 204
 205          !------------------------------------------------------
 206          !  3.1     Extract satellite id and scan position.
 207
 208          if (nint(bfr1bhdr(bufr_satellite_id)) == bufsat_dmsp16) then
 209             satellite_id = 16
 210          else if (nint(bfr1bhdr(bufr_satellite_id)) == bufsat_dmsp17) then
 211             satellite_id = 17
 212          else if (nint(bfr1bhdr(bufr_satellite_id)) == bufsat_dmsp18) then
 213             satellite_id = 18
 214          end if
 215
 216          ! 3.3 Find wrfvar instrument index from RTTOV instrument triplet
 217          !     go to next data if id is not in the lists
 218
 219          inst = 0
 220          do i = 1, rtminit_nsensor
 221             if (platform_id  == rtminit_platform(i)      &
 222                .and. satellite_id == rtminit_satid(i)    &
 223                .and. sensor_id    == rtminit_sensor(i)) then
 224                inst = i
 225                exit
 226             end if
 227          end do
 228          if (inst == 0) cycle read_loop
 229
 230          !  3.1     Extract scan number and scan position.
 231
 232          slnm = nint(bfr1bhdr(11))
 233          ifov = nint(bfr1bhdr(bufr_ifov))
 234
 235          !  3.2     Extract date information.
 236
 237          year   = bfr1bhdr(bufr_year)
 238          month  = bfr1bhdr(bufr_month)
 239          day    = bfr1bhdr(bufr_day)
 240          hour   = bfr1bhdr(bufr_hour)
 241          minute = bfr1bhdr(bufr_minute)
 242          second = bfr1bhdr(bufr_second)
 243
 244          write(unit=info%date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 245             year, '-', month, '-', day, '_', hour, ':', minute, ':', second
 246
 247          !  QC3: time consistency check with the background date
 248
 249          if (year <= 99) then
 250             if (year < 78) then
 251                year = year + 2000
 252             else
 253                year = year + 1900
 254             end if
 255          end if
 256
 257          call da_get_julian_time(year,month,day,hour,minute,obs_time)
 258
 259          if (obs_time < time_slots(0) .or.  &
 260             obs_time >= time_slots(num_fgat_time)) cycle read_loop
 261
 262          ! 3.2.1 determine FGAT index ifgat
 263
 264          do ifgat=1,num_fgat_time
 265             if (obs_time >= time_slots(ifgat-1) .and.  &
 266                 obs_time  < time_slots(ifgat)) exit
 267          end do
 268
 269          num_ssmis_global = num_ssmis_global + 1
 270          ptotal(ifgat) = ptotal(ifgat) + 1
 271
 272          if (outside) cycle ! No good for this PE
 273
 274          num_ssmis_local = num_ssmis_local + 1
 275
 276          !  Make Thinning
 277          !  Map obs to thinning grid
 278          !-------------------------------------------------------------------
 279          if (thinning) then
 280             dlat_earth = info%lat
 281             dlon_earth = info%lon
 282             if (dlon_earth<zero) dlon_earth = dlon_earth+r360
 283             if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 284             dlat_earth = dlat_earth*deg2rad
 285             dlon_earth = dlon_earth*deg2rad
 286             timedif = 0. !2.0_r_kind*abs(tdiff)        ! range:  0 to 6
 287             terrain = 0.01_r_kind*abs(bfr1bhdr(13))
 288             crit = 1. !0.01_r_kind+terrain + timedif !+ 10._r_kind*float(iskip)
 289             call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 290             if (.not. iuse) then
 291                num_ssmis_thinned=num_ssmis_thinned+1
 292                cycle
 293             end if
 294          end if
 295
 296          num_ssmis_used = num_ssmis_used + 1
 297          nread(inst) = nread(inst) + 1
 298
 299          if (num_ssmis_used > max_ssmis_input) then
 300             write(unit=message(1),fmt='(A,I10,A)') &
 301                "Max number of ssmis",max_ssmis_input," reached"
 302             call da_warning(__FILE__,__LINE__,message(1:1))
 303             num_ssmis_used = num_ssmis_used - 1
 304             exit read_loop
 305          end if
 306
 307          ! 3.4 extract satellite and solar angle
 308
 309          ! 3.5 extract surface information
 310
 311          ! 3.6 extract channel bright temperature
 312
 313          tb_inv(1:nchan) = missing_r
 314
 315          do jc = 1, nchan
 316             bch = nint(bufrtbb(1,jc))     !ch index from bufr
 317             tb_inv(jc) = bufrtbb(2,jc)
 318             if (tb_inv(jc) < tbmin .or. tb_inv(jc) > tbmax .or. bch /= jc) then
 319                 tb_inv(jc) = missing_r
 320             end if
 321          end do
 322
 323          if ( maxval(tb_inv(:)) > missing_r ) then
 324
 325             !  4.0   assign information to sequential radiance structure
 326             !--------------------------------------------------------------------------
 327             allocate (p % tb_inv (1:nchan))
 328             p%info                = info
 329             p%loc                 = loc
 330             p%landsea_mask        = landsea_mask
 331             p%scanline            = slnm
 332             p%scanpos             = ifov
 333             p%satzen              = incangl
 334             p%satazi              = 0.0     ! dummy value
 335             p%solzen              = 0.0     ! dummy value
 336             p%tb_inv(1:nchan)     = tb_inv(1:nchan)
 337             p%sensor_index        = inst
 338             p%ifgat               = ifgat
 339             p%rain_flag           = rain_flag
 340
 341             allocate (p%next)   ! add next data
 342             p => p%next
 343             nullify (p%next)
 344
 345          else
 346
 347             num_ssmis_local = num_ssmis_local - 1
 348             num_ssmis_used  = num_ssmis_used - 1
 349
 350          end if
 351
 352       end do read_loop
 353
 354    end do subset_loop
 355
 356   call closbf(lnbufr)
 357   close(lnbufr)
 358
 359 end do bufrfile
 360
 361    if (thinning .and. num_ssmis_global > 0 ) then
 362
 363 #ifdef DM_PARALLEL
 364
 365       ! Get minimum crit and associated processor index.
 366       j = 0
 367       do ifgat = 1, num_fgat_time
 368          do n = 1, iv%num_inst
 369             j = j + thinning_grid(n,ifgat)%itxmax
 370          end do
 371       end do
 372
 373       allocate ( in  (j) )
 374       allocate ( out (j) )
 375       j = 0
 376       do ifgat = 1, num_fgat_time
 377          do n = 1, iv%num_inst
 378             do i = 1, thinning_grid(n,ifgat)%itxmax
 379                j = j + 1
 380                in(j) = thinning_grid(n,ifgat)%score_crit(i)
 381             end do
 382          end do
 383       end do
 384       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 385
 386       call wrf_dm_bcast_real (out, j)
 387
 388       j = 0
 389       do ifgat = 1, num_fgat_time
 390          do n = 1, iv%num_inst
 391             do i = 1, thinning_grid(n,ifgat)%itxmax
 392                j = j + 1
 393                if ( ABS(out(j)-thinning_grid(n,ifgat)%score_crit(i)) > 1.0E-10 ) thinning_grid(n,ifgat)%ibest_obs(i) = 0
 394             end do
 395          end do
 396       end do
 397
 398       deallocate( in  )
 399       deallocate( out )
 400
 401 #endif
 402
 403       ! Delete the nodes which being thinning out
 404       p => head
 405       prev => head
 406       head_found = .false.
 407       num_ssmis_used_tmp = num_ssmis_used
 408       do j = 1, num_ssmis_used_tmp
 409          n = p%sensor_index
 410          ifgat = p%ifgat
 411          found = .false.
 412
 413          do i = 1, thinning_grid(n,ifgat)%itxmax
 414             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 415                found = .true.
 416                exit
 417             end if
 418          end do
 419
 420          ! free current data
 421          if ( .not. found ) then
 422             current => p
 423             p => p%next
 424             if ( head_found ) then
 425                prev%next => p
 426             else
 427                head => p
 428                prev => p
 429             end if
 430             deallocate ( current % tb_inv )
 431             deallocate ( current )
 432             num_ssmis_thinned = num_ssmis_thinned + 1
 433             num_ssmis_used = num_ssmis_used - 1
 434             nread(n) = nread(n) - 1
 435             continue
 436          end if
 437
 438          if ( found .and. head_found ) then
 439             prev => p
 440             p => p%next
 441             continue
 442          end if
 443
 444          if ( found .and. .not. head_found ) then
 445             head_found = .true.
 446             head => p
 447             prev => p
 448             p => p%next
 449          end if
 450
 451       end do
 452
 453    end if  ! End of thinning
 454
 455    iv%total_rad_pixel   = iv%total_rad_pixel + num_ssmis_used
 456    iv%total_rad_channel = iv%total_rad_channel + num_ssmis_used*nchan
 457
 458    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_ssmis_used
 459    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_ssmis_global
 460
 461    do i = 1, num_fgat_time
 462       ptotal(i) = ptotal(i) + ptotal(i-1)
 463       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 464    end do
 465    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 466       write(unit=message(1),fmt='(A,I10,A,I10)') &
 467           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 468       call da_warning(__FILE__,__LINE__,message(1:1))
 469    endif
 470
 471    write(unit=stdout,fmt='(a)') 'num_ssmis_file, num_ssmis_global, num_ssmis_local, num_ssmis_used, num_ssmis_thinned'
 472    write(stdout,*) num_ssmis_file, num_ssmis_global, num_ssmis_local, num_ssmis_used, num_ssmis_thinned
 473
 474    deallocate(tb_inv)
 475
 476    !  5.0 allocate innovation radiance structure
 477    !----------------------------------------------------------------
 478    do i = 1, iv%num_inst
 479       if (nread(i) < 1) cycle
 480       iv%instid(i)%num_rad  = nread(i)
 481       iv%instid(i)%info%nlocal = nread(i)
 482       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 483         'Allocating space for radiance innov structure', &
 484          i, iv%instid(i)%rttovid_string, iv%instid(i)%num_rad
 485
 486       call da_allocate_rad_iv (i, nchan, iv)
 487
 488    end do
 489
 490    !  6.0 assign sequential structure to innovation structure
 491    !-------------------------------------------------------------
 492    nread(1:rtminit_nsensor) = 0
 493    p => head
 494
 495    do n = 1, num_ssmis_used
 496       i = p%sensor_index
 497       nread(i) = nread(i) + 1
 498       call da_initialize_rad_iv (i, nread(i), iv, p)
 499
 500       iv%instid(i)%rain_flag(n) = p%rain_flag
 501
 502       current => p
 503       p => p%next
 504
 505       ! free current data
 506       deallocate ( current % tb_inv )
 507       deallocate ( current )
 508
 509    end do
 510
 511    deallocate ( p )
 512    deallocate (nread)
 513    deallocate (ptotal)
 514
 515    call closbf(lnbufr)
 516    close(lnbufr)
 517 !   call da_free_unit(lnbufr)
 518
 519    call da_trace_exit("da_read_obs_bufrssmis")
 520 #else
 521    call da_error(__FILE__,__LINE__,(/"Needs to be compiled with a BUFR library"/))
 522 #endif
 523
 524
 525 end subroutine da_read_obs_bufrssmis
 526