da_scan_bufr_obs.inc
References to this file elsewhere.
1 subroutine da_scan_bufr_obs (ob, xp, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Scan BUFR observation file for input to wrfvar
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (xpose_type), intent(in) :: xp ! Domain decomposition vars.
10 type (ob_type), intent(inout) :: ob
11
12 character(*), intent (in) :: filename
13
14 #ifdef BUFR
15
16 type (multi_level_type) :: platform
17 logical :: outside
18 integer :: total_obs, &
19 num_sound, &
20 num_synop, &
21 num_pilot, &
22 num_geoamv, &
23 num_polaramv, &
24 num_satem, &
25 num_airep, &
26 num_metar, &
27 num_ships, &
28 num_gpspw, &
29 num_gpsref, &
30 num_ssmi_retrieval, &
31 num_ssmi_tb , &
32 num_ssmt1, &
33 num_ssmt2, &
34 num_pseudo, num_qscat, &
35 num_profiler, num_buoy, num_bogus
36 character(len=40) :: hdstr
37 character(len=10) :: date
38 character(len=8) :: subset
39 real,dimension(7) :: hdr
40
41 integer :: iost, ndup, n, report
42 integer :: year,month, day,hour
43 integer :: iret,idate, kx, nrecs,miscd, t29
44 integer :: iunit
45
46 ! default value
47 ob%ob_numb(ob%current_ob_time)%sound = ob%num_sound
48 ob%ob_numb(ob%current_ob_time)%synop = ob%num_synop
49 ob%ob_numb(ob%current_ob_time)%pilot = ob%num_pilot
50 ob%ob_numb(ob%current_ob_time)%satem = ob%num_satem
51 ob%ob_numb(ob%current_ob_time)%geoamv = ob%num_geoamv
52 ob%ob_numb(ob%current_ob_time)%polaramv = ob%num_polaramv
53 ob%ob_numb(ob%current_ob_time)%airep = ob%num_airep
54 ob%ob_numb(ob%current_ob_time)%gpspw = ob%num_gpspw
55 ob%ob_numb(ob%current_ob_time)%gpsref = ob%num_gpsref
56 ob%ob_numb(ob%current_ob_time)%metar = ob%num_metar
57 ob%ob_numb(ob%current_ob_time)%ships = ob%num_ships
58 ob%ob_numb(ob%current_ob_time)%qscat = ob%num_qscat
59 ob%ob_numb(ob%current_ob_time)%buoy = ob%num_buoy
60 ob%ob_numb(ob%current_ob_time)%profiler = ob%num_profiler
61
62 ob%ob_numb(ob%current_ob_time)%ssmt1 = ob%num_ssmt1
63 ob%ob_numb(ob%current_ob_time)%ssmt2 = ob%num_ssmt2
64 ! ob%ob_numb(ob%current_ob_time)%ssmi_tb = ob%num_ssmi_tb
65 ! ob%ob_numb(ob%current_ob_time)%ssmi_retrieval = ob%num_ssmi_retrieval
66
67 ! open FILE
68 ! ---------
69
70 call da_get_unit(iunit)
71 open(unit = iunit, FILE = trim(filename), &
72 iostat = iost, form = 'unformatted', STATUS = 'OLD')
73
74 if (iost /= 0) then
75 write(unit=message(1), fmt='(A,I3,A,A)') &
76 'ERROR in OBS inPUT FILE unit ',iunit, &
77 'OBS FILENAME:', trim(filename)
78 message(2)="BUFR GTS OBSERVATIONS CANNOT BE FOUND OR CANNOT BE openED"
79 call da_warning(__FILE__,__LINE__,message(1:2))
80 call da_free_unit(iunit)
81 return
82 end if
83
84 !--------------------------------------------------------------------------
85 ! Initialise the obs-counter
86 !--------------------------------------------------------------------------
87
88 total_obs = 0
89 num_sound = 0
90 num_synop = 0
91 num_pilot = 0
92 num_geoamv = 0
93 num_polaramv = 0
94 num_satem = 0
95 num_airep = 0
96 num_metar = 0
97 num_ships = 0
98 num_gpspw = 0
99 num_gpsref = 0
100 num_ssmi_retrieval = 0
101 num_ssmi_tb = 0
102 num_ssmt1 = 0
103 num_ssmt2 = 0
104 num_pseudo = 0
105 num_qscat = 0
106 num_profiler = 0
107 num_buoy = 0
108 num_bogus = 0
109
110 hdstr='SID XOB YOB DHR TYP ELV T29 '
111 nrecs = 0
112 miscd = 0
113
114 !--------------------------------
115 ! open bufr file then check date
116 !--------------------------------
117
118 call datelen(10)
119
120 call openbf(iunit,'in',iunit)
121 call readmg(iunit,subset,idate,iret)
122
123 if (iret/=0) then
124 message(1) = "No BUFR observations"
125 write(unit=message(2), fmt='(a, i4)') &
126 'return code from readmg:', iret, &
127 'Reach the end of obs unit: ', iunit
128 call da_error(__FILE__,__LINE__,message(1:2))
129 end if
130
131 write(unit=date,fmt='(i10)') idate
132 read (unit=date,fmt='(i4,3i2)') year,month, day,hour
133 write(unit=message(1),fmt=*) 'Bufr file date is ',year,month, day,hour
134 call da_message(message(1:1))
135
136 report = 0 ! report number in file
137
138 reports: do
139
140 report = report+1
141
142 call readsb(iunit,iret)
143
144 if (iret/=0) then
145 call readmg(iunit,subset,idate,iret)
146
147 if (iret/=0) then
148 write (unit=message(1), fmt='(a, i4)') &
149 'return code from readmg:', iret, &
150 'Reach the end of obs unit: ', iunit
151 call da_warning(__FILE__,__LINE__,message(1:1))
152 exit reports
153 end if
154
155 cycle reports
156 end if
157
158 nrecs=nrecs+1
159 call ufbint(iunit,hdr,7,1,iret,hdstr)
160
161 platform % info % name(1:8) = subset
162 platform % info % id = hdstr(1:5)
163
164 if (hdr(2) > 180.) hdr(2)=hdr(2)-360.
165
166 ! Put a check on Lat
167 hdr(3) = MAX(hdr(3), -89.95)
168 hdr(3) = Min(hdr(3), 89.95)
169 platform%info%lon = hdr(2)
170 platform%info%lat = hdr(3)
171
172 if (platform%info%lon == 180.0 ) platform%info%lon =-180.000
173 ! Fix funny wind direction at Poles
174 if (platform%info%lat <= -89.95 .or. platform%info%lat >= 89.95) then
175 platform%info%lon = 0.0
176 end if
177
178 ! Restrict to a range of reports, useful for debugging
179
180 if (report < report_start) then
181 cycle
182 end if
183
184 if (report > report_end) then
185 exit
186 end if
187
188 call da_ll_to_xy (platform%info, platform%loc, &
189 xp, outside)
190
191 if (outside) cycle reports
192
193 total_obs = total_obs + 1
194
195 t29 = int(0.1 + hdr(7))
196 kx=int(0.1+hdr(5))
197
198 if (kx .eq. 183) then ! reset kx based on t29
199 if (t29 .eq. 511) kx = 181
200 if (t29 .eq. 512) kx = 187
201 if (t29 .eq. 522) kx = 180
202 if (t29 .eq. 523) kx = 180
203 if (t29 .eq. 531) kx = 180
204 if (t29 .eq. 561) kx = 180
205 if (t29 .eq. 562) kx = 180
206 end if
207
208 ! if ((kx >= 160) .and. (kx <= 179)) then ! bypass satellite data
209 ! if (t29 .eq. 61 .or. t29 .eq. 63 .or. t29 .ge. 571) then
210 ! cycle reports
211 ! end if
212
213 ! Loop over duplicating obs for global
214 ndup = 1
215 if (global .and. &
216 (platform%loc%i < xp%ids .or. platform%loc%i >= xp%ide)) ndup= 2
217
218 if (Testing_WRFVAR) ndup = 1
219
220 do n = 1, ndup
221
222 select case(t29)
223
224 case (11, 12, 13, 22, 23, 31)
225 select case(kx)
226 case (120, 122, 132, 220, 222, 232) ; ! Sound
227 if (.not.use_SoundObs) cycle reports
228 num_sound = num_sound + 1
229
230 case (221) ; ! Pilot
231 if (.not.Use_PilotObs) cycle reports
232 num_pilot = num_pilot + 1
233 end select
234
235 case (41)
236 ! case (130:131, 133, 230:231, 233) ; ! Airep
237 if (.not.use_AirepObs) cycle reports
238 num_airep = num_airep + 1
239
240 ! case (180, 182, 280, 282) ; ! Ships and buoys
241 case (522, 523) ;
242 if (.not.use_ShipsObs) cycle reports
243 num_ships = num_ships + 1
244
245 case (531, 561, 562) ; ! Buoy
246 if (.not.use_BuoyObs) cycle reports
247 num_buoy = num_buoy + 1
248
249 case (511)
250 ! case (181, 281) ; ! Synop
251 if (.not.use_SynopObs) cycle reports
252 num_synop = num_synop + 1
253
254 case (512)
255 ! case (187 , 287) ; ! Metar
256 if (.not.use_MetarObs) cycle reports
257 num_metar = num_metar + 1
258
259 case (63)
260 ! case (242:246, 252:253, 255); ! Geo. CMVs
261 if (.not.use_GeoAMVObs) cycle reports
262 num_geoamv = num_geoamv + 1
263
264 case (71, 72) ; ! Profiler & VADWND - NEXRAD winds
265 if (.not.use_ProfilerObs) cycle reports
266 num_profiler = num_profiler + 1
267
268 case (582) ; ! Quikscat
269 if (.not.use_qscatobs) cycle reports
270 num_qscat = num_qscat + 1
271
272 case (583) ; ! GPS IPW
273 if (.not.use_GpspwObs) cycle reports
274 num_gpspw = num_gpspw + 1
275
276 case (584) ; ! GPS REF
277 if (.not.use_GpsrefObs) cycle reports
278 num_gpsref = num_gpsref + 1
279
280 case default ;
281 select case (kx)
282 case (111 , 210) ; ! Tropical Cyclone Bogus
283 ! Note Tropical cyclone Bougus is given type 135 in Obs-ascii
284 if (.not.use_BogusObs) cycle reports
285 num_bogus = num_bogus + 1
286
287 case default
288 miscd = miscd + 1
289 message(1)='unsaved obs found:'
290 write(unit=message(2), fmt='(2a)') &
291 'platform%info%platform=', platform%info%platform
292 call da_warning(__FILE__,__LINE__,message(1:2))
293 end select
294 end select
295 end do ! loop over duplicate
296 end do reports
297
298 rewind iunit
299 ! close(iunit)
300
301 write (unit=message(1),fmt=*) 'total bufr records =',nrecs
302 write (unit=message(2),fmt=*) 'number in domain =',total_obs
303 write (unit=message(3),fmt=*) 'number Unhandled =',miscd
304 call da_message(message(1:3))
305
306 !------------------------------------------------------------------------
307 ! Check the numbers again, make sure we have the right number.
308 !------------------------------------------------------------------------
309
310 ob%num_sound = ob%num_sound + num_sound
311 ob%num_synop = ob%num_synop + num_synop
312 ob%num_pilot = ob%num_pilot + num_pilot
313 ob%num_satem = ob%num_satem + num_satem
314 ob%num_geoamv = ob%num_geoamv + num_geoamv
315 ob%num_polaramv = ob%num_polaramv + num_polaramv
316 ob%num_airep = ob%num_airep + num_airep
317 ob%num_gpspw = ob%num_gpspw + num_gpspw
318 ob%num_gpsref = ob%num_gpsref + num_gpsref
319 ob%num_metar = ob%num_metar + num_metar
320 ob%num_ships = ob%num_ships + num_ships
321 ob%num_qscat = ob%num_qscat + num_qscat
322 ob%num_buoy = ob%num_buoy + num_buoy
323 ob%num_profiler = ob%num_profiler + num_profiler
324 ob%num_bogus = ob%num_bogus + num_bogus
325
326 ob%num_ssmt1 = ob%num_ssmt1 + num_ssmt1
327 ob%num_ssmt2 = ob%num_ssmt2 + num_ssmt2
328
329 ! ob%num_ssmi_tb = ob%num_ssmi_tb + num_ssmi_tb
330 ! ob%num_ssmi_retrieval = ob%num_ssmi_retrieval + num_ssmi_retrieval
331
332 ob%ob_numb(ob%current_ob_time)%sound = ob%num_sound
333 ob%ob_numb(ob%current_ob_time)%synop = ob%num_synop
334 ob%ob_numb(ob%current_ob_time)%pilot = ob%num_pilot
335 ob%ob_numb(ob%current_ob_time)%satem = ob%num_satem
336 ob%ob_numb(ob%current_ob_time)%geoamv = ob%num_geoamv
337 ob%ob_numb(ob%current_ob_time)%polaramv = ob%num_polaramv
338 ob%ob_numb(ob%current_ob_time)%airep = ob%num_airep
339 ob%ob_numb(ob%current_ob_time)%gpspw = ob%num_gpspw
340 ob%ob_numb(ob%current_ob_time)%gpsref = ob%num_gpsref
341 ob%ob_numb(ob%current_ob_time)%metar = ob%num_metar
342 ob%ob_numb(ob%current_ob_time)%ships = ob%num_ships
343 ob%ob_numb(ob%current_ob_time)%qscat = ob%num_qscat
344 ob%ob_numb(ob%current_ob_time)%buoy = ob%num_buoy
345 ob%ob_numb(ob%current_ob_time)%profiler = ob%num_profiler
346 ob%ob_numb(ob%current_ob_time)%bogus = ob%num_bogus
347
348 ob%ob_numb(ob%current_ob_time)%ssmt1 = ob%num_ssmt1
349 ob%ob_numb(ob%current_ob_time)%ssmt2 = ob%num_ssmt2
350
351 ! ob%ob_numb(ob%current_ob_time)%ssmi_tb = ob%num_ssmi_tb
352 ! ob%ob_numb(ob%current_ob_time)%ssmi_retrieval = ob%num_ssmi_retrieval
353
354 ! PRinT out
355 ! =============
356
357 if (print_detail_obs) then
358 write(unit=stdout, fmt='(a)') ' '
359
360 write(unit=stdout, fmt='(5x,a,i6,a)') &
361 'Scan: ', ob%current_ob_time, ' ob time:', &
362 'Scan: ', total_obs, ' Total reports.', &
363 'Scan: ', num_sound, ' SOUND reports,', &
364 'Scan: ', num_synop, ' SYNOP reports,', &
365 'Scan: ', num_pilot, ' PILOT reports,', &
366 'Scan: ', num_satem, ' SATEM reports,', &
367 'Scan: ', num_geoamv, ' Geo. AMVs reports,', &
368 'Scan: ', num_polaramv, ' Polar AMVs reports,', &
369 'Scan: ', num_airep, ' AIREP reports,', &
370 'Scan: ', num_gpspw, ' GPSPW reports,', &
371 'Scan: ', num_gpsref, ' GPSRF reports,', &
372 'Scan: ', num_metar, ' METAR reports,', &
373 'Scan: ', num_ships, ' SHIP reports,', &
374 'Scan: ', num_qscat, ' QSCAT reports,', &
375 'Scan: ', num_profiler, ' Profiler reports,', &
376 'Scan: ', num_buoy, ' Buoy reports,', &
377 'Scan: ', num_bogus, ' TCBGS reports,', &
378 'Scan: ', num_ssmt1, ' SSMT1 reports,', &
379 'Scan: ', num_ssmt2, ' SSMT2 reports.'
380 write(unit=stdout,fmt=*) ' '
381
382 ! write(unit=stdout, fmt='(5x,a,i6,a)') &
383 ! 'Scan: ', num_ssmi_retrieval , ' SSMI_RETRIEVAL reports,', &
384 ! 'Scan: ', num_ssmi_tb , ' SSMI_TB reports.'
385
386 end if
387
388 close(iunit)
389 call da_free_unit(iunit)
390 #endif
391
392 end subroutine da_scan_bufr_obs
393
394