da_read_obs_ascii.inc
References to this file elsewhere.
1 subroutine da_read_obs_ascii (iv, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Read a GTS observation file
5 !
6 ! Modifications: 10/26/2004 Syed RH Rizvi
7 !
8 ! Fix Obs-Long as -180 if it is +180
9 !
10 ! Date: 03/04/2005 - Syed RH Rizvi
11 !
12 ! (a) AMVs from Geostationary and Polar orbiting satellite are
13 ! seperated & used as profile. Currently there is a provision
14 ! to use MODIS winds only.
15 ! (b) MODIS obs error are currently assigned through namelist
16 ! parameter (modis_cmv_error)
17 !
18 ! 03/30/2005 Syed RH Rizvi
19 ! For global option duplicate obs at East/West boundary
20 !
21 ! 08/30/2005 Syed RH Rizvi
22 ! Writing original errors & thicknesses
23 ! desired for outputting QC obs with NTMAX = 0
24 !
25 !---------------------------------------------------------------------------
26
27 implicit none
28
29 type (iv_type), intent(inout) :: iv
30 character(len=*), intent(in) :: filename
31
32 character (len = 10) :: fmt_name
33
34 character (len = 160) :: info_string
35
36 character (len = 160) :: fmt_info, fmt_loc, fmt_each
37
38 integer :: i, j, iost, nlevels, old_nlevels, fm,iunit
39
40 type (multi_level_type) :: platform
41
42 logical :: outside
43 logical :: outside_all
44 integer :: surface_level
45 real :: height_error, u_comp, v_comp
46 integer :: nlocal(num_ob_indexes)
47 integer :: ntotal(num_ob_indexes)
48
49 integer :: ndup, n, report, obs_index
50
51 if (trace_use) call da_trace_entry("da_read_obs_ascii")
52
53 ! Local counts of obs, initialise from end of last file
54
55 nlocal(:) = iv%info(:)%plocal(iv%time-1)
56 ntotal(:) = iv%info(:)%ptotal(iv%time-1)
57
58 ! open file
59 ! =========
60
61 call da_get_unit(iunit)
62 open(unit = iunit, &
63 FILE = trim(filename), &
64 FORM = 'FORMATTED', &
65 ACCESS = 'SEQUENTIAL', &
66 iostat = iost, &
67 STATUS = 'OLD')
68
69 if (iost /= 0) then
70 write(unit=message(1),fmt='(A,I5,A)') &
71 "Error",iost," opening gts obs file "//trim(filename)
72 call da_warning(__FILE__,__LINE__,message(1:1))
73 call da_free_unit(iunit)
74 if (trace_use) call da_trace_exit("da_read_obs_ascii")
75 return
76 end if
77
78 ! read header
79 ! ===========
80
81 do
82 read (unit = iunit, fmt = '(A)', iostat = iost) info_string
83 if (iost /= 0) then
84 call da_warning(__FILE__,__LINE__, &
85 (/"Problem reading gts obs header, skipping file"/))
86 if (trace_use) call da_trace_exit("da_read_obs_ascii")
87 return
88 end if
89
90 if (info_string(1:6) == 'EACH ') exit
91 end do
92
93 ! read formats
94 ! ------------
95
96 read (iunit, fmt = '(A,1X,A)', iostat = iost) &
97 fmt_name, fmt_info, &
98 fmt_name, fmt_loc, &
99 fmt_name, fmt_each
100
101 if (iost /= 0) then
102 call da_warning(__FILE__,__LINE__, &
103 (/"Problem reading gts obs file, skipping"/))
104 if (trace_use) call da_trace_exit("da_read_obs_ascii")
105 return
106 end if
107
108 if (print_detail_obs) then
109 write(unit=stdout, fmt='(2a)') &
110 'fmt_info=', fmt_info, &
111 'fmt_loc =', fmt_loc, &
112 'fmt_each=', fmt_each
113 end if
114
115 ! skip line
116 ! ----------
117
118 read (iunit, fmt = '(a)') fmt_name
119
120 ! loop over records
121 ! -----------------
122
123 report = 0 ! report number in file
124
125 reports: &
126 do
127
128 report = report+1
129
130 ! read station general info
131 ! =============================
132
133 read (iunit, fmt = fmt_info, iostat = iost) &
134 platform%info%platform, &
135 platform%info%date_char, &
136 platform%info%name, &
137 platform%info%levels, &
138 platform%info%lat, &
139 platform%info%lon, &
140 platform%info%elv, &
141 platform%info%id
142
143 if (iost /= 0) then
144 ! FIX? This is expected, but its unclear how we handle failure
145 ! here without assuming the fortran2003 convention on
146 ! error statuses
147 exit reports
148 end if
149
150 if (print_detail_obs) then
151 write(unit=stdout, fmt=fmt_info) &
152 platform%info%platform, &
153 platform%info%date_char, &
154 platform%info%name, &
155 platform%info%levels, &
156 platform%info%lat, &
157 platform%info%lon, &
158 platform%info%elv, &
159 platform%info%id
160 end if
161
162 if (platform%info%lon == 180.0 ) platform%info%lon =-180.000
163 ! Fix funny wind direction at Poles
164 if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
165 platform%info%lon = 0.0
166 end if
167
168 if (platform%info%platform(6:6) == ' ') then
169 read(platform%info%platform(4:5), '(I2)') fm
170 else
171 read(platform%info%platform(4:6), '(I3)') fm
172 end if
173
174 ! read model location
175 ! =========================
176
177 read (iunit, fmt = fmt_loc) platform%loc%slp, platform%loc%pw
178
179 ! levels < 1 and not GPSPW, go back to reports
180
181 if ((platform%info%levels < 1) .AND. (index(platform%info%platform, 'GPSPW') <= 0)) then
182 cycle reports
183 end if
184
185 ! read each level
186 ! ---------------
187
188 do i = 1, platform%info%levels
189 platform%each (i) = each_level_type(missing_r, missing, -1.0, & ! height
190 field_type(missing_r, missing, missing_r), & ! u
191 field_type(missing_r, missing, missing_r), & ! v
192 field_type(missing_r, missing, missing_r), & ! p
193 field_type(missing_r, missing, missing_r), & ! t
194 field_type(missing_r, missing, missing_r), & ! q
195 field_type(missing_r, missing, missing_r), & ! rh
196 field_type(missing_r, missing, missing_r), & ! td
197 field_type(missing_r, missing, missing_r)) ! speed
198
199 read (unit = iunit, fmt = trim (fmt_each)) &
200 platform%each (i)%p, &
201 platform%each (i)%speed, &
202 ! Here the 'direction' is stored in platform%each (i)%v:
203 platform%each (i)%v, &
204 platform%each (i)%height, &
205 platform%each (i)%height_qc, &
206 height_error, &
207 platform%each (i)%t, &
208 platform%each (i)%td, &
209 platform%each (i)%rh
210
211 if (print_detail_obs) then
212 write(unit=stdout, fmt = '(a, i6)') 'i=', i
213 write(unit=stdout, fmt = trim (fmt_each)) &
214 platform%each (i)%p, &
215 platform%each (i)%speed, &
216 platform%each (i)%v, &
217 platform%each (i)%height, &
218 platform%each (i)%height_qc, &
219 height_error, &
220 platform%each (i)%t, &
221 platform%each (i)%td, &
222 platform%each (i)%rh
223 end if
224
225 ! Uncomment the following if errors in reported heights (test later):
226 ! if (fm == 97 .or. fm == 96 .or. fm == 42) then
227 ! platform%each (i)%height_qc = -5 ! Aircraft heights are not above surface.
228 ! end if
229
230 ! To convert the wind speed and direction to
231 ! the model wind components: u, v
232
233 if (platform%each (i)%speed%qc /= missing_data .and. &
234 platform%each (i)%v%qc /= missing_data) then
235 call da_ffdduv (platform%each (i)%speed%inv, platform%each (i)%v%inv, &
236 U_comp, V_comp, platform%info%lon, convert_fd2uv)
237 platform%each (i)%u = platform%each (i)%speed
238 platform%each (i)%v = platform%each (i)%speed
239 platform%each (i)%u%inv = U_comp
240 platform%each (i)%v%inv = V_comp
241 else
242 platform%each (i)%u%inv = missing_r
243 platform%each (i)%v%inv = missing_r
244 platform%each (i)%u%error = missing_r
245 platform%each (i)%v%error = missing_r
246 platform%each (i)%u%qc = missing_data
247 platform%each (i)%v%qc = missing_data
248 end if
249 end do
250
251 ! Restrict to a range of reports, useful for debugging
252
253 if (report < report_start) then
254 cycle
255 end if
256
257 if (report > report_end) then
258 exit
259 end if
260
261 call da_llxy (platform%info, platform%loc, outside, outside_all)
262
263 if (outside_all) then
264 cycle reports
265 end if
266
267 if (print_detail_obs) then
268 ! Simplistic approach, could be improved to get it all done on PE 0
269 if (.NOT. outside) then
270 write(unit=stdout,fmt='(A,I5,A,F8.2,A,F8.2,A,I3,A,2F8.2)') &
271 "Report",report," at",platform%info%lon,"E",platform%info%lat, &
272 "N on processor", myproc,", position", platform%loc%x,platform%loc%y
273 end if
274 end if
275
276 nlevels = platform%info%levels
277 if (nlevels > max_ob_levels) then
278 nlevels = max_ob_levels
279 write(unit=message(1), fmt='(/4(2a,2x),a,2f8.2,2x,2(a,f9.2,2x),2(a,i4,2x)/)') &
280 'Station: ', trim(platform%info%name), &
281 'ID: ', trim(platform%info%id), &
282 'Platform: ', trim(platform%info%platform), &
283 'Date: ', trim(platform%info%date_char), &
284 'At Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
285 'At elv: ', platform%info%elv, &
286 'with pstar: ', platform%info%pstar, &
287 'Has level: ', platform%info%levels, &
288 'which is great than max_ob_levels: ', max_ob_levels
289
290 write (unit=message(2), fmt = '(A,1X,A,1X,A,1X,I4,2f9.3,f9.1,1X,A)') &
291 platform%info%name, &
292 platform%info%platform, &
293 platform%info%date_char, &
294 platform%info%levels, &
295 platform%info%lat, &
296 platform%info%lon, &
297 platform%info%elv, &
298 platform%info%id
299 call da_warning(__FILE__,__LINE__,message(1:2))
300
301 platform%info%levels = nlevels
302 else if (nlevels < 1) then
303 ! Not GPSPW and GPSZD data:
304 if (fm /= 111 .and. fm /= 114) then
305 cycle reports
306 end if
307 end if
308
309 if (fm /= 86) call da_obs_proc_station(platform)
310
311 nlevels = platform%info%levels
312 ! Loop over duplicating obs for global
313 ndup = 1
314 if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
315
316 ! It is possible that logic for counting obs is incorrect for the
317 ! global case with >1 MPI tasks due to obs duplication, halo, etc.
318 ! TBH: 20050913
319
320 if (.not.outside) then
321 if (print_detail_obs .and. ndup > 1) then
322 write(unit=stdout, fmt = fmt_info) &
323 platform%info%platform, &
324 platform%info%date_char, &
325 platform%info%name, &
326 platform%info%levels, &
327 platform%info%lat, &
328 platform%info%lon, &
329 platform%info%elv, &
330 platform%info%id
331
332 write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
333 ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
334 platform%loc%i, platform%loc%j, &
335 platform%loc%dx, platform%loc%dxm, &
336 platform%loc%dy, platform%loc%dym
337 end if
338 end if
339
340 obs_index = fm_index(fm)
341
342 dup_loop: do n = 1, ndup
343 select case(fm)
344
345 case (12) ;
346 if (.not. use_synopobs) cycle reports
347 if (n==1) ntotal(synop) = ntotal(synop)+1
348 if (outside) cycle reports
349 nlocal(synop) = nlocal(synop) + 1
350 if (nlocal(synop) > iv%info(synop)%nlocal) cycle reports
351
352 iv%synop(nlocal(synop))%h = platform%each(1)%height
353 iv%synop(nlocal(synop))%u = platform%each(1)%u
354 iv%synop(nlocal(synop))%v = platform%each(1)%v
355 iv%synop(nlocal(synop))%t = platform%each(1)%t
356 iv%synop(nlocal(synop))%q = platform%each(1)%q
357 iv%synop(nlocal(synop))%p = platform%each(1)%p
358
359 case (13, 17) ; ! ships
360 if (.not. use_shipsobs) cycle reports
361 if (n==1) ntotal(ships) = ntotal(ships) + 1
362 if (outside) cycle reports
363 nlocal(ships) = nlocal(ships) + 1
364 if (nlocal(ships) > iv%info(ships)%nlocal) cycle reports
365
366 iv%ships(nlocal(ships))%h = platform%each(1)%height
367 iv%ships(nlocal(ships))%u = platform%each(1)%u
368 iv%ships(nlocal(ships))%v = platform%each(1)%v
369 iv%ships(nlocal(ships))%t = platform%each(1)%t
370 iv%ships(nlocal(ships))%p = platform%each(1)%p
371 iv%ships(nlocal(ships))%q = platform%each(1)%q
372
373 case (15:16) ;
374 if (.not. use_metarobs) cycle reports
375 if (n==1) ntotal(metar) = ntotal(metar) + 1
376 if (outside) cycle reports
377 nlocal(metar) = nlocal(metar) + 1
378 if (nlocal(metar) > iv%info(metar)%nlocal) cycle reports
379
380 iv%metar(nlocal(metar))%h = platform%each(1)%height
381 iv%metar(nlocal(metar))%u = platform%each(1)%u
382 iv%metar(nlocal(metar))%v = platform%each(1)%v
383 iv%metar(nlocal(metar))%t = platform%each(1)%t
384 iv%metar(nlocal(metar))%p = platform%each(1)%p
385 iv%metar(nlocal(metar))%q = platform%each(1)%q
386
387 case (32:34) ;
388 if (.not. use_pilotobs) cycle reports
389 if (n==1) ntotal(pilot) = ntotal(pilot) + 1
390 if (outside) cycle reports
391 nlocal(pilot) = nlocal(pilot) + 1
392 if (nlocal(pilot) > iv%info(pilot)%nlocal) cycle reports
393
394 allocate (iv%pilot(nlocal(pilot))%p (1:nlevels))
395 allocate (iv%pilot(nlocal(pilot))%u (1:nlevels))
396 allocate (iv%pilot(nlocal(pilot))%v (1:nlevels))
397
398 do i = 1, nlevels
399 iv%pilot(nlocal(pilot))%p(i) = platform%each(i)%p%inv
400 iv%pilot(nlocal(pilot))%u(i) = platform%each(i)%u
401 iv%pilot(nlocal(pilot))%v(i) = platform%each(i)%v
402 end do
403
404 case (35:38) ;
405 if (.not. use_soundobs) cycle reports
406 if (n==1) ntotal(sound) = ntotal(sound) + 1
407 if (n==1) ntotal(sonde_sfc) = ntotal(sonde_sfc) + 1
408 if (outside) cycle reports
409 nlocal(sound) = nlocal(sound) + 1
410 nlocal(sonde_sfc) = nlocal(sonde_sfc) + 1
411 if (nlocal(sound) > iv%info(sound)%nlocal) cycle reports
412
413 ! Search to see if we have surface obs.
414
415 surface_level = 0
416
417 do i = 1, nlevels
418 ! if (elevation and height are the same, it is surface.)
419 if (abs(platform%info%elv - platform%each(i)%height) < 0.1) then
420 surface_level = i
421
422 ! Save surface pressure.
423 iv%sonde_sfc(nlocal(sonde_sfc))%h = platform%each(i)%height
424 iv%sonde_sfc(nlocal(sonde_sfc))%u = platform%each(i)%u
425 iv%sonde_sfc(nlocal(sonde_sfc))%v = platform%each(i)%v
426 iv%sonde_sfc(nlocal(sonde_sfc))%t = platform%each(i)%t
427 iv%sonde_sfc(nlocal(sonde_sfc))%q = platform%each(i)%q
428 iv%sonde_sfc(nlocal(sonde_sfc))%p = platform%each(i)%p
429
430 exit
431 end if
432 end do
433
434 ! processing the sonde_sfc data:
435
436 old_nlevels = nlevels
437
438 if (surface_level > 0) then
439 nlevels = nlevels - 1
440 else
441 iv%sonde_sfc(nlocal(sonde_sfc))%h = missing_r
442 iv%sonde_sfc(nlocal(sonde_sfc))%u%inv = missing_r
443 iv%sonde_sfc(nlocal(sonde_sfc))%u%qc = missing
444 iv%sonde_sfc(nlocal(sonde_sfc))%u%error = abs(missing_r)
445 iv%sonde_sfc(nlocal(sonde_sfc))%v = iv%sonde_sfc(nlocal(sonde_sfc))%u
446 iv%sonde_sfc(nlocal(sonde_sfc))%t = iv%sonde_sfc(nlocal(sonde_sfc))%u
447 iv%sonde_sfc(nlocal(sonde_sfc))%p = iv%sonde_sfc(nlocal(sonde_sfc))%u
448 iv%sonde_sfc(nlocal(sonde_sfc))%q = iv%sonde_sfc(nlocal(sonde_sfc))%u
449 end if
450
451 if (nlevels > 0) then
452
453 allocate (iv%sound(nlocal(sound))%h (1:nlevels))
454 allocate (iv%sound(nlocal(sound))%p (1:nlevels))
455 allocate (iv%sound(nlocal(sound))%u (1:nlevels))
456 allocate (iv%sound(nlocal(sound))%v (1:nlevels))
457 allocate (iv%sound(nlocal(sound))%t (1:nlevels))
458 allocate (iv%sound(nlocal(sound))%q (1:nlevels))
459
460 j = 0
461 ! do i = 1, iv%sound(nlocal(sound))%info%levels
462 do i = 1, old_nlevels
463 if (i == surface_level) cycle
464
465 j=j+1
466
467 iv%sound(nlocal(sound))%h(j) = platform%each(i)%height
468 iv%sound(nlocal(sound))%p(j) = platform%each(i)%p%inv
469 iv%sound(nlocal(sound))%u(j) = platform%each(i)%u
470 iv%sound(nlocal(sound))%v(j) = platform%each(i)%v
471 iv%sound(nlocal(sound))%t(j) = platform%each(i)%t
472 iv%sound(nlocal(sound))%q(j) = platform%each(i)%q
473 end do
474 end if
475
476 ! iv%sound(nlocal(sound))%info%levels = nlevels
477
478
479 case (86) ;
480 ! Reject cloudy satem obs.
481 if (.not. use_satemobs .or. platform%loc%pw%inv > 10.0) cycle reports
482 if (n==1) ntotal(satem) = ntotal(satem) + 1
483 if (outside) cycle reports
484 nlocal(satem) = nlocal(satem) + 1
485 if (nlocal(satem) > iv%info(satem)%nlocal) cycle reports
486
487 ! The satem ref_p is put in the SLP position in OBSPROC data stream.
488
489 iv%satem(nlocal(satem))%ref_p= platform%loc%slp%inv
490
491 allocate (iv%satem(nlocal(satem))%p (1:nlevels))
492 allocate (iv%satem(nlocal(satem))%thickness(1:nlevels))
493 allocate (iv%satem(nlocal(satem))%org_thickness(1:nlevels))
494
495 iv%satem(nlocal(satem))%p(1) = platform%each(1)%p%inv
496 iv%satem(nlocal(satem))%thickness(1) = platform%each(1)%t
497 ! write original thickness errors for filtered obs
498 if (anal_type_qcobs) then
499 do i = 1, nlevels
500 iv%satem(nlocal(satem))%org_thickness(i) = platform%each(i)%t
501 end do
502 end if
503
504 ! Splitting the reported satem data into smaller layers.
505
506 do i = 2, nlevels
507 iv%satem(nlocal(satem))%p(i) = platform%each(i)%p%inv
508 iv%satem(nlocal(satem))%thickness(i) = platform%each(i)%t
509 if (platform%each(i)%t%qc /= missing_data .and. &
510 platform%each(i-1)%t%qc /= missing_data) then
511 iv%satem(nlocal(satem))%thickness(i)%inv = &
512 platform%each(i)%t%inv - platform%each(i-1)%t%inv
513 else
514 iv%satem(nlocal(satem))%thickness(i)%inv = missing_r
515 iv%satem(nlocal(satem))%thickness(i)%qc = missing_data
516 end if
517 end do
518
519 ! Thickness error (m):
520
521 do i = nlevels, 2, -1
522 iv%satem(nlocal(satem))%thickness(i)%error = &
523 sqrt(iv%satem(nlocal(satem))%thickness(i )%error ** 2 + &
524 iv%satem(nlocal(satem))%thickness(i-1)%error ** 2) / gravity
525 end do
526
527 iv%satem(nlocal(satem))%thickness(1)%error = &
528 sqrt(iv%satem(nlocal(satem))%thickness(1)%error ** 2 + &
529 platform%loc%pw%error ** 2) / gravity
530
531 ! Geostationary ot Polar orbitting Satellite AMVs:
532
533 case (88) ;
534 if (index(platform%info%name, 'MODIS') > 0 .or. &
535 index(platform%info%name, 'modis') > 0) then
536 if (.not. use_polaramvobs) cycle reports
537 if (n==1) ntotal(polaramv) = ntotal(polaramv) + 1
538 if (outside) cycle reports
539 nlocal(polaramv) = nlocal(polaramv) + 1
540 if (nlocal(polaramv) > iv%info(polaramv)%nlocal) cycle reports
541
542 allocate (iv%polaramv(nlocal(polaramv))%p (1:nlevels))
543 allocate (iv%polaramv(nlocal(polaramv))%u (1:nlevels))
544 allocate (iv%polaramv(nlocal(polaramv))%v (1:nlevels))
545
546 do i = 1, nlevels
547 iv%polaramv(nlocal(polaramv))%p(i) = platform%each(i)%p%inv
548 iv%polaramv(nlocal(polaramv))%u(i) = platform%each(i)%u
549 iv%polaramv(nlocal(polaramv))%v(i) = platform%each(i)%v
550 end do
551 obs_index = polaramv ! geoamv is the fm_index value for 88
552 else
553 if (.not. use_geoamvobs) cycle reports
554 if (n==1) ntotal(geoamv) = ntotal(geoamv) + 1
555 if (outside) cycle reports
556 nlocal(geoamv) = nlocal(geoamv) + 1
557 if (nlocal(geoamv) > iv%info(geoamv)%nlocal) cycle reports
558
559 allocate (iv%geoamv(nlocal(geoamv))%p (1:nlevels))
560 allocate (iv%geoamv(nlocal(geoamv))%u (1:nlevels))
561 allocate (iv%geoamv(nlocal(geoamv))%v (1:nlevels))
562
563 do i = 1, nlevels
564 iv%geoamv(nlocal(geoamv))%p(i) = platform%each(i)%p%inv
565 iv%geoamv(nlocal(geoamv))%u(i) = platform%each(i)%u
566 iv%geoamv(nlocal(geoamv))%v(i) = platform%each(i)%v
567 end do
568 end if
569
570 case (42,96:97) ;
571 if (.not. use_airepobs) cycle reports
572 if (n==1) ntotal(airep) = ntotal(airep) + 1
573 if (outside) cycle reports
574 nlocal(airep) = nlocal(airep) + 1
575 if (nlocal(airep) > iv%info(airep)%nlocal) cycle reports
576
577 allocate (iv%airep(nlocal(airep))%h (1:nlevels))
578 allocate (iv%airep(nlocal(airep))%p (1:nlevels))
579 allocate (iv%airep(nlocal(airep))%u (1:nlevels))
580 allocate (iv%airep(nlocal(airep))%v (1:nlevels))
581 allocate (iv%airep(nlocal(airep))%t (1:nlevels))
582
583 do i = 1, nlevels
584 iv%airep(nlocal(airep))%h(i) = platform%each(i)%height
585 iv%airep(nlocal(airep))%p(i) = platform%each(i)%p%inv
586 iv%airep(nlocal(airep))%u(i) = platform%each(i)%u
587 iv%airep(nlocal(airep))%v(i) = platform%each(i)%v
588 iv%airep(nlocal(airep))%t(i) = platform%each(i)%t
589 end do
590
591 case (111, 114) ;
592 if (.not. use_gpspwobs) cycle reports
593 if (n==1) ntotal(gpspw) = ntotal(gpspw) + 1
594 if (outside) cycle reports
595 nlocal(gpspw) = nlocal(gpspw) + 1
596 if (nlocal(gpspw) > iv%info(gpspw)%nlocal) cycle reports
597
598 iv%gpspw(nlocal(gpspw))%tpw = platform%loc%pw
599
600 case (116) ;
601 if (.not. use_gpsrefobs) cycle reports
602 if (n==1) ntotal(gpsref) = ntotal(gpsref) + 1
603 if (outside) cycle reports
604 nlocal(gpsref) = nlocal(gpsref) + 1
605 if (nlocal(gpsref) > iv%info(gpsref)%nlocal) cycle reports
606
607 allocate (iv%gpsref(nlocal(gpsref))%h (1:nlevels))
608 allocate (iv%gpsref(nlocal(gpsref))%ref(1:nlevels))
609 allocate (iv%gpsref(nlocal(gpsref))% p(1:nlevels))
610 allocate (iv%gpsref(nlocal(gpsref))% t(1:nlevels))
611 allocate (iv%gpsref(nlocal(gpsref))% q(1:nlevels))
612
613 do i = 1, nlevels
614 iv%gpsref(nlocal(gpsref))%h(i) = platform%each(i)%height
615
616 ! In OBSPROC, use "td" field to store "gpsref"
617 iv%gpsref(nlocal(gpsref))%ref(i) = platform%each(i)%td
618 ! Keep the retrieved p and t (and q) as "field_type":
619 iv%gpsref(nlocal(gpsref))%p(i) = platform%each(i)%p
620 iv%gpsref(nlocal(gpsref))%t(i) = platform%each(i)%t
621 iv%gpsref(nlocal(gpsref))%q(i) = platform%each(i)%q
622 end do
623
624 case (121) ;
625 if (.not. use_ssmt1obs) cycle reports
626 if (n==1) ntotal(ssmt1) = ntotal(ssmt1) + 1
627 if (outside) cycle reports
628 nlocal(ssmt1) = nlocal(ssmt1) + 1
629 if (nlocal(ssmt1) > iv%info(ssmt2)%nlocal) cycle reports
630
631 allocate (iv%ssmt1(nlocal(ssmt1))%h (1:nlevels))
632 allocate (iv%ssmt1(nlocal(ssmt1))%p (1:nlevels))
633 allocate (iv%ssmt1(nlocal(ssmt1))%t (1:nlevels))
634
635 do i = 1, nlevels
636 iv%ssmt1(nlocal(ssmt1))%h(i) = platform%each(i)%height
637 iv%ssmt1(nlocal(ssmt1))%p(i) = platform%each(i)%p%inv
638 iv%ssmt1(nlocal(ssmt1))%t(i) = platform%each(i)%t
639 end do
640
641 case (122) ;
642 if (.not. use_ssmt2obs) cycle reports
643 if (n==1) ntotal(ssmt2) = ntotal(ssmt2) + 1
644 if (outside) cycle reports
645 nlocal(ssmt2) = nlocal(ssmt2) + 1
646 if (nlocal(ssmt2) > iv%info(ssmt2)%nlocal) cycle reports
647
648 allocate (iv%ssmt2(nlocal(ssmt2))%h (1:nlevels))
649 allocate (iv%ssmt2(nlocal(ssmt2))%p (1:nlevels))
650 allocate (iv%ssmt2(nlocal(ssmt2))%rh(1:nlevels))
651
652 do i = 1, nlevels
653 iv%ssmt2(nlocal(ssmt2))% h(i) = platform%each(i)%height
654 iv%ssmt2(nlocal(ssmt2))% p(i) = platform%each(i)%p%inv
655 iv%ssmt2(nlocal(ssmt2))%rh(i) = platform%each(i)%rh
656 end do
657
658 case (281) ;
659 if (.not. use_qscatobs) cycle reports
660 if (n==1) ntotal(qscat) = ntotal(qscat) + 1
661 if (outside) cycle reports
662 nlocal(qscat) = nlocal(qscat) + 1
663 if (nlocal(qscat) > iv%info(qscat)%nlocal) cycle reports
664
665 iv%qscat(nlocal(qscat))%h = platform%each(1)%height
666 iv%qscat(nlocal(qscat))%u = platform%each(1)%u
667 iv%qscat(nlocal(qscat))%v = platform%each(1)%v
668
669 ! Impose minimum observation error = 1.0m/s for Quikscat data:
670 iv%qscat(nlocal(qscat))%u%error = max(platform%each(1)%u%error,1.0)
671 iv%qscat(nlocal(qscat))%v%error = max(platform%each(1)%v%error,1.0)
672
673 case (132) ; ! profiler
674 if (.not. use_profilerobs) cycle reports
675 if (n==1) ntotal(profiler) = ntotal(profiler) + 1
676 if (outside) cycle reports
677 nlocal(profiler) = nlocal(profiler) + 1
678 if (nlocal(profiler) > iv%info(profiler)%nlocal) cycle reports
679
680 allocate (iv%profiler(nlocal(profiler))%p (1:nlevels))
681 allocate (iv%profiler(nlocal(profiler))%u (1:nlevels))
682 allocate (iv%profiler(nlocal(profiler))%v (1:nlevels))
683
684 do i = 1, nlevels
685 iv%profiler(nlocal(profiler))%p(i) = platform%each(i)%p%inv
686 iv%profiler(nlocal(profiler))%u(i) = platform%each(i)%u
687 iv%profiler(nlocal(profiler))%v(i) = platform%each(i)%v
688 end do
689
690 case (135) ; ! Bogus
691 if (.not. use_bogusobs) cycle reports
692 if (n==1) ntotal(bogus) = ntotal(bogus) + 1
693 if (outside) cycle reports
694 nlocal(bogus) = nlocal(bogus) + 1
695 if (nlocal(bogus) > iv%info(bogus)%nlocal) cycle reports
696
697 if (nlocal(bogus) > max_bogus_input) then
698 write(unit=message(1),fmt='(A,I6,A,I6)') &
699 'Bogus #=', nlocal(bogus), ' > max_bogus_input=', max_bogus_input
700 call da_error(__FILE__,__LINE__,message(1:1))
701 end if
702
703 allocate (iv%bogus(nlocal(bogus))%h (1:nlevels))
704 allocate (iv%bogus(nlocal(bogus))%p (1:nlevels))
705 allocate (iv%bogus(nlocal(bogus))%u (1:nlevels))
706 allocate (iv%bogus(nlocal(bogus))%v (1:nlevels))
707 allocate (iv%bogus(nlocal(bogus))%t (1:nlevels))
708 allocate (iv%bogus(nlocal(bogus))%q (1:nlevels))
709
710 do i = 1, nlevels
711 iv%bogus(nlocal(bogus))%h(i) = platform%each(i)%height
712 iv%bogus(nlocal(bogus))%p(i) = platform%each(i)%p%inv
713 iv%bogus(nlocal(bogus))%u(i) = platform%each(i)%u
714 iv%bogus(nlocal(bogus))%v(i) = platform%each(i)%v
715 iv%bogus(nlocal(bogus))%t(i) = platform%each(i)%t
716 iv%bogus(nlocal(bogus))%q(i) = platform%each(i)%q
717 end do
718
719 iv%bogus(nlocal(bogus))%slp = platform%loc%slp
720
721 if (print_detail_obs) then
722 write(unit=stdout,fmt=*) 'nlevels=', nlevels
723 write(unit=stdout,fmt=*) 'iv%info(bogus)%nlocal,slp', iv%info(bogus)%nlocal, &
724 iv % bogus (nlocal(bogus)) % slp
725 do i=1,nlevels
726 write(unit=stdout,fmt=*) 'nlocal(bogus), i ', nlocal(bogus),i
727 write(unit=stdout,fmt=*) 'iv%bogus(nlocal(bogus))%u,v=', &
728 iv%bogus(nlocal(bogus))%u(i),iv%bogus(nlocal(bogus))%v(i)
729 write(unit=stdout,fmt=*) 'iv%bogus(nlocal(bogus))%t,q=', &
730 iv%bogus(nlocal(bogus))%t(i),iv%bogus(nlocal(bogus))%q(i)
731 end do
732 write(unit=stdout,fmt='(2(a,i4))') 'nlocal(bogus)=',nlocal(bogus), &
733 'nlevels=',nlevels
734 end if
735
736 case (18,19) ; ! buoy
737 if (.not. use_buoyobs) cycle reports
738 if (n==1) ntotal(buoy) = ntotal(buoy) + 1
739 if (outside) cycle reports
740 nlocal(buoy) = nlocal(buoy) + 1
741 if (nlocal(buoy) > iv%info(buoy)%nlocal) cycle reports
742
743 iv%buoy(nlocal(buoy))%h = platform%each(1)%height
744 iv%buoy(nlocal(buoy))%u = platform%each(1)%u
745 iv%buoy(nlocal(buoy))%v = platform%each(1)%v
746 iv%buoy(nlocal(buoy))%t = platform%each(1)%t
747 iv%buoy(nlocal(buoy))%p = platform%each(1)%p
748 iv%buoy(nlocal(buoy))%q = platform%each(1)%q
749
750 case (133) ; ! AIRS retrievals
751 if (.not. use_airsretobs) cycle reports
752 if (n==1) ntotal(airsr) = ntotal(airsr) + 1
753 if (outside) cycle reports
754 nlocal(airsr) = nlocal(airsr) + 1
755 if (nlocal(airsr) > iv%info(airsr)%nlocal) cycle reports
756
757 allocate (iv%airsr(nlocal(airsr))%h (1:nlevels))
758 allocate (iv%airsr(nlocal(airsr))%p (1:nlevels))
759 allocate (iv%airsr(nlocal(airsr))%t (1:nlevels))
760 allocate (iv%airsr(nlocal(airsr))%q (1:nlevels))
761 do i = 1, nlevels
762 iv%airsr(nlocal(airsr))%h(i) = platform%each(i)%height
763 iv%airsr(nlocal(airsr))%p(i) = platform%each(i)%p%inv
764 iv%airsr(nlocal(airsr))%t(i) = platform%each(i)%t
765 iv%airsr(nlocal(airsr))%q(i) = platform%each(i)%q
766 end do
767
768 case default;
769
770 write(unit=message(1), fmt='(a)') 'unsaved obs found:'
771 write(unit=message(2), fmt='(2a)') &
772 'platform%info%platform=', platform%info%platform
773 write(unit=message(3), fmt='(a, i3)') &
774 'platform%info%levels=', platform%info%levels
775 call da_warning(__FILE__,__LINE__,message(1:3))
776
777 end select
778
779 iv%info(obs_index)%name(nlocal(obs_index)) = platform%info%name
780 iv%info(obs_index)%platform(nlocal(obs_index)) = platform%info%platform
781 iv%info(obs_index)%id(nlocal(obs_index)) = platform%info%id
782 iv%info(obs_index)%date_char(nlocal(obs_index)) = platform%info%date_char
783 ! nlevels adjusted for some obs types so use that
784 iv%info(obs_index)%levels(nlocal(obs_index)) = nlevels
785 iv%info(obs_index)%lat(:,nlocal(obs_index)) = platform%info%lat
786 iv%info(obs_index)%lon(:,nlocal(obs_index)) = platform%info%lon
787 iv%info(obs_index)%elv(nlocal(obs_index)) = platform%info%elv
788 iv%info(obs_index)%pstar(nlocal(obs_index)) = platform%info%pstar
789
790 iv%info(obs_index)%slp(nlocal(obs_index)) = platform%loc%slp
791 iv%info(obs_index)%pw(nlocal(obs_index)) = platform%loc%pw
792 iv%info(obs_index)%x(:,nlocal(obs_index)) = platform%loc%x
793 iv%info(obs_index)%y(:,nlocal(obs_index)) = platform%loc%y
794 iv%info(obs_index)%i(:,nlocal(obs_index)) = platform%loc%i
795 iv%info(obs_index)%j(:,nlocal(obs_index)) = platform%loc%j
796 iv%info(obs_index)%dx(:,nlocal(obs_index)) = platform%loc%dx
797 iv%info(obs_index)%dxm(:,nlocal(obs_index)) = platform%loc%dxm
798 iv%info(obs_index)%dy(:,nlocal(obs_index)) = platform%loc%dy
799 iv%info(obs_index)%dym(:,nlocal(obs_index)) = platform%loc%dym
800 iv%info(obs_index)%proc_domain(:,nlocal(obs_index)) = platform%loc%proc_domain
801
802 iv%info(obs_index)%obs_global_index(nlocal(obs_index)) = ntotal(obs_index)
803
804 ! special case for sonde_sfc, duplicate sound info
805 if (obs_index == sound) then
806 iv%info(sonde_sfc)%name(nlocal(sonde_sfc)) = platform%info%name
807 iv%info(sonde_sfc)%platform(nlocal(sonde_sfc)) = platform%info%platform
808 iv%info(sonde_sfc)%id(nlocal(sonde_sfc)) = platform%info%id
809 iv%info(sonde_sfc)%date_char(nlocal(sonde_sfc)) = platform%info%date_char
810 iv%info(sonde_sfc)%levels(nlocal(sonde_sfc)) = 1
811 iv%info(sonde_sfc)%lat(:,nlocal(sonde_sfc)) = platform%info%lat
812 iv%info(sonde_sfc)%lon(:,nlocal(sonde_sfc)) = platform%info%lon
813 iv%info(sonde_sfc)%elv(nlocal(sonde_sfc)) = platform%info%elv
814 iv%info(sonde_sfc)%pstar(nlocal(sonde_sfc)) = platform%info%pstar
815
816 iv%info(sonde_sfc)%slp(nlocal(sonde_sfc)) = platform%loc%slp
817 iv%info(sonde_sfc)%pw(nlocal(sonde_sfc)) = platform%loc%pw
818 iv%info(sonde_sfc)%x(:,nlocal(sonde_sfc)) = platform%loc%x
819 iv%info(sonde_sfc)%y(:,nlocal(sonde_sfc)) = platform%loc%y
820 iv%info(sonde_sfc)%i(:,nlocal(sonde_sfc)) = platform%loc%i
821 iv%info(sonde_sfc)%j(:,nlocal(sonde_sfc)) = platform%loc%j
822 iv%info(sonde_sfc)%dx(:,nlocal(sonde_sfc)) = platform%loc%dx
823 iv%info(sonde_sfc)%dxm(:,nlocal(sonde_sfc)) = platform%loc%dxm
824 iv%info(sonde_sfc)%dy(:,nlocal(sonde_sfc)) = platform%loc%dy
825 iv%info(sonde_sfc)%dym(:,nlocal(sonde_sfc)) = platform%loc%dym
826 iv%info(sonde_sfc)%proc_domain(:,nlocal(sonde_sfc)) = platform%loc%proc_domain
827
828 iv%info(sonde_sfc)%obs_global_index(nlocal(sonde_sfc)) = ntotal(sonde_sfc)
829 end if
830
831 if (global .and. n < 2) then
832 if (test_wrfvar) exit dup_loop
833 if (platform%loc % i >= ide) then
834 platform%loc%i = platform%loc % i - ide
835 else if (platform%loc % i < ids) then
836 platform%loc%i = platform%loc % i + ide
837 end if
838
839 platform%loc%proc_domain = .not. platform%loc%proc_domain
840 end if
841 end do dup_loop
842 end do reports
843
844 close(iunit)
845 call da_free_unit(iunit)
846
847 if (trace_use) call da_trace_exit("da_read_obs_ascii")
848
849 end subroutine da_read_obs_ascii
850
851