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