da_scan_obs_bufr.inc
References to this file elsewhere.
1 subroutine da_scan_obs_bufr (iv, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Read BUFR observation file for input to wrfvar
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (iv_type), intent(inout) :: iv
10 character(len=*), optional, intent(in) :: filename
11
12 #ifdef BUFR
13
14 type (multi_level_type) :: platform
15 logical :: outside, outside_all
16 integer :: nlocal(num_ob_indexes)
17 integer :: ntotal(num_ob_indexes)
18
19 character(len=40) :: obstr,hdstr
20 character(len=8) :: subset
21 real :: hdr(7)
22 real :: pmo(2,1)
23 ! real :: drf9,255)
24 real :: obs(9,255),qms(9,255),oes(9,255)
25 real :: time,woe,qob,toe,qoe,poe,pob,tob,zob,rob,roe
26
27 integer :: iost, ndup, n, i, j, k, surface_level, report
28 integer :: iret, idate, kx, nlevels, nrecs,miscd, t29, jx, ix
29 ! integer :: cat,zqm,pqm,qqm,tqm,wqm,pwq,pmq,rqm
30 integer :: iunit, fm, obs_index
31
32 if (trace_use) call da_trace_entry("da_scan_obs_bufr")
33
34 ! open file
35 ! ---------
36 call da_get_unit(iunit)
37 if (present(filename)) then
38 open(unit = iunit, FILE = trim(filename), &
39 iostat = iost, form = 'unformatted', STATUS = 'OLD')
40 if (iost /= 0) then
41 write(unit=stdout, fmt='(/A,I3,3(2X,A)/)') &
42 'error in obs input file unit ',iunit, &
43 'obs filename:', trim(filename), &
44 'for gts observations cannot be found or cannot be opened'
45 return
46 end if
47 end if
48
49 hdstr='SID XOB YOB DHR TYP ELV T29 '
50 obstr='POB QOB TOB ZOB UOB VOB PWO ROB CAT '
51
52 nrecs = 0
53 miscd = 0
54
55 !--------------------------------
56 ! open bufr file then check date
57 !--------------------------------
58
59 call datelen(10)
60
61 call openbf(iunit,'IN',iunit)
62
63 call readmg(iunit,subset,idate,iret)
64
65 if (iret/=0) then
66 write(unit=stdout, fmt='(a, i4)') &
67 'return code from readmg:', iret, &
68 'Reach the end of obs unit: ', iunit
69
70 call da_error(__FILE__,__LINE__, &
71 (/"No BUFR observations"/))
72 end if
73
74 report = 0 ! report number in file
75 reports: do
76 report = report+1
77
78 call readsb(iunit,iret)
79
80 if (iret/=0) then
81 call readmg(iunit,subset,idate,iret)
82
83 if (iret/=0) then
84 write(unit=stdout, fmt='(a, i4)') &
85 'return code from readmg:', iret, &
86 'Reach the end of prepbufr obs unit: ', iunit
87
88 exit reports
89 end if
90
91 cycle reports
92 end if
93
94 nrecs=nrecs+1
95
96 call ufbint(iunit,hdr,7,1,iret,hdstr)
97
98 platform % info % name(1:8) = subset
99 platform % info % id = hdstr(1:5)
100
101 if (hdr(2) > 180.0) hdr(2)=hdr(2)-360.0
102
103 ! Put a check on Lat
104
105 hdr(3) = MAX(hdr(3), -89.95)
106 hdr(3) = Min(hdr(3), 89.95)
107 platform%info%lon = hdr(2)
108 platform%info%lat = hdr(3)
109
110 ! Put a check on Lat
111
112 if (platform%info%lon == 180.0 ) platform%info%lon =-180.000
113 ! Fix funny wind direction at Poles
114 if (platform%info%lat <= -89.95 .or. platform%info%lat >= 89.95) then
115 platform%info%lon = 0.0
116 end if
117
118 ! Restrict to a range of reports, useful for debugging
119
120 if (report < report_start) cycle
121
122 if (report > report_end) exit
123
124 call da_llxy (platform%info, platform%loc,outside, outside_all)
125
126 if (outside_all) cycle reports
127
128 t29 = int(0.1 + hdr(7))
129 kx=int(0.1+hdr(5))
130
131 if (kx == 183) then ! reset kx based on t29
132 if (t29 .eq. 511) kx = 181
133 if (t29 .eq. 512) kx = 187
134 if (t29 .eq. 522) kx = 180
135 if (t29 .eq. 523) kx = 180
136 if (t29 .eq. 531) kx = 180
137 if (t29 .eq. 561) kx = 180
138 if (t29 .eq. 562) kx = 180
139 end if
140
141 ! WHY?
142 ! if ((kx >= 160) .and. (kx <= 179)) then ! bypass satellite data
143 ! if (t29 .eq. 61 .or. t29 .eq. 63 .or. t29 .ge. 571) then
144 ! cycle reports
145 ! end if
146
147 ! Conventional data
148
149 call ufbint(iunit,obs,9,255,nlevels,obstr)
150 if (nlevels > max_ob_levels) nlevels = max_ob_levels
151
152 if ((nlevels < 1) .and. ((kx /= 164) .or. (kx /= 174))) cycle reports
153
154 !---------------------------------------------------------------------------
155 ! This is basically converting rh to q i
156 ! Method :
157 ! if rh, temp and pr all available computes Qs otherwise sets Qs= missing
158 ! if rh > 100 sets q = qs otherwise q = rh*Qs/100.0
159 ! Note: Currently da_obs_proc_station is active only for ob_format_ascii
160 ! call da_obs_proc_station(platform)
161 !---------------------------------------------------------------------------
162
163 ! Loop over duplicating obs for global
164 ndup = 1
165 if (global .and. &
166 (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
167 if (test_transforms) ndup = 1
168
169 ! It is possible that logic for counting obs is incorrect for the
170 ! global case with >1 MPI tasks due to obs duplication, halo, etc.
171 ! TBH: 20050913
172 dup_loop: do n = 1, ndup
173 select case(t29)
174 case (11, 12, 13, 22, 23, 31)
175 select case (kx)
176 case (120, 122, 132, 220, 222, 232) ; ! Sound
177 if (.not.use_soundobs .or. iv%info(sound)%ntotal == max_sound_input) cycle reports
178 if (n==1) iv%info(sound)%ntotal = iv%info(sound)%ntotal + 1
179 if (n==1) iv%info(sonde_sfc)%ntotal = iv%info(sonde_sfc)%ntotal + 1
180 if (outside) cycle reports
181 iv%info(sound)%nlocal = iv%info(sound)%nlocal + 1
182 iv%info(sonde_sfc)%nlocal = iv%info(sonde_sfc)%nlocal + 1
183 fm = 35
184 case (221) ; ! Pilot
185 if (.not.use_pilotobs .or. iv%info(pilot)%ntotal == max_pilot_input) cycle reports
186 if (n==1) iv%info(pilot)%ntotal = iv%info(pilot)%ntotal + 1
187 if (outside) cycle reports
188 iv%info(pilot)%nlocal = iv%info(pilot)%nlocal + 1
189 fm = 32
190 end select
191
192 case (41)
193 ! case (130:131, 133, 230:231, 233) ; ! Airep
194 if (.not.use_airepobs .or. iv%info(airep)%ntotal == max_airep_input) cycle reports
195 if (n==1) iv%info(airep)%ntotal = iv%info(airep)%ntotal + 1
196 if (outside) cycle reports
197 iv%info(airep)%nlocal = iv%info(airep)%nlocal + 1
198 fm = 42
199
200 case (522, 523); ! Ships
201 if (.not.use_shipsobs .or. iv%info(ships)%ntotal == max_ships_input) cycle reports
202 if (n==1) iv%info(ships)%ntotal = iv%info(ships)%ntotal + 1
203 if (outside) cycle reports
204 iv%info(ships)%nlocal = iv%info(ships)%nlocal + 1
205 fm = 13
206
207 case (531, 561, 562) ; ! Buoy
208 if (.not.use_buoyobs .or. iv%info(buoy)%ntotal == max_buoy_input) cycle reports
209 if (n==1) iv%info(buoy)%ntotal = iv%info(buoy)%ntotal + 1
210 if (outside) cycle reports
211 iv%info(buoy)%nlocal = iv%info(buoy)%nlocal + 1
212 fm = 18
213
214 case (511)
215 ! case (181, 281) ; ! Synop
216 if (.not.use_synopobs .or. iv%info(synop)%ntotal == max_synop_input) cycle reports
217 if (n==1) iv%info(synop)%ntotal = iv%info(synop)%ntotal + 1
218 if (outside) cycle reports
219 iv%info(synop)%nlocal = iv%info(synop)%nlocal + 1
220 fm = 12
221
222 case (512)
223 ! case (187, 287) ; ! Metar
224 if (.not.use_metarobs .or. iv%info(metar)%ntotal == max_metar_input) cycle reports
225 if (n==1) iv%info(metar)%ntotal = iv%info(metar)%ntotal + 1
226 if (outside) cycle reports
227 iv%info(metar)%nlocal = iv%info(metar)%nlocal + 1
228 fm = 15
229
230 case (63)
231 ! case (242:246, 252:253, 255) ; ! Geo. CMVs
232 if (.not.use_geoamvobs .or. iv%info(geoamv)%ntotal == max_geoamv_input) cycle reports
233 if (n==1) iv%info(geoamv)%ntotal = iv%info(geoamv)%ntotal + 1
234 if (outside) cycle reports
235 iv%info(geoamv)%nlocal = iv%info(geoamv)%nlocal + 1
236 fm = 88
237
238 case (582)
239 if (.not.use_qscatobs .or. iv%info(qscat)%ntotal == max_qscat_input) cycle reports
240 if (n==1) iv%info(qscat)%ntotal = iv%info(qscat)%ntotal + 1
241 if (outside) cycle reports
242 iv%info(qscat)%nlocal = iv%info(qscat)%nlocal + 1
243 fm = 281
244
245 case (583) ! GPS PW
246 if (.not.use_gpspwobs .or. iv%info(gpspw)%ntotal == max_gpspw_input) cycle reports
247 if (n==1) iv%info(gpspw)%ntotal = iv%info(gpspw)%ntotal + 1
248 if (outside) cycle reports
249 iv%info(gpspw)%nlocal = iv%info(gpspw)%nlocal + 1
250 fm = 111
251
252 case (584) ! GPS REF
253 if (.not.use_gpsrefobs .or. iv%info(gpsref)%ntotal == max_gpsref_input) cycle reports
254 if (n==1) iv%info(gpsref)%ntotal = iv%info(gpsref)%ntotal + 1
255 if (outside) cycle reports
256 iv%info(gpsref)%nlocal = iv%info(gpsref)%nlocal + 1
257 fm = 116
258
259 case (71, 72)
260 if (.not.use_profilerobs .or. iv%info(profiler)%ntotal == max_profiler_input) cycle reports
261 if (n==1) iv%info(profiler)%ntotal = iv%info(profiler)%ntotal + 1
262 if (outside) cycle reports
263 iv%info(profiler)%nlocal = iv%info(profiler)%nlocal + 1
264 fm = 132
265 case (571, 65)
266 if (.not. use_ssmiretrievalobs .or. iv%info(ssmi_rv)%ntotal == max_ssmi_rv_input) cycle reports
267 if (n==1) iv%info(ssmi_rv)%ntotal = iv%info(ssmi_rv)%ntotal + 1
268 if (outside) cycle reports
269 iv%info(ssmi_rv)%nlocal = iv%info(ssmi_rv)%nlocal + 1
270 fm = 125 ! ssmi wind speed & tpw
271 case default
272 select case (kx)
273 case (111 , 210) ; ! Tropical Cyclone Bogus
274 ! Note Tropical cyclone Bougus is given type 135 in Obs-ascii
275 if (.not.use_bogusobs .or. iv%info(bogus)%ntotal == max_bogus_input) cycle reports
276 if (n==1) iv%info(bogus)%ntotal = iv%info(bogus)%ntotal + 1
277 if (outside) cycle reports
278 iv%info(bogus)%nlocal = iv%info(bogus)%nlocal + 1
279 fm = 135
280
281 case default
282 write(unit=message(1), fmt='(a, 2i12)') &
283 'unsaved obs found with kx & t29= ',kx,t29
284 call da_warning(__FILE__,__LINE__,message(1:1))
285 end select
286 end select
287 end do dup_loop
288 obs_index=fm_index(fm)
289 iv%info(obs_index)%max_lev = max(iv%info(obs_index)%max_lev, nlevels)
290 end do reports
291 iv%info(sonde_sfc)%max_lev=1
292
293 close(iunit)
294 call da_free_unit(iunit)
295
296 if (trace_use) call da_trace_exit("da_scan_obs_bufr")
297 #endif
298
299 end subroutine da_scan_obs_bufr
300