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, old_nlevels,nlevels, nrecs,miscd, t29, jx, ix
29 integer :: cat,zqm,pqm,qqm,tqm,wqm,pwq,pmq,rqm
30 integer :: iunit, fm , obs_index
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
176 ! WHY?
177 ! no need for balloon drift
178 ! select case(kx)
179 ! case (120, 132, 220:221, 232) ; ! read drift
180 ! call ufbint(iunit,drf,9,255,iret,drift)
181 ! do k=1,nlevels
182 ! if (drf(1,k) .ge. 360.0)drf(1,k)=drf(1,k)-360.
183 ! if (drf(1,k) .lt. 0.0)drf(1,k)=drf(1,k)+360.
184 ! if (abs(drf(2,k)) .gt. 1000.0 .or. abs(drf(1,k)) .gt. 1000.0) then
185 ! drf(2,k)=hdr(3)
186 ! drf(1,k)=hdr(2)
187 ! end if
188 ! if (abs(drf(3,k)) .gt. 3.0)then
189 ! drf(3,k)=time
190 ! end if
191 ! end do
192 ! case default ; ! Do nothing
193 ! end select
194
195 platform % info % levels =nlevels
196 platform % info % elv =hdr(6)
197
198 platform % loc % slp %inv = missing_r
199 platform % loc % slp %qc = missing_data
200 platform % loc % slp %error=200.0
201 platform % loc % pw %inv = missing_r
202 platform % loc % pw %qc = missing_data
203 platform % loc % pw %error=0.2
204
205 do i=1,max_ob_levels
206 platform % each (i) % height = missing_r
207 platform % each (i) % height_qc = missing_data
208
209 platform % each (i) % zk = missing_r
210
211 platform % each (i) % u % inv = missing_r
212 platform % each (i) % u % qc = missing_data
213 platform % each (i) % u % error =10.0
214
215 platform % each (i) % v = platform % each (i) % u
216
217 platform % each (i) % t % inv = missing_r
218 platform % each (i) % t % qc = missing_data
219 platform % each (i) % t % error =5.0
220
221 platform % each (i) % p % inv = missing_r
222 platform % each (i) % p % qc = missing_data
223 platform % each (i) % p % error =200.0
224
225 platform % each (i) % q % inv = missing_r
226 platform % each (i) % q % qc = missing_data
227 platform % each (i) % q % error =0.1
228 end do
229
230 do k = 1, platform % info % levels
231 pob=obs(1,k)
232 ! scale q and compute t from tv and set units to Kelvin, if they aren't missing
233
234 if (obs(2,k) > 0.0 .and. obs(2,k) < 1e9) then
235 obs(2,k) = obs(2,k)*1e-6
236 if (obs(3,k) > -200.0 .and. obs(3,k) < 300.0) then
237 obs(3,k) = (obs(3,k)+t_kelvin) / (1.0 + 0.608 * obs(2,k))
238 end if
239 else
240 if (obs(3,k) > -200.0 .and. obs(3,k) < 300.0) then
241 obs(3,k) = obs(3,k) + t_kelvin
242 end if
243 end if
244 qob=obs(2,k)
245 tob=obs(3,k)
246 zob=obs(4,k)
247 rob=obs(8,k)
248 pqm = qflag
249 if (abs(qms(1,k)) < qflag) pqm=nint(qms(1,k))
250 qqm = qflag
251 if (abs(qms(2,k)) < qflag) qqm=nint(qms(2,k))
252 tqm = qflag
253 if (abs(qms(3,k)) < qflag) tqm=nint(qms(3,k))
254 zqm = qflag
255 if (abs(qms(4,k)) < qflag) zqm=nint(qms(4,k))
256 wqm = qflag
257 if (abs(qms(5,k)) < qflag) wqm=nint(qms(5,k))
258 pwq = qflag
259 if (abs(qms(7,k)) < qflag) pwq=nint(qms(7,k))
260 rqm = qflag
261 if (abs(qms(9,k)) < qflag) rqm=nint(qms(9,k))
262 pmq = qflag
263 if (abs(pmo(2,1)) < qflag) pmq=nint(pmo(2,1))
264 cat=nint(obs(9,k))
265
266 if (pmq < qflag) then
267 platform % loc % slp % inv =pmo(1,1)*100.0
268 platform % loc % slp % qc =pmq
269 platform % loc % slp % error = 200.0
270 end if
271
272 if (pwq < qflag) then
273 platform % loc % pw % inv = obs(7,k) * 0.1 ! convert to cm
274 platform % loc % pw % qc = pwq
275 platform % loc % pw % error = 0.2 ! hardwired to 0.2 cm
276 end if
277
278 poe=100.0
279
280 if (tqm<qflag) then
281 toe=oes(3,k)
282 platform % each (k) % t % inv =obs(3,k)
283 platform % each (k) % t % qc =tqm
284 platform % each (k) % t % error =toe
285 end if
286
287 if (wqm < qflag) then
288 woe=oes(5,k)
289 ! WHY?
290 ! if (wqm==0) woe=woe*0.9
291 ! if (wqm==3) woe=woe*1.2
292 platform % each (k) % u % inv =obs(5,k)
293 platform % each (k) % v % inv =obs(6,k)
294 platform % each (k) % u % qc =wqm
295 platform % each (k) % u % error =woe
296 platform % each (k) % v % qc =wqm
297 platform % each (k) % v % error =woe
298
299 ! Convert earth wind to model wind.
300 ! if( kx /= 283 .and. kx /= 152) & ! skip SSMI speed or TPW
301 call da_earth_2_model_wind(obs(5,k), obs(6,k), &
302 platform % each (k) % u % inv, &
303 platform % each (k) % v % inv, &
304 platform%info%lon)
305 if (platform % each (k) % u % inv == 0.0 .and. platform % each (k) % v % inv == 0.0) then
306 platform % each (k) % u % inv = missing_r
307 platform % each (k) % v % inv = missing_r
308 platform % each (k) % u % qc = missing_data
309 platform % each (k) % v % qc = missing_data
310 end if
311 end if
312
313 if (qqm<qflag .and. obs(2,k)<1.0e9) then
314 qoe=oes(2,k)
315 platform % each (k) % q % inv =obs(2,k)
316 platform % each (k) % q % qc =qqm
317 platform % each (k) % q % error =qoe
318 end if
319
320 if (zqm<qflag)then
321 platform % each (k) % height =zob
322 platform % each (k) % height_qc =zqm
323 end if
324
325 if (rqm<qflag)then
326 roe = oes(8,k)
327 platform % each (k) % u % inv =rob
328 platform % each (k) % u % qc =rqm
329 platform % each (k) % u % error =roe
330 end if
331
332 if (pqm<qflag)then
333 platform % each (k) % p % inv =obs(1,k)*100.0
334 platform % each (k) % p % qc =pqm
335 platform % each (k) % p % error =poe
336 end if
337 end do
338
339 nlevels = platform%info%levels
340
341 if (nlevels > max_ob_levels) then
342 nlevels = max_ob_levels
343
344 write(unit=stderr, fmt='(/a/)') 'Warning: Too many levels.'
345
346 write(unit=stderr, fmt='(/2a/2a/2x,a,2f8.2,2(2x,a,f9.2)/2(2x,a,i4)/)') &
347 'Subset: ', platform%info%name(1:8), &
348 'Platfrom: ', trim(platform%info%platform), &
349 'Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
350 'elv: ', platform%info%elv, &
351 'pstar: ', platform%info%pstar, &
352 'level: ', platform%info%levels, &
353 'kx: ', kx
354 else if ((nlevels < 1) .and. ((kx /= 164) .or. (kx /= 174))) then
355 write(unit=stderr, fmt='(/a/)') &
356 'Warning: Too few levels.'
357
358 write(unit=stderr, fmt='(/2a/2a/2x,a,2f8.2,2(2x,a,f9.2)/2(2x,a,i4)/)') &
359 'Subset: ', platform%info%name(1:8), &
360 'Platfrom: ', trim(platform%info%platform), &
361 'Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
362 'elv: ', platform%info%elv, &
363 'pstar: ', platform%info%pstar, &
364 'level: ', platform%info%levels, &
365 'kx: ', kx
366
367 cycle reports
368 end if
369
370 !---------------------------------------------------------------------------
371 ! This is basically converting rh to q i
372 ! Method :
373 ! if rh, temp and pr all available computes Qs otherwise sets Qs= missing
374 ! if rh > 100 sets q = qs otherwise q = rh*Qs/100.0
375 ! Note: Currently da_obs_proc_station is active only for ob_format_ascii
376 ! call da_obs_proc_station(platform)
377 !---------------------------------------------------------------------------
378
379 ! Loop over duplicating obs for global
380 ndup = 1
381 if (global .and. &
382 (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
383
384 ! It is possible that logic for counting obs is incorrect for the
385 ! global case with >1 MPI tasks due to obs duplication, halo, etc.
386 ! TBH: 20050913
387 if (.not.outside) then
388 if (print_detail_obs .and. ndup > 1) then
389 write(unit=stdout, fmt = '(A12,1X,A19,1X,A40,1X,I6,3(F12.3,11X),6X,A5)') &
390 platform%info%platform, &
391 platform%info%date_char, &
392 platform%info%name, &
393 platform%info%levels, &
394 platform%info%lat, &
395 platform%info%lon, &
396 platform%info%elv, &
397 platform%info%id
398
399 write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
400 ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
401 platform%loc%i, platform%loc%j, &
402 platform%loc%dx, platform%loc%dxm, &
403 platform%loc%dy, platform%loc%dym
404 end if
405 end if
406 dup_loop: do n = 1, ndup
407 select case(t29)
408 case (11, 12, 13, 22, 23, 31)
409 select case (kx)
410 case (120, 122, 132, 220, 222, 232) ; ! Sound
411 if (.not.use_soundobs) cycle reports
412 if( n==1 ) ntotal(sound) = ntotal(sound) + 1
413 if( n==1 ) ntotal(sonde_sfc) = ntotal(sonde_sfc) + 1
414 if (outside) then
415 cycle reports
416 end if
417
418
419 platform % info % platform ='FM-35 TEMP'
420 nlocal(sound) = nlocal(sound) + 1
421 nlocal(sonde_sfc) = nlocal(sonde_sfc) + 1
422 if (nlocal(sound) > iv%info(sound)%nlocal) cycle reports
423 fm = 35
424 old_nlevels = nlevels
425 ! Track serial obs index for reassembly of obs during
426 ! bit-for-bit tests with different numbers of MPI tasks.
427
428 ! Search to see if we have surface obs.
429
430 surface_level = 0
431
432 do i = 1, nlevels
433 ! if (elevation and height are the same, it is surface)
434 if (abs(platform%info%elv - &
435 platform%each(i)%height) < 0.1) then
436 surface_level = i
437
438 ! Save surface pressure.
439 iv%sonde_sfc(nlocal(sonde_sfc))%h = platform%each(i)%height
440 iv%sonde_sfc(nlocal(sonde_sfc))%u = platform%each(i)%u
441 iv%sonde_sfc(nlocal(sonde_sfc))%v = platform%each(i)%v
442 iv%sonde_sfc(nlocal(sonde_sfc))%t = platform%each(i)%t
443 iv%sonde_sfc(nlocal(sonde_sfc))%q = platform%each(i)%q
444 iv%sonde_sfc(nlocal(sonde_sfc))%p = platform%each(i)%p
445 exit
446 end if
447 end do
448
449 ! processing the sound_sfc data:
450
451 if (surface_level > 0) then
452 nlevels = nlevels - 1
453 else
454 iv%sonde_sfc(nlocal(sonde_sfc))%h = missing_r
455 iv%sonde_sfc(nlocal(sonde_sfc))%u%inv = missing_r
456 iv%sonde_sfc(nlocal(sonde_sfc))%u%qc = missing
457 iv%sonde_sfc(nlocal(sonde_sfc))%u%error = abs(missing_r)
458 iv%sonde_sfc(nlocal(sonde_sfc))%v = iv%sonde_sfc(nlocal(sonde_sfc))%u
459 iv%sonde_sfc(nlocal(sonde_sfc))%t = iv%sonde_sfc(nlocal(sonde_sfc))%u
460 iv%sonde_sfc(nlocal(sonde_sfc))%p = iv%sonde_sfc(nlocal(sonde_sfc))%u
461 iv%sonde_sfc(nlocal(sonde_sfc))%q = iv%sonde_sfc(nlocal(sonde_sfc))%u
462 end if
463
464 if (nlevels < 1) cycle reports
465
466 allocate (iv%sound(nlocal(sound))%h (1:nlevels))
467 allocate (iv%sound(nlocal(sound))%p (1:nlevels))
468 allocate (iv%sound(nlocal(sound))%u (1:nlevels))
469 allocate (iv%sound(nlocal(sound))%v (1:nlevels))
470 allocate (iv%sound(nlocal(sound))%t (1:nlevels))
471 allocate (iv%sound(nlocal(sound))%q (1:nlevels))
472
473 j = 0
474 do i = 1, old_nlevels
475 if (i == surface_level) cycle
476
477 j=j+1
478
479 iv%sound(nlocal(sound))%h(j) = platform%each(i)%height
480 iv%sound(nlocal(sound))%p(j) = platform%each(i)%p%inv
481 iv%sound(nlocal(sound))%u(j) = platform%each(i)%u
482 iv%sound(nlocal(sound))%v(j) = platform%each(i)%v
483 iv%sound(nlocal(sound))%t(j) = platform%each(i)%t
484 iv%sound(nlocal(sound))%q(j) = platform%each(i)%q
485 end do
486
487 case (221) ; ! Pilot
488 if (.not.use_pilotobs) cycle reports
489 if( n==1 ) ntotal(pilot) = ntotal(pilot) + 1
490 if (outside) then
491 cycle reports
492 end if
493
494 nlocal(pilot) = nlocal(pilot) + 1
495 if (nlocal(pilot) > iv%info(pilot)%nlocal) cycle reports
496 fm = 32
497 ! Track serial obs index for reassembly of obs during
498 ! bit-for-bit tests with different numbers of MPI tasks.
499
500 allocate (iv%pilot(nlocal(pilot))%p (1:nlevels))
501 allocate (iv%pilot(nlocal(pilot))%u (1:nlevels))
502 allocate (iv%pilot(nlocal(pilot))%v (1:nlevels))
503
504 do i = 1, nlevels
505 iv%pilot(nlocal(pilot))%p(i) = platform%each(i)%p%inv
506 iv%pilot(nlocal(pilot))%u(i) = platform%each(i)%u
507 iv%pilot(nlocal(pilot))%v(i) = platform%each(i)%v
508 end do
509 end select
510
511 case (41)
512 ! case (130:131, 133, 230:231, 233) ; ! Airep
513 if (.not.use_airepobs) cycle reports
514 if( n==1 ) ntotal(airep) = ntotal(airep) + 1
515 if (outside) then
516 cycle reports
517 end if
518
519
520 platform % info % platform ='FM-97 AIREP'
521 nlocal(airep) = nlocal(airep) + 1
522 if (nlocal(airep) > iv%info(airep)%nlocal) cycle reports
523 fm = 42
524 ! Track serial obs index for reassembly of obs during
525 ! bit-for-bit tests with different numbers of MPI tasks.
526
527 allocate (iv % airep (nlocal(airep)) % h (1:nlevels))
528 allocate (iv % airep (nlocal(airep)) % p (1:nlevels))
529 allocate (iv % airep (nlocal(airep)) % u (1:nlevels))
530 allocate (iv % airep (nlocal(airep)) % v (1:nlevels))
531 allocate (iv % airep (nlocal(airep)) % t (1:nlevels))
532
533 do i = 1, nlevels
534 iv % airep (nlocal(airep)) % h(i) = platform % each(i) % height
535 iv % airep (nlocal(airep)) % p(i) = platform % each(i) % p % inv
536 iv % airep (nlocal(airep)) % u(i) = platform % each(i) % u
537 iv % airep (nlocal(airep)) % v(i) = platform % each(i) % v
538 iv % airep (nlocal(airep)) % t(i) = platform % each(i) % t
539 end do
540
541 ! case (180, 182, 280, 282) ; ! Ships and buoys
542
543 case (522, 523); ! Ships
544 if (.not.use_shipsobs) cycle reports
545 if( n==1 ) ntotal(ships) = ntotal(ships) + 1
546 if (outside) then
547 cycle reports
548 end if
549
550
551 platform % info % platform ='FM-13 SHIP '
552 nlocal(ships) = nlocal(ships) + 1
553 if (nlocal(ships) > iv%info(ships)%nlocal) cycle reports
554 fm = 13
555 ! Track serial obs index for reassembly of obs during
556 ! bit-for-bit tests with different numbers of MPI tasks.
557
558 iv % ships (nlocal(ships)) % h = platform % each(1) % height
559 iv % ships (nlocal(ships)) % u = platform % each(1) % u
560 iv % ships (nlocal(ships)) % v = platform % each(1) % v
561 iv % ships (nlocal(ships)) % t = platform % each(1) % t
562 ! WHY?
563 ! iv % ships (nlocal(ships)) % p = platform % loc % slp
564 iv % ships (nlocal(ships)) % p = platform % each(1) % p
565 iv % ships (nlocal(ships)) % q = platform % each(1) % q
566
567 case (531, 561, 562) ; ! Buoy
568 if (.not.use_buoyobs) cycle reports
569 if( n==1 ) ntotal(buoy) = ntotal(buoy) + 1
570 if (outside) then
571 cycle reports
572 end if
573
574 platform % info % platform ='FM-18 BUOY '
575 nlocal(buoy) = nlocal(buoy) + 1
576 if (nlocal(buoy) > iv%info(buoy)%nlocal) cycle reports
577 fm = 18
578
579 iv%buoy(nlocal(buoy))%h = platform%each(1)%height
580 iv%buoy(nlocal(buoy))%u = platform%each(1)%u
581 iv%buoy(nlocal(buoy))%v = platform%each(1)%v
582 iv%buoy(nlocal(buoy))%t = platform%each(1)%t
583 iv%buoy(nlocal(buoy))%p = platform%each(1)%p
584 iv%buoy(nlocal(buoy))%q = platform%each(1)%q
585
586 case (511)
587 if (.not.use_synopobs) cycle reports
588 ! case (181, 281) ; ! Synop
589
590 if( n==1 ) ntotal(synop) = ntotal(synop) + 1
591 if (outside) then
592 cycle reports
593 end if
594
595
596 platform % info % platform ='FM-12 SYNOP'
597
598 nlocal(synop) = nlocal(synop) + 1
599 if (nlocal(synop) > iv%info(synop)%nlocal) cycle reports
600 fm = 12
601
602 iv % synop (nlocal(synop)) % h = platform % each(1) % height
603 iv % synop (nlocal(synop)) % u = platform % each(1) % u
604 iv % synop (nlocal(synop)) % v = platform % each(1) % v
605 iv % synop (nlocal(synop)) % t = platform % each(1) % t
606 ! WHY?
607 ! iv % synop (nlocal(synop)) % p = platform % loc % slp
608 iv % synop (nlocal(synop)) % p = platform % each(1) % p
609 iv % synop (nlocal(synop)) % q = platform % each(1) % q
610
611 if (iv % synop(nlocal(synop)) % h < platform % info % elv) then
612 iv % synop(nlocal(synop)) % h = platform % info % elv
613 end if
614
615 case (512)
616 if (.not.use_metarobs) cycle reports
617 ! case (187, 287) ; ! Metar
618 if( n==1 ) ntotal(metar) = ntotal(metar) + 1
619 if (outside) then
620 cycle reports
621 end if
622
623
624 platform % info % platform ='FM-15 METAR'
625 nlocal(metar) = nlocal(metar) + 1
626 if (nlocal(metar) > iv%info(metar)%nlocal) cycle reports
627 fm = 15
628
629 iv % metar (nlocal(metar)) % h = platform % each(1) % height
630 iv % metar (nlocal(metar)) % u = platform % each(1) % u
631 iv % metar (nlocal(metar)) % v = platform % each(1) % v
632 iv % metar (nlocal(metar)) % t = platform % each(1) % t
633 iv % metar (nlocal(metar)) % p = platform % each(1) % p
634 iv % metar (nlocal(metar)) % q = platform % each(1) % q
635
636 case (63)
637 if (.not.use_geoamvobs) cycle reports
638 ! case (242:246, 252:253, 255) ; ! Geo. CMVs
639 if( n==1 ) ntotal(geoamv) = ntotal(geoamv) + 1
640 if (outside) then
641 cycle reports
642 end if
643
644 platform % info % platform ='FM-88 SATOB'
645 nlocal(geoamv) = nlocal(geoamv) + 1
646 if (nlocal(geoamv) > iv%info(geoamv)%nlocal) cycle reports
647 fm = 88
648
649 allocate (iv%geoamv(nlocal(geoamv))%p (1:nlevels))
650 allocate (iv%geoamv(nlocal(geoamv))%u (1:nlevels))
651 allocate (iv%geoamv(nlocal(geoamv))%v (1:nlevels))
652
653 do i = 1, nlevels
654 iv % geoamv (nlocal(geoamv)) % p(i) = platform % each(i) % p % inv
655 iv % geoamv (nlocal(geoamv)) % u(i) = platform % each(i) % u
656 iv % geoamv (nlocal(geoamv)) % v(i) = platform % each(i) % v
657 iv % info(geoamv)%zk(i,nlocal(geoamv)) = platform % each(i) % zk
658 end do
659
660 case (582)
661 if (.not.use_qscatobs) cycle reports
662 if( n==1 ) ntotal(qscat) = ntotal(qscat) + 1
663 if (outside) then
664 cycle reports
665 end if
666
667 platform % info % platform ='FM-281 Quiks'
668 nlocal(qscat) = nlocal(qscat) + 1
669 if (nlocal(qscat) > iv%info(qscat)%nlocal) cycle reports
670 fm = 281
671
672 ! WHY?
673 ! iv%qscat(nlocal(qscat))%h = platform%each(1)%height
674 ! prepbufr uses pressure not height, so hardwire height to
675 ! 0 (sea-level)
676 iv%qscat(nlocal(qscat))%h = 0.0
677 iv%qscat(nlocal(qscat))%u = platform%each(1)%u
678 iv%qscat(nlocal(qscat))%v = platform%each(1)%v
679 iv%qscat(nlocal(qscat))%u%error = max(platform%each(1)%u%error,1.0)
680 iv%qscat(nlocal(qscat))%v%error = max(platform%each(1)%v%error,1.0)
681
682 case (583) ! GPS PW
683 if (.not.use_gpspwobs) cycle reports
684 if( n==1 ) ntotal(gpspw) = ntotal(gpspw) + 1
685 if (outside) then
686 cycle reports
687 end if
688
689 platform % info % platform ='FM-111 GPSPW'
690 nlocal(gpspw) = nlocal(gpspw) + 1
691 if (nlocal(gpspw) > iv%info(gpspw)%nlocal) cycle reports
692 fm = 111
693
694 iv%gpspw(nlocal(gpspw))%tpw = platform%loc%pw
695
696 case (584) ! GPS REF
697 if (.not.use_gpsrefobs) cycle reports
698 if( n==1 ) ntotal(gpsref) = ntotal(gpsref) + 1
699 if (outside) then
700 cycle reports
701 end if
702
703 platform % info % platform ='FM-116 GPSRF'
704 nlocal(gpsref) = nlocal(gpsref) + 1
705 if (nlocal(gpsref) > iv%info(gpsref)%nlocal) cycle reports
706 fm = 116
707
708 allocate (iv%gpsref(nlocal(gpsref))%h (1:nlevels))
709 allocate (iv%gpsref(nlocal(gpsref))%ref(1:nlevels))
710 allocate (iv%gpsref(nlocal(gpsref))% p(1:nlevels))
711 allocate (iv%gpsref(nlocal(gpsref))% t(1:nlevels))
712 allocate (iv%gpsref(nlocal(gpsref))% q(1:nlevels))
713
714 do i = 1, nlevels
715 iv%gpsref(nlocal(gpsref))%h(i) = platform%each(i)%height
716 ! u is used to store ref in this routine (ascii version uses td)
717 iv%gpsref(nlocal(gpsref))%ref(i) = platform%each(i)%u
718 ! Keep the retrieved p and t (and q) as "field_type":
719 iv%gpsref(nlocal(gpsref))%p(i) = platform%each(i)%p
720 iv%gpsref(nlocal(gpsref))%t(i) = platform%each(i)%t
721 iv%gpsref(nlocal(gpsref))%q(i) = platform%each(i)%q
722 end do
723
724 case (71, 72)
725 ! case (223, 224 ) ; ! Profiler & VADWND - NEXRAD winds
726 if (.not.use_profilerobs) cycle reports
727 if( n==1 ) ntotal(profiler) = ntotal(profiler)
728 if (outside) then
729 cycle reports
730 end if
731
732 platform % info % platform ='FM-132 PRFLR'
733 nlocal(profiler) = nlocal(profiler) + 1
734 if (nlocal(profiler) > iv%info(profiler)%nlocal) cycle reports
735 fm = 132
736 ! Track serial obs index for reassembly of obs during bit-for-bit
737 ! tests with different numbers of MPI tasks.
738
739 allocate (iv%profiler(nlocal(profiler))%p (1:nlevels))
740 allocate (iv%profiler(nlocal(profiler))%u (1:nlevels))
741 allocate (iv%profiler(nlocal(profiler))%v (1:nlevels))
742
743 do i = 1, nlevels
744 iv%profiler(nlocal(profiler))%p(i) = platform%each(i)%p%inv
745 iv%profiler(nlocal(profiler))%u(i) = platform%each(i)%u
746 iv%profiler(nlocal(profiler))%v(i) = platform%each(i)%v
747 end do
748
749 case (571, 65) ! SSM/I wind speed & TPW
750 if (.not. use_ssmiretrievalobs) cycle reports
751 if( n==1 ) ntotal(ssmi_rv) = ntotal(ssmi_rv) + 1
752 if (outside) then
753 cycle reports
754 end if
755
756 nlocal(ssmi_rv) = nlocal(ssmi_rv) + 1
757 if (nlocal(ssmi_rv) > iv%info(ssmi_rv)%nlocal) cycle reports
758 platform % info % platform ='FM-125 SSMI'
759 fm = 125
760 select case (kx)
761 case ( 283)
762 !------------------
763 ! write(unit=stdout,fmt=*) 'report, ',c8id(1:5),' typ = ',kx, &
764 ! ' nlevels = ',nlevels
765 ! do ix = 1,nlevels
766 ! write(unit=stdout,fmt='(a5,2x,6(f10.4,2x))')c8id(1:5), (hdr(jx),jx=2,7)
767 ! write(unit=stdout,fmt='(9(f10.4,2x))') (obs(jx,ix),jx=1,9)
768 ! write(unit=stdout,fmt='(9(f10.4,2x))') (qms(jx,ix),jx=1,9)
769 ! write(unit=stdout,fmt=*) 'PMO = ',pmo(1,1)
770 ! end do
771 !------------------
772 do i = 1, nlevels
773 iv%ssmi_rv(nlocal(ssmi_rv))%speed = platform%each(i)%u
774 iv%ssmi_rv(nlocal(ssmi_rv))%speed%inv = sqrt ( &
775 platform%each(i)%u%inv * platform%each(i)%u%inv + &
776 platform%each(i)%v%inv * platform%each(i)%v%inv )
777 iv%ssmi_rv(nlocal(ssmi_rv))%tpw = platform%loc%pw
778 end do
779 case ( 152 )
780 !------------------
781 ! write(unit=stdout,fmt=*) 'report, ',c8id(1:5),' typ = ',kx, &
782 ! ' nlevels = ',nlevels
783 ! do ix = 1,nlevels
784 ! write(unit=stdout,fmt='(a5,2x,6(f10.4,2x))')c8id(1:5), (hdr(jx),jx=2,7)
785 ! write(unit=stdout,fmt='(9(f10.4,2x))') (obs(jx,ix),jx=1,9)
786 ! write(unit=stdout,fmt='(9(f10.4,2x))') (qms(jx,ix),jx=1,9)
787 ! write(unit=stdout,fmt=*) 'PMO = ',pmo(1,1)
788 ! end do
789 !------------------
790 do i = 1, nlevels
791 iv%ssmi_rv(nlocal(ssmi_rv))%speed = platform%each(i)%u
792 iv%ssmi_rv(nlocal(ssmi_rv))%tpw = platform%loc%pw
793 end do
794 case default
795 exit dup_loop
796 end select
797
798 case default
799 select case (kx)
800 case (111 , 210) ; ! Tropical Cyclone Bogus
801 ! Note Tropical cyclone Bougus is given type 135 in Obs-ascii
802 if( n==1 ) ntotal(bogus) = ntotal(bogus) + 1
803 if (outside) then
804 cycle reports
805 end if
806
807 if (.not.use_bogusobs) cycle reports
808 platform % info % platform ='FM-135 TCBOG'
809 nlocal(bogus) = nlocal(bogus) + 1
810 fm = 135
811 ! Track serial obs index for reassembly of obs during bit-for-bit
812 ! tests with different numbers of MPI tasks.
813
814 if (nlocal(bogus) > max_bogus_input) then
815 write(unit=message(1),fmt='(A,I6,A,I6)') &
816 'Bogus #=', nlocal(bogus), ' > max_bogus_input=', max_bogus_input
817 call da_error(__FILE__,__LINE__,message(1:1))
818 end if
819
820 allocate (iv%bogus(nlocal(bogus))%h (1:nlevels))
821 allocate (iv%bogus(nlocal(bogus))%p (1:nlevels))
822 allocate (iv%bogus(nlocal(bogus))%u (1:nlevels))
823 allocate (iv%bogus(nlocal(bogus))%v (1:nlevels))
824 allocate (iv%bogus(nlocal(bogus))%t (1:nlevels))
825 allocate (iv%bogus(nlocal(bogus))%q (1:nlevels))
826
827 do i = 1, nlevels
828 iv%bogus(nlocal(bogus))%h(i) = platform%each(i)%height
829 iv%bogus(nlocal(bogus))%p(i) = platform%each(i)%p%inv
830 iv%bogus(nlocal(bogus))%u(i) = platform%each(i)%u
831 iv%bogus(nlocal(bogus))%v(i) = platform%each(i)%v
832 iv%bogus(nlocal(bogus))%t(i) = platform%each(i)%t
833 iv%bogus(nlocal(bogus))%q(i) = platform%each(i)%q
834 end do
835
836 iv%bogus(nlocal(bogus))%slp = platform%loc%slp
837
838 case default
839 exit dup_loop
840 end select
841 end select
842
843 ! rizvi's add start
844 obs_index = fm_index(fm)
845 iv%info(obs_index)%name(nlocal(obs_index)) = platform%info%name
846 iv%info(obs_index)%platform(nlocal(obs_index)) = platform%info%platform
847 iv%info(obs_index)%id(nlocal(obs_index)) = platform%info%id
848 iv%info(obs_index)%date_char(nlocal(obs_index)) = platform%info%date_char
849 ! nlevels adjusted for some obs types so use that
850 iv%info(obs_index)%levels(nlocal(obs_index)) = nlevels
851 iv%info(obs_index)%lat(:,nlocal(obs_index)) = platform%info%lat
852 iv%info(obs_index)%lon(:,nlocal(obs_index)) = platform%info%lon
853 iv%info(obs_index)%elv(nlocal(obs_index)) = platform%info%elv
854 iv%info(obs_index)%pstar(nlocal(obs_index)) = platform%info%pstar
855
856 iv%info(obs_index)%slp(nlocal(obs_index)) = platform%loc%slp
857 iv%info(obs_index)%pw(nlocal(obs_index)) = platform%loc%pw
858 iv%info(obs_index)%x(:,nlocal(obs_index)) = platform%loc%x
859 iv%info(obs_index)%y(:,nlocal(obs_index)) = platform%loc%y
860 iv%info(obs_index)%i(:,nlocal(obs_index)) = platform%loc%i
861 iv%info(obs_index)%j(:,nlocal(obs_index)) = platform%loc%j
862 iv%info(obs_index)%dx(:,nlocal(obs_index)) = platform%loc%dx
863 iv%info(obs_index)%dxm(:,nlocal(obs_index)) = platform%loc%dxm
864 iv%info(obs_index)%dy(:,nlocal(obs_index)) = platform%loc%dy
865 iv%info(obs_index)%dym(:,nlocal(obs_index)) = platform%loc%dym
866 iv%info(obs_index)%proc_domain(:,nlocal(obs_index)) = platform%loc%proc_domain
867
868 iv%info(obs_index)%obs_global_index(nlocal(obs_index)) = ntotal(obs_index)
869 ! special case for sonde_sfc, duplicate sound info
870 if (obs_index == sound) then
871
872 iv%info(sonde_sfc)%name(nlocal(sonde_sfc)) = platform%info%name
873 iv%info(sonde_sfc)%platform(nlocal(sonde_sfc)) = platform%info%platform
874 iv%info(sonde_sfc)%id(nlocal(sonde_sfc)) = platform%info%id
875 iv%info(sonde_sfc)%date_char(nlocal(sonde_sfc)) = platform%info%date_char
876 iv%info(sonde_sfc)%levels(nlocal(sonde_sfc)) = 1
877 iv%info(sonde_sfc)%lat(:,nlocal(sonde_sfc)) = platform%info%lat
878 iv%info(sonde_sfc)%lon(:,nlocal(sonde_sfc)) = platform%info%lon
879 iv%info(sonde_sfc)%elv(nlocal(sonde_sfc)) = platform%info%elv
880 iv%info(sonde_sfc)%pstar(nlocal(sonde_sfc)) = platform%info%pstar
881
882 iv%info(sonde_sfc)%slp(nlocal(sonde_sfc)) = platform%loc%slp
883 iv%info(sonde_sfc)%pw(nlocal(sonde_sfc)) = platform%loc%pw
884 iv%info(sonde_sfc)%x(:,nlocal(sonde_sfc)) = platform%loc%x
885 iv%info(sonde_sfc)%y(:,nlocal(sonde_sfc)) = platform%loc%y
886 iv%info(sonde_sfc)%i(:,nlocal(sonde_sfc)) = platform%loc%i
887 iv%info(sonde_sfc)%j(:,nlocal(sonde_sfc)) = platform%loc%j
888 iv%info(sonde_sfc)%dx(:,nlocal(sonde_sfc)) = platform%loc%dx
889 iv%info(sonde_sfc)%dxm(:,nlocal(sonde_sfc)) = platform%loc%dxm
890 iv%info(sonde_sfc)%dy(:,nlocal(sonde_sfc)) = platform%loc%dy
891 iv%info(sonde_sfc)%dym(:,nlocal(sonde_sfc)) = platform%loc%dym
892 iv%info(sonde_sfc)%proc_domain(:,nlocal(sonde_sfc)) = platform%loc%proc_domain
893
894 iv%info(sonde_sfc)%obs_global_index(nlocal(sonde_sfc)) = ntotal(sonde_sfc)
895 end if
896 !rizvi add ends
897
898 if (global .and. n < 2) then
899 if (test_transforms) exit dup_loop
900 if (platform%loc % i >= ide) then
901 platform%loc%i = platform%loc % i - ide
902 else if (platform%loc % i < ids) then
903 platform%loc%i = platform%loc % i + ide
904 end if
905 platform%loc%proc_domain = .not. platform%loc%proc_domain
906 end if
907 end do dup_loop
908 end do reports
909
910 call closbf(iunit)
911 close(iunit)
912 call da_free_unit(iunit)
913
914 if (trace_use) call da_trace_exit("da_read_obs_bufr")
915 #else
916 call da_error(__FILE__,__LINE__,(/"must compile with BUFR library"/))
917 #endif
918
919 end subroutine da_read_obs_bufr