da_read_bufrairs.inc

References to this file elsewhere.
1 subroutine da_read_bufrairs(obstype,iv,xp,infile)
2 !--------------------------------------------------------
3 !  Purpose: read in NCEP bufr eos AIRS/AMSUA/HSB 1b data 
4 !            to innovation structure
5 !
6 !   METHOD: use F90 sequantial data structure to avoid read file twice
7 !            so that da_scan_bufrairs is not necessary any more.
8 !            1. read file radiance data in sequential data structure
9 !            2. do gross QC check
10 !            3. assign sequential data structure to innovation structure
11 !                  and deallocate sequential data structure
12 !
13 !  HISTORY: 2006/01/03 - Creation            Zhiquan Liu
14 !
15 !------------------------------------------------------------------------------
16 
17   implicit none
18 
19   character(9)      ,  intent (in)  :: obstype
20   character(100)    ,  intent (in)  :: infile
21   type (xpose_type) ,  intent (in)  :: xp
22   type (ob_type)    ,intent (inout) :: iv
23 
24 #ifdef BUFR
25 
26 ! Number of channels for sensors in BUFR
27   integer(i_kind),parameter :: N_AIRSCHAN = 281	 !- 281 subset ch out of 2378 ch for AIRS
28   integer(i_kind),parameter :: N_AMSUCHAN =  15  
29   integer(i_kind),parameter :: N_HSBCHAN  =   4
30   integer(i_kind),parameter :: N_MAXCHAN  = 350
31   integer(i_kind),parameter :: maxinfo    =  12
32 
33 
34 ! BUFR format for AQUASPOT 
35   integer(i_kind),parameter :: N_AQUASPOT_LIST = 25
36   type aquaspot_list
37      sequence
38      real(r_double) :: said   ! Satellite identifier
39      real(r_double) :: orbn   ! Orbit number
40      real(r_double) :: slnm   ! Scan line number 
41      real(r_double) :: mjfc   ! Major frame count
42      real(r_double) :: selv   ! Height of station
43      real(r_double) :: soza   ! Solar zenith angle
44      real(r_double) :: solazi ! Solar azimuth angle
45      real(r_double) :: intms(2,9) ! SATELLITE inSTRUMENT TEMPERATURES
46   end type aquaspot_list
47 
48 
49 ! BUFR format for AIRSSPOT 
50   integer(i_kind),parameter :: N_AIRSSPOT_LIST = 12
51   type airsspot_list
52      sequence
53      real(r_double) :: siid  ! Satellite instruments
54      real(r_double) :: year
55      real(r_double) :: mnth
56      real(r_double) :: days
57      real(r_double) :: hour
58      real(r_double) :: minu
59      real(r_double) :: seco
60      real(r_double) :: clath ! Latitude (high accuracy)
61      real(r_double) :: clonh ! Longitude (high accuracy)
62      real(r_double) :: saza  ! Satellite zenith angle 
63      real(r_double) :: bearaz ! Bearing or azimuth
64      real(r_double) :: fovn  ! Field of view number
65   end type airsspot_list
66 
67 
68 ! BUFR format for AIRSCHAN 
69   integer(i_kind),parameter :: N_AIRSCHAN_LIST = 4
70   type airschan_list
71      sequence
72      real(r_double) :: chnm    ! Channel number
73      real(r_double) :: logrcw  ! Log-10 of (Temperature-radiance central wavenumber
74      real(r_double) :: acqf    ! Channel quality flags for ATOVS
75      real(r_double) :: tmbrst  ! Brightness temperature
76   end type airschan_list
77   
78 
79 ! BUFR talble file sequencial number
80   character(len=512)  :: table_file
81 
82 
83 ! Variables for BUFR IO    
84   type(aquaspot_list) :: aquaspot
85   type(airsspot_list) :: airsspot
86   type(airschan_list) :: airschan(N_MAXCHAN)
87   
88   real(r_kind)      :: step, start
89   real(r_kind)      :: airdata(N_AIRSCHAN+maxinfo)
90   character(len=8)  :: subset
91   character(len=4)  :: senname
92   character(len=8)  :: spotname
93   character(len=8)  :: channame
94   integer(i_kind)   :: nchan,nchanr
95   integer(i_kind)   :: iret
96 
97 
98 ! Work variables for time
99   integer(i_kind)   :: idate, ifgat
100   integer(i_kind)   :: idate5(6)
101   character(len=10) :: date
102   integer(i_kind)   :: nmind
103   integer(i_kind)   :: iy, im, idd, ihh
104   real              :: obs_time
105 
106 
107 ! Other work variables
108   integer(i_kind)  :: nreal, nobs, ityp,ityp2
109   integer(i_kind)  :: itx, k, itt, iobsout
110   integer(i_kind),dimension(19)::icount
111   real(r_kind)     :: rlat, rlon, dx, dy, dx1, dy1, sstx, dlon, dlat
112   real(r_kind)     :: sza, timedif, pred, crit1
113   integer(i_kind)  :: klat1, klon1, klatp1, klonp1
114   integer(i_kind)  :: ifov,size
115   integer(i_kind)  :: inst,platform_id,satellite_id,sensor_id
116   logical          :: iflag,outside,inside_halo
117   integer(i_kind)  :: i, l, n, error, airs_table_unit
118   integer          :: iost, lnbufr
119   real(r_kind),allocatable,dimension(:,:):: airdata_all
120 
121 ! Set standard parameters
122   real(r_kind)     :: POinT001 =   0.001_r_kind
123   real(r_kind)     :: POinT01  =   0.01_r_kind
124   real(r_kind)     :: TEN      =  10._r_kind
125   real(r_kind)     :: R45      =  45._r_kind
126   real(r_kind)     :: R60      =  60._r_kind
127   real(r_kind)     :: R90      =  90._r_kind
128   real(r_kind)     :: R180     = 180._r_kind
129   real(r_kind)     :: R360     = 360._r_kind
130 
131   logical           :: airs, eos_amsua, hsb, airstab
132   type(info_type)            :: info
133   type(model_loc_type)       :: loc
134   type (datalink_type), pointer   :: head, p, current
135 
136    if (trace_use) call da_trace_entry("da_read_bufrairs")
137 
138 !  0.0  Initialize variables
139 !-----------------------------------
140   platform_id  = 9   ! eos series
141   satellite_id = 2   ! eos-2
142   nreal  = maxinfo
143   nobs   = 0
144   airs=      obstype == 'airs     '
145   eos_amsua= obstype == 'eos_amsua'
146   hsb=       obstype == 'hsb      '
147 
148   icount=0
149   if(airs)then
150      sensor_id = 11
151      step   = 1.1_r_kind
152      start = -48.9_r_kind
153      senname = 'AIRS'
154      nchan  = N_AIRSCHAN
155      nchanr = N_AIRSCHAN
156   else if(eos_amsua)then
157      sensor_id = 3
158      step   = three + one/three
159      start  = -48.33_r_kind
160      senname = 'AMSU'
161      nchan  = N_AMSUCHAN
162      nchanr = N_AMSUCHAN
163   else if(hsb)then
164      sensor_id = 12
165      step   = 1.1_r_kind
166      start  = -48.95_r_kind
167      senname = 'HSB'
168      nchan  = N_HSBCHAN
169      nchanr = N_HSBCHAN+1
170   end if
171   spotname = trim(senname)//'SPOT'
172   channame = trim(senname)//'CHAN'
173  
174       do inst = 1, rtminit_nsensor
175         if (    platform_id  == rtminit_platform(inst) &
176           .and. satellite_id == rtminit_satid(inst)    &
177           .and. sensor_id    == rtminit_sensor(inst)    ) then
178             exit
179         end if
180       end do
181 
182       if ( inst == rtminit_nsensor .and.           &
183            platform_id  /= rtminit_platform(inst)  &
184           .or. satellite_id /= rtminit_satid(inst) &
185           .or. sensor_id /= rtminit_sensor(inst)  ) return
186 
187 !    1.0  Open BUFR table and BUFR file
188 !--------------------------------------------------------------
189   table_file = 'gmao_airs_bufr.tbl'      ! make table file name
190   inquire(file=table_file,exist=airstab)
191   if (airstab) then
192       if (print_detail_radiance) then
193          write(unit=message(1),fmt=*) &
194             'Reading BUFR Table A file: ',trim(table_file)
195          call da_warning(__FILE__,__LINE__,message(1:1))
196       end if
197       call da_get_unit(airs_table_unit)
198       open(unit=airs_table_unit,file=table_file,iostat = iost)
199       if (iost /= 0) then
200          call da_error(__FILE__,__LINE__, &
201             (/"Cannot open file "//trim(table_file)/))
202       end if
203   end if
204 
205 ! Open BUFR file
206   call da_get_unit(lnbufr)
207   open(unit=lnbufr,file=trim(infile),form='unformatted',iostat = iost)
208   if (iost /= 0) then
209      call da_warning(__FILE__,__LINE__, &
210         (/"Cannot open file "//trim(infile)/))
211      return
212   end if
213   if ( airstab ) then
214      call openbf(lnbufr,'IN',airs_table_unit)
215   else
216      call openbf(lnbufr,'IN',lnbufr)
217   end if
218   call datelen(10)
219 
220 
221 !   2.0  Read header
222 !---------------------------
223   call readmg(lnbufr,subset,idate,iret)
224 
225   iy = 0
226   im = 0
227   idd = 0
228   ihh = 0
229   if( iret /= 0 ) goto 1000     ! no data?
230 
231   write(unit=date,fmt='( i10)') idate
232   read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
233    if (print_detail_radiance) then
234       write(unit=message(1),fmt=*) 'Bufr file date is ', &
235          iy,im,idd,ihh
236       call da_message(message(1:1))
237    end if
238   
239 !   3.0 Loop over observations
240 !----------------------------
241 
242   allocate ( head )
243 !  allocate ( head % tb (1:nchan) )
244   nullify  ( head % next )
245   p => head
246   size = 0
247 
248   loop_obspoints: do
249 
250 !   3.1 Read headder
251 !-------------------------------
252      call readsb(lnbufr,iret)
253 
254      if( iret /=0 )then
255         call readmg(lnbufr,subset,idate,iret)
256         if( iret /= 0 ) exit loop_obspoints     ! end of file
257         cycle loop_obspoints
258      end if
259   
260 !   3.2 Read AQUASPOT
261 !------------------------
262      call ufbseq(lnbufr,aquaspot,N_AQUASPOT_LIST,1,iret,'AQUASPOT')
263 
264 !   3.3 Read AIRSSPOT or AMSUSPOT or HSBSPOT
265 !-------------------------------------------------
266      call ufbseq(lnbufr,airsspot,N_AIRSSPOT_LIST,1,iret,spotname)
267 
268 !   3.4 Read AIRSCHAN or AMSUCHAN or HSBCHAN
269 !-------------------------------------------
270      call ufbseq(lnbufr,airschan,N_AIRSCHAN_LIST,N_MAXCHAN,iret,channame)
271 
272      if (iret /= nchanr) then
273         write(unit=message(1),fmt=*) &
274             'Cannot read ', channame, &
275             ' bufr data:', &
276             iret, ' ch data is read instead of', nchanr
277         call da_warning(__FILE__,__LINE__,message(1:1))
278         cycle loop_obspoints
279      end if
280 
281 !   4.0  Check observing position (lat/lon)
282 !      QC1:  juge if data is in the domain, 
283 !            read next record if not
284 !------------------------------------------
285      if( abs(airsspot%clath) > R90  .or. &
286           abs(airsspot%clonh) > R360 .or. &
287           (abs(airsspot%clath) == R90 .and. airsspot%clonh /= ZERO) )then
288         cycle loop_obspoints
289      end if
290 
291 !    Retrieve observing position
292      if(airsspot%clonh >= R360) then
293         airsspot%clonh = airsspot%clonh - R360
294 !     else if(airsspot%clonh < ZERO) then
295 !        airsspot%clonh = airsspot%clonh + R360
296      end if
297 
298         info%lat  = airsspot%clath
299         info%lon  = airsspot%clonh 
300         call da_ll_to_xy (info, loc, xp, outside, inside_halo )
301         if (outside) cycle LOOP_OBSPOinTS
302 
303 !  4.1  Check obs time
304 !-------------------------------------
305      idate5(1) = airsspot%year ! year
306      idate5(2) = airsspot%mnth ! month
307      idate5(3) = airsspot%days ! day
308      idate5(4) = airsspot%hour ! hour
309      idate5(5) = airsspot%minu ! minute
310      idate5(6) = airsspot%seco ! second
311 
312      if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. &
313           idate5(2) <    1 .or. idate5(2) >   12 .or. &
314           idate5(3) <    1 .or. idate5(3) >   31 .or. &
315           idate5(4) <    0 .or. idate5(4) >   24 .or. &
316           idate5(5) <    0 .or. idate5(5) >   60 .or. &
317           idate5(6) <    0 .or. idate5(6) >   60 ) then
318         cycle loop_obspoints
319      end if
320 
321 !  QC3: time consistency check with the background date
322 
323       if (idate5(1) .LE. 99) then
324         if (idate5(1) .LT. 78) then
325           idate5(1) = idate5(1) + 2000
326         else
327           idate5(1) = idate5(1) + 1900
328         end if
329       end if
330 
331       call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
332 
333       if ( obs_time < time_slots(0) .or.  &
334            obs_time >= time_slots(num_fgat_time) ) cycle
335 
336 !  3.2.1   determine FGAT index ifgat
337 !
338        do ifgat=1,num_fgat_time
339            if ( obs_time >= time_slots(ifgat-1) .and.  &
340                 obs_time  < time_slots(ifgat) ) exit
341        end do
342 
343        write(unit=info%date_char, &
344          fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
345          idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
346          ':', idate5(5), ':', idate5(6)
347 
348        info%elv = 0.  !aquaspot%selv
349 
350 !  4.2  Check observational info
351 !-------------------------------------------------------
352      if( airsspot%fovn <    0._r_kind .or. airsspot%fovn > 100._r_kind .or. &
353           airsspot%saza < -360._r_kind .or. airsspot%saza > 360._r_kind .or. &
354           aquaspot%soza < -180._r_kind .or. aquaspot%soza > 180._r_kind )then
355          write(unit=message(1),fmt=*) &
356             'Cannot read ', channame, ' bufr data:', &
357             ' strange obs info(fov,saza,soza):', &
358             airsspot%fovn, airsspot%saza, aquaspot%soza
359         call da_warning(__FILE__,__LINE__,message(1:1))
360         cycle loop_obspoints
361      end if
362 
363 !    Retrieve observing info
364      ifov = int( airsspot%fovn + POinT001 )
365      sza  = abs(airsspot%saza)
366 !     if( ((airs .or. hsb) .and. ifov <= 45) .or. &
367 !          ( eos_amsua     .and. ifov <= 15) )then
368 !        sza = - sza
369 !     end if
370 
371 !     airdata(6) = (start + float(ifov-1)*step)  ! look angle (deg)
372 !     airdata(9) = ZERO                          ! surface height
373 !     airdata(10)= POinT001                      ! land sea mask
374 
375 
376 !   4.3 Retrieve Tb
377 !-----------------------
378      iflag = .false.
379   
380      do l=1,nchan
381         airdata(l+nreal) = airschan(l)%tmbrst	         ! brightness temerature
382         if( airdata(l+nreal) > 0._r_kind .and. airdata(l+nreal) < 500._r_kind )then
383            iflag = .true.
384         else
385            airdata(l+nreal) = missing_r
386         end if
387      end do
388 
389      if ( .not. iflag )then
390         write(unit=message(1),fmt=*) &
391           'Error in reading ', channame, ' bufr data:', &
392           ' all tb data is missing'
393         call da_warning(__FILE__,__LINE__,message(1:1))
394         cycle loop_obspoints
395      end if
396 
397      nobs=nobs + 1
398 
399 !  4.0   assign information to sequential radiance structure
400 !--------------------------------------------------------------------------
401    allocate ( p % tb_inv (1:nchan) )
402    p%info             = info
403    p%loc              = loc
404    p%landsea_mask     = POinT001
405    p%scanline         = int(aquaspot%slnm + POinT001)
406    p%scanpos          = ifov
407    p%satzen           = sza
408    p%satazi           = (start + float(ifov-1)*step)  ! look angle (deg) ! airsspot%bearaz
409    p%solzen           = aquaspot%soza
410    p%solazi           = aquaspot%solazi
411    p%tb_inv(1:nchan) = airdata(nreal+1:nreal+nchan)
412    p%sensor_index         = inst
413    p%ifgat                = ifgat
414 
415    size = size + 1
416    allocate ( p%next, stat=error)   ! add next data
417    if (error /= 0 ) then
418       call da_error(__FILE__,__LINE__, &
419           (/"Cannot allocate radiance structure"/))
420    end if
421 
422    p => p%next
423    nullify (p%next)
424 
425   end do loop_obspoints
426 
427    iv%total_obs = iv%total_obs + size
428    iv%total_rad_pixel   = iv%total_rad_pixel + size
429    iv%total_rad_channel = iv%total_rad_channel + size*nchan
430 
431 !  5.0 allocate innovation radiance structure
432 !----------------------------------------------------------------
433 !  do i = 1, iv%num_inst
434    if ( nobs > 0 ) then
435       iv%instid(inst)%num_rad = nobs
436       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
437         'Allocating space for radiance innov structure', &
438          inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
439       if (rtm_option == rtm_option_crtm) then
440          iv%instid(inst)%nlevels = xp%kme-xp%kms+1
441          iv%instid(inst)%nchannels=nchan
442       end if
443    end if
444 !  end do
445 
446 !   6.0 assign sequential structure to innovation structure
447 !-------------------------------------------------------------
448 #ifdef RTTOV
449       if (rtm_option == rtm_option_rttov) then
450          call rttov_setupchan(1, nchan, coefs(inst), &   ! in
451             iv%instid(inst)%nfrequencies,iv%instid(inst)%nchannels, &
452             iv%instid(inst)%nbtout)      ! out
453       end if
454 #endif
455 
456   n = 0
457   p => head
458 
459       allocate (iv%instid(inst)%info(iv%instid(inst)%num_rad))
460       allocate (iv%instid(inst)%loc(iv%instid(inst)%num_rad))
461       allocate (iv%instid(inst)%loc_i(iv%instid(inst)%num_rad))
462       allocate (iv%instid(inst)%loc_j(iv%instid(inst)%num_rad))
463       allocate (iv%instid(inst)%loc_k(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
464       allocate (iv%instid(inst)%loc_dz(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
465       allocate (iv%instid(inst)%loc_dzm(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
466       allocate (iv%instid(inst)%zk(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
467       allocate (iv%instid(inst)%t(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
468       allocate (iv%instid(inst)%mr(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
469       allocate (iv%instid(inst)%loc_dx(iv%instid(inst)%num_rad))
470       allocate (iv%instid(inst)%loc_dy(iv%instid(inst)%num_rad))
471       allocate (iv%instid(inst)%loc_dxm(iv%instid(inst)%num_rad))
472       allocate (iv%instid(inst)%loc_dym(iv%instid(inst)%num_rad))
473       allocate (iv%instid(inst)%tm(xp%kms:xp%kme,iv%instid(inst)%num_rad))
474       allocate (iv%instid(inst)%qm(xp%kms:xp%kme,iv%instid(inst)%num_rad))
475       allocate (iv%instid(inst)%qrn(xp%kms:xp%kme,iv%instid(inst)%num_rad))
476       allocate (iv%instid(inst)%qcw(xp%kms:xp%kme,iv%instid(inst)%num_rad))
477       allocate (iv%instid(inst)%qci(xp%kms:xp%kme,iv%instid(inst)%num_rad))
478       allocate (iv%instid(inst)%qsn(xp%kms:xp%kme,iv%instid(inst)%num_rad))
479       allocate (iv%instid(inst)%qgr(xp%kms:xp%kme,iv%instid(inst)%num_rad))
480       allocate (iv%instid(inst)%pm(xp%kms:xp%kme,iv%instid(inst)%num_rad))
481       allocate (iv%instid(inst)%pf(0:xp%kme,iv%instid(inst)%num_rad))
482       allocate (iv%instid(inst)%u10(iv%instid(inst)%num_rad))
483       allocate (iv%instid(inst)%v10(iv%instid(inst)%num_rad))
484       allocate (iv%instid(inst)%t2m(iv%instid(inst)%num_rad))
485       allocate (iv%instid(inst)%q2m(iv%instid(inst)%num_rad))
486       allocate (iv%instid(inst)%mr2m(iv%instid(inst)%num_rad))
487       allocate (iv%instid(inst)%psfc(iv%instid(inst)%num_rad))
488       allocate (iv%instid(inst)%ts(iv%instid(inst)%num_rad))
489       allocate (iv%instid(inst)%smois(iv%instid(inst)%num_rad))
490       allocate (iv%instid(inst)%tslb(iv%instid(inst)%num_rad))
491       allocate (iv%instid(inst)%snowh(iv%instid(inst)%num_rad))
492       allocate (iv%instid(inst)%isflg(iv%instid(inst)%num_rad))
493       allocate (iv%instid(inst)%landsea_mask(iv%instid(inst)%num_rad))
494       allocate (iv%instid(inst)%elevation(iv%instid(inst)%num_rad))
495       allocate (iv%instid(inst)%vegfra(iv%instid(inst)%num_rad))
496       allocate (iv%instid(inst)%vegtyp(iv%instid(inst)%num_rad))
497       allocate (iv%instid(inst)%soiltyp(iv%instid(inst)%num_rad))
498       allocate (iv%instid(inst)%clwp(iv%instid(inst)%num_rad))
499       allocate (iv%instid(inst)%ps(iv%instid(inst)%num_rad))
500       allocate (iv%instid(inst)%tb_xb(nchan,iv%instid(inst)%num_rad))
501       allocate (iv%instid(inst)%tb_qc(nchan,iv%instid(inst)%num_rad))
502       allocate (iv%instid(inst)%tb_inv(nchan,iv%instid(inst)%num_rad))
503       allocate (iv%instid(inst)%tb_error(nchan,iv%instid(inst)%num_rad))
504       allocate (iv%instid(inst)%emiss(iv%instid(inst)%nchannels,iv%instid(inst)%num_rad))
505       allocate(iv%instid(inst)%scanpos(iv%instid(inst)%num_rad))
506       allocate(iv%instid(inst)%scanline(iv%instid(inst)%num_rad))
507       allocate(iv%instid(inst)%ifgat(iv%instid(inst)%num_rad))
508       allocate(iv%instid(inst)%cloud_flag(nchan,iv%instid(inst)%num_rad))
509       allocate(iv%instid(inst)%satzen(iv%instid(inst)%num_rad))
510       allocate(iv%instid(inst)%satazi(iv%instid(inst)%num_rad))
511       allocate(iv%instid(inst)%solzen(iv%instid(inst)%num_rad))
512       allocate(iv%instid(inst)%solazi(iv%instid(inst)%num_rad))
513       allocate(iv%instid(inst)%proc_domain(iv%instid(inst)%num_rad))
514       if (rtm_option == rtm_option_crtm) then
515         allocate(iv%instid(inst)%water_coverage(iv%instid(inst)%num_rad))
516         allocate(iv%instid(inst)%land_coverage(iv%instid(inst)%num_rad))
517         allocate(iv%instid(inst)%ice_coverage(iv%instid(inst)%num_rad))
518         allocate(iv%instid(inst)%snow_coverage(iv%instid(inst)%num_rad))
519         allocate(iv%instid(inst)%ps_jacobian(nchan,iv%instid(inst)%num_rad))
520         allocate(iv%instid(inst)%t_jacobian(nchan,xp%kte,iv%instid(inst)%num_rad))
521         allocate(iv%instid(inst)%q_jacobian(nchan,xp%kte,iv%instid(inst)%num_rad))
522       end if
523 
524   do i = 1, size
525 !   inst = p%sensor_index
526    n = n + 1
527 
528    iv%instid(inst)%info(n)   =  p%info
529    iv%instid(inst)%loc_i(n)        = p%loc%i
530    iv%instid(inst)%loc_j(n)        = p%loc%j
531    iv%instid(inst)%loc_k(:,n)      = 0
532    iv%instid(inst)%loc_dx(n)       = p%loc%dx
533    iv%instid(inst)%loc_dy(n)       = p%loc%dy
534    iv%instid(inst)%loc_dz(:,n)     = 0.0
535    iv%instid(inst)%loc_dxm(n)      = p%loc%dxm
536    iv%instid(inst)%loc_dym(n)      = p%loc%dym
537    iv%instid(inst)%loc_dzm(:,n)    = 0.0
538 
539       ! z done in da_get_innov_vector_rad
540       iv%instid(inst)%t(:,n)          = 0.0
541       iv%instid(inst)%mr(:,n)         = 0.0
542       iv%instid(inst)%tm(:,n)         = 0.0
543       iv%instid(inst)%qm(:,n)         = 0.0
544       iv%instid(inst)%qrn(:,n)        = 0.0
545       iv%instid(inst)%qcw(:,n)        = 0.0
546       iv%instid(inst)%qci(:,n)        = 0.0
547       iv%instid(inst)%qsn(:,n)        = 0.0
548       iv%instid(inst)%qgr(:,n)        = 0.0
549       iv%instid(inst)%pm(:,n)         = 0.0
550       iv%instid(inst)%u10(n)          = 0.0
551       iv%instid(inst)%v10(n)          = 0.0
552       iv%instid(inst)%t2m(n)          = 0.0
553       iv%instid(inst)%q2m(n)          = 0.0
554       iv%instid(inst)%mr2m(n)         = 0.0
555       iv%instid(inst)%psfc(n)         = 0.0
556       iv%instid(inst)%ts(n)           = 0.0
557       iv%instid(inst)%smois(n)        = 0.0
558       iv%instid(inst)%tslb(n)         = 0.0
559       iv%instid(inst)%snowh(n)        = 0.0
560       iv%instid(inst)%isflg(n)        = 0
561       iv%instid(inst)%soiltyp(n)      = 0.0
562       iv%instid(inst)%landsea_mask(n) = p%landsea_mask
563       iv%instid(inst)%elevation(n)    = 0.0
564       iv%instid(inst)%vegfra(n)       = 0.0
565       iv%instid(inst)%vegtyp(n)       = 0.0
566       iv%instid(inst)%clwp(n)         = 0.0
567       iv%instid(inst)%ps(n)           = 0.0
568       iv%instid(inst)%tb_xb(:,n)      = 0
569       iv%instid(inst)%tb_inv(:,n)     = p%tb_inv(:)
570       iv%instid(inst)%tb_qc(:,n)      = 0
571       iv%instid(inst)%tb_error(:,n)   = 500.
572       iv%instid(inst)%emiss(:,n)      = 0.0
573       iv%instid(inst)%scanpos(n)      = p%scanpos
574       ! iv%instid(inst)%scanline(n)    = p%scanline
575       iv%instid(inst)%scanline(n)     = 0
576       iv%instid(inst)%ifgat(n)        = p%ifgat
577       iv%instid(inst)%cloud_flag(:,n) = 1        ! no cloud
578       iv%instid(inst)%satzen(n)       = p%satzen
579       iv%instid(inst)%satazi(n)       = p%satazi
580       iv%instid(inst)%solzen(n)       = p%solzen
581       ! iv%instid(inst)%solazi(n)     = p%solazi
582       iv%instid(inst)%solazi(n)       = 0.0
583       iv%instid(inst)%proc_domain(n)  = .false.
584 
585    current => p
586    p => p%next
587 
588 ! free current data
589    deallocate ( current % tb_inv )
590    deallocate ( current )
591 
592  end do
593 
594 1000 continue
595    call closbf(lnbufr)
596    close(lnbufr)
597    call da_free_unit(lnbufr)
598    if (airstab) then
599       close(airs_table_unit)
600       call da_free_unit(airs_table_unit)
601    end if
602 
603    if (trace_use) call da_trace_exit("da_read_bufrairs")
604 #endif
605 
606 end subroutine da_read_bufrairs