da_scan_ssmi.inc
References to this file elsewhere.
1 subroutine da_scan_ssmi (ob, xp, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Read the header of a 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 logical :: isfilter,ipass
29 logical :: outside
30 integer :: irain, iprecip
31
32 isfilter = .true. ! filter out rain points
33 irain = 0
34
35 ob%ob_numb(ob%current_ob_time)%ssmi_retrieval = ob%num_ssmi_retrieval
36 ob%ob_numb(ob%current_ob_time)%ssmi_tb = ob%num_ssmi_tb
37
38 ! open FILE
39 call da_get_unit(iunit)
40 open(unit = iunit, &
41 FILE = trim(filename), &
42 FORM = 'FORMATTED', &
43 ACCESS = 'SEQUENTIAL', &
44 iostat = iost, &
45 STATUS = 'OLD')
46
47 if (iost /= 0) then
48 call da_warning(__FILE__,__LINE__, &
49 (/"Cannot open SSMI file "//trim(filename)/))
50 return
51 end if
52
53 rewind (unit = iunit)
54
55 ! 2. read HEADER
56 ! ===============
57
58 ! 2.1 read FIRST LINE
59 ! ---------------
60
61 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
62
63 ! 2.2 PROCESS ERROR
64 ! -------------
65
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_scan_ssmi.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 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)') &
99 fmt_name, fmt_info, &
100 fmt_name, fmt_loc
101
102 ! SKIP 1 line
103 ! -----------
104
105 read (unit=iunit, fmt = '(A)') fmt_name
106
107 ! write (unit = stdout, fmt = '(/,2(A,F3.0,/))') &
108 ! 'Error max test ratio for SWS = ',max_error_uv, &
109 ! 'Error max test ratio for PW = ',max_error_pw, &
110 ! 'Error max test ratio for PW = ',max_error_tb
111
112
113 ! loop over records
114 ! -----------------
115
116 reports: do
117
118 ! read station general info
119 ! =========================
120
121 read (unit=iunit, fmt = fmt_info, iostat = iost) &
122 info % platform, &
123 info % date_char, &
124 info % name, &
125 info % levels, &
126 info % lat, &
127 info % lon, &
128 info % elv, &
129 info % id
130
131 read(unit=info % platform (4:6),fmt='(I3)') fm
132
133 if (iost /= 0) then
134 exit reports
135 end if
136
137 ob % total_obs = ob % total_obs + 1
138
139 select case(fm)
140 case (125) ;
141 ! read SURFACE WinD SPEED AND PRECIPITABLE WATER
142 read (unit=iunit, fmt = fmt_loc) speed, tpw
143 case (126) ;
144 read (unit=iunit, fmt = fmt_loc) &
145 tb19v, tb19h, tb22v, tb37v, tb37h, tb85v, tb85h
146 case default;
147 write(unit=stdout, fmt='(/a/)') &
148 'warning warning warning warning warning :'
149 write(unit=stdout, fmt='(a, i6)') 'unsaved ssmi obs found, fm=', fm
150 write(unit=stdout, fmt='(a, 2f12.6)') &
151 'info%(lon,lat)=', info%lon, info%lat
152 end select
153
154 ! check if obs is in horizontal domain
155 ! ====================================
156
157 ! Compute the model horizontal coordinate x, y
158 ! Check if obs is wihin horizontal domain
159
160 call da_ll_to_xy (info, loc, xp, outside)
161
162 if (outside) cycle reports
163
164 select case(fm)
165 case (125) ;
166 if (.not. Use_SsmiRetrievalObs) cycle reports
167
168 ! Check if at least one field is present
169 if ((tpw % qc == missing) .AND. (speed % qc == missing)) then
170 cycle reports
171 end if
172
173 ob % num_ssmi_retrieval = ob % num_ssmi_retrieval + 1
174
175 case (126) ;
176 if (.not. Use_SsmiTbObs) cycle reports
177
178 ! Check if at least one field is present
179
180 if ((tb19v % qc == missing) .AND. (tb19h % qc == missing) .AND. &
181 (tb22v % qc == missing) .AND. &
182 (tb37v % qc == missing) .AND. (tb37h % qc == missing) .AND. &
183 (tb85v % qc == missing) .AND. (tb85h % qc == missing)) then
184 cycle reports
185 end if
186
187 ! FILTER RAin PIXELS
188 ! ====================================
189
190 if (isfilter) then
191
192 ipass = .false.
193 iprecip = 0
194 call filter(tb19v%inv, tb19h%inv, tb22v%inv, tb37v%inv, &
195 tb37h%inv, tb85v%inv, tb85h%inv, ipass, iprecip, &
196 info%lat)
197 if (iprecip == 1) then
198 irain = irain + 1
199 cycle reports
200 end if
201 end if
202
203 ob % num_ssmi_tb = ob % num_ssmi_tb + 1
204
205 case default;
206 ! Do nothing.
207 end select
208
209 end do reports
210
211 close(unit=iunit)
212 call da_free_unit(iunit)
213
214 write(unit=stdout, fmt='(5x,a,i6,a)') &
215 'Read: ', ob % num_ssmi_retrieval , ' SSMI_RETRIEVAL reports,', &
216 'Read: ', ob % num_ssmi_tb , ' SSMI_Tb reports,', &
217 'Read: ', ob % total_obs, ' Total Observations.'
218
219 write(unit=stdout, fmt='(/,5x,a,i6/)') &
220 '** Rain contaminated SSMI_Tb =', irain
221
222 ob%ob_numb(ob%current_ob_time)%ssmi_retrieval = ob%num_ssmi_retrieval
223 ob%ob_numb(ob%current_ob_time)%ssmi_tb = ob%num_ssmi_tb
224
225 end subroutine da_scan_ssmi
226
227