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
33 integer :: irain, iprecip
34
35 count_obs_ssmir = count_obs_number_type(0, 0, 0, 0)
36 count_obs_ssmit = count_obs_number_type(0, 0, 0, 0)
37
38 isfilter = .true. ! filter out rain points
39 irain = 0
40
41 ! open file
42
43 call da_get_unit(iunit)
44 open(unit = iunit, &
45 FILE = trim(filename), &
46 FORM = 'FORMATTED', &
47 ACCESS = 'SEQUENTIAL', &
48 iostat = iost, &
49 STATUS = 'OLD')
50
51 if (iost /= 0) then
52 call da_warning(__FILE__,__LINE__, &
53 (/"Cannot open SSMI file "//trim(filename)/))
54 return
55 end if
56
57 rewind (unit = iunit)
58
59 ! 2. read header
60 ! ===============
61
62 ! 2.1 read first line
63 ! ---------------
64
65 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
66 if (iost /= 0) then
67 call da_error(__FILE__,__LINE__, &
68 (/"Cannot read SSMI file"//trim(filename)/))
69 else
70 write(unit=stdout, fmt='(/2a/)') &
71 'in da_read_ssmi_info.inc, char_ned=', trim(char_ned)
72 end if
73
74 ! 2.3 read number of reports
75 ! ---------------------
76
77 do
78 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
79
80 if (iost /= 0) then
81 call da_error(__FILE__,__LINE__, &
82 (/"Cannot read SSMI file"//trim(filename)/))
83 end if
84
85 if (char_ned(1:6) == 'NESTIX') exit
86
87 end do
88
89 do
90
91 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
92 if (char_ned(1:6) == 'INFO ') exit
93
94 end do
95
96 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
97
98 ! read formats
99 ! ------------
100
101 read (unit=iunit, fmt = '(A,1X,A)') &
102 fmt_name, fmt_info, &
103 fmt_name, fmt_loc
104
105 ! SKIP 1 line
106 ! -----------
107
108 read (unit=iunit, fmt = '(A)') fmt_name
109
110 ! loop over records
111 ! -----------------
112
113 reports: do
114
115 ! read station general info
116 ! =========================
117
118 read (unit=iunit, fmt = fmt_info, iostat = iost) &
119 info % platform, &
120 info % date_char, &
121 info % name, &
122 info % levels, &
123 info % lat, &
124 info % lon, &
125 info % elv, &
126 info % id
127
128 read(unit=info % platform (4:6),fmt='(I3)') fm
129
130 if (iost /= 0) then
131 exit reports
132 end if
133
134 ob % total_obs = ob % total_obs + 1
135
136 select case(fm)
137 case (125) ;
138 ! read surface wind speed and precipitable water
139 read (unit=iunit, fmt = fmt_loc) speed, tpw
140 case (126) ;
141 read (unit=iunit, fmt = fmt_loc) &
142 tb19v, tb19h, tb22v, tb37v, tb37h, tb85v, tb85h
143
144 tb19v % error = tb19v % error + 2.
145 tb19h % error = tb19h % error + 2.
146 tb22v % error = tb22v % error + 2.
147 tb37h % error = tb37h % error + 2.
148 tb37v % error = tb37v % error + 2.
149 tb85h % error = tb85h % error + 2.
150 tb85v % error = tb85v % error + 2.
151
152 case default;
153 write(unit=message(1), fmt='(a, i6)') 'unsaved ssmi obs found, fm=', fm
154 write(unit=message(2), fmt='(a, 2f12.6)') &
155 'info%(lon,lat)=', info%lon, info%lat
156 call da_warning(__FILE__,__LINE__,message(1:2))
157 end select
158
159 ! check if obs is in horizontal domain
160 ! ====================================
161
162 ! Compute the model horizontal coordinate x, y
163 ! Check if obs is wihin horizontal domain
164
165 call da_ll_to_xy (info, loc, xp, outside)
166
167 if (outside) cycle reports
168
169 loc % pw % inv = missing_r
170 loc % pw % qc = missing
171 loc % slp = loc % pw
172
173 select case(fm)
174 case (125) ;
175 if (.not. Use_SsmiRetrievalObs) cycle reports
176
177 ! Check if at least one field is present
178 if ((tpw % qc == missing) .AND. (speed % qc == missing)) then
179 count_obs_ssmir % num_missing = &
180 count_obs_ssmir % num_missing + 1
181 cycle reports
182 end if
183
184 ! fill permanent structure
185 ! ========================
186
187 ob % num_ssmi_retrieval = ob % num_ssmi_retrieval + 1
188
189 ! One more data used
190
191 count_obs_ssmir % num_used = count_obs_ssmir % num_used + 1
192
193 ! Initialize other non read fields
194
195 ob % ssmi_retrieval (ob%num_ssmi_retrieval) % loc = loc
196
197 ! Initialize with read fields
198
199 ob % ssmi_retrieval (ob%num_ssmi_retrieval) % info = info
200
201 ob % ssmi_retrieval (ob%num_ssmi_retrieval) % speed = speed
202 ob % ssmi_retrieval (ob%num_ssmi_retrieval) % tpw = tpw
203
204 case (126) ;
205 if (.not. Use_SsmiTbObs) cycle reports
206
207 ! Check if at least one field is present
208
209 if ((tb19v % qc == missing) .AND. (tb19h % qc == missing) .AND. &
210 (tb22v % qc == missing) .AND. &
211 (tb37v % qc == missing) .AND. (tb37h % qc == missing) .AND. &
212 (tb85v % qc == missing) .AND. (tb85h % qc == missing)) then
213 count_obs_ssmit % num_missing = &
214 count_obs_ssmit % num_missing + 1
215 ! write (unit=stdout,fmt=*) 'missing data'
216 cycle reports
217 end if
218
219 ! filter rain pixels
220 ! ====================================
221
222 if (isfilter) then
223 ipass = .false.
224 iprecip = 0
225 call filter(tb19v%inv, tb19h%inv, tb22v%inv, tb37v%inv, &
226 tb37h%inv, tb85v%inv, tb85h%inv, ipass, iprecip, &
227 info%lat)
228 if (iprecip .eq. 1) then
229 tb19v % qc = -88.
230 tb19h % qc = -88.
231 tb22v % qc = -88.
232 tb37h % qc = -88.
233 tb37v % qc = -88.
234 tb85h % qc = -88.
235 tb85v % qc = -88.
236 irain = irain + 1
237 cycle reports
238 end if
239 end if
240
241 ! fill permanent structure
242 ! ========================
243
244 ! One more data read in
245
246 ob % num_ssmi_tb = ob % num_ssmi_tb + 1
247
248 ! One more data used
249
250 count_obs_ssmit % num_used = count_obs_ssmit % num_used + 1
251
252 ! Initialize other non read fields
253
254 ob % ssmi_tb (ob%num_ssmi_tb) % loc = loc
255
256 ! Initialize with read fields
257
258 ob % ssmi_tb (ob % num_ssmi_tb) % info = info
259
260 ob % ssmi_tb (ob % num_ssmi_tb) % tb19v = tb19v
261 ob % ssmi_tb (ob % num_ssmi_tb) % tb19h = tb19h
262 ob % ssmi_tb (ob % num_ssmi_tb) % tb22v = tb22v
263 ob % ssmi_tb (ob % num_ssmi_tb) % tb37v = tb37v
264 ob % ssmi_tb (ob % num_ssmi_tb) % tb37h = tb37h
265 ob % ssmi_tb (ob % num_ssmi_tb) % tb85v = tb85v
266 ob % ssmi_tb (ob % num_ssmi_tb) % tb85h = tb85h
267
268 case default;
269 ! Do nothing.
270 end select
271 end do reports
272
273 close(iunit)
274 call da_free_unit(iunit)
275
276 write(unit=stdout, fmt='(5x,a,i6,a)') &
277 'Read: ', ob % num_ssmi_retrieval , ' SSMI_RETRIEVAL reports,', &
278 'Read: ', ob % num_ssmi_tb , ' SSMI_Tb reports,', &
279 'Read: ', ob % total_obs, ' Total reports.', &
280 'There are ', ob % total_obs &
281 - ob % num_ssmi_retrieval &
282 - ob % num_ssmi_tb, &
283 ' reports unsaved.'
284
285 write(unit=stdout, fmt='(/,25x,a)') &
286 ' used outdomain max_er_chk missing'
287 write(unit=stdout, fmt='(4x,a,4i11)') &
288 'SSMI_RETRIEVAL_diag: ', count_obs_ssmir
289 write(unit=stdout, fmt='(4x,a,4i11)') &
290 'SSMI_Tb_diag: ', count_obs_ssmit
291 write(unit=stdout, fmt='(/,a)') &
292 "------------------------------------------------------------------------------"
293
294 write(unit=stdout, fmt='(/,5x,a,i6/)') '** Rain contaminated SSMI_Tb =', irain
295
296 end subroutine da_read_ssmi
297
298