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 call da_free_unit(iunit)
51 return
52 end if
53
54 rewind (unit = iunit)
55
56 ! 2. read HEADER
57 ! ===============
58
59 ! 2.1 read FIRST LINE
60 ! ---------------
61
62 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
63
64 ! 2.2 PROCESS ERROR
65 ! -------------
66
67 if (iost /= 0) then
68 call da_error(__FILE__,__LINE__, &
69 (/"Cannot read SSMI file "//trim(filename)/))
70 else
71 write(unit=stdout, fmt='(/2a/)') &
72 'in da_scan_ssmi.inc, char_ned=', trim(char_ned)
73 end if
74
75 ! 2.3 read NUMBER OF REPORTS
76 ! ---------------------
77
78 do
79 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
80
81 if (iost /= 0) then
82 call da_error(__FILE__,__LINE__, &
83 (/"Cannot read SSMI file "//trim(filename)/))
84 end if
85
86 if (char_ned(1:6) == 'NESTIX') exit
87 end do
88
89 do
90 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
91 if (char_ned(1:6) == 'INFO ') exit
92 end do
93
94 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
95
96 ! read FORMATS
97 ! ------------
98
99 read (unit=iunit, fmt = '(A,1X,A)') &
100 fmt_name, fmt_info, &
101 fmt_name, fmt_loc
102
103 ! SKIP 1 line
104 ! -----------
105
106 read (unit=iunit, fmt = '(A)') fmt_name
107
108 ! write (unit = stdout, fmt = '(/,2(A,F3.0,/))') &
109 ! 'Error max test ratio for SWS = ',max_error_uv, &
110 ! 'Error max test ratio for PW = ',max_error_pw, &
111 ! 'Error max test ratio for PW = ',max_error_tb
112
113
114 ! loop over records
115 ! -----------------
116
117 reports: do
118
119 ! read station general info
120 ! =========================
121
122 read (unit=iunit, fmt = fmt_info, iostat = iost) &
123 info % platform, &
124 info % date_char, &
125 info % name, &
126 info % levels, &
127 info % lat, &
128 info % lon, &
129 info % elv, &
130 info % id
131
132 read(unit=info % platform (4:6),fmt='(I3)') fm
133
134 if (iost /= 0) then
135 exit reports
136 end if
137
138 ob % total_obs = ob % total_obs + 1
139
140 select case(fm)
141 case (125) ;
142 ! read SURFACE WinD SPEED AND PRECIPITABLE WATER
143 read (unit=iunit, fmt = fmt_loc) speed, tpw
144 case (126) ;
145 read (unit=iunit, fmt = fmt_loc) &
146 tb19v, tb19h, tb22v, tb37v, tb37h, tb85v, tb85h
147 case default;
148 write(unit=stdout, fmt='(/a/)') &
149 'warning warning warning warning warning :'
150 write(unit=stdout, fmt='(a, i6)') 'unsaved ssmi obs found, fm=', fm
151 write(unit=stdout, fmt='(a, 2f12.6)') &
152 'info%(lon,lat)=', info%lon, info%lat
153 end select
154
155 ! check if obs is in horizontal domain
156 ! ====================================
157
158 ! Compute the model horizontal coordinate x, y
159 ! Check if obs is wihin horizontal domain
160
161 call da_ll_to_xy (info, loc, xp, outside)
162
163 if (outside) cycle reports
164
165 select case(fm)
166 case (125) ;
167 if (.not. Use_SsmiRetrievalObs) cycle reports
168
169 ! Check if at least one field is present
170 if ((tpw % qc == missing) .AND. (speed % qc == missing)) then
171 cycle reports
172 end if
173
174 ob % num_ssmi_retrieval = ob % num_ssmi_retrieval + 1
175
176 case (126) ;
177 if (.not. Use_SsmiTbObs) cycle reports
178
179 ! Check if at least one field is present
180
181 if ((tb19v % qc == missing) .AND. (tb19h % qc == missing) .AND. &
182 (tb22v % qc == missing) .AND. &
183 (tb37v % qc == missing) .AND. (tb37h % qc == missing) .AND. &
184 (tb85v % qc == missing) .AND. (tb85h % qc == missing)) then
185 cycle reports
186 end if
187
188 ! FILTER RAin PIXELS
189 ! ====================================
190
191 if (isfilter) then
192
193 ipass = .false.
194 iprecip = 0
195 call filter(tb19v%inv, tb19h%inv, tb22v%inv, tb37v%inv, &
196 tb37h%inv, tb85v%inv, tb85h%inv, ipass, iprecip, &
197 info%lat)
198 if (iprecip == 1) then
199 irain = irain + 1
200 cycle reports
201 end if
202 end if
203
204 ob % num_ssmi_tb = ob % num_ssmi_tb + 1
205
206 case default;
207 ! Do nothing.
208 end select
209
210 end do reports
211
212 close(unit=iunit)
213 call da_free_unit(iunit)
214
215 write(unit=stdout, fmt='(5x,a,i6,a)') &
216 'Read: ', ob % num_ssmi_retrieval , ' SSMI_RETRIEVAL reports,', &
217 'Read: ', ob % num_ssmi_tb , ' SSMI_Tb reports,', &
218 'Read: ', ob % total_obs, ' Total Observations.'
219
220 write(unit=stdout, fmt='(/,5x,a,i6/)') &
221 '** Rain contaminated SSMI_Tb =', irain
222
223 ob%ob_numb(ob%current_ob_time)%ssmi_retrieval = ob%num_ssmi_retrieval
224 ob%ob_numb(ob%current_ob_time)%ssmi_tb = ob%num_ssmi_tb
225
226 end subroutine da_scan_ssmi
227
228