da_scan_obs_ascii.inc
References to this file elsewhere.
1 subroutine da_scan_obs_ascii (iv, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Scan WRFVAR GTS observation file
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (iv_type), intent(inout) :: iv
10 character(*), intent(in) :: filename
11
12 character (len = 10) :: fmt_name
13 character (len = 160) :: info_string
14 character (len = 160) :: fmt_info
15 character (len = 160) :: fmt_loc
16 character (len = 160) :: fmt_each
17
18 integer :: i, iost, fm, report, iunit
19 type (multi_level_type) :: platform
20 logical :: outside, outside_all
21 real :: height_error
22 integer :: ndup, n, obs_index
23
24 if (trace_use) call da_trace_entry("da_scan_obs_ascii")
25
26 ! open file
27 ! ---------
28 call da_get_unit(iunit)
29 open(unit = iunit, &
30 FILE = trim(filename), &
31 FORM = 'FORMATTED', &
32 ACCESS = 'SEQUENTIAL', &
33 iostat = iost, &
34 STATUS = 'OLD')
35
36 if (iost /= 0) then
37 write(unit=message(1),fmt='(A,I5,A)') &
38 "Error",iost," opening gts obs file "//trim(filename)
39 call da_warning(__FILE__,__LINE__,message(1:1))
40 call da_free_unit(iunit)
41 if (trace_use) call da_trace_exit("da_scan_obs_ascii")
42 return
43 end if
44
45 ! read header
46
47 head_info: do
48 read (unit = iunit, fmt = '(A)', iostat = iost) info_string
49 if (iost /= 0) then
50 write(unit=message(1),fmt='(A,I3,A,I3)') &
51 "Error",iost,"reading gts obs header on unit",iunit
52 call da_warning(__FILE__,__LINE__,message(1:1))
53 if (trace_use) call da_trace_exit("da_scan_obs_ascii")
54 return
55 end if
56 if (info_string(1:6) == 'EACH ') exit
57 end do head_info
58
59 ! read formats
60
61 read (iunit, fmt = '(A,1X,A)', iostat = iost) &
62 fmt_name, fmt_info, &
63 fmt_name, fmt_loc, &
64 fmt_name, fmt_each
65
66 if (iost /= 0) then
67 write(unit=message(1),fmt='(A,I3,A,I3)') &
68 "Error",iost,"reading gts obs formats on unit",iunit
69 call da_warning(__FILE__,__LINE__,message(1:1))
70 if (trace_use) call da_trace_exit("da_scan_obs_ascii")
71 return
72 end if
73
74 ! skip units
75 read (iunit, fmt = '(A)') fmt_name
76
77 ! loop over records
78
79 report = 0 ! report number in file
80
81 reports: do
82 report = report+1
83
84 ! read station general info
85
86 read (iunit, fmt = fmt_info, iostat = iost) &
87 platform%info%platform, &
88 platform%info%date_char, &
89 platform%info%name, &
90 platform%info%levels, &
91 platform%info%lat, &
92 platform%info%lon, &
93 platform%info%elv, &
94 platform%info%id
95
96 if (iost /= 0) then
97 ! FIX? This is expected, but its unclear how we handle failure
98 ! here without assuming the fortran2003 convention on
99 ! error statuses
100 exit reports
101 end if
102
103 if (print_detail_obs) then
104 write(unit=stdout, fmt = fmt_info) &
105 platform%info%platform, &
106 platform%info%date_char, &
107 platform%info%name, &
108 platform%info%levels, &
109 platform%info%lat, &
110 platform%info%lon, &
111 platform%info%elv, &
112 platform%info%id
113 end if
114
115 if (platform%info%lon == 180.0) platform%info%lon =-180.000
116 ! WHY?
117 ! Fix funny wind direction at South Poles
118 ! if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
119 ! platform%info%lon = 0.0
120 ! end if
121
122 read (platform%info%platform(4:6), '(I3)') fm
123
124 ! read model location
125 read (iunit, fmt = fmt_loc) platform%loc%slp, platform%loc%pw
126
127 ! levels < 1 and not GPSPW, go back to reports
128
129 if ((platform%info%levels < 1) .AND. &
130 (index(platform%info%platform, 'GPSPW') <= 0)) then
131 cycle reports
132 end if
133
134 ! read each level
135
136 do i = 1, platform%info%levels
137 platform%each (i) = each_level_type(missing_r, missing, -1.0, & ! height
138 field_type(missing_r, missing, missing_r), & ! u
139 field_type(missing_r, missing, missing_r), & ! v
140 field_type(missing_r, missing, missing_r), & ! p
141 field_type(missing_r, missing, missing_r), & ! t
142 field_type(missing_r, missing, missing_r), & ! q
143 field_type(missing_r, missing, missing_r), & ! rh
144 field_type(missing_r, missing, missing_r), & ! td
145 field_type(missing_r, missing, missing_r)) ! speed
146
147 read (unit = iunit, fmt = trim (fmt_each)) &
148 platform%each (i)%p, &
149 platform%each (i)%speed, &
150 ! Here the 'direction' is stored in platform%each (i)%v:
151 platform%each (i)%v, &
152 platform%each (i)%height, &
153 platform%each (i)%height_qc, &
154 height_error, &
155 platform%each (i)%t, &
156 platform%each (i)%td, &
157 platform%each (i)%rh
158 end do
159
160 ! Restrict to a range of reports, useful for debugging
161
162 if (report < report_start) cycle
163 if (report > report_end) exit
164
165 call da_llxy (platform%info, platform%loc, outside, outside_all)
166
167 if (platform%info%levels < 1) then
168 if (fm /= 111) then
169 cycle reports
170 end if
171 end if
172
173 ! Loop over duplicating obs for global
174 ndup = 1
175 if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
176
177 if (test_wrfvar) ndup = 1
178
179 obs_index = fm_index(fm)
180
181 do n = 1, ndup
182 select case(fm)
183
184 case (12) ;
185 if (.not.use_synopobs .or. iv%info(synop)%ntotal == max_synop_input) cycle reports
186 if (n==1) iv%info(synop)%ntotal = iv%info(synop)%ntotal + 1
187 if (outside) cycle reports
188 iv%info(synop)%nlocal = iv%info(synop)%nlocal + 1
189
190 case (13, 17) ;
191 if (.not.use_shipsobs .or. iv%info(ships)%ntotal == max_ships_input) cycle reports
192 if (n==1) iv%info(ships)%ntotal = iv%info(ships)%ntotal + 1
193 if (outside) cycle reports
194 iv%info(ships)%nlocal = iv%info(ships)%nlocal + 1
195
196 case (15:16) ;
197 if (.not.use_metarobs .or. iv%info(metar)%ntotal == max_metar_input) cycle reports
198 if (n==1) iv%info(metar)%ntotal = iv%info(metar)%ntotal + 1
199 if (outside) cycle reports
200 iv%info(metar)%nlocal = iv%info(metar)%nlocal + 1
201
202 case (32:34) ;
203 if (.not.use_pilotobs .or. iv%info(pilot)%ntotal == max_pilot_input) cycle reports
204 if (n==1) iv%info(pilot)%ntotal = iv%info(pilot)%ntotal + 1
205 if (outside) cycle reports
206 iv%info(pilot)%nlocal = iv%info(pilot)%nlocal + 1
207
208 case (35:38) ;
209 if (.not.use_soundobs .or. iv%info(sound)%ntotal == max_sound_input) cycle reports
210 if (n==1) iv%info(sound)%ntotal = iv%info(sound)%ntotal + 1
211 if (n==1) iv%info(sonde_sfc)%ntotal = iv%info(sonde_sfc)%ntotal + 1
212 if (outside) cycle reports
213 iv%info(sound)%nlocal = iv%info(sound)%nlocal + 1
214 iv%info(sonde_sfc)%nlocal = iv%info(sonde_sfc)%nlocal + 1
215
216 case (86) ;
217 if (.not.use_satemobs .or. iv%info(satem)%ntotal == max_satem_input) cycle reports
218 ! Reject cloudy satem obs.
219 if (platform%loc%pw%inv > 10.0) then
220 cycle reports
221 end if
222 if (n==1) iv%info(satem)%ntotal = iv%info(satem)%ntotal + 1
223 if (outside) cycle reports
224 iv%info(satem)%nlocal = iv%info(satem)%nlocal + 1
225
226 case (88) ;
227 ! Geostationary or Polar orbitting Satellite AMVs:
228 if (index(platform%info%name, 'MODIS') > 0 .or. &
229 index(platform%info%name, 'modis') > 0) then
230 if (.not.use_polaramvobs .or. iv%info(polaramv)%ntotal == max_polaramv_input) cycle reports
231 if (n==1) iv%info(polaramv)%ntotal = iv%info(polaramv)%ntotal + 1
232 if (outside) cycle reports
233 iv%info(polaramv)%nlocal = iv%info(polaramv)%nlocal + 1
234 obs_index = polaramv ! geoamv is the fm_index value for 88
235 else
236 if (.not.use_geoamvobs .or. iv%info(geoamv)%ntotal == max_geoamv_input) cycle reports
237 if (n==1) iv%info(geoamv)%ntotal = iv%info(geoamv)%ntotal + 1
238 if (outside) cycle reports
239 iv%info(geoamv)%nlocal = iv%info(geoamv)%nlocal + 1
240 end if
241
242 case (42,96:97) ;
243 if (.not.use_airepobs .or. iv%info(airep)%ntotal == max_airep_input) cycle reports
244 if (n==1) iv%info(airep)%ntotal = iv%info(airep)%ntotal + 1
245 if (outside) cycle reports
246 iv%info(airep)%nlocal = iv%info(airep)%nlocal + 1
247
248 case (111) ;
249 if (.not.use_gpspwobs .or. iv%info(gpspw)%ntotal == max_gpspw_input) cycle reports
250 if (n==1) iv%info(gpspw)%ntotal = iv%info(gpspw)%ntotal + 1
251 if (outside) cycle reports
252 iv%info(gpspw)%nlocal = iv%info(gpspw)%nlocal + 1
253
254 case (116) ;
255 if (.not.use_gpsrefobs .or. iv%info(gpsref)%ntotal == max_gpsref_input) cycle reports
256 if (n==1) iv%info(gpsref)%ntotal = iv%info(gpsref)%ntotal + 1
257 if (outside) cycle reports
258 iv%info(gpsref)%nlocal = iv%info(gpsref)%nlocal + 1
259
260 case (121) ;
261 ! SSM/T1 temperatures
262 if (.not.use_ssmt1obs .or. iv%info(ssmt1)%ntotal == max_ssmt1_input) cycle reports
263 if (n==1) iv%info(ssmt1)%ntotal = iv%info(ssmt1)%ntotal + 1
264 if (outside) cycle reports
265 iv%info(ssmt2)%nlocal = iv%info(ssmt2)%nlocal + 1
266
267 case (122) ;
268 ! SSM/T2 relative humidities:
269 if (.not.use_ssmt2obs .or. iv%info(ssmt2)%ntotal == max_ssmt2_input) cycle reports
270 if (n==1) iv%info(ssmt2)%ntotal = iv%info(ssmt2)%ntotal + 1
271 if (outside) cycle reports
272 iv%info(ssmt2)%nlocal = iv%info(ssmt2)%nlocal + 1
273
274 case (281) ;
275 ! Scatterometer:
276 if (.not.use_qscatobs .or. iv%info(qscat)%ntotal == max_qscat_input) cycle reports
277 if (n==1) iv%info(qscat)%ntotal = iv%info(qscat)%ntotal + 1
278 if (outside) cycle reports
279 iv%info(qscat)%nlocal = iv%info(qscat)%nlocal + 1
280
281 case (132) ;
282 if (.not.use_profilerobs .or. iv%info(profiler)%ntotal == max_profiler_input) cycle reports
283 if (n==1) iv%info(profiler)%ntotal = iv%info(profiler)%ntotal + 1
284 if (outside) cycle reports
285 iv%info(profiler)%nlocal = iv%info(profiler)%nlocal + 1
286
287 case (135) ;
288 if (.not.use_bogusobs .or. iv%info(bogus)%ntotal == max_bogus_input) cycle reports
289 if (n==1) iv%info(bogus)%ntotal = iv%info(bogus)%ntotal + 1
290 if (outside) cycle reports
291 iv%info(bogus)%nlocal = iv%info(bogus)%nlocal + 1
292
293 case (18,19) ; ! buoy
294 if (.not.use_buoyobs .or. iv%info(buoy)%ntotal == max_buoy_input) cycle reports
295 if (n==1) iv%info(buoy)%ntotal = iv%info(buoy)%ntotal + 1
296 if (outside) cycle reports
297 iv%info(buoy)%nlocal = iv%info(buoy)%nlocal + 1
298
299 case (133) ; ! AIRS retrievals
300 if (.not.use_airsretobs .or. iv%info(airsr)%ntotal == max_airsr_input) cycle reports
301 if (n==1) iv%info(airsr)%ntotal = iv%info(airsr)%ntotal + 1
302 if (outside) cycle reports
303 iv%info(airsr)%nlocal = iv%info(airsr)%nlocal + 1
304
305 case default;
306 write(unit=message(1), fmt='(a)') 'unsaved obs found:'
307 write(unit=message(2), fmt='(2a)') &
308 'platform%info%platform=', platform%info%platform
309 write(unit=message(3), fmt='(a, i3)') &
310 'platform%info%levels=', platform%info%levels
311 call da_warning(__FILE__,__LINE__,message(1:3))
312 end select
313 iv%info(obs_index)%max_lev = max(iv%info(obs_index)%max_lev, platform%info%levels)
314 end do ! loop over duplicate
315 end do reports
316
317 iv%info(sonde_sfc)%max_lev=1
318
319 close(iunit)
320
321 call da_free_unit(iunit)
322
323 if (trace_use) call da_trace_exit("da_scan_obs_ascii")
324
325 end subroutine da_scan_obs_ascii
326
327