da_read_obs_bufrtovs.inc
References to this file elsewhere.
1 subroutine da_read_obs_bufrtovs (obstype,iv, infile)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: read in NCEP bufr tovs 1b data to innovation structure
5 !
6 ! METHOD: use F90 sequential data structure to avoid reading file twice
7 ! so that da_scan_bufrtovs 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
14 implicit none
15
16 character(5) , intent (in) :: obstype
17 character(20) , intent (in) :: infile
18 type (iv_type) , intent (inout) :: iv
19
20 #ifdef BUFR
21
22 integer :: iost
23 integer(i_kind), allocatable :: nread(:)
24
25 integer(i_kind),parameter:: n1bhdr=14
26 integer(i_kind),parameter:: maxinfo=12
27 integer(i_kind),parameter:: maxchanl=100
28
29 logical hirs2,hirs3,hirs4,msu,amsua,amsub,mhs
30 logical outside, outside_all, iuse
31 integer :: inst
32
33 character(10) date
34 character(8) subset,subfgn
35 character(80) hdr1b
36 integer(i_kind) ihh,i,k,ifov,idd,ireadmg,ireadsb
37 integer(i_kind) iret,idate,im,iy,nchan
38
39 ! thinning variables
40 integer(i_kind) itt,itx,iobs,iout
41 real(r_kind) terrain,timedif,crit,dist
42 real(r_kind) dlon_earth,dlat_earth
43
44 real(r_kind) tbmin,tbmax, tbbad
45 real(r_kind) panglr,rato
46 ! real(r_kind) rmask
47 real(r_kind) step,start
48
49 real(r_double),dimension(maxinfo+maxchanl):: data1b8
50 real(r_double),dimension(n1bhdr):: bfr1bhdr
51
52 ! Instrument triplet, follow the convension of RTTOV
53 integer :: platform_id, satellite_id, sensor_id
54
55 ! pixel information
56 integer :: year,month,day,hour,minute,second ! observation time
57 real :: obs_time
58 ! real :: rlat, rlon ! lat/lon in degrees for Anfovs
59 real :: satzen, satazi, solzen ! scan angles for Anfovs
60 integer :: landsea_mask
61 real :: srf_height
62 ! channels' bright temperature
63 real , allocatable :: tb_inv(:) ! bright temperatures
64 ! end type bright_temperature
65
66 type (datalink_type), pointer :: head, p, current
67
68 integer :: ifgat
69 type(info_type) :: info
70 type(model_loc_type) :: loc
71
72 data hdr1b /'SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SAZA SOZA HOLS LSQL'/
73 ! data hdr1b /'FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SAZA SOZA HOLS LSQL SLNM BEARAZ'/
74
75 data tbmin,tbmax,tbbad / 50.0_r_kind, 550.0_r_kind, -9.99e11_r_kind /
76 integer :: num_tovs_local, num_tovs_file, num_tovs_global, num_tovs_selected
77 integer :: num_tovs_thinned, num_tovs_used
78 integer :: lnbufr
79 integer :: n
80 integer(i_kind), allocatable :: ptotal(:)
81
82 call da_trace_entry("da_read_obs_bufrtovs")
83
84 ! Initialize variables
85
86 nchan = 20
87 allocate(nread(1:rtminit_nsensor))
88 allocate(ptotal(0:num_fgat_time))
89 nread(1:rtminit_nsensor) = 0
90 ptotal(0:num_fgat_time) = 0
91
92 ! Set various variables depending on type of data to be read
93
94 platform_id = 1 !! for NOAA series
95
96 hirs2= obstype == 'hirs2'
97 hirs3= obstype == 'hirs3'
98 hirs4= obstype == 'hirs4'
99 msu= obstype == 'msu '
100 amsua= obstype == 'amsua'
101 amsub= obstype == 'amsub'
102 mhs= obstype == 'mhs '
103
104 if (hirs2) then
105 sensor_id = 0
106 step = 1.80_r_kind
107 start = -49.5_r_kind
108 nchan=nchan_hirs2
109 subfgn='NC021021'
110 rato=1.1363987_r_kind
111 else if (hirs3) then
112 sensor_id = 0
113 step = 1.80_r_kind
114 start = -49.5_r_kind
115 nchan=nchan_hirs3
116 subfgn='NC021025'
117 else if (hirs4) then
118 sensor_id = 0
119 step = 1.80_r_kind
120 start = -49.5_r_kind
121 nchan=nchan_hirs4
122 subfgn='NC021028'
123 else if (mhs) then
124 sensor_id = 15
125 step = 10.0_r_kind/9.0_r_kind
126 start = -445.0_r_kind/9.0_r_kind
127 nchan=nchan_mhs
128 subfgn='NC021027'
129 else if (msu) then
130 sensor_id = 1
131 step = 9.474_r_kind
132 start = -47.37_r_kind
133 nchan=nchan_msu
134 subfgn='NC021022'
135 rato=1.1363987_r_kind
136 else if (amsua) then
137 sensor_id = 3
138 step = three + one/three
139 start = -48.33_r_kind
140 nchan=nchan_amsua
141 subfgn='NC021023'
142 else if (amsub) then
143 sensor_id = 4
144 step = 1.1_r_kind
145 start = -48.95_r_kind
146 nchan=nchan_amsub
147 subfgn='NC021024'
148 end if
149
150 allocate (tb_inv(nchan))
151
152 ! 0.0 Open unit to satellite bufr file and read file header
153 !--------------------------------------------------------------
154
155 call da_get_unit(lnbufr)
156 open(unit=lnbufr,file=trim(infile),form='unformatted', &
157 iostat = iost, status = 'old')
158 if (iost /= 0) then
159 call da_warning(__FILE__,__LINE__, &
160 (/"Cannot open file "//trim(infile)/))
161 call da_free_unit(lnbufr)
162 call da_trace_exit("da_read_obs_bufrtovs")
163 return
164 end if
165
166 call openbf(lnbufr,'IN',lnbufr)
167 call datelen(10)
168 call readmg(lnbufr,subset,idate,iret)
169 if (subset /= subfgn) then
170 message(1)='The file title does not match the data subset'
171 write(unit=message(2),fmt=*) &
172 'infile=', lnbufr, infile,' subset=', subset, ' subfgn=',subfgn
173 call da_error(__FILE__,__LINE__,message(1:2))
174 end if
175
176 iy=0
177 im=0
178 idd=0
179 ihh=0
180 write(unit=date,fmt='( i10)') idate
181 read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
182 write(unit=stdout,fmt=*) &
183 'Bufr file date is ',iy,im,idd,ihh,infile
184
185 ! Loop to read bufr file and assign information to a sequential structure
186 !-------------------------------------------------------------------------
187
188 allocate (head)
189 ! allocate ( head % tb_inv (1:nchan) )
190 nullify ( head % next )
191 p => head
192
193 num_tovs_file = 0 ! number of obs in file
194 num_tovs_global = 0 ! number of obs within whole domain
195 num_tovs_local = 0 ! number of obs within tile
196 num_tovs_thinned = 0 ! number of obs rejected by thinning
197 num_tovs_used = 0 ! number of obs entered into innovation computation
198 num_tovs_selected = 0 ! number of obs limited for debuging
199
200 if (tovs_start > 1) then
201 write (unit=stdout,fmt='(A,I6)') " Skipping tovs obs before", tovs_start
202 end if
203
204 obs: do while (ireadmg(lnbufr,subset,idate)==0 .and. subset==subfgn)
205 do while (ireadsb(lnbufr)==0)
206
207 ! 1.0 Read header record and data record
208
209 call ufbint(lnbufr,bfr1bhdr,n1bhdr,1,iret,hdr1b)
210 call ufbrep(lnbufr,data1b8,1,nchan,iret,'TMBR')
211 ! call ufbrep(lnbufr,data1b8,1,1,iret,'BEARAZ')
212
213 ! check if observation outside range
214
215 num_tovs_file = num_tovs_file + 1
216
217 ! 2.0 Extract observation location and other required information
218 ! QC1: judge if data is in the domain, read next record if not
219 !------------------------------------------------------------------------
220 ! rlat = bfr1bhdr(bufr_lat)
221 ! rlon = bfr1bhdr(bufr_lat)
222 ! if (rlon < 0.0) rlon = rlon+360.0
223
224 info%lat = bfr1bhdr(bufr_lat)
225 info%lon = bfr1bhdr(bufr_lon)
226 call da_llxy (info, loc, outside, outside_all)
227
228 if (outside_all) cycle
229
230 ! 3.0 Extract other information
231 !------------------------------------------------------
232 ! 3.1 Extract satellite id and scan position.
233
234 satellite_id = nint(bfr1bhdr(bufr_satellite_id))-191
235 ifov = nint(bfr1bhdr(bufr_ifov))
236
237 ! QC2: limb pixel rejected (not implemented)
238
239 ! 3.2 Extract date information.
240
241 year = bfr1bhdr(bufr_year)
242 month = bfr1bhdr(bufr_month)
243 day = bfr1bhdr(bufr_day)
244 hour = bfr1bhdr(bufr_hour)
245 minute = bfr1bhdr(bufr_minute)
246 second = bfr1bhdr(bufr_second)
247
248 write(unit=info%date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
249 year, '-', month, '-', day, '_', hour, ':', minute, ':', second
250
251 ! QC3: time consistency check with the background date
252
253 if (year <= 99) then
254 if (year < 78) then
255 year = year + 2000
256 else
257 year = year + 1900
258 end if
259 end if
260
261 call da_get_julian_time(year,month,day,hour,minute,obs_time)
262
263 if (obs_time < time_slots(0) .or. &
264 obs_time >= time_slots(num_fgat_time)) cycle
265
266 ! 3.2.1 determine FGAT index ifgat
267
268 do ifgat=1,num_fgat_time
269 if (obs_time >= time_slots(ifgat-1) .and. &
270 obs_time < time_slots(ifgat)) exit
271 end do
272
273 ! 3.3 Find wrfvar instrument index from RTTOV instrument triplet
274 ! go to next data if id is not in the lists
275
276 inst = 0
277 do i = 1, rtminit_nsensor
278 if (platform_id == rtminit_platform(i) &
279 .and. satellite_id == rtminit_satid(i) &
280 .and. sensor_id == rtminit_sensor(i)) then
281 inst = i
282 exit
283 end if
284 end do
285 if (inst == 0) cycle
286
287 num_tovs_global = num_tovs_global + 1
288 ptotal(ifgat) = ptotal(ifgat) + 1
289
290 if (num_tovs_global < tovs_start) then
291 cycle
292 end if
293
294 if (num_tovs_global > tovs_end) then
295 write (unit=stdout,fmt='(A,I6)') " Skipping radiance obs after", tovs_end
296 exit obs
297 end if
298
299 num_tovs_selected = num_tovs_selected + 1
300
301 if (num_tovs_selected > max_tovs_input) then
302 write(unit=message(1),fmt='(A,I10,A)') &
303 "Max number of tovs",max_tovs_input," reached"
304 call da_warning(__FILE__,__LINE__,message(1:1))
305 num_tovs_selected = num_tovs_selected-1
306 num_tovs_global = num_tovs_global-1
307 ptotal(ifgat) = ptotal(ifgat) - 1
308 exit obs
309 end if
310
311 if (outside) cycle ! No good for this PE
312 num_tovs_local = num_tovs_local + 1
313
314 ! Make Thinning
315 ! Map obs to thinning grid
316 !-------------------------------------------------------------------
317 if (thinning) then
318 dlat_earth = info%lat
319 dlon_earth = info%lon
320 if (dlon_earth<zero) dlon_earth = dlon_earth+r360
321 if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
322 dlat_earth = dlat_earth*deg2rad
323 dlon_earth = dlon_earth*deg2rad
324 timedif = 0.0 !2.0_r_kind*abs(tdiff) ! range: 0 to 6
325 terrain = 0.01_r_kind*abs(bfr1bhdr(13))
326 crit = 1.0 !0.01_r_kind+terrain + timedif !+ 10.0_r_kind*float(iskip)
327 iobs=0 ! argument is inout
328 call map2grids(inst,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
329 if (.not. iuse) then
330 num_tovs_thinned=num_tovs_thinned+1
331 cycle
332 end if
333 end if
334
335 num_tovs_used = num_tovs_used + 1
336 nread(inst) = nread(inst) + 1
337
338 ! 3.4 extract satellite and solar angle
339
340 panglr=(start+float(ifov-1)*step)*deg2rad
341 if (hirs2 .or. msu) then
342 satzen = asin(rato*sin(panglr))*rad2deg
343 satzen = abs(satzen)
344 else
345 satzen = bfr1bhdr(bufr_satzen) !*deg2rad ! local zenith angle
346 satzen = abs(satzen)
347 ! if (amsua .and. ifov .le. 15) satzen=-satzen
348 ! if (amsub .and. ifov .le. 45) satzen=-satzen
349 ! if (hirs3 .and. ifov .le. 28) satzen=-satzen
350 end if
351 satazi = panglr*rad2deg ! look angle
352 ! if (satazi<0.0) satazi = satazi+360.0
353 solzen = bfr1bhdr(bufr_solzen) ! solar zenith angle
354
355 ! 3.5 extract surface information
356
357 srf_height = bfr1bhdr(bufr_station_height) ! station height
358 landsea_mask = nint(bfr1bhdr(bufr_landsea_mask)) ! 0:land ; 1:sea (same as RTTOV)
359 ! rmask=one ! land
360 ! if (nint(bfr1bhdr(bufr_landsea_mask))==1) rmask=0.0_r_kind ! reverse the land/sea mask in bufr
361 ! landsea_mask = rmask+.001_r_kind ! land sea mask
362
363 info%elv = srf_height
364
365 ! 3.6 extract channel bright temperature
366
367 tb_inv(1:nchan) = data1b8(1:nchan)
368 do k = 1, nchan
369 if ( tb_inv(k) < tbmin .or. tb_inv(k) > tbmax) &
370 tb_inv(k) = missing_r
371 end do
372
373 ! 4.0 assign information to sequential radiance structure
374 !--------------------------------------------------------------------------
375 allocate (p % tb_inv (1:nchan))
376 p%info = info
377 p%loc = loc
378 p%landsea_mask = landsea_mask
379 p%scanpos = ifov
380 p%satzen = satzen
381 p%satazi = satazi
382 p%solzen = solzen
383 p%tb_inv(1:nchan) = tb_inv(1:nchan)
384 p%sensor_index = inst
385 p%ifgat = ifgat
386
387 allocate (p%next) ! add next data
388
389 p => p%next
390 nullify (p%next)
391 end do
392 end do obs
393
394 iv%total_rad_pixel = iv%total_rad_pixel + num_tovs_used
395 iv%total_rad_channel = iv%total_rad_channel + num_tovs_used*nchan
396
397 iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_tovs_used
398 iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_tovs_global
399
400 do i = 1, num_fgat_time
401 ptotal(i) = ptotal(i) + ptotal(i-1)
402 iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
403 end do
404 if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
405 write(unit=message(1),fmt='(A,I10,A,I10)') &
406 "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
407 call da_warning(__FILE__,__LINE__,message(1:1))
408 endif
409
410 write(unit=stdout,fmt='(a)') 'num_tovs_file num_tovs_global num_tovs_local num_tovs_used num_tovs_thinned'
411 write(unit=stdout,fmt='(5i10)') num_tovs_file,num_tovs_global, num_tovs_local,num_tovs_used,num_tovs_thinned
412
413 deallocate(tb_inv)
414 call closbf(lnbufr)
415 call da_free_unit(lnbufr)
416
417 ! 5.0 allocate innovation radiance structure
418 !----------------------------------------------------------------
419
420 do i = 1, iv%num_inst
421 if (nread(i) < 1) cycle
422 iv%instid(i)%num_rad = nread(i)
423 iv%instid(i)%info%nlocal = nread(i)
424 write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
425 'Allocating space for radiance innov structure', &
426 i, iv%instid(i)%rttovid_string, iv%instid(i)%num_rad
427 if (rtm_option == rtm_option_crtm) then
428 iv%instid(i)%nlevels = kme-kms+1
429 iv%instid(i)%nchannels=nchan
430 end if
431
432 #ifdef RTTOV
433 if (rtm_option == rtm_option_rttov) then
434 call rttov_setupchan(1, nchan, coefs(i), & ! in
435 iv%instid(i)%nfrequencies,iv%instid(i)%nchannels, &
436 iv%instid(i)%nbtout) ! out
437 end if
438 #endif
439
440 call da_allocate_rad_iv(i,nchan,iv)
441
442 end do
443
444 ! 6.0 assign sequential structure to innovation structure
445 !-------------------------------------------------------------
446 nread(1:rtminit_nsensor) = 0
447 p => head
448 ! do while ( associated(p) )
449
450 do n = 1, num_tovs_used
451 i = p%sensor_index
452 nread(i) = nread(i) + 1
453
454 call da_initialize_rad_iv (i, nread(i), iv, p)
455
456 current => p
457 p => p%next
458
459 ! free current data
460 deallocate ( current % tb_inv )
461 deallocate ( current )
462 end do
463
464 deallocate ( p )
465
466 deallocate (nread)
467 deallocate (ptotal)
468
469 ! check if sequential structure has been freed
470 !
471 ! p => head
472 ! do i = 1, num_rad_selected
473 ! write (unit=stdout,fmt=*) i, p%tb_inv(1:nchan)
474 ! p => p%next
475 ! end do
476
477 call da_trace_exit("da_read_obs_bufrtovs")
478 #else
479 call da_error(__FILE__,__LINE__,(/"Needs to be compiled with a BUFR library"/))
480 #endif
481
482 end subroutine da_read_obs_bufrtovs
483
484