da_scan_obs.inc
References to this file elsewhere.
1 subroutine da_scan_obs (ob, xp, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Scan WRFVAR GTS observation file
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (xpose_type), intent(in) :: xp ! Domain decomposition vars.
10 type (ob_type), intent(inout) :: ob
11
12 character(*), intent (in) :: filename
13
14 character (LEN = 10) :: fmt_name
15
16 character (LEN = 160) :: info_string
17
18 character (LEN = 160) :: fmt_info, &
19 fmt_loc, &
20 fmt_each
21
22 integer :: i, iost, fm, report,iunit
23
24 type (multi_level_type) :: platform
25
26 logical :: outside, outside_all
27
28 real :: height_error
29
30
31 integer :: total_obs, &
32 num_sound, &
33 num_synop, &
34 num_pilot, &
35 num_geoamv, &
36 num_polaramv, &
37 num_satem, &
38 num_airep, &
39 num_metar, &
40 num_ships, &
41 num_gpspw, &
42 num_gpsref, &
43 num_ssmi_retrieval, &
44 num_ssmi_tb , &
45 num_ssmt1, &
46 num_ssmt2, num_airsr, &
47 num_qscat, &
48 num_profiler, num_buoy, num_bogus
49
50 integer :: ndup, n
51
52 ! default value
53
54 ob%ob_numb(ob%current_ob_time)%sound = ob%num_sound
55 ob%ob_numb(ob%current_ob_time)%synop = ob%num_synop
56 ob%ob_numb(ob%current_ob_time)%pilot = ob%num_pilot
57 ob%ob_numb(ob%current_ob_time)%satem = ob%num_satem
58 ob%ob_numb(ob%current_ob_time)%geoamv = ob%num_geoamv
59 ob%ob_numb(ob%current_ob_time)%polaramv = ob%num_polaramv
60 ob%ob_numb(ob%current_ob_time)%airep = ob%num_airep
61 ob%ob_numb(ob%current_ob_time)%gpspw = ob%num_gpspw
62 ob%ob_numb(ob%current_ob_time)%gpsref = ob%num_gpsref
63 ob%ob_numb(ob%current_ob_time)%metar = ob%num_metar
64 ob%ob_numb(ob%current_ob_time)%ships = ob%num_ships
65 ob%ob_numb(ob%current_ob_time)%qscat = ob%num_qscat
66 ob%ob_numb(ob%current_ob_time)%buoy = ob%num_buoy
67 ob%ob_numb(ob%current_ob_time)%profiler = ob%num_profiler
68
69 ob%ob_numb(ob%current_ob_time)%ssmt1 = ob%num_ssmt1
70 ob%ob_numb(ob%current_ob_time)%ssmt2 = ob%num_ssmt2
71 ob%ob_numb(ob%current_ob_time)%airsr = ob%num_airsr
72 ! WHY
73 ! ob%ob_numb(ob%current_ob_time)%ssmi_tb = ob%num_ssmi_tb
74 ! ob%ob_numb(ob%current_ob_time)%ssmi_retrieval = ob%num_ssmi_retrieval
75
76 ! open file
77 ! ---------
78 call da_get_unit(iunit)
79 open(unit = iunit, &
80 FILE = trim(filename), &
81 FORM = 'FORMATTED', &
82 ACCESS = 'SEQUENTIAL', &
83 iostat = iost, &
84 STATUS = 'OLD')
85
86 if (iost /= 0) then
87 write(unit=message(1),fmt='(A,I5,A)') &
88 "Error",iost," opening gts obs file "//trim(filename)
89 call da_warning(__FILE__,__LINE__,message(1:1))
90 call da_free_unit(iunit)
91 return
92 end if
93
94 total_obs = 0
95 num_sound = 0
96 num_synop = 0
97 num_pilot = 0
98 num_geoamv = 0
99 num_polaramv = 0
100 num_satem = 0
101 num_airep = 0
102 num_metar = 0
103 num_ships = 0
104 num_gpspw = 0
105 num_gpsref = 0
106 num_ssmi_retrieval = 0
107 num_ssmi_tb = 0
108 num_ssmt1 = 0
109 num_ssmt2 = 0
110 num_qscat = 0
111 num_profiler = 0
112 num_buoy = 0
113 num_bogus = 0
114 num_airsr = 0
115
116 ! read header
117
118 head_info: do
119 read (unit = iunit, fmt = '(A)', iostat = iost) info_string
120 if (iost /= 0) then
121 write(unit=message(1),fmt='(A,I3,A,I3)') &
122 "Error",iost,"reading gts obs header on unit",iunit
123 call da_warning(__FILE__,__LINE__,message(1:1))
124 return
125 end if
126 if (info_string(1:6) == 'EACH ') exit
127 end do head_info
128
129 ! read formats
130
131 read (iunit, fmt = '(A,1X,A)', iostat = iost) &
132 fmt_name, fmt_info, &
133 fmt_name, fmt_loc, &
134 fmt_name, fmt_each
135
136 if (iost /= 0) then
137 write(unit=message(1),fmt='(A,I3,A,I3)') &
138 "Error",iost,"reading gts obs formats on unit",iunit
139 call da_warning(__FILE__,__LINE__,message(1:1))
140 return
141 end if
142
143 ! skip units
144 read (iunit, fmt = '(A)') fmt_name
145
146 ! loop over records
147
148 report = 0 ! report number in file
149
150 reports: &
151 do
152
153 report = report+1
154
155 ! read station general info
156
157 read (iunit, fmt = fmt_info, iostat = iost) &
158 platform%info%platform, &
159 platform%info%date_char, &
160 platform%info%name, &
161 platform%info%levels, &
162 platform%info%lat, &
163 platform%info%lon, &
164 platform%info%elv, &
165 platform%info%id
166
167 if (iost /= 0) then
168 ! JRB This is expected, but its unclear how we handle failure
169 ! here without assuming the fortran2003 convention on
170 ! error statuses
171 exit reports
172 end if
173
174 if (print_detail_obs) then
175 write(unit=stdout, fmt = fmt_info) &
176 platform%info%platform, &
177 platform%info%date_char, &
178 platform%info%name, &
179 platform%info%levels, &
180 platform%info%lat, &
181 platform%info%lon, &
182 platform%info%elv, &
183 platform%info%id
184 end if
185
186 if (platform%info%lon == 180.0) platform%info%lon =-180.000
187 ! WHY
188 ! Fix funny wind direction at South Poles
189 ! if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
190 ! platform%info%lon = 0.0
191 ! end if
192
193 read (platform%info%platform(4:6), '(I3)') fm
194
195 ! read model location
196 read (iunit, fmt = fmt_loc) platform%loc%slp, platform%loc%pw
197
198 total_obs = total_obs + 1
199
200 ! levels < 1 and not GPSPW, go back to reports
201
202 if ((platform%info%levels < 1) .AND. &
203 (index(platform%info%platform, 'GPSPW') <= 0)) then
204 cycle reports
205 end if
206
207 ! read each level
208
209 do i = 1, platform%info%levels
210 platform%each (i) = each_level_type(missing_r, missing, -1.0, & ! height
211 field_type(missing_r, missing, missing_r), & ! u
212 field_type(missing_r, missing, missing_r), & ! v
213 field_type(missing_r, missing, missing_r), & ! p
214 field_type(missing_r, missing, missing_r), & ! t
215 field_type(missing_r, missing, missing_r), & ! q
216 field_type(missing_r, missing, missing_r), & ! rh
217 field_type(missing_r, missing, missing_r), & ! td
218 field_type(missing_r, missing, missing_r)) ! speed
219
220 read (unit = iunit, fmt = trim (fmt_each)) &
221 platform%each (i)%p, &
222 platform%each (i)%speed, &
223 ! Here the 'direction' is stored in platform%each (i)%v:
224 platform%each (i)%v, &
225 platform%each (i)%height, &
226 platform%each (i)%height_qc, &
227 height_error, &
228 platform%each (i)%t, &
229 platform%each (i)%td, &
230 platform%each (i)%rh
231 end do
232
233 ! Restrict to a range of reports, useful for debugging
234
235 if (report < report_start) then
236 cycle
237 end if
238
239 if (report > report_end) then
240 exit
241 end if
242
243 call da_ll_to_xy (platform%info, platform%loc, &
244 xp, outside, outside_all)
245
246 if (outside) then
247 cycle reports
248 end if
249
250 if (platform%info%levels < 1) then
251 if (fm /= 111) then
252 cycle reports
253 end if
254 end if
255
256 ! Loop over duplicating obs for global
257 ndup = 1
258 if (global .and. &
259 (platform%loc%i < xp%ids .or. platform%loc%i >= xp%ide)) ndup= 2
260
261 if (Testing_WRFVAR) ndup = 1
262
263 do n = 1, ndup
264 select case(fm)
265
266 case (12) ;
267 if (.not.use_SynopObs) cycle reports
268 num_synop = num_synop + 1
269
270 case (13, 17) ; ! ships
271 if (.not.use_ShipsObs) cycle reports
272 num_ships = num_ships + 1
273
274 case (15:16) ;
275 if (.not.use_MetarObs) cycle reports
276 num_metar = num_metar + 1
277
278 case (32:34) ;
279 if (.not.use_PilotObs) cycle reports
280 num_pilot = num_pilot + 1
281
282 case (35:38) ;
283 if (.not.use_SoundObs) cycle reports
284 num_sound = num_sound + 1
285
286 case (86) ;
287 if (.not.use_SatemObs) cycle reports
288
289 ! Reject cloudy Satem obs.
290 if (platform%loc%pw%inv > 10.) then
291 cycle reports
292 end if
293
294 num_satem = num_satem + 1
295
296 case (88) ;
297 ! Geostationary or Polar orbitting Satellite AMVs:
298 if (index(platform%info%name, 'MODIS') > 0 .or. &
299 index(platform%info%name, 'modis') > 0) then
300 if (.not.use_PolarAMVObs) cycle reports
301 num_polaramv = num_polaramv + 1
302 else
303 if (.not.use_GeoAMVObs) cycle reports
304 num_geoamv = num_geoamv + 1
305 end if
306
307 case (42,96:97) ;
308 if (.not.use_AirepObs) cycle reports
309 num_airep = num_airep + 1
310
311 case (111) ;
312 if (.not.use_GpspwObs) cycle reports
313 num_gpspw = num_gpspw + 1
314
315 case (116) ;
316 if (.not.use_GpsrefObs) cycle reports
317 num_gpsref = num_gpsref + 1
318
319 case (121) ;
320 ! SSM/T1 temperatures
321 if (.not.use_ssmt1obs) cycle reports
322 num_ssmt1 = num_ssmt1 + 1
323
324 case (122) ;
325 ! SSM/T2 relative humidities:
326 if (.not.use_ssmt2obs) cycle reports
327 num_ssmt2 = num_ssmt2 + 1
328
329 case (281) ;
330 ! Scatterometer:
331 if (.not.use_qscatobs) cycle reports
332 num_qscat = num_qscat + 1
333
334 case (132) ;
335 if (.not.use_ProfilerObs) cycle reports
336 num_profiler = num_profiler + 1
337
338 case (135) ;
339 if (.not.use_BogusObs) cycle reports
340 num_bogus = num_bogus + 1
341
342 case (18,19) ; ! buoy
343 if (.not.use_BuoyObs) cycle reports
344 num_buoy = num_buoy + 1
345
346 case (133) ; ! AIRS retrievals
347 if (.not.use_AIRSRETObs) cycle reports
348 num_airsr = num_airsr + 1
349
350 case default;
351 write(unit=message(1), fmt='(a)') 'unsaved obs found:'
352 write(unit=message(2), fmt='(2a)') &
353 'platform%info%platform=', platform%info%platform
354 write(unit=message(3), fmt='(a, i3)') &
355 'platform%info%levels=', platform%info%levels
356 call da_warning(__FILE__,__LINE__,message(1:3))
357 end select
358 end do ! loop over duplicate
359 end do reports
360
361 close(iunit)
362 call da_free_unit(iunit)
363
364 !------------------------------------------------------------------------
365 ! Check the numbers again, make sure we have the right number.
366 !------------------------------------------------------------------------
367
368 ob%num_sound = ob%num_sound + num_sound
369 ob%num_synop = ob%num_synop + num_synop
370 ob%num_pilot = ob%num_pilot + num_pilot
371 ob%num_satem = ob%num_satem + num_satem
372 ob%num_geoamv = ob%num_geoamv + num_geoamv
373 ob%num_polaramv = ob%num_polaramv + num_polaramv
374 ob%num_airep = ob%num_airep + num_airep
375 ob%num_gpspw = ob%num_gpspw + num_gpspw
376 ob%num_gpsref = ob%num_gpsref + num_gpsref
377 ob%num_metar = ob%num_metar + num_metar
378 ob%num_ships = ob%num_ships + num_ships
379 ob%num_qscat = ob%num_qscat + num_qscat
380 ob%num_buoy = ob%num_buoy + num_buoy
381 ob%num_profiler = ob%num_profiler + num_profiler
382 ob%num_bogus = ob%num_bogus + num_bogus
383
384 ob%num_ssmt1 = ob%num_ssmt1 + num_ssmt1
385 ob%num_ssmt2 = ob%num_ssmt2 + num_ssmt2
386 ob%num_airsr = ob%num_airsr + num_airsr
387
388 ! ob%num_ssmi_tb = ob%num_ssmi_tb + num_ssmi_tb
389 ! ob%num_ssmi_retrieval = ob%num_ssmi_retrieval + num_ssmi_retrieval
390
391 ob%ob_numb(ob%current_ob_time)%sound = ob%num_sound
392 ob%ob_numb(ob%current_ob_time)%synop = ob%num_synop
393 ob%ob_numb(ob%current_ob_time)%pilot = ob%num_pilot
394 ob%ob_numb(ob%current_ob_time)%satem = ob%num_satem
395 ob%ob_numb(ob%current_ob_time)%geoamv = ob%num_geoamv
396 ob%ob_numb(ob%current_ob_time)%polaramv = ob%num_polaramv
397 ob%ob_numb(ob%current_ob_time)%airep = ob%num_airep
398 ob%ob_numb(ob%current_ob_time)%gpspw = ob%num_gpspw
399 ob%ob_numb(ob%current_ob_time)%gpsref = ob%num_gpsref
400 ob%ob_numb(ob%current_ob_time)%metar = ob%num_metar
401 ob%ob_numb(ob%current_ob_time)%ships = ob%num_ships
402 ob%ob_numb(ob%current_ob_time)%qscat = ob%num_qscat
403 ob%ob_numb(ob%current_ob_time)%buoy = ob%num_buoy
404 ob%ob_numb(ob%current_ob_time)%profiler = ob%num_profiler
405 ob%ob_numb(ob%current_ob_time)%bogus = ob%num_bogus
406
407 ob%ob_numb(ob%current_ob_time)%ssmt1 = ob%num_ssmt1
408 ob%ob_numb(ob%current_ob_time)%ssmt2 = ob%num_ssmt2
409 ob%ob_numb(ob%current_ob_time)%airsr = ob%num_airsr
410
411 ! WHY
412 ! ob%ob_numb(ob%current_ob_time)%ssmi_tb = ob%num_ssmi_tb
413 ! ob%ob_numb(ob%current_ob_time)%ssmi_retrieval = ob%num_ssmi_retrieval
414
415 ! print out
416
417 if (print_detail_obs) then
418 write(unit=stdout, fmt='(a)') ' '
419 write(unit=stdout, fmt='(5x,a,i6,a)') &
420 'Scan: ', ob%current_ob_time, ' ob time:', &
421 'Scan: ', total_obs, ' Total reports.', &
422 'Scan: ', num_sound, ' SOUND reports,', &
423 'Scan: ', num_synop, ' SYNOP reports,', &
424 'Scan: ', num_pilot, ' PILOT reports,', &
425 'Scan: ', num_satem, ' SATEM reports,', &
426 'Scan: ', num_geoamv,' Geo. AMVs rports,', &
427 'Scan: ', num_polaramv,' Polar AMVs rports,', &
428 'Scan: ', num_airep, ' AIREP reports,', &
429 'Scan: ', num_gpspw, ' GPSPW reports,', &
430 'Scan: ', num_gpsref,' GPSRF reports,', &
431 'Scan: ', num_metar, ' METAR reports,', &
432 'Scan: ', num_ships, ' SHIP reports,', &
433 'Scan: ', num_qscat, ' QSCAT reports,', &
434 'Scan: ', num_profiler, ' Profiler reports,', &
435 'Scan: ', num_buoy, ' Buoy reports,', &
436 'Scan: ', num_bogus, ' TCBGS reports,', &
437 'Scan: ', num_ssmt1, ' SSMT1 reports,', &
438 'Scan: ', num_ssmt2, ' SSMT2 reports.', &
439 'Scan: ', num_airsr, ' AIRS retrieval reports'
440
441 ! WHY
442 ! write(unit=stdout, fmt='(5x,a,i6,a)') &
443 ! 'Scan: ', num_ssmi_retrieval , ' SSMI_RETRIEVAL reports,', &
444 ! 'Scan: ', num_ssmi_tb , ' SSMI_TB reports.'
445
446 end if
447
448 end subroutine da_scan_obs
449
450