da_read_obs_bufr.inc
References to this file elsewhere.
1 subroutine da_read_obs_bufr (iv, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Read BUFR observation file for input to wrfvar
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (iv_type), intent(inout) :: iv
10 character(len=*), optional, intent(in) :: filename
11
12 #ifdef BUFR
13
14 type (multi_level_type) :: platform
15 logical :: outside, outside_all
16 integer :: nlocal(num_ob_indexes)
17 integer :: ntotal(num_ob_indexes)
18
19 character(len=40) :: obstr,drift,hdstr,qmstr,oestr
20 character(len=8) :: subset, c8id
21 real :: hdr(7)
22 real :: pmo(2,1)
23 ! real :: drf9,255)
24 real :: obs(9,255),qms(9,255),oes(9,255)
25 real :: time,woe,qob,toe,qoe,poe,pob,tob,zob,rob,roe
26
27 integer :: iost, ndup, n, i, j, k, surface_level, report
28 integer :: iret, idate, kx, nlevels, nrecs,miscd, t29, jx, ix
29 integer :: cat,zqm,pqm,qqm,tqm,wqm,pwq,pmq,rqm
30 integer :: iunit
31
32 integer, parameter :: qflag=4 ! Flag for retaining data
33 equivalence (hdr(1), c8id)
34
35 if (trace_use) call da_trace_entry("da_read_obs_bufr")
36
37 nlocal(:) = iv%info(:)%plocal(iv%time-1)
38 ntotal(:) = iv%info(:)%ptotal(iv%time-1)
39
40 ! open file
41 ! ---------
42 call da_get_unit(iunit)
43 if (present(filename)) then
44 open(unit = iunit, FILE = trim(filename), &
45 iostat = iost, form = 'unformatted', STATUS = 'OLD')
46 if (iost /= 0) then
47 write(unit=stdout, fmt='(/A,I3,3(2X,A)/)') &
48 'error in obs input file unit ',iunit, &
49 'obs filename:', trim(filename), &
50 'for gts observations cannot be found or cannot be opened'
51 return
52 end if
53 end if
54
55 hdstr='SID XOB YOB DHR TYP ELV T29 '
56 obstr='POB QOB TOB ZOB UOB VOB PWO ROB CAT '
57 drift='XDR YDR HRDR '
58 qmstr='PQM QQM TQM ZQM WQM NUL PWQ PMQ RQM '
59 oestr='POE QOE TOE NUL WOE NUL PWE ROE '
60
61 nrecs = 0
62 miscd = 0
63
64 ! read bufr observation file
65 ! call openbf(iunit,'in',iunit)
66 call readmg(iunit,subset,idate,iret)
67
68 if (iret/=0) then
69 write(unit=stdout, fmt='(a, i4)') &
70 'return code from readmg:', iret, &
71 'Reach the end of obs unit: ', iunit
72
73 call da_error(__FILE__,__LINE__, &
74 (/"No BUFR observations"/))
75 end if
76
77 report = 0 ! report number in file
78 reports: do
79 report = report+1
80
81 call readsb(iunit,iret)
82
83 if (iret/=0) then
84 call readmg(iunit,subset,idate,iret)
85
86 if (iret/=0) then
87 write(unit=stdout, fmt='(a, i4)') &
88 'return code from readmg:', iret, &
89 'Reach the end of prepbufr obs unit: ', iunit
90
91 exit reports
92 end if
93
94 cycle reports
95 end if
96
97 nrecs=nrecs+1
98
99 call ufbint(iunit,hdr,7,1,iret,hdstr)
100
101 platform % info % name(1:8) = subset
102
103 platform % info % id = c8id(1:5)
104 if (hdr(2) > 180.0) hdr(2)=hdr(2)-360.0
105
106 ! Put a check on Lat
107
108 hdr(3) = MAX(hdr(3), -89.95)
109 hdr(3) = Min(hdr(3), 89.95)
110 platform%info%lon = hdr(2)
111 platform%info%lat = hdr(3)
112
113 ! Put a check on Lat
114
115 if (platform%info%lon == 180.0 ) platform%info%lon =-180.000
116 ! Fix funny wind direction at Poles
117 if (platform%info%lat <= -89.95 .or. platform%info%lat >= 89.95) then
118 platform%info%lon = 0.0
119 end if
120
121 ! Restrict to a range of reports, useful for debugging
122
123 if (report < report_start) then
124 cycle
125 end if
126
127 if (report > report_end) then
128 exit
129 end if
130
131 call da_llxy (platform%info, platform%loc,outside, outside_all)
132
133 if (outside_all) cycle reports
134
135 if (print_detail_obs) then
136 ! Simplistic approach, could be improved to get it all done on PE 0
137 if (.NOT. outside) then
138 write(unit=stdout,fmt='(A,I4,A,F8.2,A,F8.2,A,I3,A,2I3)') &
139 "Report",report,"at",platform%info%lon,"E",platform%info%lat,"N", &
140 "on processor", myproc,"position", platform%loc%x,platform%loc%y
141 end if
142 end if
143
144 time=hdr(4)
145 if (time > 3.0) time=3.0
146 if (time < -3.0) time=-3.0
147
148 t29 = int(0.1 + hdr(7))
149 kx=int(0.1+hdr(5))
150
151 if (kx == 183) then ! reset kx based on t29
152 if (t29 .eq. 511) kx = 181
153 if (t29 .eq. 512) kx = 187
154 if (t29 .eq. 522) kx = 180
155 if (t29 .eq. 523) kx = 180
156 if (t29 .eq. 531) kx = 180
157 if (t29 .eq. 561) kx = 180
158 if (t29 .eq. 562) kx = 180
159 end if
160
161 ! WHY?
162 ! if ((kx >= 160) .and. (kx <= 179)) then ! bypass satellite data
163 ! if (t29 .eq. 61 .or. t29 .eq. 63 .or. t29 .ge. 571) then
164 ! cycle reports
165 ! end if
166
167 ! Conventional data
168
169 call ufbint(iunit,pmo,2,1,nlevels,'PMO PMQ')
170 ! write(unit=stdout,fmt=*) 'PMO = ',pmo(1,1)
171
172 call ufbint(iunit,obs,9,255,nlevels,obstr)
173 call ufbint(iunit,qms,9,255,nlevels,qmstr)
174 call ufbint(iunit,oes,9,255,nlevels,oestr)
175 if (c8id(1:5) .eq. 'RJOR ') then
176 ! if (kx .eq. 180) then
177 write(unit=stdout,fmt=*) 'report, ',c8id(1:5),' typ = ',kx, &
178 ' nlevels = ',nlevels
179 do ix = 1,nlevels
180 write(unit=stdout,fmt='(9(f10.4,2x))') (obs(jx,ix),jx=1,9)
181 write(unit=stdout,fmt='(9(f10.4,2x))') (qms(jx,ix),jx=1,9)
182 write(unit=stdout,fmt=*) 'PMO = ',pmo(1,1)
183 end do
184 ! kx = 180
185 ! end if
186 end if
187
188 ! WHY?
189 ! no need for balloon drift
190 ! select case(kx)
191 ! case (120, 132, 220:221, 232) ; ! read drift
192 ! call ufbint(iunit,drf,9,255,iret,drift)
193 ! do k=1,nlevels
194 ! if (drf(1,k) .ge. 360.0)drf(1,k)=drf(1,k)-360.
195 ! if (drf(1,k) .lt. 0.0)drf(1,k)=drf(1,k)+360.
196 ! if (abs(drf(2,k)) .gt. 1000.0 .or. abs(drf(1,k)) .gt. 1000.0) then
197 ! drf(2,k)=hdr(3)
198 ! drf(1,k)=hdr(2)
199 ! end if
200 ! if (abs(drf(3,k)) .gt. 3.0)then
201 ! drf(3,k)=time
202 ! end if
203 ! end do
204 ! case default ; ! Do nothing
205 ! end select
206
207 platform % info % levels =nlevels
208 platform % info % elv =hdr(6)
209
210 platform % loc % slp %inv = missing_r
211 platform % loc % slp %qc = missing_data
212 platform % loc % slp %error=200.0
213 platform % loc % pw %inv = missing_r
214 platform % loc % pw %qc = missing_data
215 platform % loc % pw %error=0.2
216
217 do i=1,max_ob_levels
218 platform % each (i) % height = missing_r
219 platform % each (i) % height_qc = missing_data
220
221 platform % each (i) % zk = missing_r
222
223 platform % each (i) % u % inv = missing_r
224 platform % each (i) % u % qc = missing_data
225 platform % each (i) % u % error =10.0
226
227 platform % each (i) % v = platform % each (i) % u
228
229 platform % each (i) % t % inv = missing_r
230 platform % each (i) % t % qc = missing_data
231 platform % each (i) % t % error =5.0
232
233 platform % each (i) % p % inv = missing_r
234 platform % each (i) % p % qc = missing_data
235 platform % each (i) % p % error =200.0
236
237 platform % each (i) % q % inv = missing_r
238 platform % each (i) % q % qc = missing_data
239 platform % each (i) % q % error =0.1
240 end do
241
242 do k = 1, platform % info % levels
243 pob=obs(1,k)
244 ! scale q and compute t from tv and set units to Kelvin, if they aren't missing
245
246 if (obs(2,k) > 0.0 .and. obs(2,k) < 1e9) then
247 obs(2,k) = obs(2,k)*1e-6
248 if (obs(3,k) > -200.0 .and. obs(3,k) < 300.0) then
249 obs(3,k) = (obs(3,k)+t_kelvin) / (1.0 + 0.608 * obs(2,k))
250 end if
251 else
252 if (obs(3,k) > -200.0 .and. obs(3,k) < 300.0) then
253 obs(3,k) = obs(3,k) + t_kelvin
254 end if
255 end if
256 qob=obs(2,k)
257 tob=obs(3,k)
258 zob=obs(4,k)
259 rob=obs(8,k)
260 pqm = qflag
261 if (abs(qms(1,k)) < qflag) pqm=nint(qms(1,k))
262 qqm = qflag
263 if (abs(qms(2,k)) < qflag) qqm=nint(qms(2,k))
264 tqm = qflag
265 if (abs(qms(3,k)) < qflag) tqm=nint(qms(3,k))
266 zqm = qflag
267 if (abs(qms(4,k)) < qflag) zqm=nint(qms(4,k))
268 wqm = qflag
269 if (abs(qms(5,k)) < qflag) wqm=nint(qms(5,k))
270 pwq = qflag
271 if (abs(qms(7,k)) < qflag) pwq=nint(qms(7,k))
272 rqm = qflag
273 if (abs(qms(9,k)) < qflag) rqm=nint(qms(9,k))
274 pmq = qflag
275 if (abs(pmo(2,1)) < qflag) pmq=nint(pmo(2,1))
276 cat=nint(obs(9,k))
277
278 if (pmq < qflag) then
279 platform % loc % slp % inv =pmo(1,1)*100.0
280 platform % loc % slp % qc =pmq
281 platform % loc % slp % error = 200.0
282 end if
283
284 if (pwq < qflag) then
285 platform % loc % pw % inv = obs(7,k) * 0.1 ! convert to cm
286 platform % loc % pw % qc = pwq
287 platform % loc % pw % error = 0.2 ! hardwired to 0.2 cm
288 end if
289
290 poe=100.0
291
292 if (tqm<qflag) then
293 toe=oes(3,k)
294 platform % each (k) % t % inv =obs(3,k)
295 platform % each (k) % t % qc =tqm
296 platform % each (k) % t % error =toe
297 end if
298
299 if (wqm < qflag) then
300 woe=oes(5,k)
301 ! WHY?
302 ! if (wqm==0) woe=woe*0.9
303 ! if (wqm==3) woe=woe*1.2
304 platform % each (k) % u % inv =obs(5,k)
305 platform % each (k) % v % inv =obs(6,k)
306 platform % each (k) % u % qc =wqm
307 platform % each (k) % u % error =woe
308 platform % each (k) % v % qc =wqm
309 platform % each (k) % v % error =woe
310
311 ! Convert earth wind to model wind.
312 call da_earth_2_model_wind(obs(5,k), obs(6,k), &
313 platform % each (k) % u % inv, &
314 platform % each (k) % v % inv, &
315 platform%info%lon)
316 if (platform % each (k) % u % inv == 0.0 .and. platform % each (k) % v % inv == 0.0) then
317 platform % each (k) % u % inv = missing_r
318 platform % each (k) % v % inv = missing_r
319 platform % each (k) % u % qc = missing_data
320 platform % each (k) % v % qc = missing_data
321 end if
322 end if
323
324 if (qqm<qflag .and. obs(2,k)<1.0e9) then
325 qoe=oes(2,k)
326 platform % each (k) % q % inv =obs(2,k)
327 platform % each (k) % q % qc =qqm
328 platform % each (k) % q % error =qoe
329 end if
330
331 if (zqm<qflag)then
332 platform % each (k) % height =zob
333 platform % each (k) % height_qc =zqm
334 end if
335
336 if (rqm<qflag)then
337 roe = oes(8,k)
338 platform % each (k) % u % inv =rob
339 platform % each (k) % u % qc =rqm
340 platform % each (k) % u % error =roe
341 end if
342
343 if (pqm<qflag)then
344 platform % each (k) % p % inv =obs(1,k)*100.0
345 platform % each (k) % p % qc =pqm
346 platform % each (k) % p % error =poe
347 end if
348 end do
349
350 nlevels = platform%info%levels
351
352 if (nlevels > max_ob_levels) then
353 nlevels = max_ob_levels
354
355 write(unit=stderr, fmt='(/a/)') 'Warning: Too many levels.'
356
357 write(unit=stderr, fmt='(/2a/2a/2x,a,2f8.2,2(2x,a,f9.2)/2(2x,a,i4)/)') &
358 'Subset: ', platform%info%name(1:8), &
359 'Platfrom: ', trim(platform%info%platform), &
360 'Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
361 'elv: ', platform%info%elv, &
362 'pstar: ', platform%info%pstar, &
363 'level: ', platform%info%levels, &
364 'kx: ', kx
365 else if ((nlevels < 1) .and. ((kx /= 164) .or. (kx /= 174))) then
366 write(unit=stderr, fmt='(/a/)') &
367 'Warning: Too few levels.'
368
369 write(unit=stderr, fmt='(/2a/2a/2x,a,2f8.2,2(2x,a,f9.2)/2(2x,a,i4)/)') &
370 'Subset: ', platform%info%name(1:8), &
371 'Platfrom: ', trim(platform%info%platform), &
372 'Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
373 'elv: ', platform%info%elv, &
374 'pstar: ', platform%info%pstar, &
375 'level: ', platform%info%levels, &
376 'kx: ', kx
377
378 cycle reports
379 end if
380
381 !---------------------------------------------------------------------------
382 ! This is basically converting rh to q i
383 ! Method :
384 ! if rh, temp and pr all available computes Qs otherwise sets Qs= missing
385 ! if rh > 100 sets q = qs otherwise q = rh*Qs/100.0
386 ! Note: Currently da_obs_proc_station is active only for ob_format_ascii
387 ! call da_obs_proc_station(platform)
388 !---------------------------------------------------------------------------
389
390 ! Loop over duplicating obs for global
391 ndup = 1
392 if (global .and. &
393 (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
394
395 ! It is possible that logic for counting obs is incorrect for the
396 ! global case with >1 MPI tasks due to obs duplication, halo, etc.
397 ! TBH: 20050913
398 if (.not.outside) then
399 if (print_detail_obs .and. ndup > 1) then
400 write(unit=stdout, fmt = '(A12,1X,A19,1X,A40,1X,I6,3(F12.3,11X),6X,A5)') &
401 platform%info%platform, &
402 platform%info%date_char, &
403 platform%info%name, &
404 platform%info%levels, &
405 platform%info%lat, &
406 platform%info%lon, &
407 platform%info%elv, &
408 platform%info%id
409
410 write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
411 ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
412 platform%loc%i, platform%loc%j, &
413 platform%loc%dx, platform%loc%dxm, &
414 platform%loc%dy, platform%loc%dym
415 end if
416 end if
417
418 dup_loop: do n = 1, ndup
419 select case(t29)
420 case (11, 12, 13, 22, 23, 31)
421 select case (kx)
422 case (120, 122, 132, 220, 222, 232) ; ! Sound
423 ntotal(sound) = ntotal(sound) + 1
424 ntotal(sonde_sfc) = ntotal(sonde_sfc) + 1
425 if (outside) then
426 cycle reports
427 end if
428
429 if (.not.use_soundobs) cycle reports
430
431 platform % info % platform ='FM-35 TEMP'
432 nlocal(sound) = nlocal(sound) + 1
433 nlocal(sonde_sfc) = nlocal(sonde_sfc) + 1
434 ! Track serial obs index for reassembly of obs during
435 ! bit-for-bit tests with different numbers of MPI tasks.
436
437 ! Search to see if we have surface obs.
438
439 surface_level = 0
440
441 do i = 1, nlevels
442 ! if (elevation and height are the same, it is surface)
443 if (abs(platform%info%elv - &
444 platform%each(i)%height) < 0.1) then
445 surface_level = i
446
447 ! Save surface pressure.
448 iv%sonde_sfc(nlocal(sonde_sfc))%h = platform%each(i)%height
449 iv%sonde_sfc(nlocal(sonde_sfc))%u = platform%each(i)%u
450 iv%sonde_sfc(nlocal(sonde_sfc))%v = platform%each(i)%v
451 iv%sonde_sfc(nlocal(sonde_sfc))%t = platform%each(i)%t
452 iv%sonde_sfc(nlocal(sonde_sfc))%q = platform%each(i)%q
453 iv%sonde_sfc(nlocal(sonde_sfc))%p = platform%each(i)%p
454 exit
455 end if
456 end do
457
458 ! processing the sound_sfc data:
459
460 if (surface_level > 0) then
461 nlevels = nlevels - 1
462 else
463 iv%sonde_sfc(nlocal(sonde_sfc))%h = missing_r
464 iv%sonde_sfc(nlocal(sonde_sfc))%u%inv = missing_r
465 iv%sonde_sfc(nlocal(sonde_sfc))%u%qc = missing
466 iv%sonde_sfc(nlocal(sonde_sfc))%u%error = abs(missing_r)
467 iv%sonde_sfc(nlocal(sonde_sfc))%v = iv%sonde_sfc(nlocal(sonde_sfc))%u
468 iv%sonde_sfc(nlocal(sonde_sfc))%t = iv%sonde_sfc(nlocal(sonde_sfc))%u
469 iv%sonde_sfc(nlocal(sonde_sfc))%p = iv%sonde_sfc(nlocal(sonde_sfc))%u
470 iv%sonde_sfc(nlocal(sonde_sfc))%q = iv%sonde_sfc(nlocal(sonde_sfc))%u
471 end if
472
473 if (nlevels < 1) cycle reports
474
475 allocate (iv%sound(nlocal(sound))%h (1:nlevels))
476 allocate (iv%sound(nlocal(sound))%p (1:nlevels))
477 allocate (iv%sound(nlocal(sound))%u (1:nlevels))
478 allocate (iv%sound(nlocal(sound))%v (1:nlevels))
479 allocate (iv%sound(nlocal(sound))%t (1:nlevels))
480 allocate (iv%sound(nlocal(sound))%q (1:nlevels))
481
482 j = 0
483 do i = 1, nlevels
484 if (i == surface_level) cycle
485
486 j=j+1
487
488 iv%sound(nlocal(sound))%h(j) = platform%each(i)%height
489 iv%sound(nlocal(sound))%p(j) = platform%each(i)%p%inv
490 iv%sound(nlocal(sound))%u(j) = platform%each(i)%u
491 iv%sound(nlocal(sound))%v(j) = platform%each(i)%v
492 iv%sound(nlocal(sound))%t(j) = platform%each(i)%t
493 iv%sound(nlocal(sound))%q(j) = platform%each(i)%q
494 end do
495
496 case (221) ; ! Pilot
497 ntotal(pilot) = ntotal(pilot) + 1
498 if (outside) then
499 cycle reports
500 end if
501
502 if (.not.use_pilotobs) cycle reports
503 nlocal(pilot) = nlocal(pilot) + 1
504 ! Track serial obs index for reassembly of obs during
505 ! bit-for-bit tests with different numbers of MPI tasks.
506
507 allocate (iv%pilot(nlocal(pilot))%p (1:nlevels))
508 allocate (iv%pilot(nlocal(pilot))%u (1:nlevels))
509 allocate (iv%pilot(nlocal(pilot))%v (1:nlevels))
510
511 do i = 1, nlevels
512 iv%pilot(nlocal(pilot))%p(i) = platform%each(i)%p%inv
513 iv%pilot(nlocal(pilot))%u(i) = platform%each(i)%u
514 iv%pilot(nlocal(pilot))%v(i) = platform%each(i)%v
515 end do
516 end select
517
518 case (41)
519 ! case (130:131, 133, 230:231, 233) ; ! Airep
520 ntotal(airep) = ntotal(airep) + 1
521 if (outside) then
522 cycle reports
523 end if
524
525 if (.not.use_airepobs) cycle reports
526
527 platform % info % platform ='FM-97 AIREP'
528 nlocal(airep) = nlocal(airep) + 1
529 ! Track serial obs index for reassembly of obs during
530 ! bit-for-bit tests with different numbers of MPI tasks.
531 platform%loc%obs_global_index = ntotal(airep)
532
533 allocate (iv % airep (nlocal(airep)) % h (1:nlevels))
534 allocate (iv % airep (nlocal(airep)) % p (1:nlevels))
535 allocate (iv % airep (nlocal(airep)) % u (1:nlevels))
536 allocate (iv % airep (nlocal(airep)) % v (1:nlevels))
537 allocate (iv % airep (nlocal(airep)) % t (1:nlevels))
538
539 do i = 1, nlevels
540 iv % airep (nlocal(airep)) % h(i) = platform % each(i) % height
541 iv % airep (nlocal(airep)) % p(i) = platform % each(i) % p % inv
542 iv % airep (nlocal(airep)) % u(i) = platform % each(i) % u
543 iv % airep (nlocal(airep)) % v(i) = platform % each(i) % v
544 iv % airep (nlocal(airep)) % t(i) = platform % each(i) % t
545 end do
546
547 ! case (180, 182, 280, 282) ; ! Ships and buoys
548
549 case (522, 523); ! Ships
550 ntotal(ships) = ntotal(ships) + 1
551 if (outside) then
552 cycle reports
553 end if
554
555 if (.not.use_shipsobs) cycle reports
556
557 platform % info % platform ='FM-13 SHIP '
558 nlocal(ships) = nlocal(ships) + 1
559 ! Track serial obs index for reassembly of obs during
560 ! bit-for-bit tests with different numbers of MPI tasks.
561 platform%loc%obs_global_index = ntotal(ships)
562
563 iv % ships (nlocal(ships)) % h = platform % each(1) % height
564 iv % ships (nlocal(ships)) % u = platform % each(1) % u
565 iv % ships (nlocal(ships)) % v = platform % each(1) % v
566 iv % ships (nlocal(ships)) % t = platform % each(1) % t
567 ! WHY?
568 ! iv % ships (nlocal(ships)) % p = platform % loc % slp
569 iv % ships (nlocal(ships)) % p = platform % each(1) % p
570 iv % ships (nlocal(ships)) % q = platform % each(1) % q
571
572 case (531, 561, 562) ; ! Buoy
573 ntotal(buoy) = ntotal(buoy) + 1
574 if (outside) then
575 cycle reports
576 end if
577
578 if (.not.use_buoyobs) cycle reports
579 platform % info % platform ='FM-18 BUOY '
580 nlocal(buoy) = nlocal(buoy) + 1
581 platform%loc%obs_global_index = ntotal(buoy)
582
583 iv%buoy(nlocal(buoy))%h = platform%each(1)%height
584 iv%buoy(nlocal(buoy))%u = platform%each(1)%u
585 iv%buoy(nlocal(buoy))%v = platform%each(1)%v
586 iv%buoy(nlocal(buoy))%t = platform%each(1)%t
587 iv%buoy(nlocal(buoy))%p = platform%each(1)%p
588 iv%buoy(nlocal(buoy))%q = platform%each(1)%q
589
590 case (511)
591 ! case (181, 281) ; ! Synop
592
593 ntotal(synop) = ntotal(synop) + 1
594 if (outside) then
595 cycle reports
596 end if
597
598 if (.not.use_synopobs) cycle reports
599
600 platform % info % platform ='FM-12 SYNOP'
601
602 nlocal(synop) = nlocal(synop) + 1
603 platform%loc%obs_global_index = ntotal(synop)
604
605 iv % synop (nlocal(synop)) % h = platform % each(1) % height
606 iv % synop (nlocal(synop)) % u = platform % each(1) % u
607 iv % synop (nlocal(synop)) % v = platform % each(1) % v
608 iv % synop (nlocal(synop)) % t = platform % each(1) % t
609 ! WHY?
610 ! iv % synop (nlocal(synop)) % p = platform % loc % slp
611 iv % synop (nlocal(synop)) % p = platform % each(1) % p
612 iv % synop (nlocal(synop)) % q = platform % each(1) % q
613
614 if (iv % synop(nlocal(synop)) % h < platform % info % elv) then
615 iv % synop(nlocal(synop)) % h = platform % info % elv
616 end if
617
618 case (512)
619 ! case (187, 287) ; ! Metar
620 ntotal(metar) = ntotal(metar) + 1
621 if (outside) then
622 cycle reports
623 end if
624
625 if (.not.use_metarobs) cycle reports
626
627 platform % info % platform ='FM-15 METAR'
628 nlocal(metar) = nlocal(metar) + 1
629 platform%loc%obs_global_index = ntotal(metar)
630
631 iv % metar (nlocal(metar)) % h = platform % each(1) % height
632 iv % metar (nlocal(metar)) % u = platform % each(1) % u
633 iv % metar (nlocal(metar)) % v = platform % each(1) % v
634 iv % metar (nlocal(metar)) % t = platform % each(1) % t
635 iv % metar (nlocal(metar)) % p = platform % each(1) % p
636 iv % metar (nlocal(metar)) % q = platform % each(1) % q
637
638 case (63)
639 ! case (242:246, 252:253, 255) ; ! Geo. CMVs
640 ntotal(geoamv) = ntotal(geoamv) + 1
641 if (outside) then
642 cycle reports
643 end if
644
645 if (.not.use_geoamvobs) cycle reports
646 platform % info % platform ='FM-88 SATOB'
647 nlocal(geoamv) = nlocal(geoamv) + 1
648 platform%loc%obs_global_index = ntotal(geoamv)
649
650 allocate (iv%geoamv(nlocal(geoamv))%p (1:nlevels))
651 allocate (iv%geoamv(nlocal(geoamv))%u (1:nlevels))
652 allocate (iv%geoamv(nlocal(geoamv))%v (1:nlevels))
653
654 do i = 1, nlevels
655 iv % geoamv (nlocal(geoamv)) % p(i) = platform % each(i) % p % inv
656 iv % geoamv (nlocal(geoamv)) % u(i) = platform % each(i) % u
657 iv % geoamv (nlocal(geoamv)) % v(i) = platform % each(i) % v
658 iv % info(geoamv)%zk(i,nlocal(geoamv)) = platform % each(i) % zk
659 end do
660
661 case (582)
662 ntotal(qscat) = ntotal(qscat) + 1
663 if (outside) then
664 cycle reports
665 end if
666
667 if (.not.use_qscatobs) cycle reports
668 platform % info % platform ='FM-281 Quiks'
669 nlocal(qscat) = nlocal(qscat) + 1
670
671 ! WHY?
672 ! iv%qscat(nlocal(qscat))%h = platform%each(1)%height
673 ! prepbufr uses pressure not height, so hardwire height to
674 ! 0 (sea-level)
675 iv%qscat(nlocal(qscat))%h = 0.0
676 iv%qscat(nlocal(qscat))%u = platform%each(1)%u
677 iv%qscat(nlocal(qscat))%v = platform%each(1)%v
678 iv%qscat(nlocal(qscat))%u%error = max(platform%each(1)%u%error,1.0)
679 iv%qscat(nlocal(qscat))%v%error = max(platform%each(1)%v%error,1.0)
680
681 case (583) ! GPS PW
682 ntotal(gpspw) = ntotal(gpspw) + 1
683 if (outside) then
684 cycle reports
685 end if
686
687 if (.not.use_gpspwobs) cycle reports
688 platform % info % platform ='FM-111 GPSPW'
689 nlocal(gpspw) = nlocal(gpspw) + 1
690
691 iv%gpspw(nlocal(gpspw))%tpw = platform%loc%pw
692
693 case (584) ! GPS REF
694 ntotal(gpsref) = ntotal(gpsref) + 1
695 if (outside) then
696 cycle reports
697 end if
698
699 if (.not.use_gpsrefobs) cycle reports
700 platform % info % platform ='FM-116 GPSRF'
701 nlocal(gpsref) = nlocal(gpsref) + 1
702
703 allocate (iv%gpsref(nlocal(gpsref))%h (1:nlevels))
704 allocate (iv%gpsref(nlocal(gpsref))%ref(1:nlevels))
705 allocate (iv%gpsref(nlocal(gpsref))% p(1:nlevels))
706 allocate (iv%gpsref(nlocal(gpsref))% t(1:nlevels))
707 allocate (iv%gpsref(nlocal(gpsref))% q(1:nlevels))
708
709 do i = 1, nlevels
710 iv%gpsref(nlocal(gpsref))%h(i) = platform%each(i)%height
711 ! u is used to store ref in this routine (ascii version uses td)
712 iv%gpsref(nlocal(gpsref))%ref(i) = platform%each(i)%u
713 ! Keep the retrieved p and t (and q) as "field_type":
714 iv%gpsref(nlocal(gpsref))%p(i) = platform%each(i)%p
715 iv%gpsref(nlocal(gpsref))%t(i) = platform%each(i)%t
716 iv%gpsref(nlocal(gpsref))%q(i) = platform%each(i)%q
717 end do
718
719 case (71, 72)
720 ! case (223, 224 ) ; ! Profiler & VADWND - NEXRAD winds
721 ntotal(profiler) = ntotal(profiler)
722 if (outside) then
723 cycle reports
724 end if
725
726 if (.not.use_profilerobs) cycle reports
727 platform % info % platform ='FM-132 PRFLR'
728 nlocal(profiler) = nlocal(profiler) + 1
729 ! Track serial obs index for reassembly of obs during bit-for-bit
730 ! tests with different numbers of MPI tasks.
731
732 allocate (iv%profiler(nlocal(profiler))%p (1:nlevels))
733 allocate (iv%profiler(nlocal(profiler))%u (1:nlevels))
734 allocate (iv%profiler(nlocal(profiler))%v (1:nlevels))
735
736 do i = 1, nlevels
737 iv%profiler(nlocal(profiler))%p(i) = platform%each(i)%p%inv
738 iv%profiler(nlocal(profiler))%u(i) = platform%each(i)%u
739 iv%profiler(nlocal(profiler))%v(i) = platform%each(i)%v
740 end do
741
742 case default
743 select case (kx)
744 case (111 , 210) ; ! Tropical Cyclone Bogus
745 ! Note Tropical cyclone Bougus is given type 135 in Obs-ascii
746 ntotal(bogus) = ntotal(bogus) + 1
747 if (outside) then
748 cycle reports
749 end if
750
751 if (.not.use_bogusobs) cycle reports
752 platform % info % platform ='FM-135 TCBOG'
753 nlocal(bogus) = nlocal(bogus) + 1
754 ! Track serial obs index for reassembly of obs during bit-for-bit
755 ! tests with different numbers of MPI tasks.
756 platform%loc%obs_global_index = ntotal(bogus)
757
758 if (nlocal(bogus) > max_bogus_input) then
759 write(unit=message(1),fmt='(A,I6,A,I6)') &
760 'Bogus #=', nlocal(bogus), ' > max_bogus_input=', max_bogus_input
761 call da_error(__FILE__,__LINE__,message(1:1))
762 end if
763
764 allocate (iv%bogus(nlocal(bogus))%h (1:nlevels))
765 allocate (iv%bogus(nlocal(bogus))%p (1:nlevels))
766 allocate (iv%bogus(nlocal(bogus))%u (1:nlevels))
767 allocate (iv%bogus(nlocal(bogus))%v (1:nlevels))
768 allocate (iv%bogus(nlocal(bogus))%t (1:nlevels))
769 allocate (iv%bogus(nlocal(bogus))%q (1:nlevels))
770
771 do i = 1, nlevels
772 iv%bogus(nlocal(bogus))%h(i) = platform%each(i)%height
773 iv%bogus(nlocal(bogus))%p(i) = platform%each(i)%p%inv
774 iv%bogus(nlocal(bogus))%u(i) = platform%each(i)%u
775 iv%bogus(nlocal(bogus))%v(i) = platform%each(i)%v
776 iv%bogus(nlocal(bogus))%t(i) = platform%each(i)%t
777 iv%bogus(nlocal(bogus))%q(i) = platform%each(i)%q
778 end do
779
780 iv%bogus(nlocal(bogus))%slp = platform%loc%slp
781
782 case default
783 write(unit=message(1), fmt='(a, i12)') &
784 'unsaved obs found with type: ',kx
785 call da_warning(__FILE__,__LINE__,message(1:1))
786 end select
787 end select
788
789 if (global .and. n < 2) then
790 if (test_wrfvar) exit dup_loop
791 if (platform%loc % i >= ide) then
792 platform%loc%i = platform%loc % i - ide
793 else if (platform%loc % i < ids) then
794 platform%loc%i = platform%loc % i + ide
795 end if
796 platform%loc%proc_domain = .not. platform%loc%proc_domain
797 end if
798 end do dup_loop
799 end do reports
800
801 call closbf(iunit)
802 close(iunit)
803 call da_free_unit(iunit)
804
805 if (trace_use) call da_trace_exit("da_read_obs_bufr")
806 #else
807 call da_error(__FILE__,__LINE__,(/"must compile with BUFR library"/))
808 #endif
809
810 end subroutine da_read_obs_bufr
811
812