da_scan_obs_bufr.inc
References to this file elsewhere.
1 subroutine da_scan_obs_bufr (iv, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Scan BUFR observation file for input to wrfvar
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (iv_type), intent(inout) :: iv
10 character(len=*), intent(in) :: filename
11
12 #ifdef BUFR
13
14 type (multi_level_type) :: platform
15 logical :: outside
16 logical :: outside_all
17
18 character(len=40) :: hdstr
19 character(len=10) :: date
20 character(len=8) :: subset
21 real :: hdr(7)
22
23 integer :: iost, ndup, n, report
24 integer :: year,month, day,hour
25 integer :: iret,idate, kx, nrecs,miscd, t29
26 integer :: iunit
27
28 if (trace_use) call da_trace_entry("da_scan_obs_bufr")
29
30 ! open FILE
31 ! ---------
32
33 call da_get_unit(iunit)
34 open(unit = iunit, FILE = trim(filename), &
35 iostat = iost, form = 'unformatted', STATUS = 'OLD')
36
37 if (iost /= 0) then
38 write(unit=message(1), fmt='(A,I3,A,A)') &
39 'error in obs input file unit ',iunit, &
40 'obs filename:', trim(filename)
41 message(2)="bufr gts observations cannot be found or cannot be opened"
42 call da_warning(__FILE__,__LINE__,message(1:2))
43 call da_free_unit(iunit)
44 return
45 end if
46
47 !--------------------------------------------------------------------------
48 ! Initialise the obs-counter
49 !--------------------------------------------------------------------------
50
51 hdstr='SID XOB YOB DHR TYP ELV T29 '
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 call readmg(iunit,subset,idate,iret)
63
64 if (iret/=0) then
65 message(1) = "No bufr observations"
66 write(unit=message(2), fmt='(a, i4)') &
67 'return code from readmg:', iret, &
68 'Reach the end of obs unit: ', iunit
69 call da_error(__FILE__,__LINE__,message(1:2))
70 end if
71
72 write(unit=date,fmt='(i10)') idate
73 read (unit=date,fmt='(i4,3i2)') year,month, day,hour
74 write(unit=message(1),fmt=*) 'Bufr file date is ',year,month, day,hour
75 call da_message(message(1:1))
76
77 report = 0 ! report number in file
78
79 reports: do
80 report = report+1
81
82 call readsb(iunit,iret)
83
84 if (iret/=0) then
85 call readmg(iunit,subset,idate,iret)
86 if (iret/=0) then
87 write (unit=message(1), fmt='(a, i4)') &
88 'return code from readmg:', iret, &
89 'Reach the end of obs unit: ', iunit
90 call da_warning(__FILE__,__LINE__,message(1:1))
91 exit reports
92 end if
93
94 cycle reports
95 end if
96
97 nrecs=nrecs+1
98 call ufbint(iunit,hdr,7,1,iret,hdstr)
99
100 platform % info % name(1:8) = subset
101 platform % info % id = hdstr(1:5)
102
103 if (hdr(2) > 180.0) hdr(2)=hdr(2)-360.0
104
105 ! Put a check on Lat
106 hdr(3) = MAX(hdr(3), -89.95)
107 hdr(3) = Min(hdr(3), 89.95)
108 platform%info%lon = hdr(2)
109 platform%info%lat = hdr(3)
110
111 if (platform%info%lon == 180.0 ) platform%info%lon =-180.000
112 ! Fix funny wind direction at Poles
113 if (platform%info%lat <= -89.95 .or. platform%info%lat >= 89.95) then
114 platform%info%lon = 0.0
115 end if
116
117 ! Restrict to a range of reports, useful for debugging
118
119 if (report < report_start) cycle
120 if (report > report_end) exit
121
122 call da_llxy (platform%info, platform%loc, outside, outside_all)
123
124 if (outside_all) cycle reports
125
126 t29 = int(0.1 + hdr(7))
127 kx=int(0.1+hdr(5))
128
129 if (kx .eq. 183) then ! reset kx based on t29
130 if (t29 .eq. 511) kx = 181
131 if (t29 .eq. 512) kx = 187
132 if (t29 .eq. 522) kx = 180
133 if (t29 .eq. 523) kx = 180
134 if (t29 .eq. 531) kx = 180
135 if (t29 .eq. 561) kx = 180
136 if (t29 .eq. 562) kx = 180
137 end if
138
139 ! if ((kx >= 160) .and. (kx <= 179)) then ! bypass satellite data
140 ! if (t29 .eq. 61 .or. t29 .eq. 63 .or. t29 .ge. 571) then
141 ! cycle reports
142 ! end if
143
144 ! Loop over duplicating obs for global
145 ndup = 1
146 if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) ndup= 2
147
148 if (test_wrfvar) ndup = 1
149
150 do n = 1, ndup
151 select case(t29)
152
153 case (11, 12, 13, 22, 23, 31)
154 select case(kx)
155 case (120, 122, 132, 220, 222, 232) ; ! Sound
156 if (.not.use_soundobs .or. iv%info(sound)%ntotal == max_sound_input) cycle reports
157 if (n==1) iv%info(sound)%ntotal = iv%info(sound)%ntotal + 1
158 if (n==1) iv%info(sonde_sfc)%ntotal = iv%info(sonde_sfc)%ntotal + 1
159 if (outside) cycle reports
160 iv%info(sound)%nlocal = iv%info(sound)%nlocal + 1
161 iv%info(sonde_sfc)%nlocal = iv%info(sonde_sfc)%nlocal + 1
162 platform%loc%obs_global_index = iv%info(sound)%ntotal
163
164 case (221) ; ! Pilot
165 if (.not.use_pilotobs .or. iv%info(pilot)%ntotal == max_pilot_input) cycle reports
166 if (n==1) iv%info(pilot)%ntotal = iv%info(pilot)%ntotal + 1
167 if (outside) cycle reports
168 iv%info(pilot)%nlocal = iv%info(pilot)%nlocal + 1
169 platform%loc%obs_global_index = iv%info(pilot)%ntotal
170 end select
171
172 case (41)
173 ! case (130:131, 133, 230:231, 233) ; ! Airep
174 if (.not.use_airepobs .or. iv%info(airep)%ntotal == max_airep_input) cycle reports
175 if (n==1) iv%info(airep)%ntotal = iv%info(airep)%ntotal + 1
176 if (outside) cycle reports
177 iv%info(airep)%nlocal = iv%info(airep)%nlocal + 1
178 platform%loc%obs_global_index = iv%info(airep)%ntotal
179
180 ! case (180, 182, 280, 282) ; ! Ships
181 case (522, 523) ;
182 if (.not.use_shipsobs .or. iv%info(ships)%ntotal == max_ships_input) cycle reports
183 if (n==1) iv%info(ships)%ntotal = iv%info(ships)%ntotal + 1
184 if (outside) cycle reports
185 iv%info(ships)%nlocal = iv%info(ships)%nlocal + 1
186 platform%loc%obs_global_index = iv%info(ships)%ntotal
187
188 case (531, 561, 562) ; ! Buoy
189 if (.not.use_synopobs .or. iv%info(synop)%ntotal == max_synop_input) cycle reports
190 if (n==1) iv%info(synop)%ntotal = iv%info(synop)%ntotal + 1
191 if (outside) cycle reports
192 iv%info(synop)%nlocal = iv%info(synop)%nlocal + 1
193 platform%loc%obs_global_index = iv%info(synop)%ntotal
194
195 case (511)
196 ! case (181, 281) ; ! Synop
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 platform%loc%obs_global_index = iv%info(metar)%ntotal
202
203 case (512)
204 ! case (187 , 287) ; ! Metar
205 if (.not.use_metarobs .or. iv%info(metar)%ntotal == max_metar_input) cycle reports
206 if (n==1) iv%info(metar)%ntotal = iv%info(metar)%ntotal + 1
207 if (outside) cycle reports
208 iv%info(metar)%nlocal = iv%info(metar)%nlocal + 1
209 platform%loc%obs_global_index = iv%info(metar)%ntotal
210
211 case (63)
212 ! case (242:246, 252:253, 255); ! Geo. AMVs
213 if (.not.use_geoamvobs .or. iv%info(geoamv)%ntotal == max_geoamv_input) cycle reports
214 if (n==1) iv%info(geoamv)%ntotal = iv%info(geoamv)%ntotal + 1
215 if (outside) cycle reports
216 iv%info(geoamv)%nlocal = iv%info(geoamv)%nlocal + 1
217 platform%loc%obs_global_index = iv%info(geoamv)%ntotal
218
219 case (71, 72) ; ! Profiler & VADWND - NEXRAD winds
220 if (.not.use_polaramvobs .or. iv%info(polaramv)%ntotal == max_polaramv_input) cycle reports
221 if (n==1) iv%info(polaramv)%ntotal = iv%info(polaramv)%ntotal + 1
222 if (outside) cycle reports
223 iv%info(polaramv)%nlocal = iv%info(polaramv)%nlocal + 1
224 platform%loc%obs_global_index = iv%info(polaramv)%ntotal
225
226 case (582) ; ! Quikscat
227 if (.not.use_qscatobs .or. iv%info(qscat)%ntotal == max_qscat_input) cycle reports
228 if (n==1) iv%info(qscat)%ntotal = iv%info(qscat)%ntotal + 1
229 if (outside) cycle reports
230 iv%info(qscat)%nlocal = iv%info(qscat)%nlocal + 1
231 platform%loc%obs_global_index = iv%info(qscat)%ntotal
232
233 case (583) ; ! GPS IPW
234 if (.not.use_gpspwobs .or. iv%info(gpspw)%ntotal == max_gpspw_input) cycle reports
235 if (n==1) iv%info(gpspw)%ntotal = iv%info(gpspw)%ntotal + 1
236 if (outside) cycle reports
237 iv%info(gpspw)%nlocal = iv%info(gpspw)%nlocal + 1
238 platform%loc%obs_global_index = iv%info(gpspw)%ntotal
239
240 case (584) ; ! GPS REF
241 if (.not.use_gpsrefobs .or. iv%info(gpsref)%ntotal == max_gpsref_input) cycle reports
242 if (n==1) iv%info(gpsref)%ntotal = iv%info(gpsref)%ntotal + 1
243 if (outside) cycle reports
244 iv%info(gpsref)%nlocal = iv%info(gpsref)%nlocal + 1
245 platform%loc%obs_global_index = iv%info(gpsref)%ntotal
246
247 case default ;
248 select case (kx)
249 case (111 , 210) ; ! Tropical Cyclone Bogus
250 ! Note Tropical cyclone Bougus is given type 135 in Obs-ascii
251 if (.not.use_bogusobs) cycle reports
252 iv%info(bogus)%nlocal = iv%info(bogus)%nlocal + 1
253
254 case default
255 miscd = miscd + 1
256 message(1)='unsaved obs found:'
257 write(unit=message(2), fmt='(2a)') &
258 'platform%info%platform=', platform%info%platform
259 call da_warning(__FILE__,__LINE__,message(1:2))
260 end select
261 end select
262 end do ! loop over duplicate
263 end do reports
264
265 write (unit=message(1),fmt=*) 'total bufr records =',nrecs
266 write (unit=message(3),fmt=*) 'number unhandled =',miscd
267 call da_message(message(1:3))
268
269 close(iunit)
270 call da_free_unit(iunit)
271
272 if (trace_use) call da_trace_exit("da_scan_obs_bufr")
273 #endif
274
275 end subroutine da_scan_obs_bufr
276
277