da_read_obs_bufrairs.inc
References to this file elsewhere.
1 subroutine da_read_obs_bufrairs(obstype,iv,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 (iv_type) ,intent (inout) :: iv
22
23 #ifdef BUFR
24
25 ! Number of channels for sensors in BUFR
26 integer(i_kind),parameter :: N_AIRSCHAN = 281 !- 281 subset ch out of 2378 ch for AIRS
27 integer(i_kind),parameter :: N_AMSUCHAN = 15
28 integer(i_kind),parameter :: N_HSBCHAN = 4
29 integer(i_kind),parameter :: N_MAXCHAN = 350
30 integer(i_kind),parameter :: maxinfo = 12
31
32
33 ! BUFR format for AQUASPOT
34 integer(i_kind),parameter :: N_AQUASPOT_LIST = 25
35 type aquaspot_list
36 sequence
37 real(r_double) :: said ! Satellite identifier
38 real(r_double) :: orbn ! Orbit number
39 real(r_double) :: slnm ! Scan line number
40 real(r_double) :: mjfc ! Major frame count
41 real(r_double) :: selv ! Height of station
42 real(r_double) :: soza ! Solar zenith angle
43 real(r_double) :: solazi ! Solar azimuth angle
44 real(r_double) :: intms(2,9) ! SATELLITE inSTRUMENT TEMPERATURES
45 end type aquaspot_list
46
47
48 ! BUFR format for AIRSSPOT
49 integer(i_kind),parameter :: N_AIRSSPOT_LIST = 12
50 type airsspot_list
51 sequence
52 real(r_double) :: siid ! Satellite instruments
53 real(r_double) :: year
54 real(r_double) :: mnth
55 real(r_double) :: days
56 real(r_double) :: hour
57 real(r_double) :: minu
58 real(r_double) :: seco
59 real(r_double) :: clath ! Latitude (high accuracy)
60 real(r_double) :: clonh ! Longitude (high accuracy)
61 real(r_double) :: saza ! Satellite zenith angle
62 real(r_double) :: bearaz ! Bearing or azimuth
63 real(r_double) :: fovn ! Field of view number
64 end type airsspot_list
65
66
67 ! BUFR format for AIRSCHAN
68 integer(i_kind),parameter :: N_AIRSCHAN_LIST = 4
69 type airschan_list
70 sequence
71 real(r_double) :: chnm ! Channel number
72 real(r_double) :: logrcw ! Log-10 of (Temperature-radiance central wavenumber
73 real(r_double) :: acqf ! Channel quality flags for ATOVS
74 real(r_double) :: tmbrst ! Brightness temperature
75 end type airschan_list
76
77
78 ! BUFR talble file sequencial number
79 character(len=512) :: table_file
80
81
82 ! Variables for BUFR IO
83 type(aquaspot_list) :: aquaspot
84 type(airsspot_list) :: airsspot
85 type(airschan_list) :: airschan(N_MAXCHAN)
86
87 real(r_kind) :: step, start
88 real(r_kind) :: airdata(N_AIRSCHAN+maxinfo)
89 character(len=8) :: subset
90 character(len=4) :: senname
91 character(len=8) :: spotname
92 character(len=8) :: channame
93 integer(i_kind) :: nchan,nchanr
94 integer(i_kind) :: iret
95
96
97 ! Work variables for time
98 integer(i_kind) :: idate, ifgat
99 integer(i_kind) :: idate5(6)
100 character(len=10) :: date
101 integer(i_kind) :: nmind
102 integer(i_kind) :: iy, im, idd, ihh
103 real :: obs_time
104
105
106 ! Other work variables
107 integer(i_kind) :: nreal, nobs, ityp,ityp2
108 integer(i_kind) :: itx, k, itt, iobsout
109 integer(i_kind),dimension(19)::icount
110 real(r_kind) :: rlat, rlon, dx, dy, dx1, dy1, sstx, dlon, dlat
111 real(r_kind) :: sza, timedif, pred, crit1
112 integer(i_kind) :: klat1, klon1, klatp1, klonp1
113 integer(i_kind) :: ifov,size
114 integer(i_kind) :: inst,platform_id,satellite_id,sensor_id
115 logical :: iflag,outside,inside_halo
116 integer(i_kind) :: i, l, n, error, airs_table_unit
117 integer :: iost, lnbufr
118 real(r_kind),allocatable,dimension(:,:):: airdata_all
119
120 ! Set standard parameters
121 real(r_kind) :: POinT001 = 0.001_r_kind
122 real(r_kind) :: POinT01 = 0.01_r_kind
123 real(r_kind) :: TEN = 10.0_r_kind
124 real(r_kind) :: R45 = 45.0_r_kind
125 real(r_kind) :: R60 = 60.0_r_kind
126 real(r_kind) :: R90 = 90.0_r_kind
127 real(r_kind) :: R180 = 180.0_r_kind
128 real(r_kind) :: R360 = 360.0_r_kind
129
130 logical :: airs, eos_amsua, hsb, airstab
131 type(info_type) :: info
132 type(model_loc_type) :: loc
133 type (datalink_type), pointer :: head, p, current
134
135 if (trace_use_dull) call da_trace_entry("da_read_obs_bufrairs")
136
137 ! 0.0 Initialize variables
138 !-----------------------------------
139 platform_id = 9 ! eos series
140 satellite_id = 2 ! eos-2
141 nreal = maxinfo
142 nobs = 0
143 airs= obstype == 'airs '
144 eos_amsua= obstype == 'eos_amsua'
145 hsb= obstype == 'hsb '
146
147 icount=0
148 if(airs)then
149 sensor_id = 11
150 step = 1.1_r_kind
151 start = -48.9_r_kind
152 senname = 'AIRS'
153 nchan = N_AIRSCHAN
154 nchanr = N_AIRSCHAN
155 else if(eos_amsua)then
156 sensor_id = 3
157 step = three + one/three
158 start = -48.33_r_kind
159 senname = 'AMSU'
160 nchan = N_AMSUCHAN
161 nchanr = N_AMSUCHAN
162 else if(hsb)then
163 sensor_id = 12
164 step = 1.1_r_kind
165 start = -48.95_r_kind
166 senname = 'HSB'
167 nchan = N_HSBCHAN
168 nchanr = N_HSBCHAN+1
169 end if
170 spotname = trim(senname)//'SPOT'
171 channame = trim(senname)//'CHAN'
172
173 do inst = 1, rtminit_nsensor
174 if ( platform_id == rtminit_platform(inst) &
175 .and. satellite_id == rtminit_satid(inst) &
176 .and. sensor_id == rtminit_sensor(inst) ) then
177 exit
178 end if
179 end do
180
181 if ( inst == rtminit_nsensor .and. &
182 platform_id /= rtminit_platform(inst) &
183 .or. satellite_id /= rtminit_satid(inst) &
184 .or. sensor_id /= rtminit_sensor(inst) ) return
185
186 ! 1.0 Open BUFR table and BUFR file
187 !--------------------------------------------------------------
188 table_file = 'gmao_airs_bufr.tbl' ! make table file name
189 inquire(file=table_file,exist=airstab)
190 if (airstab) then
191 if (print_detail_rad) then
192 write(unit=message(1),fmt=*) &
193 'Reading BUFR Table A file: ',trim(table_file)
194 call da_warning(__FILE__,__LINE__,message(1:1))
195 end if
196 call da_get_unit(airs_table_unit)
197 open(unit=airs_table_unit,file=table_file,iostat = iost)
198 if (iost /= 0) then
199 call da_error(__FILE__,__LINE__, &
200 (/"Cannot open file "//trim(table_file)/))
201 end if
202 end if
203
204 ! Open BUFR file
205 call da_get_unit(lnbufr)
206 open(unit=lnbufr,file=trim(infile),form='unformatted',iostat = iost)
207 if (iost /= 0) then
208 call da_warning(__FILE__,__LINE__, &
209 (/"Cannot open file "//trim(infile)/))
210 return
211 end if
212 if ( airstab ) then
213 call openbf(lnbufr,'IN',airs_table_unit)
214 else
215 call openbf(lnbufr,'IN',lnbufr)
216 end if
217 call datelen(10)
218
219
220 ! 2.0 Read header
221 !---------------------------
222 call readmg(lnbufr,subset,idate,iret)
223
224 iy = 0
225 im = 0
226 idd = 0
227 ihh = 0
228 if( iret /= 0 ) goto 1000 ! no data?
229
230 write(unit=date,fmt='( i10)') idate
231 read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
232 if (print_detail_rad) then
233 write(unit=message(1),fmt=*) 'Bufr file date is ', &
234 iy,im,idd,ihh
235 call da_message(message(1:1))
236 end if
237
238 ! 3.0 Loop over observations
239 !----------------------------
240
241 allocate ( head )
242 ! allocate ( head % tb (1:nchan) )
243 nullify ( head % next )
244 p => head
245 size = 0
246
247 loop_obspoints: do
248
249 ! 3.1 Read headder
250 !-------------------------------
251 call readsb(lnbufr,iret)
252
253 if( iret /=0 )then
254 call readmg(lnbufr,subset,idate,iret)
255 if( iret /= 0 ) exit loop_obspoints ! end of file
256 cycle loop_obspoints
257 end if
258
259 ! 3.2 Read AQUASPOT
260 !------------------------
261 call ufbseq(lnbufr,aquaspot,N_AQUASPOT_LIST,1,iret,'AQUASPOT')
262
263 ! 3.3 Read AIRSSPOT or AMSUSPOT or HSBSPOT
264 !-------------------------------------------------
265 call ufbseq(lnbufr,airsspot,N_AIRSSPOT_LIST,1,iret,spotname)
266
267 ! 3.4 Read AIRSCHAN or AMSUCHAN or HSBCHAN
268 !-------------------------------------------
269 call ufbseq(lnbufr,airschan,N_AIRSCHAN_LIST,N_MAXCHAN,iret,channame)
270
271 if (iret /= nchanr) then
272 write(unit=message(1),fmt=*) &
273 'Cannot read ', channame, &
274 ' bufr data:', &
275 iret, ' ch data is read instead of', nchanr
276 call da_warning(__FILE__,__LINE__,message(1:1))
277 cycle loop_obspoints
278 end if
279
280 ! 4.0 Check observing position (lat/lon)
281 ! QC1: juge if data is in the domain,
282 ! read next record if not
283 !------------------------------------------
284 if( abs(airsspot%clath) > R90 .or. &
285 abs(airsspot%clonh) > R360 .or. &
286 (abs(airsspot%clath) == R90 .and. airsspot%clonh /= ZERO) )then
287 cycle loop_obspoints
288 end if
289
290 ! Retrieve observing position
291 if(airsspot%clonh >= R360) then
292 airsspot%clonh = airsspot%clonh - R360
293 ! else if(airsspot%clonh < ZERO) then
294 ! airsspot%clonh = airsspot%clonh + R360
295 end if
296
297 info%lat = airsspot%clath
298 info%lon = airsspot%clonh
299 call da_llxy (info, loc, outside, inside_halo )
300 if (outside) cycle loop_obspoints
301
302 ! 4.1 Check obs time
303 !-------------------------------------
304 idate5(1) = airsspot%year ! year
305 idate5(2) = airsspot%mnth ! month
306 idate5(3) = airsspot%days ! day
307 idate5(4) = airsspot%hour ! hour
308 idate5(5) = airsspot%minu ! minute
309 idate5(6) = airsspot%seco ! second
310
311 if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. &
312 idate5(2) < 1 .or. idate5(2) > 12 .or. &
313 idate5(3) < 1 .or. idate5(3) > 31 .or. &
314 idate5(4) < 0 .or. idate5(4) > 24 .or. &
315 idate5(5) < 0 .or. idate5(5) > 60 .or. &
316 idate5(6) < 0 .or. idate5(6) > 60 ) then
317 cycle loop_obspoints
318 end if
319
320 ! QC3: time consistency check with the background date
321
322 if (idate5(1) .LE. 99) then
323 if (idate5(1) .LT. 78) then
324 idate5(1) = idate5(1) + 2000
325 else
326 idate5(1) = idate5(1) + 1900
327 end if
328 end if
329
330 call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
331
332 if ( obs_time < time_slots(0) .or. &
333 obs_time >= time_slots(num_fgat_time) ) cycle
334
335 ! 3.2.1 determine FGAT index ifgat
336 !
337 do ifgat=1,num_fgat_time
338 if ( obs_time >= time_slots(ifgat-1) .and. &
339 obs_time < time_slots(ifgat) ) exit
340 end do
341
342 write(unit=info%date_char, &
343 fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
344 idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
345 ':', idate5(5), ':', idate5(6)
346
347 info%elv = 0.0 !aquaspot%selv
348
349 ! 4.2 Check observational info
350 !-------------------------------------------------------
351 if( airsspot%fovn < 0.0_r_kind .or. airsspot%fovn > 100.0_r_kind .or. &
352 airsspot%saza < -360.0_r_kind .or. airsspot%saza > 360.0_r_kind .or. &
353 aquaspot%soza < -180.0_r_kind .or. aquaspot%soza > 180.0_r_kind )then
354 write(unit=message(1),fmt=*) &
355 'Cannot read ', channame, ' bufr data:', &
356 ' strange obs info(fov,saza,soza):', &
357 airsspot%fovn, airsspot%saza, aquaspot%soza
358 call da_warning(__FILE__,__LINE__,message(1:1))
359 cycle loop_obspoints
360 end if
361
362 ! Retrieve observing info
363 ifov = int( airsspot%fovn + POinT001 )
364 sza = abs(airsspot%saza)
365 ! if( ((airs .or. hsb) .and. ifov <= 45) .or. &
366 ! ( eos_amsua .and. ifov <= 15) )then
367 ! sza = - sza
368 ! end if
369
370 ! airdata(6) = (start + float(ifov-1)*step) ! look angle (deg)
371 ! airdata(9) = ZERO ! surface height
372 ! airdata(10)= POinT001 ! land sea mask
373
374
375 ! 4.3 Retrieve Tb
376 !-----------------------
377 iflag = .false.
378
379 do l=1,nchan
380 airdata(l+nreal) = airschan(l)%tmbrst ! brightness temerature
381 if( airdata(l+nreal) > 0.0_r_kind .and. airdata(l+nreal) < 500.0_r_kind )then
382 iflag = .true.
383 else
384 airdata(l+nreal) = missing_r
385 end if
386 end do
387
388 if ( .not. iflag )then
389 write(unit=message(1),fmt=*) &
390 'Error in reading ', channame, ' bufr data:', &
391 ' all tb data is missing'
392 call da_warning(__FILE__,__LINE__,message(1:1))
393 cycle loop_obspoints
394 end if
395
396 nobs=nobs + 1
397
398 ! 4.0 assign information to sequential radiance structure
399 !--------------------------------------------------------------------------
400 allocate ( p % tb_inv (1:nchan) )
401 p%info = info
402 p%loc = loc
403 p%landsea_mask = POinT001
404 p%scanline = int(aquaspot%slnm + POinT001)
405 p%scanpos = ifov
406 p%satzen = sza
407 p%satazi = (start + float(ifov-1)*step) ! look angle (deg) ! airsspot%bearaz
408 p%solzen = aquaspot%soza
409 p%solazi = aquaspot%solazi
410 p%tb_inv(1:nchan) = airdata(nreal+1:nreal+nchan)
411 p%sensor_index = inst
412 p%ifgat = ifgat
413
414 size = size + 1
415 allocate ( p%next, stat=error) ! add next data
416 if (error /= 0 ) then
417 call da_error(__FILE__,__LINE__, &
418 (/"Cannot allocate radiance structure"/))
419 end if
420
421 p => p%next
422 nullify (p%next)
423
424 end do loop_obspoints
425
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 = kme-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 call da_allocate_rad_iv (inst, nchan, iv)
458
459 do i = 1, size
460 ! inst = p%sensor_index
461 n = n + 1
462
463 call da_initialize_rad_iv (inst, n, iv, p)
464
465 current => p
466 p => p%next
467
468 ! free current data
469 deallocate ( current % tb_inv )
470 deallocate ( current )
471
472 end do
473
474 1000 continue
475 call closbf(lnbufr)
476 close(lnbufr)
477 call da_free_unit(lnbufr)
478 if (airstab) then
479 close(airs_table_unit)
480 call da_free_unit(airs_table_unit)
481 end if
482
483 if (trace_use_dull) call da_trace_exit("da_read_obs_bufrairs")
484
485 #endif
486
487 end subroutine da_read_obs_bufrairs