da_read_obs_radar.inc
References to this file elsewhere.
1 subroutine da_read_obs_radar (iv, filename)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: Read the radar 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 = 120) :: char_total_radar
13 character (len = 120) :: char_ned
14
15 integer :: i, j, n, iost, nlevels, fm
16
17 type (radar_multi_level_type) :: platform
18
19 logical :: outside
20
21 integer :: total_radar
22 integer :: n_dup, ndup, iunit
23 integer :: nlocal(num_ob_indexes)
24 integer :: ntotal(num_ob_indexes)
25
26 if (trace_use) call da_trace_entry("da_read_obs_radar")
27
28 nlocal(:) = iv%info(:)%plocal(iv%time-1)
29 ntotal(:) = iv%info(:)%ptotal(iv%time-1)
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_obs_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 ! 3.5.2 read each level
150
151 do i = 1, platform % info % levels
152 ! height
153 platform%each(i) = radar_each_level_type(missing_r, missing, -1.0, &
154 field_type(missing_r, missing, missing_r), & ! rv
155 field_type(missing_r, missing, missing_r)) ! rf
156
157 read (unit = iunit, fmt = '(3X, F12.1, 2(F12.3,I4,F12.3,2X))') &
158 platform % each (i) % height, &
159 platform % each (i) % rv % inv, &
160 platform % each (i) % rv % qc, &
161 platform % each (i) % rv % error, &
162 platform % each (i) % rf % inv, &
163 platform % each (i) % rf % qc, &
164 platform % each (i) % rf % error
165
166 if (platform % each (i) % rv % error == 0.0) then
167 platform % each (i) % rv % error = 1.0
168 end if
169
170 if (platform % each (i) % rf % error == 0.0) then
171 platform % each (i) % rf % error = 1.0
172 end if
173
174 if (platform % each (i) % rv % inv == missing_r .or. &
175 platform % each (i) % rv % error == missing_r) then
176 platform % each (i) % rv % qc = missing_data
177 end if
178
179 if (platform % each (i) % rf % inv == missing_r .or. &
180 platform % each (i) % rf % error == missing_r) then
181 platform % each (i) % rf % qc = missing_data
182 end if
183 end do
184
185 call da_llxy (platform%info, platform%loc, outside)
186
187 nlevels = platform%info%levels
188
189 if (nlevels > max_ob_levels) then
190 write(unit=message(1),fmt='(A,2I8)') &
191 ' radar=> nlevels > max_ob_levels:',nlevels, max_ob_levels
192 call da_warning(__FILE__,__LINE__,message(1:1))
193 nlevels = max_ob_levels
194 platform%info%levels = nlevels
195 else if (nlevels < 1) then
196 cycle reports
197 end if
198
199 ! Loop over duplicating obs for global
200 n_dup = 1
201 if (global .and. &
202 (platform%loc%i == ids .or. platform%loc%i == ide)) n_dup= 2
203 do ndup = 1, n_dup
204 select case (fm)
205
206 case (128)
207 if (n==1) ntotal(radar) = ntotal(radar) + 1
208 if (outside) cycle reports
209 nlocal(radar) = nlocal(radar) + 1
210 if (nlocal(radar) > iv%info(radar)%nlocal) cycle reports
211
212 iv % radar (nlocal(radar)) % stn_loc % lon = platform % stn % lon
213 iv % radar (nlocal(radar)) % stn_loc % lat = platform % stn % lat
214 iv % radar (nlocal(radar)) % stn_loc % elv = platform % stn % elv
215
216 iv%info(radar)%levels(nlocal(radar)) = nlevels
217 iv%info(radar)%name(nlocal(radar)) = platform%info%name
218 iv%info(radar)%platform(nlocal(radar)) = platform%info%platform
219 iv%info(radar)%id(nlocal(radar)) = platform%info%id
220 iv%info(radar)%date_char(nlocal(radar)) = platform%info%date_char
221 iv%info(radar)%lat(:,nlocal(radar)) = platform%info%lat
222 iv%info(radar)%lon(:,nlocal(radar)) = platform%info%lon
223 iv%info(radar)%elv(nlocal(radar)) = platform%info%elv
224 iv%info(radar)%pstar(nlocal(radar)) = platform%info%pstar
225
226 iv%info(radar)%slp(nlocal(radar)) = platform%loc%slp
227 iv%info(radar)%pw(nlocal(radar)) = platform%loc%pw
228 iv%info(radar)%x(:,nlocal(radar)) = platform%loc%x
229 iv%info(radar)%y(:,nlocal(radar)) = platform%loc%y
230 iv%info(radar)%i(:,nlocal(radar)) = platform%loc%i
231 iv%info(radar)%j(:,nlocal(radar)) = platform%loc%j
232 iv%info(radar)%dx(:,nlocal(radar)) = platform%loc%dx
233 iv%info(radar)%dxm(:,nlocal(radar)) = platform%loc%dxm
234 iv%info(radar)%dy(:,nlocal(radar)) = platform%loc%dy
235 iv%info(radar)%dym(:,nlocal(radar)) = platform%loc%dym
236 iv%info(radar)%proc_domain(:,nlocal(radar)) = platform%loc%proc_domain
237
238 iv%info(radar)%obs_global_index(nlocal(radar)) = ntotal(radar)
239
240 allocate (iv % radar (nlocal(radar)) % model_p (1:nlevels))
241 allocate (iv % radar (nlocal(radar)) % model_rho(1:nlevels))
242 allocate (iv % radar (nlocal(radar)) % model_qrn(1:nlevels))
243 allocate (iv % radar (nlocal(radar)) % height (1:nlevels))
244 allocate (iv % radar (nlocal(radar)) % height_qc(1:nlevels))
245 allocate (iv % radar (nlocal(radar)) % rv (1:nlevels))
246 allocate (iv % radar (nlocal(radar)) % rf (1:nlevels))
247
248 do i = 1, nlevels
249 iv % radar (nlocal(radar)) % height(i) = platform % each(i) % height
250 iv % radar (nlocal(radar)) % height_qc(i) = platform % each(i) % height_qc
251 iv % radar (nlocal(radar)) % rv(i) = platform % each(i) % rv
252 iv % radar (nlocal(radar)) % rf(i) = platform % each(i) % rf
253 end do
254
255 case default;
256 write(unit=message(1), fmt='(a)') 'Unsaved obs found:'
257 write(unit=message(2), fmt='(2a)') &
258 'platform % info % platform=', platform % info % platform
259 write(unit=message(3), fmt='(a, i3)') &
260 'platform % info % levels=', platform % info % levels
261 call da_warning(__FILE__,__LINE__,message(1:3))
262 end select
263
264 if (global .and. ndup == 1) then
265 if (platform%loc % i >= ide) then
266 platform%loc%i = ids
267 platform%loc%proc_domain = .false.
268 else if (platform%loc % i <= ids) then
269 platform%loc%i = ide
270 platform%loc%proc_domain = .false.
271 end if
272 end if
273 end do ! loop over duplicate
274 end do reports
275 end do
276
277 close(iunit)
278 call da_free_unit(iunit)
279
280 if (trace_use) call da_trace_exit("da_read_obs_radar")
281
282
283 end subroutine da_read_obs_radar
284
285