da_read_radar.inc
References to this file elsewhere.
1 subroutine da_read_radar (ob, xp, filename)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: Read the Radar 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
12 character(len=*), intent(in) :: filename
13
14 character (len = 120) :: char_total_radar
15 character (len = 120) :: char_ned
16
17 integer :: i, j, n, iost, nlevels, fm
18
19 type (radar_multi_level_type) :: platform
20
21 logical :: outside
22
23 integer :: total_obs, num_radar, total_radar
24 integer ::n_dup, ndup, iunit
25
26 if (trace_use) call da_trace_entry("da_read_radar")
27
28 total_obs = ob%total_obs
29 num_radar = ob%num_radar
30
31 ! 1. open file
32
33 call da_get_unit(iunit)
34 open(unit = iunit, &
35 FILE = trim(filename), &
36 FORM = 'FORMATTED', &
37 ACCESS = 'SEQUENTIAL', &
38 iostat = iost, &
39 STATUS = 'OLD')
40
41 if (iost /= 0) then
42 ! Missing file does not matter
43 call da_warning(__FILE__,__LINE__, &
44 (/"Cannot open radar file "//trim(filename)/))
45 call da_free_unit(iunit)
46 if (trace_use) call da_trace_exit("da_read_radar")
47 return
48 end if
49
50 ! 2. read total radar
51
52 ! 2.1 read first line
53
54 read (unit=iunit, fmt = '(A)', iostat = iost) char_total_radar
55
56 ! 2.2 process error
57
58 if (iost /= 0) then
59 call da_error(__FILE__,__LINE__, &
60 (/"Cannot read radar file"/))
61 end if
62
63 ! 2.3 total radar number
64
65 read (unit=char_total_radar (15:17),fmt='(I3)', iostat = iost) total_radar
66
67 if (print_detail_radar) write (unit=stdout,fmt='(/,A,I3,/)') &
68 ' TOTAL RADAR: ', total_radar
69
70 ! 2.4 skip one line
71
72 read (unit=iunit, fmt = '(A)', iostat = iost)
73
74 ! 3. read radar data
75
76 do n = 1, total_radar
77
78 ! 3.1 skip one blank line
79
80 read (unit=iunit, fmt = '(A)', iostat = iost)
81
82 ! 3.2 read header
83
84 read (unit=iunit, fmt = '(A)', iostat = iost) char_ned
85
86 ! 3.3 read header information
87
88 read (unit=char_ned (1:5),fmt='(A5)', iostat = iost) platform % stn % platform
89
90 if (print_detail_radar) write (unit=stdout,fmt='(A)') 'RADAR Observations information'
91
92 read (unit=char_ned (8:19),fmt='(A12)', iostat = iost) platform % stn % name
93
94 if (print_detail_radar) write (unit=stdout,fmt='(A,A5,A,A12)') &
95 ' Expect: ',platform % stn % platform, &
96 'data at station:', platform % stn % name
97
98 read (unit=char_ned(20:27),fmt='(F8.3)', iostat = iost) platform % stn % lon
99
100 read (unit=char_ned (30:37),fmt='(F8.3)', iostat = iost) platform % stn % lat
101
102 read (unit=char_ned (40:47),fmt='(F8.1)', iostat = iost) platform % stn % elv
103
104 if (print_detail_radar) write (unit=stdout,fmt='(A,2(F8.3,2X),F8.1)') &
105 'The station longitude, latutude, and altitude are: ', &
106 platform % stn % lon, &
107 platform % stn % lat, platform % stn % elv
108
109 read (unit=char_ned (50:68),fmt='(A19)', iostat = iost) platform % stn % date_char
110
111 if (print_detail_radar) write (unit=stdout,fmt='(A,A19)') &
112 'The observation time for this data is ', &
113 platform % stn % date_char
114
115 read (unit=char_ned (69:74),fmt='(I6)', iostat = iost) platform % stn % numObs
116
117 if (print_detail_radar) write (unit=stdout,fmt='(A,I6)') &
118 'Total number of Super-observations is ', &
119 platform % stn % numobs
120
121 read (unit=char_ned (75:80),fmt='(I6)', iostat = iost) platform % stn % levels
122
123 if (print_detail_radar) write (unit=stdout,fmt='(A,I6)') &
124 'Vertical layers for each Super-observation is ', &
125 platform % stn % levels
126
127 ! 3.4 skip two lines
128
129 read (unit=iunit, fmt = '(A)', iostat = iost)
130 read (unit=iunit, fmt = '(A)', iostat = iost)
131
132 ! 3.5 loop over records
133
134 reports: do j = 1, platform % stn % numobs
135
136 ! 3.5.1 read station general info
137
138 read (unit = iunit, iostat = iost, &
139 fmt = '(A12,3X,A19,2X,2(F12.3,2X),F8.1,2X,I6)') &
140 platform % info % platform, &
141 platform % info % date_char, &
142 platform % info % lat, &
143 platform % info % lon, &
144 platform % info % elv, &
145 platform % info % levels
146
147 read(platform % info % platform (4:6), '(I3)') fm
148
149 total_obs = total_obs + 1
150
151 ! 3.5.2 read each level
152
153 do i = 1, platform % info % levels
154 ! height
155 platform%each(i) = Radar_each_level_type(missing_r, missing, -1.0, &
156 field_type(missing_r, missing, missing_r), & ! rv
157 field_type(missing_r, missing, missing_r)) ! rf
158
159 read (unit = iunit, fmt = '(3X, F12.1, 2(F12.3,I4,F12.3,2X))') &
160 platform % each (i) % height, &
161 platform % each (i) % rv % inv, &
162 platform % each (i) % rv % qc, &
163 platform % each (i) % rv % error, &
164 platform % each (i) % rf % inv, &
165 platform % each (i) % rf % qc, &
166 platform % each (i) % rf % error
167
168 if (platform % each (i) % rv % error == 0.0) then
169 platform % each (i) % rv % error = 1.0
170 end if
171
172 if (platform % each (i) % rf % error == 0.0) then
173 platform % each (i) % rf % error = 1.0
174 end if
175
176 if (platform % each (i) % rv % inv == missing_r .or. &
177 platform % each (i) % rv % error == missing_r) then
178 platform % each (i) % rv % qc = missing_data
179 end if
180
181 if (platform % each (i) % rf % inv == missing_r .or. &
182 platform % each (i) % rf % error == missing_r) then
183 platform % each (i) % rf % qc = missing_data
184 end if
185 end do
186
187 call da_ll_to_xy (platform%info, platform%loc, xp, outside)
188
189 if (outside) then
190 cycle reports
191 end if
192
193 nlevels = platform%info%levels
194
195 if (nlevels > max_ob_levels) then
196 write(unit=message(1),fmt='(A,2I8)') &
197 ' Radar=> nlevels > max_ob_levels:',nlevels, max_ob_levels
198 call da_warning(__FILE__,__LINE__,message(1:1))
199 nlevels = max_ob_levels
200 platform%info%levels = nlevels
201 else if (nlevels < 1) then
202 cycle reports
203 end if
204
205 ! Loop over duplicating obs for global
206 n_dup = 1
207 if (global .and. &
208 (platform%loc%i == xp%ids .or. platform%loc%i == xp%ide)) n_dup= 2
209 do ndup = 1, n_dup
210 select case (fm)
211
212 case (128)
213
214 if (ndup==1) ob%num_radar_glo = ob%num_radar_glo + 1
215
216 num_radar = num_radar + 1
217
218 if (num_radar > max_radar) then
219 write(unit=message(1),fmt='(A,I6,A,I6)') &
220 "Radar #= ",ob % num_radar, " > max_radar = ", max_radar
221 call da_error(__FILE__,__LINE__,message(1:1))
222 end if
223 platform % loc%obs_global_index = ob%num_radar_glo
224
225 ob % radar (num_radar) % stn_loc % lon = platform % stn % lon
226 ob % radar (num_radar) % stn_loc % lat = platform % stn % lat
227 ob % radar (num_radar) % stn_loc % elv = platform % stn % elv
228
229 ob % radar (num_radar) % info = platform % info
230 ob % radar (num_radar) % loc = platform % loc
231
232 allocate (ob % radar (num_radar) % height (1:nlevels))
233 allocate (ob % radar (num_radar) % height_qc(1:nlevels))
234 allocate (ob % radar (num_radar) % zk (1:nlevels))
235 allocate (ob % radar (num_radar) % rv (1:nlevels))
236 allocate (ob % radar (num_radar) % rf (1:nlevels))
237
238 do i = 1, nlevels
239 ob % radar (num_radar) % height(i) = platform % each(i) % height
240 ob % radar (num_radar) % height_qc(i) = platform % each(i) % height_qc
241 ob % radar (num_radar) % zk(i) = platform % each(i) % zk
242 ob % radar (num_radar) % rv(i) = platform % each(i) % rv
243 ob % radar (num_radar) % rf(i) = platform % each(i) % rf
244 end do
245
246 case default;
247 write(unit=message(1), fmt='(a)') 'Unsaved obs found:'
248 write(unit=message(2), fmt='(2a)') &
249 'platform % info % platform=', platform % info % platform
250 write(unit=message(3), fmt='(a, i3)') &
251 'platform % info % levels=', platform % info % levels
252 call da_warning(__FILE__,__LINE__,message(1:3))
253 end select
254
255 if (global .and. ndup == 1) then
256 if (platform%loc % i >= xp % ide) then
257 platform%loc%i = xp % ids
258 platform%loc%proc_domain = .false.
259 else if (platform%loc % i <= xp % ids) then
260 platform%loc%i = xp % ide
261 platform%loc%proc_domain = .false.
262 end if
263 end if
264 end do ! loop over duplicate
265 end do reports
266 end do
267
268 ob % num_radar = num_radar
269
270 ! 4. print out
271
272 if (print_detail_radar) then
273 write(unit=stdout, fmt='(a)') ' '
274 write(unit=stdout, fmt='(5x,a,i6,a)') &
275 ' Read: ', ob % num_radar, ' Radar reports '
276 end if
277
278 if (ob%ob_numb(ob%current_ob_time)%radar /= ob%num_radar) then
279 write(unit=message(1),fmt='(a, i6, 2x, a, i6)') &
280 'Radar obs mismatch ob%ob_numb(ob%current_ob_time)%radar=', &
281 ob%ob_numb(ob%current_ob_time)%radar, &
282 'ob%num_radar=', ob%num_radar
283 call da_error(__FILE__,__LINE__,message(1:1))
284 end if
285
286 close(iunit)
287 call da_free_unit(iunit)
288
289 if (trace_use) call da_trace_exit("da_read_radar")
290
291
292 end subroutine da_read_radar
293
294