da_read_bufr_obs.inc
References to this file elsewhere.
1 subroutine da_read_bufr_obs (ob, xp, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Read 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), optional :: filename
13
14 #ifdef BUFR
15
16 type (multi_level_type) :: platform
17 logical :: outside, outside_all
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
37 character(len=40) :: obstr,drift,hdstr,qmstr,oestr
38 character(len=8) :: subset, c8id
39 real,dimension(7) :: hdr
40 real,dimension(2,1) :: pmo
41 ! real,dimension(9,255) :: drf
42 real,dimension(9,255) :: obs,qms,oes
43 real :: time,woe,qob,toe,qoe,poe,pob,tob,zob,rob,roe
44
45 integer :: iost, ndup, n, i, j, k, surface_level, report
46 integer :: iret, idate, kx, nlevels, nrecs,miscd, t29, jx, ix
47 integer :: cat,zqm,pqm,qqm,tqm,wqm,pwq,pmq,rqm
48 integer :: iunit
49
50 integer, parameter :: qflag=4 ! Flag for retaining data
51 equivalence (hdr(1), c8id)
52
53 total_obs = ob%total_obs
54 num_sound = ob%num_sound
55 num_synop = ob%num_synop
56 num_pilot = ob%num_pilot
57 num_geoamv = ob%num_geoamv
58 num_polaramv = ob%num_polaramv
59 num_satem = ob%num_satem
60 num_airep = ob%num_airep
61 num_metar = ob%num_metar
62 num_ships = ob%num_ships
63 num_gpspw = ob%num_gpspw
64 num_gpsref = ob%num_gpsref
65 num_ssmi_retrieval = ob%num_ssmi_retrieval
66 num_ssmi_tb = ob%num_ssmi_tb
67 num_ssmt1 = ob%num_ssmt1
68 num_ssmt2 = ob%num_ssmt2
69 num_qscat = ob%num_qscat
70 num_profiler = ob%num_profiler
71 num_buoy = ob%num_buoy
72 num_bogus= ob%num_bogus
73 num_pseudo = 0
74
75 ! open file
76 ! ---------
77 call da_get_unit(iunit)
78 if (present(filename)) then
79 open(unit = iunit, FILE = trim(filename), &
80 iostat = iost, form = 'unformatted', STATUS = 'OLD')
81 if (iost /= 0) then
82 write(unit=stdout, fmt='(/A,I3,3(2X,A)/)') &
83 'ERROR in OBS inPUT FILE unit ',iunit, &
84 'OBS FILENAME:', trim(filename), &
85 'FOR GTS OBSERVATIONS CANNOT BE FOUND OR CANNOT BE openED'
86 return
87 end if
88 end if
89
90 hdstr='SID XOB YOB DHR TYP ELV T29 '
91 obstr='POB QOB TOB ZOB UOB VOB PWO ROB CAT '
92 drift='XDR YDR HRDR '
93 qmstr='PQM QQM TQM ZQM WQM NUL PWQ PMQ RQM '
94 oestr='POE QOE TOE NUL WOE NUL PWE ROE '
95
96 nrecs = 0
97 miscd = 0
98
99 ! read bufr observation file
100 ! call openbf(iunit,'in',iunit)
101 call readmg(iunit,subset,idate,iret)
102
103 if (iret/=0) then
104 write(unit=stdout, fmt='(a, i4)') &
105 'return code from readmg:', iret, &
106 'Reach the end of obs unit: ', iunit
107
108 call da_error(__FILE__,__LINE__, &
109 (/"No BUFR observations"/))
110 end if
111
112 report = 0 ! report number in file
113 reports: do
114 report = report+1
115
116 call readsb(iunit,iret)
117
118 if (iret/=0) then
119 call readmg(iunit,subset,idate,iret)
120
121 if (iret/=0) then
122 write(unit=stdout, fmt='(a, i4)') &
123 'return code from readmg:', iret, &
124 'Reach the end of prepbufr obs unit: ', iunit
125
126 exit reports
127 end if
128
129 cycle reports
130 end if
131
132 nrecs=nrecs+1
133
134 call ufbint(iunit,hdr,7,1,iret,hdstr)
135
136 platform % info % name(1:8) = subset
137
138 platform % info % id = c8id(1:5)
139 if (hdr(2) > 180.0) hdr(2)=hdr(2)-360.0
140
141 ! Put a check on Lat
142
143 hdr(3) = MAX(hdr(3), -89.95)
144 hdr(3) = Min(hdr(3), 89.95)
145 platform%info%lon = hdr(2)
146 platform%info%lat = hdr(3)
147
148 ! Put a check on Lat
149
150 if (platform%info%lon == 180.0 ) platform%info%lon =-180.000
151 ! Fix funny wind direction at Poles
152 if (platform%info%lat <= -89.95 .or. platform%info%lat >= 89.95) then
153 platform%info%lon = 0.0
154 end if
155
156 ! Restrict to a range of reports, useful for debugging
157
158 if (report < report_start) then
159 cycle
160 end if
161
162 if (report > report_end) then
163 exit
164 end if
165
166 call da_ll_to_xy (platform%info, platform%loc,xp, outside, outside_all)
167
168 if (outside_all) cycle reports
169
170 if (print_detail_obs) then
171 ! Simplistic approach, could be improved to get it all done on PE 0
172 if (.NOT. outside) then
173 write(unit=stdout,fmt='(A,I4,A,F8.2,A,F8.2,A,I3,A,2I3)') &
174 "Report",report,"at",platform%info%lon,"E",platform%info%lat,"N", &
175 "on processor", myproc,"position", platform%loc%x,platform%loc%y
176 end if
177 end if
178
179 time=hdr(4)
180 if (time > 3.0) time=3.0
181 if (time < -3.0) time=-3.0
182
183 t29 = int(0.1 + hdr(7))
184 kx=int(0.1+hdr(5))
185
186 if (kx .eq. 183) then ! reset kx based on t29
187 if (t29 .eq. 511) kx = 181
188 if (t29 .eq. 512) kx = 187
189 if (t29 .eq. 522) kx = 180
190 if (t29 .eq. 523) kx = 180
191 if (t29 .eq. 531) kx = 180
192 if (t29 .eq. 561) kx = 180
193 if (t29 .eq. 562) kx = 180
194 end if
195
196 ! WHY
197 ! if ((kx >= 160) .and. (kx <= 179)) then ! bypass satellite data
198 ! if (t29 .eq. 61 .or. t29 .eq. 63 .or. t29 .ge. 571) then
199 ! cycle reports
200 ! end if
201 ob % total_obs = ob % total_obs + 1
202
203
204 ! Conventional data
205
206 call ufbint(iunit,pmo,2,1,nlevels,'PMO PMQ')
207 ! write(unit=stdout,fmt=*) 'PMO = ',pmo(1,1)
208
209 call ufbint(iunit,obs,9,255,nlevels,obstr)
210 call ufbint(iunit,qms,9,255,nlevels,qmstr)
211 call ufbint(iunit,oes,9,255,nlevels,oestr)
212 if (c8id(1:5) .eq. 'RJOR ') then
213 ! if (kx .eq. 180) then
214 write(unit=stdout,fmt=*) 'report, ',c8id(1:5),' typ = ',kx, &
215 ' nlevels = ',nlevels
216 do ix = 1,nlevels
217 write(unit=stdout,fmt='(9(f10.4,2x))') (obs(jx,ix),jx=1,9)
218 write(unit=stdout,fmt='(9(f10.4,2x))') (qms(jx,ix),jx=1,9)
219 write(unit=stdout,fmt=*) 'PMO = ',pmo(1,1)
220 end do
221 ! kx = 180
222 ! end if
223 end if
224
225 ! WHY
226 ! no need for balloon drift
227 ! select case(kx)
228 ! case (120, 132, 220:221, 232) ; ! read drift
229 ! call ufbint(iunit,drf,9,255,iret,drift)
230 ! do k=1,nlevels
231 ! if (drf(1,k) .ge. 360.)drf(1,k)=drf(1,k)-360.
232 ! if (drf(1,k) .lt. 0.)drf(1,k)=drf(1,k)+360.
233 ! if (abs(drf(2,k)) .gt. 1000. .or. abs(drf(1,k)) .gt. 1000.) then
234 ! drf(2,k)=hdr(3)
235 ! drf(1,k)=hdr(2)
236 ! end if
237 ! if (abs(drf(3,k)) .gt. 3.)then
238 ! drf(3,k)=time
239 ! end if
240 ! end do
241 ! case default ; ! Do nothing
242 ! end select
243
244 platform % info % levels =nlevels
245 platform % info % elv =hdr(6)
246
247 platform % loc % slp %inv = missing_r
248 platform % loc % slp %qc = missing_data
249 platform % loc % slp %error=200.
250 platform % loc % pw %inv = missing_r
251 platform % loc % pw %qc = missing_data
252 platform % loc % pw %error=.2
253
254 do i=1,max_ob_levels
255 platform % each (i) % height = missing_r
256 platform % each (i) % height_qc = missing_data
257
258 platform % each (i) % zk = missing_r
259
260 platform % each (i) % u % inv = missing_r
261 platform % each (i) % u % qc = missing_data
262 platform % each (i) % u % error =10.
263
264 platform % each (i) % v = platform % each (i) % u
265
266 platform % each (i) % t % inv = missing_r
267 platform % each (i) % t % qc = missing_data
268 platform % each (i) % t % error =5.
269
270 platform % each (i) % p % inv = missing_r
271 platform % each (i) % p % qc = missing_data
272 platform % each (i) % p % error =200.
273
274 platform % each (i) % q % inv = missing_r
275 platform % each (i) % q % qc = missing_data
276 platform % each (i) % q % error =.1
277 end do
278
279 do k = 1, platform % info % levels
280 pob=obs(1,k)
281 ! scale q and compute t from tv and set units to Kelvin, if they aren't missing
282
283 if (obs(2,k) .gt. 0. .and. obs(2,k) .lt. 1e9) then
284 obs(2,k) = obs(2,k)*1e-6
285 if (obs(3,k) .gt. -200. .and. obs(3,k) .lt. 300.) then
286 obs(3,k) = (obs(3,k)+t_kelvin) / (1. + 0.608 * obs(2,k))
287 end if
288 else
289 if (obs(3,k) .gt. -200. .and. obs(3,k) .lt. 300.) then
290 obs(3,k) = obs(3,k) + t_kelvin
291 end if
292 end if
293 qob=obs(2,k)
294 tob=obs(3,k)
295 zob=obs(4,k)
296 rob=obs(8,k)
297 pqm = qflag
298 if (abs(qms(1,k)) < qflag) pqm=nint(qms(1,k))
299 qqm = qflag
300 if (abs(qms(2,k)) < qflag) qqm=nint(qms(2,k))
301 tqm = qflag
302 if (abs(qms(3,k)) < qflag) tqm=nint(qms(3,k))
303 zqm = qflag
304 if (abs(qms(4,k)) < qflag) zqm=nint(qms(4,k))
305 wqm = qflag
306 if (abs(qms(5,k)) < qflag) wqm=nint(qms(5,k))
307 pwq = qflag
308 if (abs(qms(7,k)) < qflag) pwq=nint(qms(7,k))
309 rqm = qflag
310 if (abs(qms(9,k)) < qflag) rqm=nint(qms(9,k))
311 pmq = qflag
312 if (abs(pmo(2,1)) < qflag) pmq=nint(pmo(2,1))
313 cat=nint(obs(9,k))
314
315 if (pmq < qflag) then
316 platform % loc % slp % inv =pmo(1,1)*100.
317 platform % loc % slp % qc =pmq
318 platform % loc % slp % error = 200.0
319 end if
320
321 if (pwq < qflag) then
322 platform % loc % pw % inv = obs(7,k) * 0.1 ! convert to cm
323 platform % loc % pw % qc = pwq
324 platform % loc % pw % error = 0.2 ! hardwired to 0.2 cm
325 end if
326
327 poe=100.0
328
329 if (tqm<qflag) then
330 toe=oes(3,k)
331 platform % each (k) % t % inv =obs(3,k)
332 platform % each (k) % t % qc =tqm
333 platform % each (k) % t % error =toe
334 end if
335
336 if (wqm < qflag) then
337 woe=oes(5,k)
338 ! WHY
339 ! if (wqm==0) woe=woe*0.9
340 ! if (wqm==3) woe=woe*1.2
341 platform % each (k) % u % inv =obs(5,k)
342 platform % each (k) % v % inv =obs(6,k)
343 platform % each (k) % u % qc =wqm
344 platform % each (k) % u % error =woe
345 platform % each (k) % v % qc =wqm
346 platform % each (k) % v % error =woe
347
348 ! Convert earth wind to model wind.
349 call da_earth_2_model_wind(obs(5,k), obs(6,k), &
350 platform % each (k) % u % inv, &
351 platform % each (k) % v % inv, &
352 platform%info%lon)
353 if (platform % each (k) % u % inv == 0.0 .and. &
354 platform % each (k) % v % inv == 0.0) then
355 platform % each (k) % u % inv = missing_r
356 platform % each (k) % v % inv = missing_r
357 platform % each (k) % u % qc = missing_data
358 platform % each (k) % v % qc = missing_data
359 end if
360 end if
361
362 if (qqm<qflag .and. obs(2,k)<1.0e9) then
363 qoe=oes(2,k)
364 platform % each (k) % q % inv =obs(2,k)
365 platform % each (k) % q % qc =qqm
366 platform % each (k) % q % error =qoe
367 end if
368
369 if (zqm<qflag)then
370 platform % each (k) % height =zob
371 platform % each (k) % height_qc =zqm
372 end if
373
374 if (rqm<qflag)then
375 roe = oes(8,k)
376 platform % each (k) % u % inv =rob
377 platform % each (k) % u % qc =rqm
378 platform % each (k) % u % error =roe
379 end if
380
381 if (pqm<qflag)then
382 platform % each (k) % p % inv =obs(1,k)*100.
383 platform % each (k) % p % qc =pqm
384 platform % each (k) % p % error =poe
385 end if
386 end do
387
388 nlevels = platform%info%levels
389
390 if (nlevels > max_ob_levels) then
391 nlevels = max_ob_levels
392
393 write(unit=stderr, fmt='(/a/)') &
394 'Warning: Too many levels.'
395
396 write(unit=stderr, fmt='(/2a/2a/2x,a,2f8.2,2(2x,a,f9.2)/2(2x,a,i4)/)') &
397 'Subset: ', platform%info%name(1:8), &
398 'Platfrom: ', trim(platform%info%platform), &
399 'Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
400 'elv: ', platform%info%elv, &
401 'pstar: ', platform%info%pstar, &
402 'level: ', platform%info%levels, &
403 'kx: ', kx
404 else if ((nlevels < 1) .and. ((kx /= 164) .or. (kx /= 174))) then
405 write(unit=stderr, fmt='(/a/)') &
406 'Warning: Too few levels.'
407
408 write(unit=stderr, fmt='(/2a/2a/2x,a,2f8.2,2(2x,a,f9.2)/2(2x,a,i4)/)') &
409 'Subset: ', platform%info%name(1:8), &
410 'Platfrom: ', trim(platform%info%platform), &
411 'Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
412 'elv: ', platform%info%elv, &
413 'pstar: ', platform%info%pstar, &
414 'level: ', platform%info%levels, &
415 'kx: ', kx
416
417 cycle reports
418 end if
419
420 !---------------------------------------------------------------------------
421 ! This is basically converting rh to q i
422 ! Method :
423 ! if rh, temp and pr all available computes Qs otherwise sets Qs= missing
424 ! if rh > 100 sets q = qs otherwise q = rh*Qs/100.
425 ! Note: Currently da_obs_proc_station is active only for ob_format_ascii
426 ! call da_obs_proc_station(platform)
427 !---------------------------------------------------------------------------
428
429 ! Loop over duplicating obs for global
430 ndup = 1
431 if (global .and. &
432 (platform%loc%i < xp%ids .or. platform%loc%i >= xp%ide)) ndup= 2
433
434 ! It is possible that logic for counting obs is incorrect for the
435 ! global case with >1 MPI tasks due to obs duplication, halo, etc.
436 ! TBH: 20050913
437 if (.not.outside) then
438 if (print_detail_obs .and. ndup > 1) then
439 write(unit=stdout, fmt = '(A12,1X,A19,1X,A40,1X,I6,3(F12.3,11X),6X,A5)') &
440 platform%info%platform, &
441 platform%info%date_char, &
442 platform%info%name, &
443 platform%info%levels, &
444 platform%info%lat, &
445 platform%info%lon, &
446 platform%info%elv, &
447 platform%info%id
448
449 write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
450 ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
451 platform%loc%i, platform%loc%j, &
452 platform%loc%dx, platform%loc%dxm, &
453 platform%loc%dy, platform%loc%dym
454 end if
455 end if
456
457 dup_loop: do n = 1, ndup
458 select case(t29)
459 case (11, 12, 13, 22, 23, 31)
460 select case (kx)
461 case (120, 122, 132, 220, 222, 232) ; ! Sound
462 ob%num_sound_glo = ob%num_sound_glo + 1
463 if (outside) then
464 cycle reports
465 end if
466
467 if (.not.Use_SoundObs) cycle reports
468
469 platform % info % platform ='FM-35 TEMP'
470 num_sound = num_sound + 1
471 ! Track serial obs index for reassembly of obs during
472 ! bit-for-bit tests with different numbers of MPI tasks.
473 platform%loc%obs_global_index = ob%num_sound_glo
474
475 ob%sonde_sfc(num_sound)%info = platform%info
476 ob%sonde_sfc(num_sound)%loc = platform%loc
477
478 ! Search to see if we have surface obs.
479
480 surface_level = 0
481
482 do i = 1, nlevels
483 ! if (elevation and height are the same, it is surface)
484 if (abs(platform%info%elv - &
485 platform%each(i)%height) < 0.1) then
486 surface_level = i
487
488 ! Save surface pressure.
489 ob%sonde_sfc(num_sound)%h = platform%each(i)%height
490 ob%sonde_sfc(num_sound)%u = platform%each(i)%u
491 ob%sonde_sfc(num_sound)%v = platform%each(i)%v
492 ob%sonde_sfc(num_sound)%t = platform%each(i)%t
493 ob%sonde_sfc(num_sound)%q = platform%each(i)%q
494 ob%sonde_sfc(num_sound)%p = platform%each(i)%p
495 exit
496 end if
497 end do
498
499 ! processing the sound_sfc data:
500
501 ob%sound(num_sound)%info = platform%info
502 ob%sound(num_sound)%loc = platform%loc
503
504 if (surface_level > 0) then
505 nlevels = nlevels - 1
506 else
507 ob%sonde_sfc(num_sound)%h = missing_r
508 ob%sonde_sfc(num_sound)%u%inv = missing_r
509 ob%sonde_sfc(num_sound)%u%qc = missing
510 ob%sonde_sfc(num_sound)%u%error = abs(missing_r)
511 ob%sonde_sfc(num_sound)%v = ob%sonde_sfc(num_sound)%u
512 ob%sonde_sfc(num_sound)%t = ob%sonde_sfc(num_sound)%u
513 ob%sonde_sfc(num_sound)%p = ob%sonde_sfc(num_sound)%u
514 ob%sonde_sfc(num_sound)%q = ob%sonde_sfc(num_sound)%u
515 end if
516
517 if (nlevels < 1) cycle reports
518
519 allocate (ob%sound(num_sound)%h (1:nlevels))
520 allocate (ob%sound(num_sound)%p (1:nlevels))
521 allocate (ob%sound(num_sound)%zk(1:nlevels))
522 allocate (ob%sound(num_sound)%u (1:nlevels))
523 allocate (ob%sound(num_sound)%v (1:nlevels))
524 allocate (ob%sound(num_sound)%t (1:nlevels))
525 allocate (ob%sound(num_sound)%q (1:nlevels))
526
527 j = 0
528 do i = 1, ob%sound(num_sound)%info%levels
529 if (i == surface_level) cycle
530
531 j=j+1
532
533 ob%sound(num_sound)%h(j) = platform%each(i)%height
534 ob%sound(num_sound)%p(j) = platform%each(i)%p%inv
535 ob%sound(num_sound)%u(j) = platform%each(i)%u
536 ob%sound(num_sound)%v(j) = platform%each(i)%v
537 ob%sound(num_sound)%t(j) = platform%each(i)%t
538 ob%sound(num_sound)%q(j) = platform%each(i)%q
539 end do
540
541 ob%sound(num_sound)%info%levels = nlevels
542 case (221) ; ! Pilot
543 ob%num_pilot_glo = ob%num_pilot_glo + 1
544 if (outside) then
545 cycle reports
546 end if
547
548 if (.not.Use_PilotObs) cycle reports
549 num_pilot = num_pilot + 1
550 ! Track serial obs index for reassembly of obs during
551 ! bit-for-bit tests with different numbers of MPI tasks.
552 platform%loc%obs_global_index = ob%num_pilot_glo
553
554 ob%pilot(num_pilot)%info = platform%info
555 ob%pilot(num_pilot)%loc = platform%loc
556
557 allocate (ob%pilot(num_pilot)%p (1:nlevels))
558 allocate (ob%pilot(num_pilot)%zk(1:nlevels))
559 allocate (ob%pilot(num_pilot)%u (1:nlevels))
560 allocate (ob%pilot(num_pilot)%v (1:nlevels))
561
562 do i = 1, nlevels
563 ob%pilot(num_pilot)%p(i) = platform%each(i)%p%inv
564 ob%pilot(num_pilot)%u(i) = platform%each(i)%u
565 ob%pilot(num_pilot)%v(i) = platform%each(i)%v
566 end do
567 end select
568
569 case (41)
570 ! case (130:131, 133, 230:231, 233) ; ! Airep
571 ob%num_airep_glo = ob%num_airep_glo + 1
572 if (outside) then
573 cycle reports
574 end if
575
576 if (.not.Use_AirepObs) cycle reports
577
578 platform % info % platform ='FM-97 AIREP'
579 num_airep = num_airep + 1
580 ! Track serial obs index for reassembly of obs during
581 ! bit-for-bit tests with different numbers of MPI tasks.
582 platform%loc%obs_global_index = ob%num_airep_glo
583
584 ob % airep (num_airep) % info = platform % info
585 ob % airep (num_airep) % loc = platform % loc
586
587 nlevels= ob % airep (num_airep) % info % levels
588 allocate (ob % airep (num_airep) % h (1:nlevels))
589 allocate (ob % airep (num_airep) % p (1:nlevels))
590 allocate (ob % airep (num_airep) % zk(1:nlevels))
591 allocate (ob % airep (num_airep) % u (1:nlevels))
592 allocate (ob % airep (num_airep) % v (1:nlevels))
593 allocate (ob % airep (num_airep) % t (1:nlevels))
594
595 do i = 1, nlevels
596 ob % airep (num_airep) % h(i) = platform % each(i) % height
597 ob % airep (num_airep) % p(i) = platform % each(i) % p % inv
598 ob % airep (num_airep) % u(i) = platform % each(i) % u
599 ob % airep (num_airep) % v(i) = platform % each(i) % v
600 ob % airep (num_airep) % t(i) = platform % each(i) % t
601 ob % airep (num_airep) %zk(i) = platform % each(i) % zk
602 end do
603
604 ! case (180, 182, 280, 282) ; ! Ships and buoys
605
606 case (522, 523); ! Ships
607 ob%num_ships_glo = ob%num_ships_glo + 1
608 if (outside) then
609 cycle reports
610 end if
611
612 if (.not.Use_ShipsObs) cycle reports
613
614 platform % info % platform ='FM-13 SHIP '
615 num_ships = num_ships + 1
616 ! Track serial obs index for reassembly of obs during
617 ! bit-for-bit tests with different numbers of MPI tasks.
618 platform%loc%obs_global_index = ob%num_ships_glo
619
620 ob % ships (num_ships) % info = platform % info
621 ob % ships (num_ships) % loc = platform % loc
622
623 ob % ships (num_ships) % h = platform % each(1) % height
624 ob % ships (num_ships) % u = platform % each(1) % u
625 ob % ships (num_ships) % v = platform % each(1) % v
626 ob % ships (num_ships) % t = platform % each(1) % t
627 ! WHY
628 ! ob % ships (num_ships) % p = platform % loc % slp
629 ob % ships (num_ships) % p = platform % each(1) % p
630 ob % ships (num_ships) % q = platform % each(1) % q
631 ob % ships (num_ships) % zk= platform % each(1) % zk
632
633 case (531, 561, 562) ; ! Buoy
634 ob%num_buoy_glo = ob%num_buoy_glo + 1
635 if (outside) then
636 cycle reports
637 end if
638
639 if (.not.use_BuoyObs) cycle reports
640 platform % info % platform ='FM-18 BUOY '
641 num_buoy = num_buoy + 1
642 platform%loc%obs_global_index = ob%num_buoy_glo
643
644 ob%buoy(num_buoy)%info = platform%info
645 ob%buoy(num_buoy)%loc = platform%loc
646
647 ob%buoy(num_buoy)%h = platform%each(1)%height
648 ob%buoy(num_buoy)%u = platform%each(1)%u
649 ob%buoy(num_buoy)%v = platform%each(1)%v
650 ob%buoy(num_buoy)%t = platform%each(1)%t
651 ob%buoy(num_buoy)%p = platform%each(1)%p
652 ob%buoy(num_buoy)%q = platform%each(1)%q
653
654 case (511)
655 ! case (181, 281) ; ! Synop
656
657 ob%num_synop_glo = ob%num_synop_glo + 1
658 if (outside) then
659 cycle reports
660 end if
661
662 if (.not.Use_SynopObs) cycle reports
663
664 platform % info % platform ='FM-12 SYNOP'
665
666 num_synop = num_synop + 1
667 platform%loc%obs_global_index = ob%num_synop_glo
668
669 ob % synop (num_synop) % info = platform % info
670 ob % synop (num_synop) % loc = platform % loc
671
672 ob % synop (num_synop) % h = platform % each(1) % height
673 ob % synop (num_synop) % u = platform % each(1) % u
674 ob % synop (num_synop) % v = platform % each(1) % v
675 ob % synop (num_synop) % t = platform % each(1) % t
676 ! WHY
677 ! ob % synop (num_synop) % p = platform % loc % slp
678 ob % synop (num_synop) % p = platform % each(1) % p
679 ob % synop (num_synop) % q = platform % each(1) % q
680 ob % synop (num_synop) % zk= platform % each(1) % zk
681
682 if (ob % synop(num_synop) % h < ob % synop (num_synop) % info % elv) then
683 ob % synop(num_synop) % h = ob % synop (num_synop) % info % elv
684 end if
685
686 case (512)
687 ! case (187, 287) ; ! Metar
688 ob%num_metar_glo = ob%num_metar_glo + 1
689 if (outside) then
690 cycle reports
691 end if
692
693 if (.not.Use_MetarObs) cycle reports
694
695 platform % info % platform ='FM-15 METAR'
696 num_metar = num_metar + 1
697 platform%loc%obs_global_index = ob%num_metar_glo
698
699 ob % metar (num_metar) % info = platform % info
700 ob % metar (num_metar) % loc = platform % loc
701
702 ob % metar (num_metar) % h = platform % each(1) % height
703 ob % metar (num_metar) % u = platform % each(1) % u
704 ob % metar (num_metar) % v = platform % each(1) % v
705 ob % metar (num_metar) % t = platform % each(1) % t
706 ob % metar (num_metar) % p = platform % each(1) % p
707 ob % metar (num_metar) % q = platform % each(1) % q
708 ob % metar (num_metar) % zk= platform % each(1) % zk
709
710 case (63)
711 ! case (242:246, 252:253, 255) ; ! Geo. CMVs
712 ob%num_geoamv_glo = ob%num_geoamv_glo + 1
713 if (outside) then
714 cycle reports
715 end if
716
717 if (.not.Use_geoamvObs) cycle reports
718 platform % info % platform ='FM-88 SATOB'
719 num_geoamv = num_geoamv + 1
720 platform%loc%obs_global_index = ob%num_geoamv_glo
721
722 ob % geoamv (num_geoamv) % info = platform % info
723 ob % geoamv (num_geoamv) % loc = platform % loc
724 allocate (ob%geoamv(num_geoamv)%p (1:nlevels))
725 allocate (ob%geoamv(num_geoamv)%zk(1:nlevels))
726 allocate (ob%geoamv(num_geoamv)%u (1:nlevels))
727 allocate (ob%geoamv(num_geoamv)%v (1:nlevels))
728
729 do i = 1, nlevels
730 ob % geoamv (num_geoamv) % p(i) = platform % each(i) % p % inv
731 ob % geoamv (num_geoamv) % u(i) = platform % each(i) % u
732 ob % geoamv (num_geoamv) % v(i) = platform % each(i) % v
733 ob % geoamv (num_geoamv) %zk(i) = platform % each(i) % zk
734 end do
735
736 case (582)
737 ob%num_qscat_glo = ob%num_qscat_glo + 1
738 if (outside) then
739 cycle reports
740 end if
741
742 if (.not.use_qscatobs) cycle reports
743 platform % info % platform ='FM-281 Quiks'
744 num_qscat = num_qscat + 1
745 platform%loc%obs_global_index = ob%num_qscat_glo
746
747 ob%qscat(num_qscat)%info = platform%info
748 ob%qscat(num_qscat)%loc = platform%loc
749 ! WHY
750 ! ob%qscat(num_qscat)%h = platform%each(1)%height
751 ! prepbufr uses pressure not height, so hardwire height to
752 ! 0 (sea-level)
753 ob%qscat(num_qscat)%h = 0.
754 ob%qscat(num_qscat)%u = platform%each(1)%u
755 ob%qscat(num_qscat)%v = platform%each(1)%v
756 ob%qscat(num_qscat)%u%error = max(platform%each(1)%u%error,1.0)
757 ob%qscat(num_qscat)%v%error = max(platform%each(1)%v%error,1.0)
758
759 case (583) ! GPS PW
760 ob%num_gpspw_glo = ob%num_gpspw_glo + 1
761 if (outside) then
762 cycle reports
763 end if
764
765 if (.not.use_GpspwObs) cycle reports
766 platform % info % platform ='FM-111 GPSPW'
767 num_gpspw = num_gpspw + 1
768 platform%loc%obs_global_index = ob%num_gpspw_glo
769
770 ob%gpspw(num_gpspw)%info = platform%info
771 ob%gpspw(num_gpspw)%loc = platform%loc
772 ob%gpspw(num_gpspw)%tpw = platform%loc%pw
773
774 case (584) ! GPS REF
775 ob%num_gpsref_glo = ob%num_gpsref_glo + 1
776 if (outside) then
777 cycle reports
778 end if
779
780 if (.not.use_GpsrefObs) cycle reports
781 platform % info % platform ='FM-116 GPSRF'
782 num_gpsref = num_gpsref + 1
783 platform%loc%obs_global_index = ob%num_gpsref_glo
784
785 ob%gpsref(num_gpsref)%info = platform%info
786 ob%gpsref(num_gpsref)%loc = platform%loc
787
788 allocate (ob%gpsref(num_gpsref)%h (1:nlevels))
789 allocate (ob%gpsref(num_gpsref)%zk(1:nlevels))
790 allocate (ob%gpsref(num_gpsref)%ref(1:nlevels))
791 allocate (ob%gpsref(num_gpsref)% p(1:nlevels))
792 allocate (ob%gpsref(num_gpsref)% t(1:nlevels))
793 allocate (ob%gpsref(num_gpsref)% q(1:nlevels))
794
795 do i = 1, nlevels
796 ob%gpsref(num_gpsref)%h(i) = platform%each(i)%height
797 ! u is used to store ref in this routine (ascii version uses td)
798 ob%gpsref(num_gpsref)%ref(i) = platform%each(i)%u
799 ! Keep the retrieved p and t (and q) as "field_type":
800 ob%gpsref(num_gpsref)%p(i) = platform%each(i)%p
801 ob%gpsref(num_gpsref)%t(i) = platform%each(i)%t
802 ob%gpsref(num_gpsref)%q(i) = platform%each(i)%q
803 end do
804
805 case (71, 72)
806 ! case (223, 224 ) ; ! Profiler & VADWND - NEXRAD winds
807 ob%num_profiler_glo = ob%num_profiler_glo + 1
808 if (outside) then
809 cycle reports
810 end if
811
812 if (.not.use_ProfilerObs) cycle reports
813 platform % info % platform ='FM-132 PRFLR'
814 num_profiler = num_profiler + 1
815 ! Track serial obs index for reassembly of obs during bit-for-bit
816 ! tests with different numbers of MPI tasks.
817 platform%loc%obs_global_index = ob%num_profiler_glo
818
819 ob%profiler(num_profiler)%info = platform%info
820 ob%profiler(num_profiler)%loc = platform%loc
821
822 allocate (ob%profiler(num_profiler)%p (1:nlevels))
823 allocate (ob%profiler(num_profiler)%zk(1:nlevels))
824 allocate (ob%profiler(num_profiler)%u (1:nlevels))
825 allocate (ob%profiler(num_profiler)%v (1:nlevels))
826
827 do i = 1, nlevels
828 ob%profiler(num_profiler)%p(i) = platform%each(i)%p%inv
829 ob%profiler(num_profiler)%u(i) = platform%each(i)%u
830 ob%profiler(num_profiler)%v(i) = platform%each(i)%v
831 end do
832
833 case default
834 select case (kx)
835 case (111 , 210) ; ! Tropical Cyclone Bogus
836 ! Note Tropical cyclone Bougus is given type 135 in Obs-ascii
837 ob%num_bogus_glo = ob%num_bogus_glo + 1
838 if (outside) then
839 cycle reports
840 end if
841
842 if (.not.use_BogusObs) cycle reports
843 platform % info % platform ='FM-135 TCBOG'
844 num_bogus = num_bogus + 1
845 ! Track serial obs index for reassembly of obs during bit-for-bit
846 ! tests with different numbers of MPI tasks.
847 platform%loc%obs_global_index = ob%num_bogus_glo
848
849 if (num_bogus > max_bogus_input) then
850 write(unit=message(1),fmt='(A,I6,A,I6)') &
851 'Bogus #=', num_bogus, ' > max_bogus_input=', max_bogus_input
852 call da_error(__FILE__,__LINE__,message(1:1))
853 end if
854
855 ob%bogus(num_bogus)%info = platform%info
856 ob%bogus(num_bogus)%loc = platform%loc
857
858 allocate (ob%bogus(num_bogus)%h (1:nlevels))
859 allocate (ob%bogus(num_bogus)%p (1:nlevels))
860 allocate (ob%bogus(num_bogus)%zk(1:nlevels))
861 allocate (ob%bogus(num_bogus)%u (1:nlevels))
862 allocate (ob%bogus(num_bogus)%v (1:nlevels))
863 allocate (ob%bogus(num_bogus)%t (1:nlevels))
864 allocate (ob%bogus(num_bogus)%q (1:nlevels))
865
866 do i = 1, nlevels
867 ob%bogus(num_bogus)%h(i) = platform%each(i)%height
868 ob%bogus(num_bogus)%p(i) = platform%each(i)%p%inv
869 ob%bogus(num_bogus)%u(i) = platform%each(i)%u
870 ob%bogus(num_bogus)%v(i) = platform%each(i)%v
871 ob%bogus(num_bogus)%t(i) = platform%each(i)%t
872 ob%bogus(num_bogus)%q(i) = platform%each(i)%q
873 end do
874
875 ob%bogus(num_bogus)%slp = platform%loc%slp
876
877 case default
878 write(unit=message(1), fmt='(a, i12)') &
879 'unsaved obs found with type: ',kx
880 call da_warning(__FILE__,__LINE__,message(1:1))
881 end select
882 end select
883
884 if (global .and. n < 2) then
885 if (Testing_WRFVAR) exit dup_loop
886 if (platform%loc % i >= xp % ide) then
887 platform%loc%i = platform%loc % i - xp % ide
888 else if (platform%loc % i < xp % ids) then
889 platform%loc%i = platform%loc % i + xp % ide
890 end if
891 platform%loc%proc_domain = .not. platform%loc%proc_domain
892 end if
893 end do dup_loop
894 end do reports
895
896 call closbf(iunit)
897 close(iunit)
898 call da_free_unit(iunit)
899
900 ob%num_sound = num_sound
901 ob%num_synop = num_synop
902 ob%num_pilot = num_pilot
903 ob%num_satem = num_satem
904 ob%num_geoamv = num_geoamv
905 ob%num_polaramv = num_polaramv
906 ob%num_airep = num_airep
907 ob%num_gpspw = num_gpspw
908 ob%num_gpsref = num_gpsref
909 ob%num_metar = num_metar
910 ob%num_ships = num_ships
911 ob%num_qscat = num_qscat
912 ob%num_buoy = num_buoy
913 ob%num_profiler = num_profiler
914 ob%num_bogus = num_bogus
915
916 ! WHY
917 ! ob%num_ssmt1 = num_ssmt1
918 ! ob%num_ssmt2 = num_ssmt2
919 ! ob%num_ssmi_tb = num_ssmi_tb
920 ! ob%num_ssmi_retrieval = num_ssmi_retrieval
921
922 if (print_detail_obs) then
923 write(unit=stdout, fmt='(/a,i6/)') 'ob%current_ob_time=', ob%current_ob_time
924
925 write(unit=stdout, fmt='(5x,a,i6,a)') &
926 'Read: ', num_sound, ' SOUND reports,', &
927 'Read: ', num_synop, ' SYNOP reports,', &
928 'Read: ', num_pilot, ' PILOT reports,', &
929 'Read: ', num_satem, ' SATEM reports,', &
930 'Read: ', num_geoamv, ' Geo AMV reports,', &
931 'Read: ', num_polaramv, ' Polar AMV reports,', &
932 'Read: ', num_airep, ' AIREP reports,', &
933 'Read: ', num_gpspw , ' GPSPW/GPSZD reports,', &
934 'Read: ', num_gpsref, ' GPSRF reports,', &
935 'Read: ', num_metar, ' METAR reports,', &
936 'Read: ', num_ships , ' SHIP reports,', &
937 'Read: ', num_ssmi_retrieval , ' SSMI_RETRIEVAL reports,', &
938 'Read: ', num_ssmi_tb, ' SSMI_TB reports,', &
939 'Read: ', num_ssmt1, ' SSMT1 reports,', &
940 'Read: ', num_ssmt2, ' SSMT2 reports,', &
941 'Read: ', num_qscat, ' QSCAT reports,', &
942 'Read: ', num_profiler, ' Profiler reports,', &
943 'Read: ', num_buoy, ' Buoy reports,', &
944 'Read: ', num_bogus, ' Bogus reports,', &
945 'Read: ', total_obs, ' Total reports.', &
946 'There are ', total_obs - num_sound - num_synop &
947 - num_pilot - num_satem &
948 - num_geoamv - num_polaramv - num_airep &
949 - num_metar - num_ships &
950 - num_ssmi_retrieval &
951 - num_ssmi_tb - num_ssmt1 - num_ssmt2 &
952 - num_gpspw - num_gpsref - num_qscat &
953 - num_profiler - num_buoy, &
954 - num_bogus - num_bogus, &
955 ' reports unsaved.'
956 write(unit=stdout, fmt=*) ' '
957 end if
958
959 if ((ob%ob_numb(ob%current_ob_time)%sound /= ob%num_sound) .or. &
960 (ob%ob_numb(ob%current_ob_time)%synop /= ob%num_synop) .or. &
961 (ob%ob_numb(ob%current_ob_time)%pilot /= ob%num_pilot) .or. &
962 (ob%ob_numb(ob%current_ob_time)%satem /= ob%num_satem) .or. &
963 (ob%ob_numb(ob%current_ob_time)%geoamv /= ob%num_geoamv) .or. &
964 (ob%ob_numb(ob%current_ob_time)%polaramv /= ob%num_polaramv) .or. &
965 (ob%ob_numb(ob%current_ob_time)%airep /= ob%num_airep) .or. &
966 (ob%ob_numb(ob%current_ob_time)%gpspw /= ob%num_gpspw) .or. &
967 (ob%ob_numb(ob%current_ob_time)%gpsref /= ob%num_gpsref) .or. &
968 (ob%ob_numb(ob%current_ob_time)%metar /= ob%num_metar) .or. &
969 (ob%ob_numb(ob%current_ob_time)%ships /= ob%num_ships) .or. &
970 (ob%ob_numb(ob%current_ob_time)%qscat /= ob%num_qscat) .or. &
971 (ob%ob_numb(ob%current_ob_time)%buoy /= ob%num_buoy) .or. &
972 (ob%ob_numb(ob%current_ob_time)%bogus /= ob%num_bogus) .or. &
973 (ob%ob_numb(ob%current_ob_time)%ssmt2 /= ob%num_ssmt1) .or. &
974 (ob%ob_numb(ob%current_ob_time)%ssmt2 /= ob%num_ssmt2) .or. &
975 (ob%ob_numb(ob%current_ob_time)%profiler /= ob%num_profiler)) then
976
977 write(unit=stderr, fmt='(a, i6, 2x, a, i6)') &
978 'ob%ob_numb(ob%current_ob_time)%sound=', ob%ob_numb(ob%current_ob_time)%sound, &
979 'ob%num_sound=', ob%num_sound, &
980 'ob%ob_numb(ob%current_ob_time)%synop=', ob%ob_numb(ob%current_ob_time)%synop, &
981 'ob%num_synop=', ob%num_synop, &
982 'ob%ob_numb(ob%current_ob_time)%pilot=', ob%ob_numb(ob%current_ob_time)%pilot, &
983 'ob%num_pilot=', ob%num_pilot, &
984 'ob%ob_numb(ob%current_ob_time)%satem=', ob%ob_numb(ob%current_ob_time)%satem, &
985 'ob%num_satem=', ob%num_satem, &
986 'ob%ob_numb(ob%current_ob_time)%geoamv=', ob%ob_numb(ob%current_ob_time)%geoamv, &
987 'ob%num_geoamv=', ob%num_geoamv, &
988 'ob%ob_numb(ob%current_ob_time)%polaramv=', ob%ob_numb(ob%current_ob_time)%polaramv, &
989 'ob%num_polaramv=', ob%num_polaramv, &
990 'ob%ob_numb(ob%current_ob_time)%airep=', ob%ob_numb(ob%current_ob_time)%airep, &
991 'ob%num_airep=', ob%num_airep, &
992 'ob%ob_numb(ob%current_ob_time)%gpspw=', ob%ob_numb(ob%current_ob_time)%gpspw, &
993 'ob%num_gpspw=', ob%num_gpspw, &
994 'ob%ob_numb(ob%current_ob_time)%gpsref=', ob%ob_numb(ob%current_ob_time)%gpsref,&
995 'ob%num_gpsref=', ob%num_gpsref, &
996 'ob%ob_numb(ob%current_ob_time)%metar=', ob%ob_numb(ob%current_ob_time)%metar, &
997 'ob%num_metar=', ob%num_metar, &
998 'ob%ob_numb(ob%current_ob_time)%ships=', ob%ob_numb(ob%current_ob_time)%ships, &
999 'ob%num_ships=', ob%num_ships, &
1000 'ob%ob_numb(ob%current_ob_time)%qscat=', ob%ob_numb(ob%current_ob_time)%qscat, &
1001 'ob%num_qscat=', ob%num_qscat, &
1002 'ob%ob_numb(ob%current_ob_time)%buoy =', ob%ob_numb(ob%current_ob_time)%buoy , &
1003 'ob%num_buoy =', ob%num_buoy , &
1004 'ob%ob_numb(ob%current_ob_time)%bogus=', ob%ob_numb(ob%current_ob_time)%bogus, &
1005 'ob%num_bogus=', ob%num_bogus, &
1006 'ob%ob_numb(ob%current_ob_time)%ssmt1=', ob%ob_numb(ob%current_ob_time)%ssmt1, &
1007 'ob%num_ssmt1=', ob%num_ssmt1, &
1008 'ob%ob_numb(ob%current_ob_time)%ssmt2=', ob%ob_numb(ob%current_ob_time)%ssmt2, &
1009 'ob%num_ssmt2=', ob%num_ssmt2, &
1010 'ob%ob_numb(ob%current_ob_time)%profiler=', ob%ob_numb(ob%current_ob_time)%profiler, &
1011 'ob%num_profiler=', ob%num_profiler
1012
1013 call da_error(__FILE__,__LINE__,(/"Obs mismatch"/))
1014 end if
1015
1016 #endif
1017
1018 end subroutine da_read_bufr_obs
1019
1020