var/da/da_radiance/da_read_obs_fy3.inc

   1 subroutine da_read_obs_fy3 (obstype,iv, infile)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: read in fy3 1b data to innovation structure
   5    !
   6    !   Dong peiming 2012/03/09
   7    !   METHOD: use F90 sequential data structure to avoid reading file twice
   8    !            so that da_scan_bufrtovs is not necessary any more.
   9    !            1. read file radiance data in sequential data structure
  10    !            2. do gross QC check
  11    !            3. assign sequential data structure to innovation structure
  12    !               and deallocate sequential data structure
  13    !---------------------------------------------------------------------------
  14
  15    implicit none
  16
  17    character(5)      ,  intent (in)    :: obstype
  18    character(20)     ,  intent (in)    :: infile
  19    type (iv_type)    ,  intent (inout) :: iv
  20
  21 !
  22 TYPE type_rad_FY3
  23     INTEGER :: yyyy,mn,dd,hh,mm,ss
  24     INTEGER :: iscanline,iscanpos
  25     REAL*4  :: rlat,rlon !lat/lon in degrees for Anfovs
  26     INTEGER :: isurf_height, isurf_type !height/type for Anfovs
  27     REAL*4  :: satzen,satazi,solzen,solazi !scan angles for Anfovs
  28     REAL*4  :: tbb(20) !bright temperatures
  29 !   REAL*4  :: btemps(20)
  30     INTEGER :: iavhrr(13),ihirsflag
  31     INTEGER :: iprepro(5) ! values from pre-processing
  32     REAL*4  :: clfra ! Cloud cover (<1.0)
  33     REAL*4  :: ts ! Skin temperature
  34     REAL*4  :: tctop ! Cloud top temperature
  35 END TYPE type_rad_FY3
  36
  37   TYPE (type_rad_FY3)   :: rad
  38   integer :: iscan,nscan
  39 !
  40    integer          :: iost
  41    integer(i_kind), allocatable :: nread(:)
  42
  43 !dongpm   logical hirs2,hirs3,hirs4,msu,amsua,amsub,mhs
  44    logical mwts,mwhs
  45    logical outside, outside_all, iuse
  46    integer :: inst
  47
  48    integer(i_kind) i,j,k,ifov
  49    integer(i_kind) nchan
  50    integer :: num_bufr(7), numbufr, ibufr
  51    character(20) :: filename
  52
  53    ! thinning variables
  54    integer(i_kind) itt,itx,iobs,iout
  55    real(r_kind) terrain,timedif,crit,dist
  56    real(r_kind) dlon_earth,dlat_earth
  57
  58    real(r_kind) tbmin,tbmax, tbbad
  59    ! real(r_kind) rmask
  60
  61    ! Instrument triplet, follow the convension of RTTOV
  62    integer   :: platform_id, satellite_id, sensor_id
  63
  64    ! pixel information
  65    integer   ::  year,month,day,hour,minute,second  ! observation time
  66    real*8    ::  obs_time
  67    ! real      ::  rlat, rlon                         !  lat/lon in degrees   for Anfovs
  68    real      ::  satzen, satazi, solzen ,solazi       !  scan angles for Anfovs
  69    integer   ::  landsea_mask
  70    real      ::  srf_height
  71    ! channels' bright temperature
  72    real , allocatable ::   tb_inv(:)                    !  bright temperatures
  73    !  end type bright_temperature
  74
  75    type (datalink_type), pointer    :: head, p, current, prev
  76
  77    integer                        ::  ifgat
  78    type(info_type)                ::  info
  79    type(model_loc_type)           ::  loc
  80
  81    data tbmin,tbmax,tbbad / 50.0_r_kind, 550.0_r_kind, -9.99e11_r_kind /
  82    integer :: num_tovs_local, num_tovs_file, num_tovs_global, num_tovs_selected
  83    integer :: num_tovs_thinned, num_tovs_used, num_tovs_used_tmp
  84    integer :: lnbufr
  85    integer :: n
  86    integer(i_kind), allocatable :: ptotal(:)
  87    real , allocatable :: in(:), out(:)
  88    logical :: found, head_found
  89
  90    if (trace_use) call da_trace_entry("da_read_obs_fy3")
  91
  92    ! Initialize variables
  93
  94    nchan = 20
  95    allocate(nread(1:rtminit_nsensor))
  96    allocate(ptotal(0:num_fgat_time))
  97    nread(1:rtminit_nsensor) = 0
  98    ptotal(0:num_fgat_time) = 0
  99
 100    ! Set various variables depending on type of data to be read
 101
 102    ! platform_id  = 1                 !! for NOAA series
 103    ! platform_id  = 10                !! for METOP series
 104
 105 !dongpm   hirs2=     obstype == 'hirs2'
 106 !dongpm   hirs3=     obstype == 'hirs3'
 107 !dongpm   hirs4=     obstype == 'hirs4'
 108 !dongpm   msu=       obstype == 'msu  '
 109 !dongpm   amsua=     obstype == 'amsua'
 110 !dongpm   amsub=     obstype == 'amsub'
 111 !dongpm   mhs=       obstype == 'mhs  '
 112           mwts=      obstype == 'mwts '
 113           mwhs=      obstype == 'mwhs '
 114
 115 !dongpm   if (hirs2) then
 116 !dongpm      sensor_id    =  0
 117 !dongpm      step   = 1.80_r_kind
 118 !dongpm      start  = -49.5_r_kind
 119 !dongpm      nchan=nchan_hirs2
 120 !dongpm      subfgn='NC021021'
 121 !dongpm      rato=1.1363987_r_kind
 122 !dongpm   else if (hirs3) then
 123 !dongpm      sensor_id    =  0
 124 !dongpm      step   = 1.80_r_kind
 125 !dongpm      start  = -49.5_r_kind
 126 !dongpm      nchan=nchan_hirs3
 127 !dongpm      subfgn='NC021025'
 128 !dongpm   else if (hirs4) then
 129 !dongpm      sensor_id    =  0
 130 !dongpm      step   = 1.80_r_kind
 131 !dongpm      start  = -49.5_r_kind
 132 !dongpm      nchan=nchan_hirs4
 133 !dongpm      subfgn='NC021028'
 134 !dongpm   else if (mhs) then
 135 !dongpm      sensor_id    =  15
 136 !dongpm      step   = 10.0_r_kind/9.0_r_kind
 137 !dongpm      start  = -445.0_r_kind/9.0_r_kind
 138 !dongpm      nchan=nchan_mhs
 139 !dongpm      subfgn='NC021027'
 140 !dongpm   else if (msu) then
 141 !dongpm      sensor_id    =  1
 142 !dongpm      step   = 9.474_r_kind
 143 !dongpm      start  = -47.37_r_kind
 144 !dongpm      nchan=nchan_msu
 145 !dongpm      subfgn='NC021022'
 146 !dongpm      rato=1.1363987_r_kind
 147 !dongpm   else if (amsua) then
 148 !dongpm      sensor_id    =  3
 149 !dongpm      step   = three + one/three
 150 !dongpm      start  = -48.33_r_kind
 151 !dongpm      nchan=nchan_amsua
 152 !dongpm      subfgn='NC021023'
 153 !dongpm   else if (amsub)  then
 154 !dongpm      sensor_id    =  4
 155 !dongpm      step   = 1.1_r_kind
 156 !dongpm      start  = -48.95_r_kind
 157 !dongpm      nchan=nchan_amsub
 158 !dongpm      subfgn='NC021024'
 159 !dongpm   end if
 160           if (mwts) then
 161            sensor_id    =  40
 162            nchan=4
 163            nscan=15
 164           else if(mwhs) then
 165            sensor_id    =  41
 166            nchan=5
 167            nscan=98
 168           endif
 169
 170    allocate (tb_inv(nchan))
 171
 172    num_tovs_file     = 0    ! number of obs in file
 173    num_tovs_global   = 0    ! number of obs within whole domain
 174    num_tovs_local    = 0    ! number of obs within tile
 175    num_tovs_thinned  = 0    ! number of obs rejected by thinning
 176    num_tovs_used     = 0    ! number of obs entered into innovation computation
 177    num_tovs_selected = 0    ! number of obs limited for debuging
 178    iobs = 0                 ! for thinning, argument is inout
 179
 180    ! 0.0  Open unit to satellite bufr file and read file header
 181    !--------------------------------------------------------------
 182
 183    num_bufr(:)=0
 184    numbufr=0
 185    if (num_fgat_time>1) then
 186       do i=1,7
 187          call da_get_unit(lnbufr)
 188          write(filename,fmt='(A,2I1,A)') trim(infile),0,i,'.dat'
 189          open(unit   = lnbufr, FILE   = trim(filename),iostat =  iost, form = 'unformatted', STATUS = 'OLD')
 190          if (iost == 0) then
 191             numbufr=numbufr+1
 192             num_bufr(numbufr)=i
 193          else
 194             close (lnbufr)
 195          end if
 196          call da_free_unit(lnbufr)
 197       end do
 198    else
 199      numbufr=1
 200    end if
 201
 202    if (numbufr==0) numbufr=1
 203
 204 bufrfile:  do ibufr=1,numbufr
 205    if (num_fgat_time==1) then
 206       filename=trim(infile)//'.dat'
 207    else
 208       if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then
 209          filename=trim(infile)//'.dat'
 210       else
 211          write(filename,fmt='(A,2I1,A)') trim(infile),0,num_bufr(ibufr),'.dat'
 212       end if
 213    end if
 214
 215 !  We want to use specific unit number for bufr data, so we can control the endian format in environment.
 216    lnbufr = 99
 217
 218    open(unit=lnbufr,file=trim(filename),form='unformatted', &
 219       iostat = iost, status = 'old')
 220    if (iost /= 0) then
 221       call da_warning(__FILE__,__LINE__, &
 222          (/"Cannot open file "//filename/))
 223       if (trace_use) call da_trace_exit("da_read_obs_fy3")
 224       return
 225    end if
 226
 227
 228    if ( ibufr == 1 ) then
 229       allocate (head)
 230       !  allocate ( head % tb_inv (1:nchan) )
 231       nullify  ( head % next )
 232       p => head
 233    endif
 234
 235    if (tovs_start > 1) then
 236       write (unit=stdout,fmt='(A,I6)') "   Skipping tovs obs before", tovs_start
 237    end if
 238
 239    obs: do while (.true.)
 240         do iscan=1,nscan
 241
 242          ! 1.0     Read fy3 data
 243          read(lnbufr,end=1000) rad
 244
 245          num_tovs_file = num_tovs_file + 1
 246
 247          ! 2.0     Extract observation location and other required information
 248          !     QC1:  judge if data is in the domain, read next record if not
 249          !------------------------------------------------------------------------
 250          ! rlat = bfr1bhdr(bufr_lat)
 251          ! rlon = bfr1bhdr(bufr_lat)
 252          ! if (rlon < 0.0) rlon = rlon+360.0
 253
 254          info%lat  =  rad%rlat
 255
 256          info%lon  =  rad%rlon
 257          call da_llxy (info, loc, outside, outside_all)
 258
 259          if (outside_all) cycle
 260
 261          !  3.0     Extract other information
 262          !------------------------------------------------------
 263          !  3.1     Extract satellite id and scan position.
 264
 265          platform_id = 23
 266          if(infile(5:5)=='a') then
 267             satellite_id = 1
 268          elseif(infile(5:5)=='b') then
 269             satellite_id = 2
 270          else
 271             call da_warning(__FILE__,__LINE__,(/"Can not assimilate data from this instrument"/))
 272             if (trace_use) call da_trace_exit("da_read_obs_fy3")
 273             return
 274          endif
 275          ifov = rad%iscanpos
 276
 277          !  QC2:  limb pixel rejected (not implemented)
 278
 279          !  3.2     Extract date information.
 280
 281          year   = rad%yyyy
 282          month  = rad%mn
 283          day    = rad%dd
 284          hour   = rad%hh
 285          minute = rad%mm
 286          second = rad%ss
 287 !dongpm for test
 288 !          year   = 2008
 289 !          month  = 8
 290 !          day    = 5
 291 !          hour   = 18
 292 !          minute = 0
 293 !          second = 0
 294
 295          write(unit=info%date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
 296             year, '-', month, '-', day, '_', hour, ':', minute, ':', second
 297
 298          !  QC3: time consistency check with the background date
 299
 300          if (year <= 99) then
 301             if (year < 78) then
 302                year = year + 2000
 303             else
 304                year = year + 1900
 305             end if
 306          end if
 307
 308          call da_get_julian_time(year,month,day,hour,minute,obs_time)
 309
 310          if (obs_time < time_slots(0) .or.  &
 311             obs_time >= time_slots(num_fgat_time)) cycle
 312
 313          ! 3.2.1 determine FGAT index ifgat
 314
 315          do ifgat=1,num_fgat_time
 316             if (obs_time >= time_slots(ifgat-1) .and.  &
 317                 obs_time  < time_slots(ifgat)) exit
 318          end do
 319
 320          ! 3.3 Find wrfvar instrument index from RTTOV instrument triplet
 321          !     go to next data if id is not in the lists
 322
 323          inst = 0
 324          do i = 1, rtminit_nsensor
 325             if (platform_id  == rtminit_platform(i) &
 326                .and. satellite_id == rtminit_satid(i)    &
 327                .and. sensor_id    == rtminit_sensor(i)) then
 328                inst = i
 329                exit
 330             end if
 331          end do
 332          if (inst == 0) cycle
 333
 334          ! 3.4 extract satellite and solar angle
 335
 336             satzen = rad%satzen !*deg2rad   ! local zenith angle
 337             satzen = abs(satzen)
 338
 339             ! if (amsua .and. ifov .le. 15) satzen=-satzen
 340             ! if (amsub .and. ifov .le. 45) satzen=-satzen
 341             ! if (hirs3 .and. ifov .le. 28) satzen=-satzen
 342 !dongpm         if ( satzen > 65.0 ) cycle   ! CRTM has a satzen > 65.0 check
 343          satazi = rad%satazi*0.01           ! look angle
 344          ! if (satazi<0.0) satazi = satazi+360.0
 345          solzen = rad%solzen*0.01              ! solar zenith angle
 346          solazi = rad%solazi*0.01              !RTTOV9_3
 347
 348          num_tovs_global = num_tovs_global + 1
 349          ptotal(ifgat) = ptotal(ifgat) + 1
 350
 351          if (num_tovs_global < tovs_start) then
 352             cycle
 353          end if
 354
 355          if (num_tovs_global > tovs_end) then
 356             write (unit=stdout,fmt='(A,I6)') "   Skipping radiance obs after", tovs_end
 357             exit obs
 358          end if
 359
 360          num_tovs_selected = num_tovs_selected + 1
 361
 362          if (num_tovs_selected > max_tovs_input) then
 363             write(unit=message(1),fmt='(A,I10,A)') &
 364                "Max number of tovs",max_tovs_input," reached"
 365             call da_warning(__FILE__,__LINE__,message(1:1))
 366             num_tovs_selected = num_tovs_selected-1
 367             num_tovs_global   = num_tovs_global-1
 368             ptotal(ifgat) = ptotal(ifgat) - 1
 369             exit obs
 370          end if
 371
 372          if (outside) cycle ! No good for this PE
 373          num_tovs_local = num_tovs_local + 1
 374
 375          !  Make Thinning
 376          !  Map obs to thinning grid
 377          !-------------------------------------------------------------------
 378          if (thinning) then
 379             dlat_earth = info%lat
 380             dlon_earth = info%lon
 381             if (dlon_earth<zero) dlon_earth = dlon_earth+r360
 382             if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
 383             dlat_earth = dlat_earth*deg2rad
 384             dlon_earth = dlon_earth*deg2rad
 385             timedif = 0.0 !2.0_r_kind*abs(tdiff)        ! range:  0 to 6
 386 !dongpm            terrain = 0.01_r_kind*abs(bfr1bhdr(13))
 387             terrain = 0.01_r_kind*abs(rad%satzen)
 388             crit = 1.0 !0.01_r_kind+terrain + timedif !+ 10.0_r_kind*float(iskip)
 389             call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
 390             if (.not. iuse) then
 391                num_tovs_thinned=num_tovs_thinned+1
 392                cycle
 393             end if
 394          end if
 395
 396          num_tovs_used = num_tovs_used + 1
 397          nread(inst) = nread(inst) + 1
 398
 399          ! 3.5 extract surface information
 400          srf_height = rad%isurf_height          ! station height
 401          if (srf_height < 8888.0 .AND. srf_height > -416.0) then
 402          else
 403             srf_height = 0.0
 404          endif
 405
 406 !dongpm         landsea_mask = rad%isurf_type  ! 0:land ; 1:sea (same as RTTOV)
 407 !fy3 isurf_type is just reversed as RTTOV
 408          if(rad%isurf_type .eq. 0) then   ! sea
 409            landsea_mask = 1
 410          elseif(rad%isurf_type .eq. 1) then   !coast
 411            landsea_mask = 0
 412          elseif(rad%isurf_type .eq. 2) then   !land
 413            landsea_mask = 0
 414          else
 415            landsea_mask = rad%isurf_type
 416            write(unit=message(1),fmt='(A,I6)') 'Unknown surface type: ', landsea_mask
 417            call da_warning(__FILE__,__LINE__,message(1:1))
 418          endif
 419          ! rmask=one                          ! land
 420          ! if (nint(bfr1bhdr(bufr_landsea_mask))==1) rmask=0.0_r_kind   ! reverse the land/sea mask in bufr
 421          ! landsea_mask = rmask+.001_r_kind             ! land sea mask
 422
 423          info%elv = srf_height
 424
 425          ! 3.6 extract channel bright temperature
 426
 427          tb_inv(1:nchan) = rad%tbb(1:nchan)
 428          do k = 1, nchan
 429             if ( tb_inv(k) < tbmin .or. tb_inv(k) > tbmax) &
 430                tb_inv(k) = missing_r
 431          end do
 432          if ( all(tb_inv<0.0) ) then
 433             num_tovs_local = num_tovs_local -1
 434             num_tovs_used = num_tovs_used - 1
 435             nread(inst) = nread(inst) - 1
 436             cycle
 437          end if
 438
 439          !  4.0   assign information to sequential radiance structure
 440          !--------------------------------------------------------------------------
 441          allocate (p % tb_inv (1:nchan))
 442          p%info             = info
 443          p%loc              = loc
 444          p%landsea_mask     = landsea_mask
 445          p%scanpos          = ifov
 446          p%satzen           = satzen
 447          p%satazi           = satazi
 448          p%solzen           = solzen
 449          p%tb_inv(1:nchan)  = tb_inv(1:nchan)
 450          p%sensor_index     = inst
 451          p%ifgat            = ifgat
 452 !RTTOV9_3
 453          p%solazi           = solazi
 454  !end of RTTOV9_3
 455          allocate (p%next)   ! add next data
 456
 457          p => p%next
 458          nullify (p%next)
 459       end do
 460    end do obs
 461
 462    call closbf(lnbufr)
 463    close(lnbufr)
 464 1000  continue
 465 end do bufrfile
 466
 467    if (thinning .and. num_tovs_global > 0 ) then
 468
 469 #ifdef DM_PARALLEL
 470
 471       ! Get minimum crit and associated processor index.
 472       j = 0
 473       do ifgat = 1, num_fgat_time
 474          do n = 1, iv%num_inst
 475             j = j + thinning_grid(n,ifgat)%itxmax
 476          end do
 477       end do
 478
 479       allocate ( in  (j) )
 480       allocate ( out (j) )
 481
 482       j = 0
 483       do ifgat = 1, num_fgat_time
 484          do n = 1, iv%num_inst
 485             do i = 1, thinning_grid(n,ifgat)%itxmax
 486                j = j + 1
 487                in(j) = thinning_grid(n,ifgat)%score_crit(i)
 488             end do
 489          end do
 490       end do
 491       call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
 492
 493       call wrf_dm_bcast_real (out, j)
 494
 495       j = 0
 496       do ifgat = 1, num_fgat_time
 497          do n = 1, iv%num_inst
 498             do i = 1, thinning_grid(n,ifgat)%itxmax
 499                j = j + 1
 500                if ( ABS(out(j)-thinning_grid(n,ifgat)%score_crit(i)) > 1.0E-10 ) thinning_grid(n,ifgat)%ibest_obs(i)  = 0
 501             end do
 502          end do
 503       end do
 504
 505       deallocate( in  )
 506       deallocate( out )
 507
 508 #endif
 509
 510       ! Delete the nodes which being thinning out
 511       p => head
 512       prev => head
 513       head_found = .false.
 514       num_tovs_used_tmp = num_tovs_used
 515       do j = 1, num_tovs_used_tmp
 516          n = p%sensor_index
 517          ifgat = p%ifgat
 518          found = .false.
 519
 520          do i = 1, thinning_grid(n,ifgat)%itxmax
 521             if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
 522                found = .true.
 523                exit
 524             endif
 525          end do
 526
 527          ! free current data
 528          if ( .not. found ) then
 529             current => p
 530             p => p%next
 531             if ( head_found ) then
 532                prev%next => p
 533             else
 534                head => p
 535                prev => p
 536             endif
 537             deallocate ( current % tb_inv )
 538             deallocate ( current )
 539             num_tovs_thinned = num_tovs_thinned + 1
 540             num_tovs_used = num_tovs_used - 1
 541             nread(n) = nread(n) - 1
 542             continue
 543          endif
 544
 545          if ( found .and. head_found ) then
 546             prev => p
 547             p => p%next
 548             continue
 549          endif
 550
 551          if ( found .and. .not. head_found ) then
 552             head_found = .true.
 553             head => p
 554             prev => p
 555             p => p%next
 556          endif
 557
 558       end do
 559
 560    endif  ! End of thinning
 561
 562    iv%total_rad_pixel   = iv%total_rad_pixel   + num_tovs_used
 563    iv%total_rad_channel = iv%total_rad_channel + num_tovs_used*nchan
 564
 565    iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_tovs_used
 566    iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_tovs_global
 567
 568    do i = 1, num_fgat_time
 569       ptotal(i) = ptotal(i) + ptotal(i-1)
 570       iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
 571    end do
 572    if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
 573       write(unit=message(1),fmt='(A,I10,A,I10)') &
 574           "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
 575       call da_warning(__FILE__,__LINE__,message(1:1))
 576    endif
 577
 578    write(unit=stdout,fmt='(a)') 'num_tovs_file num_tovs_global num_tovs_local num_tovs_used num_tovs_thinned'
 579    write(unit=stdout,fmt='(5i10)') num_tovs_file,num_tovs_global, num_tovs_local,num_tovs_used,num_tovs_thinned
 580
 581    deallocate(tb_inv)
 582
 583    !  5.0 allocate innovation radiance structure
 584    !----------------------------------------------------------------
 585
 586    do i = 1, iv%num_inst
 587       if (nread(i) < 1) cycle
 588       iv%instid(i)%num_rad = nread(i)
 589       iv%instid(i)%info%nlocal = nread(i)
 590       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
 591         'Allocating space for radiance innov structure', &
 592          i, iv%instid(i)%rttovid_string, iv%instid(i)%num_rad
 593
 594       call da_allocate_rad_iv(i,nchan,iv)
 595
 596    end do
 597
 598    !  6.0 assign sequential structure to innovation structure
 599    !-------------------------------------------------------------
 600    nread(1:rtminit_nsensor) = 0
 601    p => head
 602    ! do while ( associated(p) )
 603
 604    do n = 1, num_tovs_used
 605       i = p%sensor_index
 606       nread(i) = nread(i) + 1
 607
 608       call da_initialize_rad_iv (i, nread(i), iv, p)
 609
 610       current => p
 611       p => p%next
 612
 613       ! free current data
 614       deallocate ( current % tb_inv )
 615       deallocate ( current )
 616    end do
 617
 618    deallocate ( p )
 619
 620    deallocate (nread)
 621    deallocate (ptotal)
 622
 623    ! check if sequential structure has been freed
 624    !
 625    ! p => head
 626    ! do i = 1, num_rad_selected
 627    !    write (unit=stdout,fmt=*)  i, p%tb_inv(1:nchan)
 628    !    p => p%next
 629    ! end do
 630
 631    if (trace_use) call da_trace_exit("da_read_obs_fy3")
 632
 633
 634
 635 end subroutine da_read_obs_fy3
 636
 637