da_read_bufrtovs.inc
References to this file elsewhere.
1 subroutine da_read_bufrtovs (obstype,iv,xp,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 use da_control
15
16 implicit none
17
18 character(5) , intent (in) :: obstype
19 character(20) , intent (in) :: infile
20 type (xpose_type) , intent (in) :: xp
21 type (ob_type) , intent (inout) :: iv
22
23 #ifdef BUFR
24
25 integer :: iost
26 integer(i_kind), allocatable :: nread(:)
27
28 integer(i_kind),parameter:: n1bhdr=14
29 integer(i_kind),parameter:: maxinfo=12
30 integer(i_kind),parameter:: maxchanl=100
31
32 logical hirs2,hirs3,msu,amsua,amsub
33 logical outside, outside_all
34 integer :: inst
35
36 character(10) date
37 character(8) subset,subfgn
38 character(80) hdr1b
39
40 integer(i_kind) ihh,i,k,ifov,idd,ireadmg,ireadsb
41 integer(i_kind) iret,idate,im,iy,nchan
42
43 real(r_kind) tbmin,tbmax, tbbad
44 real(r_kind) panglr,rato
45 ! real(r_kind) rmask
46 real(r_kind) step,start
47
48 real(r_double),dimension(maxinfo+maxchanl):: data1b8
49 real(r_double),dimension(n1bhdr):: bfr1bhdr
50
51 ! type bright_temperature
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_used
77 integer :: lnbufr
78 integer :: n
79
80 call da_trace_entry("da_read_bufrtovs")
81
82 ! Initialize variables
83
84 nchan = 20
85 allocate(nread(1:rtminit_nsensor))
86 nread(1:rtminit_nsensor) = 0
87
88 ! Set various variables depending on type of data to be read
89
90 call init_constants_derived
91
92 platform_id = 1 !! for NOAA series
93
94 hirs2= obstype == 'hirs2'
95 hirs3= obstype == 'hirs3'
96 msu= obstype == 'msu '
97 amsua= obstype == 'amsua'
98 amsub= obstype == 'amsub'
99
100 if (hirs2) then
101 sensor_id = 0
102 step = 1.80_r_kind
103 start = -49.5_r_kind
104 nchan=nchan_hirs2
105 subfgn='NC021021'
106 rato=1.1363987_r_kind
107 else if (hirs3) then
108 sensor_id = 0
109 step = 1.80_r_kind
110 start = -49.5_r_kind
111 nchan=nchan_hirs3
112 subfgn='NC021025'
113 else if (msu) then
114 sensor_id = 1
115 step = 9.474_r_kind
116 start = -47.37_r_kind
117 nchan=nchan_amsua
118 subfgn='NC021022'
119 rato=1.1363987_r_kind
120 else if (amsua) then
121 sensor_id = 3
122 step = three + one/three
123 start = -48.33_r_kind
124 nchan=nchan_amsua
125 subfgn='NC021023'
126 else if (amsub) then
127 sensor_id = 4
128 step = 1.1_r_kind
129 start = -48.95_r_kind
130 nchan=5
131 subfgn='NC021024'
132 end if
133
134 allocate (tb_inv(nchan))
135
136 ! 0.0 Open unit to satellite bufr file and read file header
137 !--------------------------------------------------------------
138
139 call da_get_unit(lnbufr)
140 open(unit=lnbufr,file=trim(infile),form='unformatted', &
141 iostat = iost, status = 'old')
142 if (iost /= 0) then
143 call da_warning(__FILE__,__LINE__, &
144 (/"Cannot open file "//trim(infile)/))
145 return
146 end if
147
148 call openbf(lnbufr,'IN',lnbufr)
149 call datelen(10)
150 call readmg(lnbufr,subset,idate,iret)
151 if (subset /= subfgn) then
152 message(1)='The file title does not match the data subset'
153 write(unit=message(2),fmt=*) &
154 'infile=', lnbufr, infile,' subset=', subset, ' subfgn=',subfgn
155 call da_error(__FILE__,__LINE__,message(1:2))
156 end if
157
158 iy=0
159 im=0
160 idd=0
161 ihh=0
162 write(unit=date,fmt='( i10)') idate
163 read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
164 write(unit=stdout,fmt=*) &
165 'Bufr file date is ',iy,im,idd,ihh,infile
166
167 ! Loop to read bufr file and assign information to a sequential structure
168 !-------------------------------------------------------------------------
169
170 allocate (head)
171 ! allocate ( head % tb_inv (1:nchan) )
172 nullify ( head % next )
173 p => head
174
175 num_tovs_local = 0
176 num_tovs_file = 0
177 num_tovs_global = 0
178 num_tovs_used = 0
179
180 if (tovs_start > 1) then
181 write (unit=stdout,fmt='(A,I6)') " Skipping tovs obs before", tovs_start
182 end if
183
184 obs: do while (ireadmg(lnbufr,subset,idate)==0 .and. subset==subfgn)
185 do while (ireadsb(lnbufr)==0)
186
187 ! 1.0 Read header record and data record
188
189 call ufbint(lnbufr,bfr1bhdr,n1bhdr,1,iret,hdr1b)
190 call ufbrep(lnbufr,data1b8,1,nchan,iret,'TMBR')
191 ! call ufbrep(lnbufr,data1b8,1,1,iret,'BEARAZ')
192
193 ! check if observation outside range
194
195 num_tovs_file = num_tovs_file + 1
196
197 ! 2.0 Extract observation location and other required information
198 ! QC1: judge if data is in the domain, read next record if not
199 !------------------------------------------------------------------------
200 ! rlat = bfr1bhdr(bufr_lat)
201 ! rlon = bfr1bhdr(bufr_lat)
202 ! if (rlon < 0.0) rlon = rlon+360.0
203
204 info%lat = bfr1bhdr(bufr_lat)
205 info%lon = bfr1bhdr(bufr_lon)
206 call da_ll_to_xy (info, loc, xp, outside, outside_all)
207
208 if (outside_all) cycle
209
210 ! 3.0 Extract other information
211 !------------------------------------------------------
212 ! 3.1 Extract satellite id and scan position.
213
214 satellite_id = nint(bfr1bhdr(bufr_satellite_id))-191
215 ifov = nint(bfr1bhdr(bufr_ifov))
216
217 ! QC2: limb pixel rejected (not implemented)
218
219 ! 3.2 Extract date information.
220
221 year = bfr1bhdr(bufr_year)
222 month = bfr1bhdr(bufr_month)
223 day = bfr1bhdr(bufr_day)
224 hour = bfr1bhdr(bufr_hour)
225 minute = bfr1bhdr(bufr_minute)
226 second = bfr1bhdr(bufr_second)
227
228 write(unit=info%date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
229 year, '-', month, '-', day, '_', hour, ':', minute, ':', second
230
231 ! QC3: time consistency check with the background date
232
233 if (year <= 99) then
234 if (year < 78) then
235 year = year + 2000
236 else
237 year = year + 1900
238 end if
239 end if
240
241 call da_get_julian_time(year,month,day,hour,minute,obs_time)
242
243 if (obs_time < time_slots(0) .or. &
244 obs_time >= time_slots(num_fgat_time)) cycle
245
246 ! 3.2.1 determine FGAT index ifgat
247
248 do ifgat=1,num_fgat_time
249 if (obs_time >= time_slots(ifgat-1) .and. &
250 obs_time < time_slots(ifgat)) exit
251 end do
252
253 ! 3.3 Find wrfvar instrument index from RTTOV instrument triplet
254 ! go to next data if id is not in the lists
255
256 inst = 0
257 do i = 1, rtminit_nsensor
258 if (platform_id == rtminit_platform(i) &
259 .and. satellite_id == rtminit_satid(i) &
260 .and. sensor_id == rtminit_sensor(i)) then
261 inst = i
262 exit
263 end if
264 end do
265 if (inst == 0) cycle
266
267 num_tovs_global = num_tovs_global + 1
268
269 if (num_tovs_global < tovs_start) then
270 cycle
271 end if
272
273 if (num_tovs_global > tovs_end) then
274 write (unit=stdout,fmt='(A,I6)') " Skipping radiance obs after", tovs_end
275 exit obs
276 end if
277
278 num_tovs_used = num_tovs_used + 1
279
280 if (num_tovs_used > max_tovs_input) then
281 write(unit=message(1),fmt='(A,I10,A)') &
282 "Max number of tovs",max_tovs_input," reached"
283 call da_warning(__FILE__,__LINE__,message(1:1))
284 num_tovs_used = num_tovs_used-1
285 exit obs
286 end if
287 if (outside) cycle ! No good for this PE
288
289 num_tovs_local = num_tovs_local + 1
290 nread(inst) = nread(inst) + 1
291
292 ! 3.4 extract satellite and solar angle
293
294 panglr=(start+float(ifov-1)*step)*deg2rad
295 if (hirs2 .or. msu) then
296 satzen = asin(rato*sin(panglr))*rad2deg
297 satzen = abs(satzen)
298 else
299 satzen = bfr1bhdr(bufr_satzen) !*deg2rad ! local zenith angle
300 satzen = abs(satzen)
301 ! if (amsua .and. ifov .le. 15) satzen=-satzen
302 ! if (amsub .and. ifov .le. 45) satzen=-satzen
303 ! if (hirs3 .and. ifov .le. 28) satzen=-satzen
304 end if
305 satazi = panglr*rad2deg ! look angle
306 ! if (satazi<0.0) satazi = satazi+360.0
307 solzen = bfr1bhdr(bufr_solzen) ! solar zenith angle
308
309 ! 3.5 extract surface information
310
311 srf_height = bfr1bhdr(bufr_station_height) ! station height
312 landsea_mask = nint(bfr1bhdr(bufr_landsea_mask)) ! 0:land ; 1:sea (same as RTTOV)
313 ! rmask=one ! land
314 ! if (nint(bfr1bhdr(bufr_landsea_mask))==1) rmask=0.0_r_kind ! reverse the land/sea mask in bufr
315 ! landsea_mask = rmask+.001_r_kind ! land sea mask
316
317 info%elv = srf_height
318
319 ! 3.6 extract channel bright temperature
320
321 tb_inv(1:nchan) = data1b8(1:nchan)
322 do k = 1, nchan
323 if ( tb_inv(k) < tbmin .or. tb_inv(k) > tbmax) &
324 tb_inv(k) = missing_r
325 end do
326
327 ! 4.0 assign information to sequential radiance structure
328 !--------------------------------------------------------------------------
329 allocate (p % tb_inv (1:nchan))
330 p%info = info
331 p%loc = loc
332 p%landsea_mask = landsea_mask
333 p%scanpos = ifov
334 p%satzen = satzen
335 p%satazi = satazi
336 p%solzen = solzen
337 p%tb_inv(1:nchan) = tb_inv(1:nchan)
338 p%sensor_index = inst
339 p%ifgat = ifgat
340
341 allocate (p%next) ! add next data
342
343 p => p%next
344 nullify (p%next)
345 end do
346 end do obs
347
348 iv%total_obs = iv%total_obs + num_tovs_used
349 iv%total_rad_pixel = iv%total_rad_pixel + num_tovs_used
350 iv%total_rad_channel = iv%total_rad_channel + num_tovs_used*nchan
351
352 deallocate(tb_inv)
353 call closbf(lnbufr)
354 call da_free_unit(lnbufr)
355
356 ! 5.0 allocate innovation radiance structure
357 !----------------------------------------------------------------
358
359 do i = 1, iv%num_inst
360 if (nread(i) < 1) cycle
361 iv%instid(i)%num_rad = nread(i)
362 write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
363 'Allocating space for radiance innov structure', &
364 i, iv%instid(i)%rttovid_string, iv%instid(i)%num_rad
365 if (rtm_option == rtm_option_crtm) then
366 iv%instid(i)%nlevels = xp%kme-xp%kms+1
367 iv%instid(i)%nchannels=nchan
368 end if
369
370 #ifdef RTTOV
371 if (rtm_option == rtm_option_rttov) then
372 call rttov_setupchan(1, nchan, coefs(i), & ! in
373 iv%instid(i)%nfrequencies,iv%instid(i)%nchannels, &
374 iv%instid(i)%nbtout) ! out
375 end if
376 #endif
377
378 allocate (iv%instid(i)%info(iv%instid(i)%num_rad))
379 allocate (iv%instid(i)%loc(iv%instid(i)%num_rad))
380 allocate (iv%instid(i)%loc_i(iv%instid(i)%num_rad))
381 allocate (iv%instid(i)%loc_j(iv%instid(i)%num_rad))
382 allocate (iv%instid(i)%loc_k(iv%instid(i)%nlevels,iv%instid(i)%num_rad))
383 allocate (iv%instid(i)%loc_dz(iv%instid(i)%nlevels,iv%instid(i)%num_rad))
384 allocate (iv%instid(i)%loc_dzm(iv%instid(i)%nlevels,iv%instid(i)%num_rad))
385 allocate (iv%instid(i)%zk(iv%instid(i)%nlevels,iv%instid(i)%num_rad))
386 allocate (iv%instid(i)%t(iv%instid(i)%nlevels,iv%instid(i)%num_rad))
387 allocate (iv%instid(i)%mr(iv%instid(i)%nlevels,iv%instid(i)%num_rad))
388 allocate (iv%instid(i)%loc_dx(iv%instid(i)%num_rad))
389 allocate (iv%instid(i)%loc_dy(iv%instid(i)%num_rad))
390 allocate (iv%instid(i)%loc_dxm(iv%instid(i)%num_rad))
391 allocate (iv%instid(i)%loc_dym(iv%instid(i)%num_rad))
392 allocate (iv%instid(i)%tm(xp%kms:xp%kme,iv%instid(i)%num_rad))
393 allocate (iv%instid(i)%qm(xp%kms:xp%kme,iv%instid(i)%num_rad))
394 allocate (iv%instid(i)%qrn(xp%kms:xp%kme,iv%instid(i)%num_rad))
395 allocate (iv%instid(i)%qcw(xp%kms:xp%kme,iv%instid(i)%num_rad))
396 allocate (iv%instid(i)%qci(xp%kms:xp%kme,iv%instid(i)%num_rad))
397 allocate (iv%instid(i)%qsn(xp%kms:xp%kme,iv%instid(i)%num_rad))
398 allocate (iv%instid(i)%qgr(xp%kms:xp%kme,iv%instid(i)%num_rad))
399 allocate (iv%instid(i)%pm(xp%kms:xp%kme,iv%instid(i)%num_rad))
400 allocate (iv%instid(i)%pf(0:xp%kme,iv%instid(i)%num_rad))
401 allocate (iv%instid(i)%u10(iv%instid(i)%num_rad))
402 allocate (iv%instid(i)%v10(iv%instid(i)%num_rad))
403 allocate (iv%instid(i)%t2m(iv%instid(i)%num_rad))
404 allocate (iv%instid(i)%q2m(iv%instid(i)%num_rad))
405 allocate (iv%instid(i)%mr2m(iv%instid(i)%num_rad))
406 allocate (iv%instid(i)%psfc(iv%instid(i)%num_rad))
407 allocate (iv%instid(i)%ts(iv%instid(i)%num_rad))
408 allocate (iv%instid(i)%smois(iv%instid(i)%num_rad))
409 allocate (iv%instid(i)%tslb(iv%instid(i)%num_rad))
410 allocate (iv%instid(i)%snowh(iv%instid(i)%num_rad))
411 allocate (iv%instid(i)%isflg(iv%instid(i)%num_rad))
412 allocate (iv%instid(i)%landsea_mask(iv%instid(i)%num_rad))
413 allocate (iv%instid(i)%elevation(iv%instid(i)%num_rad))
414 allocate (iv%instid(i)%vegfra(iv%instid(i)%num_rad))
415 allocate (iv%instid(i)%vegtyp(iv%instid(i)%num_rad))
416 allocate (iv%instid(i)%soiltyp(iv%instid(i)%num_rad))
417 allocate (iv%instid(i)%clwp(iv%instid(i)%num_rad))
418 allocate (iv%instid(i)%ps(iv%instid(i)%num_rad))
419 allocate (iv%instid(i)%tb_xb(nchan,iv%instid(i)%num_rad))
420 allocate (iv%instid(i)%tb_qc(nchan,iv%instid(i)%num_rad))
421 allocate (iv%instid(i)%tb_inv(nchan,iv%instid(i)%num_rad))
422 allocate (iv%instid(i)%tb_error(nchan,iv%instid(i)%num_rad))
423 allocate (iv%instid(i)%emiss(iv%instid(i)%nchannels,iv%instid(i)%num_rad))
424 allocate(iv%instid(i)%scanpos(iv%instid(i)%num_rad))
425 allocate(iv%instid(i)%scanline(iv%instid(i)%num_rad))
426 allocate(iv%instid(i)%ifgat(iv%instid(i)%num_rad))
427 allocate(iv%instid(i)%cloud_flag(nchan,iv%instid(i)%num_rad))
428 allocate(iv%instid(i)%satzen(iv%instid(i)%num_rad))
429 allocate(iv%instid(i)%satazi(iv%instid(i)%num_rad))
430 allocate(iv%instid(i)%solzen(iv%instid(i)%num_rad))
431 allocate(iv%instid(i)%solazi(iv%instid(i)%num_rad))
432 allocate(iv%instid(i)%proc_domain(iv%instid(i)%num_rad))
433 if (rtm_option == rtm_option_crtm) then
434 allocate(iv%instid(i)%water_coverage(iv%instid(i)%num_rad))
435 allocate(iv%instid(i)%land_coverage(iv%instid(i)%num_rad))
436 allocate(iv%instid(i)%ice_coverage(iv%instid(i)%num_rad))
437 allocate(iv%instid(i)%snow_coverage(iv%instid(i)%num_rad))
438 allocate(iv%instid(i)%ps_jacobian(nchan,iv%instid(i)%num_rad))
439 allocate(iv%instid(i)%t_jacobian(nchan,xp%kte,iv%instid(i)%num_rad))
440 allocate(iv%instid(i)%q_jacobian(nchan,xp%kte,iv%instid(i)%num_rad))
441 end if
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_local
451 i = p%sensor_index
452 nread(i) = nread(i) + 1
453 iv%instid(i)%info(nread(i)) = p%info
454 iv%instid(i)%loc(nread(i)) = p%loc
455 iv%instid(i)%loc_i(nread(i)) = p%loc%i
456 iv%instid(i)%loc_j(nread(i)) = p%loc%j
457 iv%instid(i)%loc_k(:,nread(i)) = 0
458 iv%instid(i)%loc_dx(nread(i)) = p%loc%dx
459 iv%instid(i)%loc_dy(nread(i)) = p%loc%dy
460 iv%instid(i)%loc_dz(:,nread(i)) = 0.0
461 iv%instid(i)%loc_dxm(nread(i)) = p%loc%dxm
462 iv%instid(i)%loc_dym(nread(i)) = p%loc%dym
463 iv%instid(i)%loc_dzm(:,nread(i)) = 0.0
464 ! z done in da_get_innov_vector_rad
465 iv%instid(i)%t(:,nread(i)) = 0.0
466 iv%instid(i)%mr(:,nread(i)) = 0.0
467 iv%instid(i)%tm(:,nread(i)) = 0.0
468 iv%instid(i)%qm(:,nread(i)) = 0.0
469 iv%instid(i)%qrn(:,nread(i)) = 0.0
470 iv%instid(i)%qcw(:,nread(i)) = 0.0
471 iv%instid(i)%qci(:,nread(i)) = 0.0
472 iv%instid(i)%qsn(:,nread(i)) = 0.0
473 iv%instid(i)%qgr(:,nread(i)) = 0.0
474 iv%instid(i)%pm(:,nread(i)) = 0.0
475 iv%instid(i)%u10(nread(i)) = 0.0
476 iv%instid(i)%v10(nread(i)) = 0.0
477 iv%instid(i)%t2m(nread(i)) = 0.0
478 iv%instid(i)%q2m(nread(i)) = 0.0
479 iv%instid(i)%mr2m(nread(i)) = 0.0
480 iv%instid(i)%psfc(nread(i)) = 0.0
481 iv%instid(i)%ts(nread(i)) = 0.0
482 iv%instid(i)%smois(nread(i)) = 0.0
483 iv%instid(i)%tslb(nread(i)) = 0.0
484 iv%instid(i)%snowh(nread(i)) = 0.0
485 iv%instid(i)%isflg(nread(i)) = 0
486 iv%instid(i)%soiltyp(nread(i)) = 0.0
487 iv%instid(i)%landsea_mask(nread(i)) = p%landsea_mask
488 iv%instid(i)%elevation(nread(i)) = 0.0
489 iv%instid(i)%vegfra(nread(i)) = 0.0
490 iv%instid(i)%vegtyp(nread(i)) = 0.0
491 iv%instid(i)%clwp(nread(i)) = 0.0
492 iv%instid(i)%ps(nread(i)) = 0.0
493 iv%instid(i)%tb_xb(:,nread(i)) = 0
494 iv%instid(i)%tb_inv(:,nread(i)) = p%tb_inv(:)
495 iv%instid(i)%tb_qc(:,nread(i)) = 0
496 iv%instid(i)%tb_error(:,nread(i)) = 500.
497 iv%instid(i)%emiss(:,nread(i)) = 0.0
498 iv%instid(i)%scanpos(nread(i)) = p%scanpos
499 ! iv%instid(i)%scanline(nread(i) = p%scanline
500 iv%instid(i)%scanline(nread(i)) = 0
501 iv%instid(i)%ifgat(nread(i)) = p%ifgat
502 iv%instid(i)%cloud_flag(:,nread(i)) = 1 ! no cloud
503 iv%instid(i)%satzen(nread(i)) = p%satzen
504 iv%instid(i)%satazi(nread(i)) = p%satazi
505 iv%instid(i)%solzen(nread(i)) = p%solzen
506 ! iv%instid(i)%solazi(nread(i)) = p%solazi
507 iv%instid(i)%solazi(nread(i)) = 0.0
508 iv%instid(i)%proc_domain(nread(i)) = .false.
509
510 current => p
511 p => p%next
512
513 ! free current data
514 deallocate ( current % tb_inv )
515 deallocate ( current )
516 end do
517
518 deallocate ( p )
519
520 deallocate (nread)
521
522 ! check if sequential structure has been freed
523 !
524 ! p => head
525 ! do i = 1, num_rad_used
526 ! write (unit=stdout,fmt=*) i, p%tb_inv(1:nchan)
527 ! p => p%next
528 ! end do
529
530 close(lnbufr)
531 call da_free_unit(lnbufr)
532
533 call da_trace_exit("da_read_bufrtovs")
534 #endif
535
536 end subroutine da_read_bufrtovs
537
538