da_read_ssmi.inc
References to this file elsewhere.
1 subroutine da_read_ssmi (ob, xp, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Read the header of a SSMI 2.0 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 character(*), intent(in) :: filename
12
13 character (LEN = 10) :: fmt_name
14
15 character (LEN = 160) :: fmt_info, &
16 fmt_loc
17 character (LEN = 120) :: char_ned
18
19 integer :: iost, fm,iunit
20
21 type (model_loc_type) :: loc
22 type (info_type) :: info
23 type (field_type) :: speed, tpw
24
25 type (field_type) :: tb19v, tb19h, tb22v
26 type (field_type) :: tb37v, tb37h, tb85v, tb85h
27
28 type (count_obs_number_type) :: count_obs_ssmir
29 type (count_obs_number_type) :: count_obs_ssmit
30
31 logical :: isfilter,ipass
32 logical :: outside, outside_all
33 integer :: irain, iprecip
34 integer :: n, ndup
35
36 count_obs_ssmir = count_obs_number_type(0, 0, 0, 0)
37 count_obs_ssmit = count_obs_number_type(0, 0, 0, 0)
38
39 isfilter = .true. ! filter out rain points
40 irain = 0
41
42 ! open file
43
44 call da_get_unit(iunit)
45 open(unit = iunit, &
46 FILE = trim(filename), &
47 FORM = 'FORMATTED', &
48 ACCESS = 'SEQUENTIAL', &
49 iostat = iost, &
50 STATUS = 'OLD')
51
52 if (iost /= 0) then
53 call da_warning(__FILE__,__LINE__, &
54 (/"Cannot open SSMI file "//trim(filename)/))
55 call da_free_unit(iunit)
56 return
57 end if
58
59 rewind (unit = iunit)
60
61 ! 2. read header
62 ! ===============
63
64 ! 2.1 read first line
65 ! ---------------
66
67 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
68 if (iost /= 0) then
69 call da_error(__FILE__,__LINE__, &
70 (/"Cannot read SSMI file"//trim(filename)/))
71 else
72 write(unit=stdout, fmt='(/2a/)') &
73 'in da_read_ssmi_info.inc, char_ned=', trim(char_ned)
74 end if
75
76 ! 2.3 read number of reports
77 ! ---------------------
78
79 do
80 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
81
82 if (iost /= 0) then
83 call da_error(__FILE__,__LINE__, &
84 (/"Cannot read SSMI file"//trim(filename)/))
85 end if
86
87 if (char_ned(1:6) == 'NESTIX') exit
88
89 end do
90
91 do
92
93 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
94 if (char_ned(1:6) == 'INFO ') exit
95
96 end do
97
98 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
99
100 ! read formats
101 ! ------------
102
103 read (unit=iunit, fmt = '(A,1X,A)') &
104 fmt_name, fmt_info, &
105 fmt_name, fmt_loc
106
107 ! SKIP 1 line
108 ! -----------
109
110 read (unit=iunit, fmt = '(A)') fmt_name
111
112 ! loop over records
113 ! -----------------
114
115 reports: do
116
117 ! read station general info
118 ! =========================
119
120 read (unit=iunit, fmt = fmt_info, iostat = iost) &
121 info % platform, &
122 info % date_char, &
123 info % name, &
124 info % levels, &
125 info % lat, &
126 info % lon, &
127 info % elv, &
128 info % id
129
130 read(unit=info % platform (4:6),fmt='(I3)') fm
131
132 if (iost /= 0) then
133 exit reports
134 end if
135
136 ob % total_obs = ob % total_obs + 1
137
138 select case(fm)
139 case (125) ;
140 ! read surface wind speed and precipitable water
141 read (unit=iunit, fmt = fmt_loc) speed, tpw
142 case (126) ;
143 read (unit=iunit, fmt = fmt_loc) &
144 tb19v, tb19h, tb22v, tb37v, tb37h, tb85v, tb85h
145
146 tb19v % error = tb19v % error + 2.
147 tb19h % error = tb19h % error + 2.
148 tb22v % error = tb22v % error + 2.
149 tb37h % error = tb37h % error + 2.
150 tb37v % error = tb37v % error + 2.
151 tb85h % error = tb85h % error + 2.
152 tb85v % error = tb85v % error + 2.
153
154 case default;
155 write(unit=message(1), fmt='(a, i6)') 'unsaved ssmi obs found, fm=', fm
156 write(unit=message(2), fmt='(a, 2f12.6)') &
157 'info%(lon,lat)=', info%lon, info%lat
158 call da_warning(__FILE__,__LINE__,message(1:2))
159 end select
160
161 ! check if obs is in horizontal domain
162 ! ====================================
163
164 ! Compute the model horizontal coordinate x, y
165 ! Check if obs is wihin horizontal domain
166
167 call da_ll_to_xy (info, loc, xp, outside, outside_all)
168
169 if (outside_all) cycle reports
170
171 loc % pw % inv = missing_r
172 loc % pw % qc = missing_data
173 loc % slp = loc % pw
174
175 ! Loop over duplicating obs for global
176 ndup = 1
177 if (global .and. &
178 (loc%i < xp%ids .or. loc%i >= xp%ide)) ndup= 2
179
180 ! It is possible that logic for counting obs is incorrect for the
181 ! global case with >1 MPI tasks due to obs duplication, halo, etc.
182 ! TBH: 20050913
183
184 if (.not.outside) then
185 if (print_detail_obs .and. ndup > 1) then
186 write(unit=stdout, fmt = fmt_info) &
187 info%platform, &
188 info%date_char, &
189 info%name, &
190 info%levels, &
191 info%lat, &
192 info%lon, &
193 info%elv, &
194 info%id
195
196 write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
197 ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
198 loc%i, loc%j, &
199 loc%dx, loc%dxm, &
200 loc%dy, loc%dym
201 end if
202 end if
203
204 dup_loop: do n = 1, ndup
205
206 select case(fm)
207 case (125) ;
208
209 if (n==1) ob%num_ssmi_retrieval_glo = ob%num_ssmi_retrieval_glo + 1
210 if (outside) then
211 cycle reports
212 end if
213
214 if (.not. Use_SsmiRetrievalObs) cycle reports
215
216 ! Check if at least one field is present
217 if ((tpw % qc == missing_data) .AND. (speed % qc == missing_data)) then
218 count_obs_ssmir % num_missing = &
219 count_obs_ssmir % num_missing + 1
220 cycle reports
221 end if
222
223 ! fill permanent structure
224 ! ========================
225
226 ob % num_ssmi_retrieval = ob % num_ssmi_retrieval + 1
227 ! Track serial obs index for reassembly of obs during bit-for-bit
228 ! tests with different numbers of MPI tasks.
229 loc%obs_global_index = ob%num_ssmi_retrieval_glo
230
231 ! One more data used
232
233 count_obs_ssmir % num_used = count_obs_ssmir % num_used + 1
234
235 ! Initialize other non read fields
236
237 ob % ssmi_retrieval (ob%num_ssmi_retrieval) % loc = loc
238
239 ! Initialize with read fields
240
241 ob % ssmi_retrieval (ob%num_ssmi_retrieval) % info = info
242
243 ob % ssmi_retrieval (ob%num_ssmi_retrieval) % speed = speed
244 ob % ssmi_retrieval (ob%num_ssmi_retrieval) % tpw = tpw
245
246 case (126) ;
247
248 if (n==1) ob%num_ssmi_tb_glo = ob%num_ssmi_tb_glo + 1
249 if (outside) then
250 cycle reports
251 end if
252
253 if (.not. Use_SsmiTbObs) cycle reports
254
255 ! Check if at least one field is present
256
257 if ((tb19v % qc == missing_data) .AND. (tb19h % qc == missing_data) .AND. &
258 (tb22v % qc == missing_data) .AND. &
259 (tb37v % qc == missing_data) .AND. (tb37h % qc == missing_data) .AND. &
260 (tb85v % qc == missing_data) .AND. (tb85h % qc == missing_data)) then
261 count_obs_ssmit % num_missing = &
262 count_obs_ssmit % num_missing + 1
263 ! write (unit=stdout,fmt=*) 'missing data'
264 cycle reports
265 end if
266
267 ! filter rain pixels
268 ! ====================================
269
270 if (isfilter) then
271 ipass = .false.
272 iprecip = 0
273 call filter(tb19v%inv, tb19h%inv, tb22v%inv, tb37v%inv, &
274 tb37h%inv, tb85v%inv, tb85h%inv, ipass, iprecip, &
275 info%lat)
276 if (iprecip .eq. 1) then
277 tb19v % qc = -88.
278 tb19h % qc = -88.
279 tb22v % qc = -88.
280 tb37h % qc = -88.
281 tb37v % qc = -88.
282 tb85h % qc = -88.
283 tb85v % qc = -88.
284 irain = irain + 1
285 cycle reports
286 end if
287 end if
288
289 ! fill permanent structure
290 ! ========================
291
292 ! One more data read in
293
294 ob % num_ssmi_tb = ob % num_ssmi_tb + 1
295 ! Track serial obs index for reassembly of obs during bit-for-bit
296 ! tests with different numbers of MPI tasks.
297 loc%obs_global_index = ob%num_ssmi_tb_glo
298
299 ! One more data used
300
301 count_obs_ssmit % num_used = count_obs_ssmit % num_used + 1
302
303 ! Initialize other non read fields
304
305 ob % ssmi_tb (ob%num_ssmi_tb) % loc = loc
306
307 ! Initialize with read fields
308
309 ob % ssmi_tb (ob % num_ssmi_tb) % info = info
310
311 ob % ssmi_tb (ob % num_ssmi_tb) % tb19v = tb19v
312 ob % ssmi_tb (ob % num_ssmi_tb) % tb19h = tb19h
313 ob % ssmi_tb (ob % num_ssmi_tb) % tb22v = tb22v
314 ob % ssmi_tb (ob % num_ssmi_tb) % tb37v = tb37v
315 ob % ssmi_tb (ob % num_ssmi_tb) % tb37h = tb37h
316 ob % ssmi_tb (ob % num_ssmi_tb) % tb85v = tb85v
317 ob % ssmi_tb (ob % num_ssmi_tb) % tb85h = tb85h
318
319 case default;
320 ! Do nothing.
321 end select
322 end do dup_loop
323 end do reports
324
325 close(iunit)
326 call da_free_unit(iunit)
327
328 write(unit=stdout, fmt='(5x,a,i6,a)') &
329 'Read: ', ob % num_ssmi_retrieval , ' SSMI_RETRIEVAL reports,', &
330 'Read: ', ob % num_ssmi_tb , ' SSMI_Tb reports,', &
331 'Read: ', ob % total_obs, ' Total reports.', &
332 'There are ', ob % total_obs &
333 - ob % num_ssmi_retrieval &
334 - ob % num_ssmi_tb, &
335 ' reports unsaved.'
336
337 write(unit=stdout, fmt='(/,25x,a)') &
338 ' used outdomain max_er_chk missing'
339 write(unit=stdout, fmt='(4x,a,4i11)') &
340 'SSMI_RETRIEVAL_diag: ', count_obs_ssmir
341 write(unit=stdout, fmt='(4x,a,4i11)') &
342 'SSMI_Tb_diag: ', count_obs_ssmit
343 write(unit=stdout, fmt='(/,a)') &
344 "------------------------------------------------------------------------------"
345
346 write(unit=stdout, fmt='(/,5x,a,i6/)') '** Rain contaminated SSMI_Tb =', irain
347
348 end subroutine da_read_ssmi
349
350