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