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