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 !  0.0  Initialize variables
137 !-----------------------------------
138   platform_id  = 9   ! eos series
139   satellite_id = 2   ! eos-2
140   nreal  = maxinfo
141   nobs   = 0
142   airs=      obstype == 'airs     '
143   eos_amsua= obstype == 'eos_amsua'
144   hsb=       obstype == 'hsb      '
145 
146   icount=0
147   if(airs)then
148      sensor_id = 11
149      step   = 1.1_r_kind
150      start = -48.9_r_kind
151      senname = 'AIRS'
152      nchan  = N_AIRSCHAN
153      nchanr = N_AIRSCHAN
154   else if(eos_amsua)then
155      sensor_id = 3
156      step   = three + one/three
157      start  = -48.33_r_kind
158      senname = 'AMSU'
159      nchan  = N_AMSUCHAN
160      nchanr = N_AMSUCHAN
161   else if(hsb)then
162      sensor_id = 12
163      step   = 1.1_r_kind
164      start  = -48.95_r_kind
165      senname = 'HSB'
166      nchan  = N_HSBCHAN
167      nchanr = N_HSBCHAN+1
168   end if
169   spotname = trim(senname)//'SPOT'
170   channame = trim(senname)//'CHAN'
171  
172       do inst = 1, rtminit_nsensor
173         if (    platform_id  == rtminit_platform(inst) &
174           .and. satellite_id == rtminit_satid(inst)    &
175           .and. sensor_id    == rtminit_sensor(inst)    ) then
176             exit
177         end if
178       end do
179 
180       if ( inst == rtminit_nsensor .and.           &
181            platform_id  /= rtminit_platform(inst)  &
182           .or. satellite_id /= rtminit_satid(inst) &
183           .or. sensor_id /= rtminit_sensor(inst)  ) return
184 
185 !    1.0  Open BUFR table and BUFR file
186 !--------------------------------------------------------------
187   table_file = 'gmao_airs_bufr.tbl'      ! make table file name
188   inquire(file=table_file,exist=airstab)
189   if (airstab) then
190       if (print_detail_rad) then
191          write(unit=message(1),fmt=*) &
192             'Reading BUFR Table A file: ',trim(table_file)
193          call da_warning(__FILE__,__LINE__,message(1:1))
194       end if
195       call da_get_unit(airs_table_unit)
196       open(unit=airs_table_unit,file=table_file,iostat = iost)
197       if (iost /= 0) then
198          call da_error(__FILE__,__LINE__, &
199             (/"Cannot open file "//trim(table_file)/))
200       end if
201   end if
202 
203 ! Open BUFR file
204   call da_get_unit(lnbufr)
205   open(unit=lnbufr,file=trim(infile),form='unformatted',iostat = iost)
206   if (iost /= 0) then
207      call da_warning(__FILE__,__LINE__, &
208         (/"Cannot open file "//trim(infile)/))
209      return
210   end if
211   if ( airstab ) then
212      call openbf(lnbufr,'IN',airs_table_unit)
213   else
214      call openbf(lnbufr,'IN',lnbufr)
215   end if
216   call datelen(10)
217 
218 
219 !   2.0  Read header
220 !---------------------------
221   call readmg(lnbufr,subset,idate,iret)
222 
223   iy = 0
224   im = 0
225   idd = 0
226   ihh = 0
227   if( iret /= 0 ) goto 1000     ! no data?
228 
229   write(unit=date,fmt='( i10)') idate
230   read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
231    if (print_detail_rad) then
232       write(unit=message(1),fmt=*) 'Bufr file date is ', &
233          iy,im,idd,ihh
234       call da_message(message(1:1))
235    end if
236   
237 !   3.0 Loop over observations
238 !----------------------------
239 
240   allocate ( head )
241 !  allocate ( head % tb (1:nchan) )
242   nullify  ( head % next )
243   p => head
244   size = 0
245 
246   loop_obspoints: do
247 
248 !   3.1 Read headder
249 !-------------------------------
250      call readsb(lnbufr,iret)
251 
252      if( iret /=0 )then
253         call readmg(lnbufr,subset,idate,iret)
254         if( iret /= 0 ) exit loop_obspoints     ! end of file
255         cycle loop_obspoints
256      end if
257   
258 !   3.2 Read AQUASPOT
259 !------------------------
260      call ufbseq(lnbufr,aquaspot,N_AQUASPOT_LIST,1,iret,'AQUASPOT')
261 
262 !   3.3 Read AIRSSPOT or AMSUSPOT or HSBSPOT
263 !-------------------------------------------------
264      call ufbseq(lnbufr,airsspot,N_AIRSSPOT_LIST,1,iret,spotname)
265 
266 !   3.4 Read AIRSCHAN or AMSUCHAN or HSBCHAN
267 !-------------------------------------------
268      call ufbseq(lnbufr,airschan,N_AIRSCHAN_LIST,N_MAXCHAN,iret,channame)
269 
270      if (iret /= nchanr) then
271         write(unit=message(1),fmt=*) &
272             'Cannot read ', channame, &
273             ' bufr data:', &
274             iret, ' ch data is read instead of', nchanr
275         call da_warning(__FILE__,__LINE__,message(1:1))
276         cycle loop_obspoints
277      end if
278 
279 !   4.0  Check observing position (lat/lon)
280 !      QC1:  juge if data is in the domain, 
281 !            read next record if not
282 !------------------------------------------
283      if( abs(airsspot%clath) > R90  .or. &
284           abs(airsspot%clonh) > R360 .or. &
285           (abs(airsspot%clath) == R90 .and. airsspot%clonh /= ZERO) )then
286         cycle loop_obspoints
287      end if
288 
289 !    Retrieve observing position
290      if(airsspot%clonh >= R360) then
291         airsspot%clonh = airsspot%clonh - R360
292 !     else if(airsspot%clonh < ZERO) then
293 !        airsspot%clonh = airsspot%clonh + R360
294      end if
295 
296         info%lat  = airsspot%clath
297         info%lon  = airsspot%clonh 
298         call da_ll_to_xy (info, loc, xp, outside, inside_halo )
299         if (outside) cycle LOOP_OBSPOinTS
300 
301 !  4.1  Check obs time
302 !-------------------------------------
303      idate5(1) = airsspot%year ! year
304      idate5(2) = airsspot%mnth ! month
305      idate5(3) = airsspot%days ! day
306      idate5(4) = airsspot%hour ! hour
307      idate5(5) = airsspot%minu ! minute
308      idate5(6) = airsspot%seco ! second
309 
310      if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. &
311           idate5(2) <    1 .or. idate5(2) >   12 .or. &
312           idate5(3) <    1 .or. idate5(3) >   31 .or. &
313           idate5(4) <    0 .or. idate5(4) >   24 .or. &
314           idate5(5) <    0 .or. idate5(5) >   60 .or. &
315           idate5(6) <    0 .or. idate5(6) >   60 ) then
316         cycle loop_obspoints
317      end if
318 
319 !  QC3: time consistency check with the background date
320 
321       if (idate5(1) .LE. 99) then
322         if (idate5(1) .LT. 78) then
323           idate5(1) = idate5(1) + 2000
324         else
325           idate5(1) = idate5(1) + 1900
326         end if
327       end if
328 
329       call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
330 
331       if ( obs_time < time_slots(0) .or.  &
332            obs_time >= time_slots(num_fgat_time) ) cycle
333 
334 !  3.2.1   determine FGAT index ifgat
335 !
336        do ifgat=1,num_fgat_time
337            if ( obs_time >= time_slots(ifgat-1) .and.  &
338                 obs_time  < time_slots(ifgat) ) exit
339        end do
340 
341        write(unit=info%date_char, &
342          fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)')  &
343          idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
344          ':', idate5(5), ':', idate5(6)
345 
346        info%elv = 0.  !aquaspot%selv
347 
348 !  4.2  Check observational info
349 !-------------------------------------------------------
350      if( airsspot%fovn <    0._r_kind .or. airsspot%fovn > 100._r_kind .or. &
351           airsspot%saza < -360._r_kind .or. airsspot%saza > 360._r_kind .or. &
352           aquaspot%soza < -180._r_kind .or. aquaspot%soza > 180._r_kind )then
353          write(unit=message(1),fmt=*) &
354             'Cannot read ', channame, ' bufr data:', &
355             ' strange obs info(fov,saza,soza):', &
356             airsspot%fovn, airsspot%saza, aquaspot%soza
357         call da_warning(__FILE__,__LINE__,message(1:1))
358         cycle loop_obspoints
359      end if
360 
361 !    Retrieve observing info
362      ifov = int( airsspot%fovn + POinT001 )
363      sza  = abs(airsspot%saza)
364 !     if( ((airs .or. hsb) .and. ifov <= 45) .or. &
365 !          ( eos_amsua     .and. ifov <= 15) )then
366 !        sza = - sza
367 !     end if
368 
369 !     airdata(6) = (start + float(ifov-1)*step)  ! look angle (deg)
370 !     airdata(9) = ZERO                          ! surface height
371 !     airdata(10)= POinT001                      ! land sea mask
372 
373 
374 !   4.3 Retrieve Tb
375 !-----------------------
376      iflag = .false.
377   
378      do l=1,nchan
379         airdata(l+nreal) = airschan(l)%tmbrst	         ! brightness temerature
380         if( airdata(l+nreal) > 0._r_kind .and. airdata(l+nreal) < 500._r_kind )then
381            iflag = .true.
382         else
383            airdata(l+nreal) = missing_r
384         end if
385      end do
386 
387      if ( .not. iflag )then
388         write(unit=message(1),fmt=*) &
389           'Error in reading ', channame, ' bufr data:', &
390           ' all tb data is missing'
391         call da_warning(__FILE__,__LINE__,message(1:1))
392         cycle loop_obspoints
393      end if
394 
395      nobs=nobs + 1
396 
397 !  4.0   assign information to sequential radiance structure
398 !--------------------------------------------------------------------------
399    allocate ( p % tb_inv (1:nchan) )
400    p%info             = info
401    p%loc              = loc
402    p%landsea_mask     = POinT001
403    p%scanline         = int(aquaspot%slnm + POinT001)
404    p%scanpos          = ifov
405    p%satzen           = sza
406    p%satazi           = (start + float(ifov-1)*step)  ! look angle (deg) ! airsspot%bearaz
407    p%solzen           = aquaspot%soza
408    p%solazi           = aquaspot%solazi
409    p%tb_inv(1:nchan) = airdata(nreal+1:nreal+nchan)
410    p%sensor_index         = inst
411    p%ifgat                = ifgat
412 
413    size = size + 1
414    allocate ( p%next, stat=error)   ! add next data
415    if (error /= 0 ) then
416       call da_error(__FILE__,__LINE__, &
417           (/"Cannot allocate radiance structure"/))
418    end if
419 
420    p => p%next
421    nullify (p%next)
422 
423   end do loop_obspoints
424 
425    iv%total_obs = iv%total_obs + size
426    iv%total_rad_pixel   = iv%total_rad_pixel + size
427    iv%total_rad_channel = iv%total_rad_channel + size*nchan
428 
429 !  5.0 allocate innovation radiance structure
430 !----------------------------------------------------------------
431 !  do i = 1, iv%num_inst
432    if ( nobs > 0 ) then
433       iv%instid(inst)%num_rad = nobs
434       write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
435         'Allocating space for radiance innov structure', &
436          inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
437       if (rtm_option == rtm_option_crtm) then
438          iv%instid(inst)%nlevels = xp%kme-xp%kms+1
439          iv%instid(inst)%nchannels=nchan
440       end if
441    end if
442 !  end do
443 
444 !   6.0 assign sequential structure to innovation structure
445 !-------------------------------------------------------------
446 #ifdef RTTOV
447       if (rtm_option == rtm_option_rttov) then
448          call rttov_setupchan(1, nchan, coefs(inst), &   ! in
449             iv%instid(inst)%nfrequencies,iv%instid(inst)%nchannels, &
450             iv%instid(inst)%nbtout)      ! out
451       end if
452 #endif
453 
454   n = 0
455   p => head
456 
457       allocate (iv%instid(inst)%info(iv%instid(inst)%num_rad))
458       allocate (iv%instid(inst)%loc(iv%instid(inst)%num_rad))
459       allocate (iv%instid(inst)%loc_i(iv%instid(inst)%num_rad))
460       allocate (iv%instid(inst)%loc_j(iv%instid(inst)%num_rad))
461 
462       allocate (iv%instid(inst)%loc_k(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
463       allocate (iv%instid(inst)%loc_dz(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
464       allocate (iv%instid(inst)%loc_dzm(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
465       allocate (iv%instid(inst)%zk(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
466       allocate (iv%instid(inst)%t(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
467       allocate (iv%instid(inst)%mr(iv%instid(inst)%nlevels,iv%instid(inst)%num_rad))
468       allocate (iv%instid(inst)%loc_dx(iv%instid(inst)%num_rad))
469       allocate (iv%instid(inst)%loc_dy(iv%instid(inst)%num_rad))
470       allocate (iv%instid(inst)%loc_dxm(iv%instid(inst)%num_rad))
471       allocate (iv%instid(inst)%loc_dym(iv%instid(inst)%num_rad))
472       allocate (iv%instid(inst)%tm(xp%kms:xp%kme,iv%instid(inst)%num_rad))
473       allocate (iv%instid(inst)%qm(xp%kms:xp%kme,iv%instid(inst)%num_rad))
474       allocate (iv%instid(inst)%qrn(xp%kms:xp%kme,iv%instid(inst)%num_rad))
475       allocate (iv%instid(inst)%qcw(xp%kms:xp%kme,iv%instid(inst)%num_rad))
476       allocate (iv%instid(inst)%qci(xp%kms:xp%kme,iv%instid(inst)%num_rad))
477       allocate (iv%instid(inst)%qsn(xp%kms:xp%kme,iv%instid(inst)%num_rad))
478       allocate (iv%instid(inst)%qgr(xp%kms:xp%kme,iv%instid(inst)%num_rad))
479       allocate (iv%instid(inst)%pm(xp%kms:xp%kme,iv%instid(inst)%num_rad))
480       allocate (iv%instid(inst)%pf(0:xp%kme,iv%instid(inst)%num_rad))
481       allocate (iv%instid(inst)%u10(iv%instid(inst)%num_rad))
482       allocate (iv%instid(inst)%v10(iv%instid(inst)%num_rad))
483       allocate (iv%instid(inst)%t2m(iv%instid(inst)%num_rad))
484       allocate (iv%instid(inst)%q2m(iv%instid(inst)%num_rad))
485       allocate (iv%instid(inst)%mr2m(iv%instid(inst)%num_rad))
486       allocate (iv%instid(inst)%psfc(iv%instid(inst)%num_rad))
487       allocate (iv%instid(inst)%ts(iv%instid(inst)%num_rad))
488       allocate (iv%instid(inst)%smois(iv%instid(inst)%num_rad))
489       allocate (iv%instid(inst)%tslb(iv%instid(inst)%num_rad))
490       allocate (iv%instid(inst)%snowh(iv%instid(inst)%num_rad))
491       allocate (iv%instid(inst)%isflg(iv%instid(inst)%num_rad))
492       allocate (iv%instid(inst)%landsea_mask(iv%instid(inst)%num_rad))
493       allocate (iv%instid(inst)%elevation(iv%instid(inst)%num_rad))
494       allocate (iv%instid(inst)%vegfra(iv%instid(inst)%num_rad))
495       allocate (iv%instid(inst)%vegtyp(iv%instid(inst)%num_rad))
496       allocate (iv%instid(inst)%soiltyp(iv%instid(inst)%num_rad))
497       allocate (iv%instid(inst)%clwp(iv%instid(inst)%num_rad))
498       allocate (iv%instid(inst)%ps(iv%instid(inst)%num_rad))
499       allocate (iv%instid(inst)%tb_xb(nchan,iv%instid(inst)%num_rad))
500       allocate (iv%instid(inst)%tb_qc(nchan,iv%instid(inst)%num_rad))
501       allocate (iv%instid(inst)%tb_inv(nchan,iv%instid(inst)%num_rad))
502       allocate (iv%instid(inst)%tb_error(nchan,iv%instid(inst)%num_rad))
503       allocate (iv%instid(inst)%emiss(iv%instid(inst)%nchannels,iv%instid(inst)%num_rad))
504       allocate(iv%instid(inst)%scanpos(iv%instid(inst)%num_rad))
505       allocate(iv%instid(inst)%scanline(iv%instid(inst)%num_rad))
506       allocate(iv%instid(inst)%ifgat(iv%instid(inst)%num_rad))
507       allocate(iv%instid(inst)%cloud_flag(nchan,iv%instid(inst)%num_rad))
508       allocate(iv%instid(inst)%satzen(iv%instid(inst)%num_rad))
509       allocate(iv%instid(inst)%satazi(iv%instid(inst)%num_rad))
510       allocate(iv%instid(inst)%solzen(iv%instid(inst)%num_rad))
511       allocate(iv%instid(inst)%solazi(iv%instid(inst)%num_rad))
512       allocate(iv%instid(inst)%proc_domain(iv%instid(inst)%num_rad))
513       if (rtm_option == rtm_option_crtm) then
514         allocate(iv%instid(inst)%water_coverage(iv%instid(inst)%num_rad))
515         allocate(iv%instid(inst)%land_coverage(iv%instid(inst)%num_rad))
516         allocate(iv%instid(inst)%ice_coverage(iv%instid(inst)%num_rad))
517         allocate(iv%instid(inst)%snow_coverage(iv%instid(inst)%num_rad))
518         !allocate(iv%instid(inst)%ps_jacobian(nchan,iv%instid(inst)%num_rad))
519         !allocate(iv%instid(inst)%t_jacobian(nchan,xp%kte,iv%instid(inst)%num_rad))
520         !allocate(iv%instid(inst)%q_jacobian(nchan,xp%kte,iv%instid(inst)%num_rad))
521       end if
522 
523   do i = 1, size
524 !   inst = p%sensor_index
525    n = n + 1
526 
527    iv%instid(inst)%info(n)   =  p%info
528    iv%instid(inst)%loc_i(n)        = p%loc%i
529    iv%instid(inst)%loc_j(n)        = p%loc%j
530    iv%instid(inst)%loc_k(:,n)      = 0
531    iv%instid(inst)%loc_dx(n)       = p%loc%dx
532    iv%instid(inst)%loc_dy(n)       = p%loc%dy
533    iv%instid(inst)%loc_dz(:,n)     = 0.0
534    iv%instid(inst)%loc_dxm(n)      = p%loc%dxm
535    iv%instid(inst)%loc_dym(n)      = p%loc%dym
536    iv%instid(inst)%loc_dzm(:,n)    = 0.0
537 
538       ! z done in da_get_innov_vector_rad
539       iv%instid(inst)%t(:,n)          = 0.0
540       iv%instid(inst)%mr(:,n)         = 0.0
541       iv%instid(inst)%tm(:,n)         = 0.0
542       iv%instid(inst)%qm(:,n)         = 0.0
543       iv%instid(inst)%qrn(:,n)        = 0.0
544       iv%instid(inst)%qcw(:,n)        = 0.0
545       iv%instid(inst)%qci(:,n)        = 0.0
546       iv%instid(inst)%qsn(:,n)        = 0.0
547       iv%instid(inst)%qgr(:,n)        = 0.0
548       iv%instid(inst)%pm(:,n)         = 0.0
549       iv%instid(inst)%u10(n)          = 0.0
550       iv%instid(inst)%v10(n)          = 0.0
551       iv%instid(inst)%t2m(n)          = 0.0
552       iv%instid(inst)%q2m(n)          = 0.0
553       iv%instid(inst)%mr2m(n)         = 0.0
554       iv%instid(inst)%psfc(n)         = 0.0
555       iv%instid(inst)%ts(n)           = 0.0
556       iv%instid(inst)%smois(n)        = 0.0
557       iv%instid(inst)%tslb(n)         = 0.0
558       iv%instid(inst)%snowh(n)        = 0.0
559       iv%instid(inst)%isflg(n)        = 0
560       iv%instid(inst)%soiltyp(n)      = 0.0
561       iv%instid(inst)%landsea_mask(n) = p%landsea_mask
562       iv%instid(inst)%elevation(n)    = 0.0
563       iv%instid(inst)%vegfra(n)       = 0.0
564       iv%instid(inst)%vegtyp(n)       = 0.0
565       iv%instid(inst)%clwp(n)         = 0.0
566       iv%instid(inst)%ps(n)           = 0.0
567       iv%instid(inst)%tb_xb(:,n)      = 0
568       iv%instid(inst)%tb_inv(:,n)     = p%tb_inv(:)
569       iv%instid(inst)%tb_qc(:,n)      = 0
570       iv%instid(inst)%tb_error(:,n)   = 500.
571       iv%instid(inst)%emiss(:,n)      = 0.0
572       iv%instid(inst)%scanpos(n)      = p%scanpos
573       ! iv%instid(inst)%scanline(n)    = p%scanline
574       iv%instid(inst)%scanline(n)     = 0
575       iv%instid(inst)%ifgat(n)        = p%ifgat
576       iv%instid(inst)%cloud_flag(:,n) = 1        ! no cloud
577       iv%instid(inst)%satzen(n)       = p%satzen
578       iv%instid(inst)%satazi(n)       = p%satazi
579       iv%instid(inst)%solzen(n)       = p%solzen
580       ! iv%instid(inst)%solazi(n)     = p%solazi
581       iv%instid(inst)%solazi(n)       = 0.0
582       iv%instid(inst)%proc_domain(n)  = .false.
583 
584    current => p
585    p => p%next
586 
587 ! free current data
588    deallocate ( current % tb_inv )
589    deallocate ( current )
590 
591  end do
592 
593 1000 continue
594    call closbf(lnbufr)
595    close(lnbufr)
596    call da_free_unit(lnbufr)
597    if (airstab) then
598       close(airs_table_unit)
599       call da_free_unit(airs_table_unit)
600    end if
601 
602 #endif
603 
604 end subroutine da_read_bufrairs