da_tune_obs_hollingsworth1.f90
References to this file elsewhere.
1 program da_tune_obs_hollingsworth1
2
3 !-----------------------------------------------------------
4 ! Author: Syed RH Rizvi org: UCAR/NCAR/MMM
5 ! Purpose: Preparing necessary diagnostic output for
6 ! Observation error tuning (Hollingsworh method)
7 ! Ref: Tellus (1986) 38, pp.111-161 (Part I & II)
8 ! Update:
9 ! 01/03/2007 GPS refractivity update Syed RH Rizvi
10 !-----------------------------------------------------------
11
12 use da_control, only : filename_len
13
14 implicit none
15
16 integer, parameter :: y_unit = 50
17 integer, parameter :: unit1 = 35
18 integer, parameter :: unit2 = 36
19 integer, parameter :: unit3 = 37
20 integer, parameter :: unit4 = 38
21 integer, parameter :: unit5 = 39
22
23 integer, parameter :: obs_qc_pointer = 0
24 real, parameter :: missing_r = -888888.0
25
26 character(len=filename_len) :: filename
27 integer :: n, k, current_time
28
29 type info_type
30 character*5 :: id
31 integer :: time
32 real :: lat
33 real :: lon
34 end type info_type
35
36 type field_type
37 real :: yo
38 real :: omb
39 integer :: qc
40 real :: err
41 real :: oma
42 end type field_type
43
44 type surfc_type
45 type (info_type) :: info
46 real :: pressure
47 type (field_type) :: u
48 type (field_type) :: v
49 type (field_type) :: t
50 type (field_type) :: p
51 type (field_type) :: q
52 end type surfc_type
53
54 type qscat_type
55 type (info_type) :: info
56 real :: height
57 type (field_type) :: u
58 type (field_type) :: v
59 end type qscat_type
60
61 type geoamv_type
62 type (info_type) :: info
63 real :: pressure
64 type (field_type) :: u
65 type (field_type) :: v
66 end type geoamv_type
67
68 type polaramv_type
69 type (info_type) :: info
70 character*2 :: channel
71 character*1 :: landmask
72 real :: pressure
73 type (field_type) :: u
74 type (field_type) :: v
75 end type polaramv_type
76
77 type gpspw_type
78 type (info_type) :: info
79 type (field_type) :: tpw
80 end type gpspw_type
81
82 type sound_type
83 type (info_type) :: info
84 integer :: numlevs
85 real, pointer :: pressure(:)
86 type (field_type), pointer :: u(:)
87 type (field_type), pointer :: v(:)
88 type (field_type), pointer :: t(:)
89 type (field_type), pointer :: q(:)
90 end type sound_type
91
92 type airsr_type
93 type (info_type) :: info
94 integer :: numlevs
95 real, pointer :: pressure(:)
96 type (field_type), pointer :: t(:)
97 type (field_type), pointer :: q(:)
98 end type airsr_type
99
100
101 type airep_type
102 type (info_type) :: info
103 integer :: numlevs
104 real, pointer :: pressure(:)
105 type (field_type), pointer :: u(:)
106 type (field_type), pointer :: v(:)
107 type (field_type), pointer :: t(:)
108 end type airep_type
109
110 type pilot_type
111 type (info_type) :: info
112 integer :: numlevs
113 real, pointer :: pressure(:)
114 type (field_type), pointer :: u(:)
115 type (field_type), pointer :: v(:)
116 end type pilot_type
117
118 type ssmir_type
119 type (info_type) :: info
120 type (field_type) :: speed
121 type (field_type) :: tpw
122 end type ssmir_type
123
124 type satem_type
125 type (info_type) :: info
126 integer :: numlevs
127 real, pointer :: pressure(:)
128 type (field_type), pointer :: thickness(:)
129 end type satem_type
130
131 type ssmt1_type
132 type (info_type) :: info
133 integer :: numlevs
134 real, pointer :: height(:)
135 type (field_type), pointer :: t(:)
136 end type ssmt1_type
137
138 type ssmt2_type
139 type (info_type) :: info
140 integer :: numlevs
141 real, pointer :: height(:)
142 type (field_type), pointer :: rh(:)
143 end type ssmt2_type
144 type bogus_type
145 type (info_type) :: info
146 integer :: numlevs
147 real, pointer :: pressure(:)
148 type (field_type), pointer :: u(:)
149 type (field_type), pointer :: v(:)
150 type (field_type), pointer :: t(:)
151 type (field_type), pointer :: q(:)
152 type (field_type) :: slp
153 end type bogus_type
154
155 type gpsref_type
156 type (info_type) :: info
157 integer :: numlevs
158 real, pointer :: pressure(:)
159 type (field_type), pointer :: ref(:)
160 end type gpsref_type
161
162 type iv_type
163 integer :: num_synop, num_metar, num_ships, &
164 num_polaramv, num_qscat, num_geoamv, num_gpspw, num_sound, &
165 num_airep, num_pilot, num_ssmir, num_airsret, &
166 num_satem, num_ssmt1, num_ssmt2, &
167 num_sonde_sfc, num_buoy, num_profiler, num_bogus, num_gpsref
168
169 type (surfc_type), pointer :: synop(:)
170 type (surfc_type), pointer :: metar(:)
171 type (surfc_type), pointer :: ships(:)
172 type (surfc_type), pointer :: sonde_sfc(:)
173 type (surfc_type), pointer :: buoy(:)
174 type (polaramv_type), pointer :: polaramv(:)
175 type (geoamv_type), pointer :: geoamv(:)
176 type (qscat_type), pointer :: qscat(:)
177 type (gpspw_type), pointer :: gpspw(:)
178 type (sound_type), pointer :: sound(:)
179 type (airsr_type), pointer :: airsr(:)
180 type (airep_type), pointer :: airep(:)
181 type (pilot_type), pointer :: pilot(:)
182 type (pilot_type), pointer :: profiler(:)
183 type (ssmir_type), pointer :: ssmir(:)
184 type (satem_type), pointer :: satem(:)
185 type (ssmt1_type), pointer :: ssmt1(:)
186 type (ssmt2_type), pointer :: ssmt2(:)
187 type (bogus_type), pointer :: bogus(:)
188 type (gpsref_type), pointer :: gpsref(:)
189 end type iv_type
190
191 type (iv_type) :: ob
192
193 !--------------------------------------------------------------------------
194 ! [1.0] Count total number of observations and allocate arrays:
195 !--------------------------------------------------------------------------
196
197 filename = 'hollingsworth1.in'
198 open( y_unit, file = filename, status = 'old' )
199
200 call da_count_obs( y_unit, ob )
201
202 !--------------------------------------------------------------------------
203 ! [2.0] Read in observation data:
204 !--------------------------------------------------------------------------
205
206 call da_read_y( y_unit, ob )
207
208 !--------------------------------------------------------------------------
209 ! [3.0] Perform basic statistics for each ob type:
210 !--------------------------------------------------------------------------
211
212 if ( ob % num_synop > 0 ) then
213 call da_calc_stats( 'synop u ', ob % num_synop, ob % synop % u )
214 call da_calc_stats( 'synop v ', ob % num_synop, ob % synop % v )
215 call da_calc_stats( 'synop t ', ob % num_synop, ob % synop % t )
216 call da_calc_stats( 'synop p ', ob % num_synop, ob % synop % p )
217 call da_calc_stats( 'synop q ', ob % num_synop, ob % synop % q )
218 write(6,*)
219 end if
220
221 if ( ob % num_buoy > 0 ) then
222 call da_calc_stats( 'buoy u ', ob % num_buoy, ob % buoy % u )
223 call da_calc_stats( 'buoy v ', ob % num_buoy, ob % buoy % v )
224 call da_calc_stats( 'buoy t ', ob % num_buoy, ob % buoy % t )
225 call da_calc_stats( 'buoy p ', ob % num_buoy, ob % buoy % p )
226 call da_calc_stats( 'buoy q ', ob % num_buoy, ob % buoy % q )
227 write(6,*)
228 end if
229
230 if ( ob % num_sonde_sfc > 0 ) then
231 call da_calc_stats( 'sonde_sfc u', ob % num_sonde_sfc, ob % sonde_sfc % u )
232 call da_calc_stats( 'sonde_sfc v', ob % num_sonde_sfc, ob % sonde_sfc % v )
233 call da_calc_stats( 'sonde_sfc t', ob % num_sonde_sfc, ob % sonde_sfc % t )
234 call da_calc_stats( 'sonde_sfc p', ob % num_sonde_sfc, ob % sonde_sfc % p )
235 call da_calc_stats( 'sonde_sfc q', ob % num_sonde_sfc, ob % sonde_sfc % q )
236 write(6,*)
237 end if
238
239 if ( ob % num_metar > 0 ) then
240 call da_calc_stats( 'metar u ', ob % num_metar, ob % metar % u )
241 call da_calc_stats( 'metar v ', ob % num_metar, ob % metar % v )
242 call da_calc_stats( 'metar t ', ob % num_metar, ob % metar % t )
243 call da_calc_stats( 'metar p ', ob % num_metar, ob % metar % p )
244 call da_calc_stats( 'metar q ', ob % num_metar, ob % metar % q )
245 write(6,*)
246 end if
247
248 if ( ob % num_ships > 0 ) then
249 call da_calc_stats( 'ships u ', ob % num_ships, ob % ships % u )
250 call da_calc_stats( 'ships v ', ob % num_ships, ob % ships % v )
251 call da_calc_stats( 'ships t ', ob % num_ships, ob % ships % t )
252 call da_calc_stats( 'ships p ', ob % num_ships, ob % ships % p )
253 call da_calc_stats( 'ships q ', ob % num_ships, ob % ships % q )
254 write(6,*)
255 end if
256
257 if ( ob % num_polaramv > 0 ) then
258 call da_calc_stats( 'polaramv u ', ob % num_polaramv, ob % polaramv % u )
259 call da_calc_stats( 'polaramv v ', ob % num_polaramv, ob % polaramv % v )
260 write(6,*)
261 end if
262
263 if ( ob % num_geoamv > 0 ) then
264 call da_calc_stats( 'geoamv u ', ob % num_geoamv, ob % geoamv % u )
265 call da_calc_stats( 'geoamv v ', ob % num_geoamv, ob % geoamv % v )
266 write(6,*)
267 end if
268
269 if ( ob % num_qscat > 0 ) then
270 call da_calc_stats( 'qscat u ', ob % num_qscat, ob % qscat % u )
271 call da_calc_stats( 'qscat v ', ob % num_qscat, ob % qscat % v )
272 write(6,*)
273 end if
274
275 if ( ob % num_gpspw > 0 ) then
276 call da_calc_stats( 'gpspw tpw', ob % num_gpspw, ob % gpspw % tpw)
277 write(6,*)
278 end if
279
280 if ( ob % num_ssmir > 0 ) then
281 call da_calc_stats( 'ssmir spd', ob % num_ssmir, ob % ssmir % speed)
282 call da_calc_stats( 'ssmir tpw', ob % num_ssmir, ob % ssmir % tpw)
283 write(6,*)
284 end if
285
286 if ( ob % num_sound > 0 ) call da_calc_stats_sound( ob % num_sound, ob % sound )
287 if ( ob % num_airsret> 0 ) call da_calc_stats_airsr( ob % num_airsret, ob % airsr )
288 if ( ob % num_airep > 0 ) call da_calc_stats_airep( ob % num_airep, ob % airep )
289 if ( ob % num_pilot > 0 ) call da_calc_stats_pilot( ob % num_pilot, ob % pilot )
290
291 if ( ob % num_profiler > 0 ) call da_calc_stats_profiler( ob % num_profiler, ob % profiler )
292 if ( ob % num_satem > 0 ) call da_calc_stats_satem( ob % num_satem, ob % satem )
293 if ( ob % num_ssmt1 > 0 ) call da_calc_stats_ssmt1( ob % num_ssmt1, ob % ssmt1 )
294 if ( ob % num_ssmt2 > 0 ) call da_calc_stats_ssmt2( ob % num_ssmt2, ob % ssmt2 )
295 if ( ob % num_bogus > 0 ) call da_calc_stats_bogus( ob % num_bogus, ob % bogus )
296 if ( ob % num_gpsref > 0 ) call da_calc_stats_gpsref( ob % num_gpsref, ob % gpsref )
297
298 !--------------------------------------------------------------------------
299 ! [4.0] Write data for post-processing:
300 !--------------------------------------------------------------------------
301 ! [4.1] Sonde O-B:
302 if ( ob % num_sound > 0 ) then
303 open( unit1, file = 'sound_u_omb.dat', status = 'unknown' )
304 open( unit2, file = 'sound_v_omb.dat', status = 'unknown' )
305 open( unit3, file = 'sound_t_omb.dat', status = 'unknown' )
306 open( unit4, file = 'sound_q_omb.dat', status = 'unknown' )
307 current_time = 1
308 do n = 1, ob % num_sound
309 do k = 1, ob % sound(n) % numlevs
310 call da_write_data( k, current_time, unit1, ob % sound(n) % info, &
311 ob % sound(n) % pressure(k), ob % sound(n) % u(k) )
312 call da_write_data( k, current_time, unit2, ob % sound(n) % info, &
313 ob % sound(n) % pressure(k), ob % sound(n) % v(k) )
314 call da_write_data( k, current_time, unit3, ob % sound(n) % info, &
315 ob % sound(n) % pressure(k), ob % sound(n) % t(k) )
316 call da_write_data( k, current_time, unit4, ob % sound(n) % info, &
317 ob % sound(n) % pressure(k), ob % sound(n) % q(k) )
318
319 end do
320 current_time = ob % sound(n) % info % time
321 end do
322
323 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
324 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
325 write(unit3,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
326 write(unit4,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
327
328 close( unit1 ); close( unit2 ); close( unit3 ); close( unit4 )
329
330 end if
331 ! [4.2] Synop O-B:
332 if ( ob % num_synop > 0 ) then
333
334 open( unit1, file = 'synop_u_omb.dat', status = 'unknown' )
335 open( unit2, file = 'synop_v_omb.dat', status = 'unknown' )
336 open( unit3, file = 'synop_t_omb.dat', status = 'unknown' )
337 open( unit4, file = 'synop_p_omb.dat', status = 'unknown' )
338 open( unit5, file = 'synop_q_omb.dat', status = 'unknown' )
339
340 current_time = 1
341 do n = 1, ob % num_synop
342 call da_write_data( 1, current_time, unit1, ob % synop(n) % info, &
343 ob % synop(n) % pressure, ob % synop(n) % u )
344 call da_write_data( 1, current_time, unit2, ob % synop(n) % info, &
345 ob % synop(n) % pressure, ob % synop(n) % v )
346 call da_write_data( 1, current_time, unit3, ob % synop(n) % info, &
347 ob % synop(n) % pressure, ob % synop(n) % t )
348 call da_write_data( 1, current_time, unit4, ob % synop(n) % info, &
349 ob % synop(n) % pressure, ob % synop(n) % p )
350 call da_write_data( 1, current_time, unit5, ob % synop(n) % info, &
351 ob % synop(n) % pressure, ob % synop(n) % q )
352 current_time = ob % synop(n) % info % time
353 end do
354
355 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
356 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
357 write(unit3,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
358 write(unit4,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
359 write(unit5,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
360
361 close( unit1 ); close( unit2 ); close( unit3 ); close( unit4 ); close( unit5 )
362 end if
363 ! [4.3] Metar O-B:
364 if ( ob % num_metar > 0 ) then
365
366 open( unit1, file = 'metar_u_omb.dat', status = 'unknown' )
367 open( unit2, file = 'metar_v_omb.dat', status = 'unknown' )
368 open( unit3, file = 'metar_t_omb.dat', status = 'unknown' )
369 open( unit4, file = 'metar_p_omb.dat', status = 'unknown' )
370 open( unit5, file = 'metar_q_omb.dat', status = 'unknown' )
371
372 current_time = 1
373 do n = 1, ob % num_metar
374 call da_write_data( 1, current_time, unit1, ob % metar(n) % info, &
375 ob % metar(n) % pressure, ob % metar(n) % u )
376 call da_write_data( 1, current_time, unit2, ob % metar(n) % info, &
377 ob % metar(n) % pressure, ob % metar(n) % v )
378 call da_write_data( 1, current_time, unit3, ob % metar(n) % info, &
379 ob % metar(n) % pressure, ob % metar(n) % t )
380 call da_write_data( 1, current_time, unit4, ob % metar(n) % info, &
381 ob % metar(n) % pressure, ob % metar(n) % p )
382 call da_write_data( 1, current_time, unit5, ob % metar(n) % info, &
383 ob % metar(n) % pressure, ob % metar(n) % q )
384 current_time = ob % metar(n) % info % time
385 end do
386
387 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
388 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
389 write(unit3,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
390 write(unit4,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
391 write(unit5,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
392
393 close( unit1 ); close( unit2 ); close( unit3 ); close( unit4 ); close( unit5 )
394 end if
395 ! [4.4] Polar AMV O-B:
396 if ( ob % num_polaramv > 0 ) then
397 open( unit1, file = 'polaramv_u_omb.dat', status = 'unknown' )
398 open( unit2, file = 'polaramv_v_omb.dat', status = 'unknown' )
399
400 current_time = 1
401 do n = 1, ob % num_polaramv
402 call da_write_data( 1, current_time, unit1, ob % polaramv(n) % info, &
403 ob % polaramv(n) % pressure, ob % polaramv(n) % u )
404 call da_write_data( 1, current_time, unit2, ob % polaramv(n) % info, &
405 ob % polaramv(n) % pressure, ob % polaramv(n) % v )
406 current_time = ob % polaramv(n) % info % time
407 end do
408
409 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
410 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
411
412 close( unit1 ); close( unit2 )
413 end if
414 ! [4.5] Geo AMV O-B:
415 if ( ob % num_geoamv > 0 ) then
416
417 open( unit1, file = 'geoamv_u_omb.dat', status = 'unknown' )
418 open( unit2, file = 'geoamvv_omb.dat', status = 'unknown' )
419
420 current_time = 1
421 do n = 1, ob % num_geoamv
422 call da_write_data( 1, current_time, unit1, ob % geoamv(n) % info, &
423 ob % geoamv(n) % pressure, ob % geoamv(n) % u )
424 call da_write_data( 1, current_time, unit2, ob % geoamv(n) % info, &
425 ob % geoamv(n) % pressure, ob % geoamv(n) % v )
426 current_time = ob % geoamv(n) % info % time
427 end do
428
429 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
430 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
431
432 close( unit1 ); close( unit2 )
433
434 end if
435 ! [4.6] Buoy O-B:
436
437 if ( ob % num_buoy > 0 ) then
438 open( unit1, file = 'buoy_u_omb.dat', status = 'unknown' )
439 open( unit2, file = 'buoy_v_omb.dat', status = 'unknown' )
440 open( unit3, file = 'buoy_t_omb.dat', status = 'unknown' )
441 open( unit4, file = 'buoy_p_omb.dat', status = 'unknown' )
442 open( unit5, file = 'buoy_q_omb.dat', status = 'unknown' )
443
444 current_time = 1
445 do n = 1, ob % num_buoy
446 call da_write_data( 1, current_time, unit1, ob % buoy(n) % info, &
447 ob % buoy(n) % pressure, ob % buoy(n) % u )
448 call da_write_data( 1, current_time, unit2, ob % buoy(n) % info, &
449 ob % buoy(n) % pressure, ob % buoy(n) % v )
450 call da_write_data( 1, current_time, unit3, ob % buoy(n) % info, &
451 ob % buoy(n) % pressure, ob % buoy(n) % t )
452 call da_write_data( 1, current_time, unit4, ob % buoy(n) % info, &
453 ob % buoy(n) % pressure, ob % buoy(n) % p )
454 call da_write_data( 1, current_time, unit5, ob % buoy(n) % info, &
455 ob % buoy(n) % pressure, ob % buoy(n) % q )
456 current_time = ob % buoy(n) % info % time
457 end do
458
459 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
460 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
461 write(unit3,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
462 write(unit4,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
463 write(unit5,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
464
465 close( unit1 ); close( unit2 ); close( unit3 ); close( unit4 ); close( unit5 )
466 end if
467
468 ! [4.7] sonde_sfc O-B:
469
470 if ( ob % num_sonde_sfc > 0 ) then
471 open( unit1, file = 'sonde_sfc_u_omb.dat', status = 'unknown' )
472 open( unit2, file = 'sonde_sfc_v_omb.dat', status = 'unknown' )
473 open( unit3, file = 'sonde_sfc_t_omb.dat', status = 'unknown' )
474 open( unit4, file = 'sonde_sfc_p_omb.dat', status = 'unknown' )
475 open( unit5, file = 'sonde_sfc_q_omb.dat', status = 'unknown' )
476
477 current_time = 1
478 do n = 1, ob % num_sonde_sfc
479 call da_write_data( 1, current_time, unit1, ob % sonde_sfc(n) % info, &
480 ob % sonde_sfc(n) % pressure, ob % sonde_sfc(n) % u )
481 call da_write_data( 1, current_time, unit2, ob % sonde_sfc(n) % info, &
482 ob % sonde_sfc(n) % pressure, ob % sonde_sfc(n) % v )
483 call da_write_data( 1, current_time, unit3, ob % sonde_sfc(n) % info, &
484 ob % sonde_sfc(n) % pressure, ob % sonde_sfc(n) % t )
485 call da_write_data( 1, current_time, unit4, ob % sonde_sfc(n) % info, &
486 ob % sonde_sfc(n) % pressure, ob % sonde_sfc(n) % p )
487 call da_write_data( 1, current_time, unit5, ob % sonde_sfc(n) % info, &
488 ob % sonde_sfc(n) % pressure, ob % sonde_sfc(n) % q )
489 current_time = ob % sonde_sfc(n) % info % time
490 end do
491
492 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
493 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
494 write(unit3,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
495 write(unit4,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
496 write(unit5,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
497
498 close( unit1 ); close( unit2 ); close( unit3 ); close( unit4 ); close( unit5 )
499 end if
500 ! [4.8] Profiler O-B:
501 if ( ob % num_profiler > 0 ) then
502 open( unit1, file = 'profiler_u_omb.dat', status = 'unknown' )
503 open( unit2, file = 'profiler_v_omb.dat', status = 'unknown' )
504
505 current_time = 1
506 do n = 1, ob % num_profiler
507 do k = 1, ob % profiler(n) % numlevs
508 call da_write_data( k, current_time, unit1, ob % profiler(n) % info, &
509 ob % profiler(n) % pressure(k), ob % profiler(n) % u(k) )
510 call da_write_data( k, current_time, unit2, ob % profiler(n) % info, &
511 ob % profiler(n) % pressure(k), ob % profiler(n) % v(k) )
512 end do
513 current_time = ob % profiler(n) % info % time
514 end do
515
516 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
517 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
518
519 close( unit1 ); close( unit2 )
520 end if
521 ! [4.9] AIRS retrievals O-B:
522 if ( ob % num_airsret > 0 ) then
523 open( unit1, file = 'airsr_t_omb.dat', status = 'unknown' )
524 open( unit2, file = 'airsr_q_omb.dat', status = 'unknown' )
525 current_time = 1
526 do n = 1, ob % num_airsret
527 do k = 1, ob % airsr(n) % numlevs
528 call da_write_data( k, current_time, unit1, ob % airsr(n) % info, &
529 ob % airsr(n) % pressure(k), ob % airsr(n) % t(k) )
530 call da_write_data( k, current_time, unit2, ob % airsr(n) % info, &
531 ob % airsr(n) % pressure(k), ob % airsr(n) % q(k) )
532
533 end do
534 current_time = ob % airsr(n) % info % time
535 end do
536
537 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
538 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
539
540 close( unit1 ); close( unit2 )
541 end if
542 ! [4.10] Pilot O-B:
543 if ( ob % num_pilot > 0 ) then
544 open( unit1, file = 'pilot_u_omb.dat', status = 'unknown' )
545 open( unit2, file = 'pilot_v_omb.dat', status = 'unknown' )
546
547 current_time = 1
548 do n = 1, ob % num_pilot
549 do k = 1, ob % pilot(n) % numlevs
550 call da_write_data( k, current_time, unit1, ob % pilot(n) % info, &
551 ob % pilot(n) % pressure(k), ob % pilot(n) % u(k) )
552 call da_write_data( k, current_time, unit2, ob % pilot(n) % info, &
553 ob % pilot(n) % pressure(k), ob % pilot(n) % v(k) )
554 end do
555 current_time = ob % pilot(n) % info % time
556 end do
557
558 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
559 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
560
561 close( unit1 ); close( unit2 )
562 end if
563 !--------------------------------------------------------------------------
564 ! [4.11] For Bogus O-B:
565 if ( ob % num_bogus > 0 ) then
566 open( unit1, file = 'bogus_u_omb.dat', status = 'unknown' )
567 open( unit2, file = 'bogus_v_omb.dat', status = 'unknown' )
568 open( unit3, file = 'bogus_t_omb.dat', status = 'unknown' )
569 open( unit4, file = 'bogus_q_omb.dat', status = 'unknown' )
570 open( unit5, file = 'bogus_slp_omb.dat', status = 'unknown' )
571 current_time = 1
572 do n = 1, ob % num_bogus
573 do k = 1, ob % bogus(n) % numlevs
574 call da_write_data( k, current_time, unit1, ob % bogus(n) % info, &
575 ob % bogus(n) % pressure(k), ob % bogus(n) % u(k) )
576 call da_write_data( k, current_time, unit2, ob % bogus(n) % info, &
577 ob % bogus(n) % pressure(k), ob % bogus(n) % v(k) )
578 call da_write_data( k, current_time, unit3, ob % bogus(n) % info, &
579 ob % bogus(n) % pressure(k), ob % bogus(n) % t(k) )
580 call da_write_data( k, current_time, unit4, ob % bogus(n) % info, &
581 ob % bogus(n) % pressure(k), ob % bogus(n) % q(k) )
582
583 end do
584 call da_write_data( 1, current_time, unit5, ob % bogus(n) % info, &
585 missing_r, ob % bogus(n) % slp )
586 current_time = ob % bogus(n) % info % time
587 end do
588
589 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
590 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
591 write(unit3,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
592 write(unit4,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
593 write(unit5,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
594
595 close( unit1 ); close( unit2 ); close( unit3 ); close( unit4 ); close( unit5 )
596
597 end if
598 !--------------------------------------------------------------------------
599 ! [4.11] For Airep O-B:
600 if ( ob % num_airep > 0 ) then
601 open( unit1, file = 'airep_u_omb.dat', status = 'unknown' )
602 open( unit2, file = 'airep_v_omb.dat', status = 'unknown' )
603 open( unit3, file = 'airep_t_omb.dat', status = 'unknown' )
604 current_time = 1
605 do n = 1, ob % num_airep
606 do k = 1, ob % airep(n) % numlevs
607 call da_write_data( k, current_time, unit1, ob % airep(n) % info, &
608 ob % airep(n) % pressure(k), ob % airep(n) % u(k) )
609 call da_write_data( k, current_time, unit2, ob % airep(n) % info, &
610 ob % airep(n) % pressure(k), ob % airep(n) % v(k) )
611 call da_write_data( k, current_time, unit3, ob % airep(n) % info, &
612 ob % airep(n) % pressure(k), ob % airep(n) % t(k) )
613 end do
614 current_time = ob % airep(n) % info % time
615 end do
616
617 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
618 write(unit2,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
619 write(unit3,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
620
621 close( unit1 ); close( unit2 ); close( unit3 )
622
623 end if
624 ! [4.12] gpsref O-B:
625 if ( ob % num_gpsref > 0 ) then
626 open( unit1, file = 'gpsref_ref_omb.dat', status = 'unknown' )
627
628 current_time = 1
629 do n = 1, ob % num_gpsref
630 do k = 1, ob % gpsref(n) % numlevs
631 call da_write_data( k, current_time, unit1, ob % gpsref(n) % info, &
632 ob % gpsref(n) % pressure(k), ob % gpsref(n) % ref(k) )
633 end do
634 current_time = ob % gpsref(n) % info % time
635 end do
636
637 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
638
639 close( unit1 )
640 end if
641
642 ! [4.13] Gpspw O-B:
643 if ( ob % num_gpspw > 0 ) then
644 open( unit1, file = 'gpspw_tpw_omb.dat', status = 'unknown' )
645
646 current_time = 1
647 do n = 1, ob % num_gpspw
648 call da_write_data( 1, current_time, unit1, ob % gpspw(n) % info, &
649 missing_r, ob % gpspw(n) % tpw )
650 current_time = ob % gpspw(n) % info % time
651 end do
652
653 write(unit1,'(a5,2f9.3,3f17.7,i8)')'*end*', 0., 0., 0., 0., 0., 0
654
655 close( unit1 )
656
657 endif
658 contains
659
660 subroutine da_write_data( lev, current_time, ounit, info, p, field )
661
662 implicit none
663
664 integer, intent(in) :: lev
665 integer, intent(in) :: current_time
666 integer, intent(in) :: ounit
667 type (info_type), intent(in) :: info
668 real, intent(in) :: p
669 type (field_type), intent(in) :: field
670
671 if ( info % time == current_time ) then ! New ob at same time:
672 write(ounit,'(a5,2f9.3,3f17.7,i8)') info % id, info % lat, info % lon, &
673 p, field % omb, field % err, field % qc
674 else
675 if ( lev == 1)write(ounit,'(a5,2f9.3,3f17.7,i8)')'*****', 0., 0., 0., 0., 0., 0
676 write(ounit,'(a5,2f9.3,3f17.7,i8)') info % id, info % lat, info % lon, &
677 p, field % omb, field % err, field % qc
678 end if
679
680 end subroutine da_write_data
681
682 !--------------------------------------------------------------------------
683
684 subroutine da_count_obs( y_unit, ob )
685
686 implicit none
687
688 integer, intent(in) :: y_unit
689 type (iv_type), intent(inout) :: ob
690
691 character*20 :: ob_name, dummy
692 integer :: num_obs, num_times, levels,k, kk
693
694 ! [1] Initialize ob numbers:
695
696 ob % num_synop = 0
697 ob % num_metar = 0
698 ob % num_ships = 0
699 ob % num_polaramv = 0
700 ob % num_geoamv = 0
701 ob % num_qscat = 0
702 ob % num_gpspw = 0
703 ob % num_sound = 0
704 ob % num_airsret = 0
705 ob % num_airep = 0
706 ob % num_pilot = 0
707 ob % num_ssmir = 0
708 ob % num_satem = 0
709 ob % num_ssmt1 = 0
710 ob % num_ssmt2 = 0
711 ob % num_sonde_sfc = 0
712 ob % num_buoy = 0
713 ob % num_profiler = 0
714 ob % num_bogus = 0
715 ob % num_gpsref = 0
716
717 num_times = 0
718
719 ! [2] Loop through input file to count number of obs:
720
721 do ! loop over entire input file: (ends in *end*)
722 do ! loop over particular time (ends in *****)
723
724 read(y_unit,'(a20,i8)')ob_name, num_obs
725
726 if ( index( ob_name,'*****') > 0 .or. index( ob_name,'*end*') > 0 ) exit
727
728 if ( index( ob_name,'synop') > 0 ) then
729 ob % num_synop = ob % num_synop + num_obs
730 else if ( index( ob_name,'metar') > 0 ) then
731 ob % num_metar = ob % num_metar + num_obs
732 else if ( index( ob_name,'ships') > 0 ) then
733 ob % num_ships = ob % num_ships + num_obs
734 else if ( index( ob_name,'polaramv') > 0 ) then
735 ob % num_polaramv = ob % num_polaramv + num_obs
736 else if ( index( ob_name,'geoamv') > 0 ) then
737 ob % num_geoamv = ob % num_geoamv + num_obs
738 else if ( index( ob_name,'qscat') > 0 ) then
739 ob % num_qscat = ob % num_qscat + num_obs
740 else if ( index( ob_name,'gpspw') > 0 ) then
741 ob % num_gpspw = ob % num_gpspw + num_obs
742 else if ( index( ob_name,'sound') > 0 ) then
743 ob % num_sound = ob % num_sound + num_obs
744 else if ( index( ob_name,'airsr') > 0 ) then
745 ob % num_airsret = ob % num_airsret + num_obs
746 else if ( index( ob_name,'airep') > 0 ) then
747 ob % num_airep = ob % num_airep + num_obs
748 else if ( index( ob_name,'pilot') > 0 ) then
749 ob % num_pilot = ob % num_pilot + num_obs
750 else if ( index( ob_name,'ssmir') > 0 ) then
751 ob % num_ssmir = ob % num_ssmir + num_obs
752 else if ( index( ob_name,'satem') > 0 ) then
753 ob % num_satem = ob % num_satem + num_obs
754 else if ( index( ob_name,'ssmt1') > 0 ) then
755 ob % num_ssmt1 = ob % num_ssmt1 + num_obs
756 else if ( index( ob_name,'ssmt2') > 0 ) then
757 ob % num_ssmt2 = ob % num_ssmt2 + num_obs
758 else if ( index( ob_name,'sonde_sfc') > 0 ) then
759 ob % num_sonde_sfc = ob % num_sonde_sfc + num_obs
760 else if ( index( ob_name,'buoy') > 0 ) then
761 ob % num_buoy = ob % num_buoy + num_obs
762 else if ( index( ob_name,'profiler') > 0 ) then
763 ob % num_profiler = ob % num_profiler + num_obs
764 else if ( index( ob_name,'bogus') > 0 ) then
765 ob % num_bogus = ob % num_bogus + num_obs
766 else if ( index( ob_name,'gpsref') > 0 ) then
767 ob % num_gpsref = ob % num_gpsref + num_obs
768 else
769 print*,' unknown obs type: ',trim(ob_name),' found in diagnostics.in file '
770 end if
771
772 do kk = 1, num_obs
773 if( index( ob_name,'bogus') > 0 ) then
774 read(y_unit,'(i8)') levels
775 read(y_unit,'(a20)') dummy
776 end if
777 read(y_unit,'(i8)') levels
778 do k=1,levels
779 read(y_unit,'(a20)') dummy
780 end do
781 end do
782 end do
783
784 if ( index( ob_name,'*end*') > 0 ) then
785 exit
786 else
787 num_times = num_times + 1
788 end if
789 end do
790
791 write(6,'(a,i8)')' Number of times read= ', num_times
792
793 ! [3] Allocate ob structures where obs exist:
794
795 if ( ob % num_synop > 0 ) then
796 allocate( ob % synop(1:ob % num_synop) )
797 write(6,'(a,i8)')' Number of synop obs = ', ob % num_synop
798 end if
799
800 if ( ob % num_metar > 0 ) then
801 allocate( ob % metar(1:ob % num_metar) )
802 write(6,'(a,i8)')' Number of metar obs = ', ob % num_metar
803 end if
804
805 if ( ob % num_ships > 0 ) then
806 allocate( ob % ships(1:ob % num_ships) )
807 write(6,'(a,i8)')' Number of ships obs = ', ob % num_ships
808 end if
809
810 if ( ob % num_polaramv > 0 ) then
811 allocate( ob % polaramv(1:ob % num_polaramv) )
812 write(6,'(a,i8)')' Number of polaramv obs = ', ob % num_polaramv
813 end if
814
815 if ( ob % num_geoamv > 0 ) then
816 allocate( ob % geoamv(1:ob % num_geoamv) )
817 write(6,'(a,i8)')' Number of geoamv obs = ', ob % num_geoamv
818 end if
819
820 if ( ob % num_qscat > 0 ) then
821 allocate( ob % qscat(1:ob % num_qscat) )
822 write(6,'(a,i8)')' Number of qscat obs = ', ob % num_qscat
823 end if
824
825 if ( ob % num_gpspw > 0 ) then
826 allocate( ob % gpspw(1:ob % num_gpspw) )
827 write(6,'(a,i8)')' Number of gpspw obs = ', ob % num_gpspw
828 end if
829
830 if ( ob % num_sound > 0 ) then
831 allocate( ob % sound(1:ob % num_sound) )
832 write(6,'(a,i8)')' Number of sound obs = ', ob % num_sound
833 end if
834
835 if ( ob % num_airsret > 0 ) then
836 allocate( ob % airsr(1:ob % num_airsret) )
837 write(6,'(a,i8)')' Number of AIRS retrievals = ', ob % num_airsret
838 end if
839
840 if ( ob % num_airep > 0 ) then
841 allocate( ob % airep(1:ob % num_airep) )
842 write(6,'(a,i8)')' Number of airep obs = ', ob % num_airep
843 end if
844
845 if ( ob % num_pilot > 0 ) then
846 allocate( ob % pilot(1:ob % num_pilot) )
847 write(6,'(a,i8)')' Number of pilot obs = ', ob % num_pilot
848 end if
849
850 if ( ob % num_ssmir > 0 ) then
851 allocate( ob % ssmir(1:ob % num_ssmir) )
852 write(6,'(a,i8)')' Number of ssmir obs = ', ob % num_ssmir
853 end if
854
855 if ( ob % num_satem > 0 ) then
856 allocate( ob % satem(1:ob % num_satem) )
857 write(6,'(a,i8)')' Number of satem obs = ', ob % num_satem
858 end if
859
860 if ( ob % num_ssmt1 > 0 ) then
861 allocate( ob % ssmt1(1:ob % num_ssmt1) )
862 write(6,'(a,i8)')' Number of ssmt1 obs = ', ob % num_ssmt1
863 end if
864
865 if ( ob % num_ssmt2 > 0 ) then
866 allocate( ob % ssmt2(1:ob % num_ssmt2) )
867 write(6,'(a,i8)')' Number of ssmt2 obs = ', ob % num_ssmt2
868 end if
869
870 if ( ob % num_buoy > 0 ) then
871 allocate( ob % buoy(1:ob % num_buoy) )
872 write(6,'(a,i8)')' Number of buoy obs = ', ob % num_buoy
873 end if
874
875 if ( ob % num_sonde_sfc > 0 ) then
876 allocate( ob % sonde_sfc(1:ob % num_sonde_sfc) )
877 write(6,'(a,i8)')' Number of sonde_sfc obs = ', ob % num_sonde_sfc
878 end if
879
880 if ( ob % num_profiler > 0 ) then
881 allocate( ob % profiler(1:ob % num_profiler) )
882 write(6,'(a,i8)')' Number of profiler obs = ', ob % num_profiler
883 end if
884 if ( ob % num_gpsref > 0 ) then
885 allocate( ob % gpsref(1:ob % num_gpsref) )
886 write(6,'(a,i8)')' Number of gpsref obs = ', ob % num_gpsref
887 end if
888
889
890 if ( ob % num_bogus > 0 ) then
891 allocate( ob % bogus(1:ob % num_bogus) )
892 write(6,'(a,i8)')' Number of bogus obs = ', ob % num_bogus
893 end if
894
895 write(6,*)
896
897 end subroutine da_count_obs
898 subroutine da_read_y( y_unit, ob )
899
900 implicit none
901
902 integer, intent(in) :: y_unit
903 type (iv_type), intent(inout) :: ob
904
905 character*20 :: ob_name, dummy
906 integer :: n, ndum, k, kdum, num_obs, num_levs
907 integer :: num_obs_sonde_sfc, num_obs_buoy, num_obs_profiler, num_obs_bogus
908 integer :: num_obs_synop, num_obs_metar, num_obs_ships, num_obs_gpsref
909 integer :: num_obs_qscat, num_obs_polaramv, num_obs_geoamv
910 integer :: num_obs_gpspw, num_obs_sound, num_obs_airsr, num_obs_airep, num_obs_pilot
911 integer :: num_obs_ssmir, num_obs_satem, num_obs_ssmt1, num_obs_ssmt2
912 integer :: synopt, metart, shipst, polaramvt, geoamvt, qscatt, gpspwt, soundt, airsrt
913 integer :: sonde_sfct, buoyt, profilert, bogust, gpsreft
914 integer :: airept, pilott, ssmirt, satemt, ssmt1t, ssmt2t
915 real :: rdum
916
917 rewind (y_unit)
918 num_obs_buoy = 0; num_obs_sonde_sfc = 0; num_obs_profiler = 0 ; num_obs_bogus = 0
919 num_obs_synop = 0; num_obs_metar = 0; num_obs_ships = 0 ; num_obs_gpsref = 0
920 num_obs_polaramv = 0; num_obs_geoamv = 0; num_obs_qscat = 0
921 num_obs_gpspw = 0; num_obs_sound = 0; num_obs_airsr = 0; num_obs_airep = 0; num_obs_pilot = 0
922 num_obs_ssmir = 0; num_obs_satem = 0; num_obs_ssmt1 = 0; num_obs_ssmt2 = 0
923 buoyt = 0; sonde_sfct = 0; profilert = 0; bogust = 0 ; gpsreft = 0
924 synopt = 0; metart = 0; shipst = 0; polaramvt = 0; geoamvt = 0; qscatt = 0
925 gpspwt = 0; soundt = 0; airsrt = 0; airept = 0; pilott = 0
926 ssmirt = 0; satemt = 0; ssmt1t = 0; ssmt2t = 0
927 !
928 do ! loop over entire input file: (ends in *end*)
929 do ! loop over particular time (ends in *****)
930
931 read(y_unit,'(a20,i8)')ob_name, num_obs
932 if ( index( ob_name,'*****') > 0 .or. index( ob_name,'*end*') > 0 ) exit
933
934 if ( index( ob_name,'synop') > 0 ) then
935 synopt = synopt + 1
936 do n = num_obs_synop + 1, num_obs_synop + num_obs
937
938 read(y_unit,'(i8)') num_levs
939 do k = 1, num_levs
940 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
941 ndum, kdum, ob % synop(n) % info % id, & ! Station
942 ob % synop(n) % info % lat, & ! Latitude
943 ob % synop(n) % info % lon, & ! Longitude
944 ob % synop(n) % pressure, & ! Obs height
945 ob % synop(n) % u, & ! O, O-B, O-A
946 ob % synop(n) % v, & ! O, O-B, O-A
947 ob % synop(n) % t, & ! O, O-B, O-A
948 ob % synop(n) % p, & ! O, O-B, O-A
949 ob % synop(n) % q ! O, O-B, O-A
950 end do
951 ob % synop(n) % info % time = synopt
952 end do
953 num_obs_synop = num_obs_synop + num_obs
954
955 elseif ( index( ob_name,'metar') > 0 ) then
956 metart = metart + 1
957 do n = num_obs_metar + 1, num_obs_metar + num_obs
958 read(y_unit,'(i8)') num_levs
959 do k = 1, num_levs
960 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
961 ndum, kdum, ob % metar(n) % info % id, & ! Station
962 ob % metar(n) % info % lat, & ! Latitude
963 ob % metar(n) % info % lon, & ! Longitude
964 ob % metar(n) % pressure, & ! Obs height
965 ob % metar(n) % u, & ! O, O-B, O-A
966 ob % metar(n) % v, & ! O, O-B, O-A
967 ob % metar(n) % t, & ! O, O-B, O-A
968 ob % metar(n) % p, & ! O, O-B, O-A
969 ob % metar(n) % q ! O, O-B, O-A
970 end do
971 ob % metar(n) % info % time = metart
972 end do
973 num_obs_metar = num_obs_metar + num_obs
974 elseif ( index( ob_name,'ships') > 0 ) then
975 shipst = shipst + 1
976 do n = num_obs_ships + 1, num_obs_ships + num_obs
977 read(y_unit,'(i8)') num_levs
978 do k = 1, num_levs
979 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
980 ndum, kdum, ob % ships(n) % info % id, & ! Station
981 ob % ships(n) % info % lat, & ! Latitude
982 ob % ships(n) % info % lon, & ! Longitude
983 ob % ships(n) % pressure, & ! Obs height
984 ob % ships(n) % u, & ! O, O-B, O-A
985 ob % ships(n) % v, & ! O, O-B, O-A
986 ob % ships(n) % t, & ! O, O-B, O-A
987 ob % ships(n) % p, & ! O, O-B, O-A
988 ob % ships(n) % q ! O, O-B, O-A
989 end do
990 ob % ships(n) % info % time = shipst
991 end do
992 num_obs_ships = num_obs_ships + num_obs
993
994 elseif ( index( ob_name,'polaramv') > 0 ) then
995 polaramvt = polaramvt + 1
996 do n = num_obs_polaramv + 1, num_obs_polaramv + num_obs
997 read(y_unit,'(i8)') num_levs
998 do k = 1, num_levs
999 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1000 ndum, kdum, ob % polaramv(n) % info % id, & ! Station
1001 ob % polaramv(n) % info % lat, & ! Latitude
1002 ob % polaramv(n) % info % lon, & ! Longitude
1003 ob % polaramv(n) % pressure, & ! Obs pressure
1004 ob % polaramv(n) % u, & ! O, O-B, O-A
1005 ob % polaramv(n) % v
1006 ob % polaramv(n) % info % time = polaramvt
1007 ob % polaramv(n) % channel = ob % polaramv(n) % info % id(3:4)
1008 ob % polaramv(n) % landmask = ob % polaramv(n) % info % id(5:5)
1009 end do
1010 end do
1011 num_obs_polaramv = num_obs_polaramv + num_obs
1012
1013 elseif ( index( ob_name,'geoamv') > 0 ) then
1014 geoamvt = geoamvt + 1
1015 do n = num_obs_geoamv + 1, num_obs_geoamv + num_obs
1016 read(y_unit,'(i8)') num_levs
1017 do k = 1, num_levs
1018 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1019 ndum, kdum, ob % geoamv(n) % info % id, & ! Station
1020 ob % geoamv(n) % info % lat, & ! Latitude
1021 ob % geoamv(n) % info % lon, & ! Longitude
1022 ob % geoamv(n) % pressure, & ! Obs pressure
1023 ob % geoamv(n) % u, & ! O, O-B, O-A
1024 ob % geoamv(n) % v
1025 end do
1026 ob % geoamv(n) % info % time = geoamvt
1027 end do
1028 num_obs_geoamv = num_obs_geoamv + num_obs
1029 elseif ( index( ob_name,'qscat') > 0 ) then
1030 qscatt = qscatt + 1
1031 do n = num_obs_qscat + 1, num_obs_qscat + num_obs
1032 read(y_unit,'(i8)') num_levs
1033 do k = 1, num_levs
1034 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1035 ndum, kdum, ob % qscat(n) % info % id, & ! Station
1036 ob % qscat(n) % info % lat, & ! Latitude
1037 ob % qscat(n) % info % lon, & ! Longitude
1038 ob % qscat(n) % height, & ! Obs height/pressure
1039 ob % qscat(n) % u, & ! O, O-B, O-A
1040 ob % qscat(n) % v
1041 ob % qscat(n) % info % time = qscatt
1042 end do
1043 end do
1044 num_obs_qscat = num_obs_qscat + num_obs
1045
1046 elseif ( index( ob_name,'gpspw') > 0 ) then
1047 gpspwt = gpspwt + 1
1048 do n = num_obs_gpspw + 1, num_obs_gpspw + num_obs
1049 read(y_unit,'(i8)') num_levs
1050 do k = 1, num_levs
1051 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1052 ndum, kdum, ob % gpspw(n) % info % id, & ! Station
1053 ob % gpspw(n) % info % lat, & ! Latitude
1054 ob % gpspw(n) % info % lon, & ! Longitude
1055 rdum, & ! Obs height
1056 ob % gpspw(n) % tpw
1057 ob % gpspw(n) % info % time = gpspwt
1058 end do
1059 end do
1060 num_obs_gpspw = num_obs_gpspw + num_obs
1061
1062 elseif ( index( ob_name,'sound') > 0 ) then
1063 soundt = soundt + 1
1064 do n = num_obs_sound + 1, num_obs_sound + num_obs
1065 read(y_unit,'(i8)')num_levs
1066 ob % sound(n) % numlevs = num_levs
1067 allocate( ob % sound(n) % pressure(1:num_levs) )
1068 allocate( ob % sound(n) % u(1:num_levs) )
1069 allocate( ob % sound(n) % v(1:num_levs) )
1070 allocate( ob % sound(n) % t(1:num_levs) )
1071 allocate( ob % sound(n) % q(1:num_levs) )
1072
1073 do k = 1, num_levs
1074 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1075 ndum, kdum, ob % sound(n) % info % id, & ! Station
1076 ob % sound(n) % info % lat, & ! Latitude
1077 ob % sound(n) % info % lon, & ! Longitude
1078 ob % sound(n) % pressure(k), & ! Obs height
1079 ob % sound(n) % u(k), & ! O, O-B, O-A
1080 ob % sound(n) % v(k), & ! O, O-B, O-A
1081 ob % sound(n) % t(k), & ! O, O-B, O-A
1082 ob % sound(n) % q(k) ! O, O-B, O-A
1083 end do
1084 ob % sound(n) % info % time = soundt
1085 end do
1086 num_obs_sound = num_obs_sound + num_obs
1087
1088 elseif ( index( ob_name,'airsr') > 0 ) then
1089 airsrt = airsrt + 1
1090 do n = num_obs_airsr + 1, num_obs_airsr + num_obs
1091 read(y_unit,'(i8)')num_levs
1092 ob % airsr(n) % numlevs = num_levs
1093 allocate( ob % airsr(n) % pressure(1:num_levs) )
1094 allocate( ob % airsr(n) % t(1:num_levs) )
1095 allocate( ob % airsr(n) % q(1:num_levs) )
1096
1097 do k = 1, num_levs
1098 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1099 ndum, kdum, ob % airsr(n) % info % id, & ! Station
1100 ob % airsr(n) % info % lat, & ! Latitude
1101 ob % airsr(n) % info % lon, & ! Longitude
1102 ob % airsr(n) % pressure(k), & ! Obs height
1103 ob % airsr(n) % t(k), & ! O, O-B, O-A
1104 ob % airsr(n) % q(k) ! O, O-B, O-A
1105 end do
1106 ob % airsr(n) % info % time = airsrt
1107 end do
1108 num_obs_airsr = num_obs_airsr + num_obs
1109
1110 elseif ( index( ob_name,'airep') > 0 ) then
1111 airept = airept + 1
1112 do n = num_obs_airep + 1, num_obs_airep + num_obs
1113 read(y_unit,'(i8)')num_levs
1114 ob % airep(n) % numlevs = num_levs
1115 allocate( ob % airep(n) % pressure(1:num_levs) )
1116 allocate( ob % airep(n) % u(1:num_levs) )
1117 allocate( ob % airep(n) % v(1:num_levs) )
1118 allocate( ob % airep(n) % t(1:num_levs) )
1119
1120 do k = 1, num_levs
1121 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1122 ndum, kdum, ob % airep(n) % info % id, & ! Station
1123 ob % airep(n) % info % lat, & ! Latitude
1124 ob % airep(n) % info % lon, & ! Longitude
1125 ob % airep(n) % pressure(k), & ! Obs height
1126 ob % airep(n) % u(k), & ! O, O-B, O-A
1127 ob % airep(n) % v(k), & ! O, O-B, O-A
1128 ob % airep(n) % t(k) ! O, O-B, O-A
1129 end do
1130 ob % airep(n) % info % time = airept
1131 end do
1132 num_obs_airep = num_obs_airep + num_obs
1133
1134 elseif ( index( ob_name,'pilot') > 0 ) then
1135 pilott = pilott + 1
1136 do n = num_obs_pilot + 1, num_obs_pilot + num_obs
1137 read(y_unit,'(i8)')num_levs
1138 ob % pilot(n) % numlevs = num_levs
1139 allocate( ob % pilot(n) % pressure(1:num_levs) )
1140 allocate( ob % pilot(n) % u(1:num_levs) )
1141 allocate( ob % pilot(n) % v(1:num_levs) )
1142
1143 do k = 1, num_levs
1144 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1145 ndum, kdum, ob % pilot(n) % info % id, & ! Station
1146 ob % pilot(n) % info % lat, & ! Latitude
1147 ob % pilot(n) % info % lon, & ! Longitude
1148 ob % pilot(n) % pressure(k), & ! Obs height
1149 ob % pilot(n) % u(k), & ! O, O-B, O-A
1150 ob % pilot(n) % v(k)
1151 end do
1152 ob % pilot(n) % info % time = pilott
1153 end do
1154 num_obs_pilot = num_obs_pilot + num_obs
1155
1156 elseif ( index( ob_name,'ssmir') > 0 ) then
1157 ssmirt = ssmirt + 1
1158 do n = num_obs_ssmir + 1, num_obs_ssmir + num_obs
1159 read(y_unit,'(i8)')num_levs
1160 do k = 1, num_levs
1161 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1162 ndum, kdum, ob % ssmir(n) % info % id, & ! Station
1163 ob % ssmir(n) % info % lat, & ! Latitude
1164 ob % ssmir(n) % info % lon, & ! Longitude
1165 rdum, & ! Obs height
1166 ob % ssmir(n) % speed, & ! O, O-B, O-A
1167 ob % ssmir(n) % tpw
1168 end do
1169 ob % ssmir(n) % info % time = ssmirt
1170 end do
1171 num_obs_ssmir = num_obs_ssmir + num_obs
1172
1173 elseif ( index( ob_name,'satem') > 0 ) then
1174 satemt = satemt + 1
1175 do n = num_obs_satem + 1, num_obs_satem + num_obs
1176 read(y_unit,'(i8)')num_levs
1177 ob % satem(n) % numlevs = num_levs
1178 allocate( ob % satem(n) % pressure(1:num_levs) )
1179 allocate( ob % satem(n) % thickness(1:num_levs) )
1180
1181 do k = 1, num_levs
1182 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1183 ndum, kdum, ob % satem(n) % info % id, & ! Station
1184 ob % satem(n) % info % lat, & ! Latitude
1185 ob % satem(n) % info % lon, & ! Longitude
1186 ob % satem(n) % pressure(k), & ! Obs height
1187 ob % satem(n) % thickness(k) ! O, O-B, O-A
1188 end do
1189 ob % satem(n) % info % time = satemt
1190 end do
1191 num_obs_satem = num_obs_satem + num_obs
1192
1193 elseif ( index( ob_name,'ssmt1') > 0 ) then
1194 ssmt1t = ssmt1t + 1
1195 do n = num_obs_ssmt1 + 1, num_obs_ssmt1 + num_obs
1196 read(y_unit,'(i8)')num_levs
1197 ob % ssmt1(n) % numlevs = num_levs
1198 allocate( ob % ssmt1(n) % height(1:num_levs) )
1199 allocate( ob % ssmt1(n) % t(1:num_levs) )
1200
1201 do k = 1, num_levs
1202 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1203 ndum, kdum, ob % ssmt1(n) % info % id, & ! Station
1204 ob % ssmt1(n) % info % lat, & ! Latitude
1205 ob % ssmt1(n) % info % lon, & ! Longitude
1206 ob % ssmt1(n) % height(k), & ! Obs height
1207 ob % ssmt1(n) % t(k) ! O, O-B, O-A
1208 end do
1209 ob % ssmt1(n) % info % time = ssmt1t
1210 end do
1211 num_obs_ssmt1 = num_obs_ssmt1 + num_obs
1212
1213 elseif ( index( ob_name,'ssmt2') > 0 ) then
1214 ssmt2t = ssmt2t + 1
1215 do n = num_obs_ssmt2 + 1, num_obs_ssmt2 + num_obs
1216 read(y_unit,'(i8)')num_levs
1217 ob % ssmt2(n) % numlevs = num_levs
1218 allocate( ob % ssmt2(n) % height(1:num_levs) )
1219 allocate( ob % ssmt2(n) % rh(1:num_levs) )
1220
1221 do k = 1, num_levs
1222 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1223 ndum, kdum, ob % ssmt2(n) % info % id, & ! Station
1224 ob % ssmt2(n) % info % lat, & ! Latitude
1225 ob % ssmt2(n) % info % lon, & ! Longitude
1226 ob % ssmt2(n) % height(k), & ! Obs height
1227 ob % ssmt2(n) % rh(k) ! O, O-B, O-A
1228 end do
1229 ob % ssmt2(n) % info % time = ssmt2t
1230 end do
1231 num_obs_ssmt2 = num_obs_ssmt2 + num_obs
1232
1233 elseif ( index( ob_name,'sonde_sfc') > 0 ) then
1234 sonde_sfct = sonde_sfct + 1
1235 do n = num_obs_sonde_sfc + 1, num_obs_sonde_sfc + num_obs
1236 read(y_unit,'(i8)') num_levs
1237 do k = 1, num_levs
1238 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1239 ndum, kdum, ob % sonde_sfc(n) % info % id, & ! Station
1240 ob % sonde_sfc(n) % info % lat, & ! Latitude
1241 ob % sonde_sfc(n) % info % lon, & ! Longitude
1242 ob % sonde_sfc(n) % pressure, & ! Obs height
1243 ob % sonde_sfc(n) % u, & ! O, O-B, O-A
1244 ob % sonde_sfc(n) % v, & ! O, O-B, O-A
1245 ob % sonde_sfc(n) % t, & ! O, O-B, O-A
1246 ob % sonde_sfc(n) % p, & ! O, O-B, O-A
1247 ob % sonde_sfc(n) % q ! O, O-B, O-A
1248 end do
1249 ob % sonde_sfc(n) % info % time = sonde_sfct
1250 end do
1251 num_obs_sonde_sfc = num_obs_sonde_sfc + num_obs
1252
1253 elseif ( index( ob_name,'buoy') > 0 ) then
1254 buoyt = buoyt + 1
1255 do n = num_obs_buoy + 1, num_obs_buoy + num_obs
1256
1257 read(y_unit,'(i8)') num_levs
1258 do k = 1, num_levs
1259 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1260 ndum, kdum, ob % buoy(n) % info % id, & ! Station
1261 ob % buoy(n) % info % lat, & ! Latitude
1262 ob % buoy(n) % info % lon, & ! Longitude
1263 ob % buoy(n) % pressure, & ! Obs height
1264 ob % buoy(n) % u, & ! O, O-B, O-A
1265 ob % buoy(n) % v, & ! O, O-B, O-A
1266 ob % buoy(n) % t, & ! O, O-B, O-A
1267 ob % buoy(n) % p, & ! O, O-B, O-A
1268 ob % buoy(n) % q ! O, O-B, O-A
1269 end do
1270 ob % buoy(n) % info % time = buoyt
1271 end do
1272 num_obs_buoy = num_obs_buoy + num_obs
1273 elseif ( index( ob_name,'profiler') > 0 ) then
1274 profilert = profilert + 1
1275 do n = num_obs_profiler + 1, num_obs_profiler + num_obs
1276 read(y_unit,'(i8)')num_levs
1277 ob % profiler(n) % numlevs = num_levs
1278 allocate( ob % profiler(n) % pressure(1:num_levs) )
1279 allocate( ob % profiler(n) % u(1:num_levs) )
1280 allocate( ob % profiler(n) % v(1:num_levs) )
1281
1282 do k = 1, num_levs
1283 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1284 ndum, kdum, ob % profiler(n) % info % id, & ! Station
1285 ob % profiler(n) % info % lat, & ! Latitude
1286 ob % profiler(n) % info % lon, & ! Longitude
1287 ob % profiler(n) % pressure(k), & ! Obs height
1288 ob % profiler(n) % u(k), & ! O, O-B, O-A
1289 ob % profiler(n) % v(k)
1290 end do
1291 ob % profiler(n) % info % time = profilert
1292 end do
1293 num_obs_profiler = num_obs_profiler + num_obs
1294
1295 elseif ( index( ob_name,'bogus') > 0 ) then
1296 bogust = bogust + 1
1297 do n = num_obs_bogus + 1, num_obs_bogus + num_obs
1298 read(y_unit,'(i8)')num_levs
1299 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1300 ndum, kdum, ob % bogus(n) % info % id, & ! Station
1301 ob % bogus(n) % info % lat, & ! Latitude
1302 ob % bogus(n) % info % lon, & ! Longitude
1303 dummy, &
1304 ob % bogus(n) % slp ! O, O-B, O-A
1305
1306
1307 read(y_unit,'(i8)')num_levs
1308 ob % bogus(n) % numlevs = num_levs
1309 allocate( ob % bogus(n) % pressure(1:num_levs) )
1310 allocate( ob % bogus(n) % u(1:num_levs) )
1311 allocate( ob % bogus(n) % v(1:num_levs) )
1312 allocate( ob % bogus(n) % t(1:num_levs) )
1313 allocate( ob % bogus(n) % q(1:num_levs) )
1314
1315 do k = 1, num_levs
1316 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1317 ndum, kdum, ob % bogus(n) % info % id, & ! Station
1318 ob % bogus(n) % info % lat, & ! Latitude
1319 ob % bogus(n) % info % lon, & ! Longitude
1320 ob % bogus(n) % pressure(k), & ! Obs height
1321 ob % bogus(n) % u(k), & ! O, O-B, O-A
1322 ob % bogus(n) % v(k), & ! O, O-B, O-A
1323 ob % bogus(n) % t(k), & ! O, O-B, O-A
1324 ob % bogus(n) % q(k) ! O, O-B, O-A
1325 end do
1326 ob % bogus(n) % info % time = bogust
1327 end do
1328 num_obs_bogus = num_obs_bogus + num_obs
1329
1330 elseif ( index( ob_name,'ssmiT') > 0 ) then
1331 do n = 1, num_obs
1332 read(y_unit,'(i8)')num_levs
1333 read(y_unit,'(a)')dummy
1334 end do
1335 elseif ( index( ob_name,'gpsref') > 0 ) then
1336 gpsreft = gpsreft + 1
1337 do n = num_obs_gpsref + 1, num_obs_gpsref + num_obs
1338 read(y_unit,'(i8)')num_levs
1339 ob % gpsref(n) % numlevs = num_levs
1340 allocate( ob % gpsref(n) % pressure(1:num_levs) )
1341 allocate( ob % gpsref(n) % ref(1:num_levs) )
1342
1343 do k = 1, num_levs
1344 read(y_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))') &
1345 ndum, kdum, ob % gpsref(n) % info % id, & ! Station
1346 ob % gpsref(n) % info % lat, & ! Latitude
1347 ob % gpsref(n) % info % lon, & ! Longitude
1348 ob % gpsref(n) % pressure(k), & ! Obs height
1349 ob % gpsref(n) % ref(k) ! O, O-B, O-A
1350 end do
1351 ob % gpsref(n) % info % time = gpsreft
1352 end do
1353 num_obs_gpsref = num_obs_gpsref + num_obs
1354
1355 else
1356 print*,' Got unknown obs_tye : ',trim(ob_name)
1357 stop
1358 end if
1359 end do ! loop over particular time (ends in *****)
1360 if ( index( ob_name,'*end*') > 0 ) exit
1361 end do ! loop over entire input file: (ends in *end*)
1362
1363 end subroutine da_read_y
1364
1365 subroutine da_calc_stats( ob_name, num_obs, field )
1366
1367
1368 implicit none
1369
1370 character*9, intent(in) :: ob_name ! Ob name
1371 integer, intent(in) :: num_obs ! Number of observations
1372 type (field_type), intent(in) :: field(:)! Obs data.
1373
1374 integer :: n, count
1375 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1376 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1377
1378 if ( num_obs > 0 ) then
1379 count = 0
1380 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1381
1382 do n = 1, num_obs
1383 call da_increment_stats( field(n), count, &
1384 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1385 mean_omb, mean_oma, stdv_omb, stdv_oma )
1386 end do
1387
1388 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')ob_name, count, '/', num_obs, &
1389 '. O-B(m,s), O-A(m,s) =', &
1390 mean_omb, stdv_omb, mean_oma, stdv_oma
1391
1392 end if
1393
1394 end subroutine da_calc_stats
1395
1396 !--------------------------------------------------------------------------
1397
1398 subroutine da_calc_stats_sound( num_obs, sound )
1399
1400 implicit none
1401
1402 integer, intent(in) :: num_obs ! Number of obs.
1403 type (sound_type), intent(in) :: sound(1:num_obs) ! Obs data.
1404
1405 integer :: n, k, count, count1
1406 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1407 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1408
1409 if ( num_obs > 0 ) then
1410
1411 count = 0; count1 = 0
1412 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1413 do n = 1, num_obs
1414 do k = 1, sound(n) % numlevs
1415 count1 = count1 + 1
1416 call da_increment_stats( sound(n) % u(k), count, &
1417 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1418 mean_omb, mean_oma, stdv_omb, stdv_oma )
1419 end do
1420 end do
1421 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'sound u ', count, '/', count1, &
1422 '. O-B(m,s), O-A(m,s) =', &
1423 mean_omb, stdv_omb, mean_oma, stdv_oma
1424
1425 count = 0; count1 = 0
1426 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1427 do n = 1, num_obs
1428 do k = 1, sound(n) % numlevs
1429 count1 = count1 + 1
1430 call da_increment_stats( sound(n) % v(k), count, &
1431 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1432 mean_omb, mean_oma, stdv_omb, stdv_oma )
1433 end do
1434 end do
1435 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'sound v ', count, '/', count1, &
1436 '. O-B(m,s), O-A(m,s) =', &
1437 mean_omb, stdv_omb, mean_oma, stdv_oma
1438
1439 count = 0; count1 = 0
1440 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1441 do n = 1, num_obs
1442 do k = 1, sound(n) % numlevs
1443 count1 = count1 + 1
1444 call da_increment_stats( sound(n) % t(k), count, &
1445 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1446 mean_omb, mean_oma, stdv_omb, stdv_oma )
1447 end do
1448 end do
1449 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'sound t ', count, '/', count1, &
1450 '. O-B(m,s), O-A(m,s) =', &
1451 mean_omb, stdv_omb, mean_oma, stdv_oma
1452
1453 count = 0; count1 = 0
1454 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1455 do n = 1, num_obs
1456 do k = 1, sound(n) % numlevs
1457 count1 = count1 + 1
1458 call da_increment_stats( sound(n) % q(k), count, &
1459 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1460 mean_omb, mean_oma, stdv_omb, stdv_oma )
1461 end do
1462 end do
1463 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'sound q ', count, '/', count1, &
1464 '. O-B(m,s), O-A(m,s) =', &
1465 mean_omb, stdv_omb, mean_oma, stdv_oma
1466 write(6,*)
1467
1468 end if
1469
1470 end subroutine da_calc_stats_sound
1471
1472 !--------------------------------------------------------------------------
1473
1474 subroutine da_calc_stats_airsr( num_obs, airsr )
1475
1476 implicit none
1477
1478 integer, intent(in) :: num_obs ! Number of obs.
1479 type (airsr_type), intent(in) :: airsr(1:num_obs) ! Obs data.
1480
1481 integer :: n, k, count, count1
1482 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1483 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1484 !--------------------------------------------------------------------------
1485
1486 if ( num_obs > 0 ) then
1487
1488
1489 count = 0; count1 = 0
1490 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1491 do n = 1, num_obs
1492 do k = 1, airsr(n) % numlevs
1493 count1 = count1 + 1
1494 call da_increment_stats( airsr(n) % t(k), count, &
1495 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1496 mean_omb, mean_oma, stdv_omb, stdv_oma )
1497 end do
1498 end do
1499 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'airsr t ', count, '/', count1, &
1500 '. O-B(m,s), O-A(m,s) =', &
1501 mean_omb, stdv_omb, mean_oma, stdv_oma
1502 count = 0; count1 = 0
1503 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1504 do n = 1, num_obs
1505 do k = 1, airsr(n) % numlevs
1506 count1 = count1 + 1
1507 call da_increment_stats( airsr(n) % q(k), count, &
1508 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1509 mean_omb, mean_oma, stdv_omb, stdv_oma )
1510 end do
1511 end do
1512 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'airsr q ', count, '/', count1, &
1513 '. O-B(m,s), O-A(m,s) =', &
1514 mean_omb, stdv_omb, mean_oma, stdv_oma
1515 write(6,*)
1516
1517 end if
1518
1519 end subroutine da_calc_stats_airsr
1520
1521 subroutine da_calc_stats_bogus( num_obs, bogus )
1522
1523 implicit none
1524
1525 integer, intent(in) :: num_obs ! Number of obs.
1526 type (bogus_type), intent(in) :: bogus(1:num_obs) ! Obs data.
1527
1528 integer :: n, k, count, count1
1529 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1530 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1531
1532 if ( num_obs > 0 ) then
1533
1534 count = 0; count1 = 0
1535 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1536 do n = 1, num_obs
1537 do k = 1, bogus(n) % numlevs
1538 count1 = count1 + 1
1539 call da_increment_stats( bogus(n) % u(k), count, &
1540 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1541 mean_omb, mean_oma, stdv_omb, stdv_oma )
1542 end do
1543 end do
1544 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'bogus u ', count, '/', count1, &
1545 '. O-B(m,s), O-A(m,s) =', &
1546 mean_omb, stdv_omb, mean_oma, stdv_oma
1547
1548 count = 0; count1 = 0
1549 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1550 do n = 1, num_obs
1551 do k = 1, bogus(n) % numlevs
1552 count1 = count1 + 1
1553 call da_increment_stats( bogus(n) % v(k), count, &
1554 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1555 mean_omb, mean_oma, stdv_omb, stdv_oma )
1556 end do
1557 end do
1558 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'bogus v ', count, '/', count1, &
1559 '. O-B(m,s), O-A(m,s) =', &
1560 mean_omb, stdv_omb, mean_oma, stdv_oma
1561
1562 count = 0; count1 = 0
1563 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1564 do n = 1, num_obs
1565 do k = 1, bogus(n) % numlevs
1566 count1 = count1 + 1
1567 call da_increment_stats( bogus(n) % t(k), count, &
1568 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1569 mean_omb, mean_oma, stdv_omb, stdv_oma )
1570 end do
1571 end do
1572 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'bogus t ', count, '/', count1, &
1573 '. O-B(m,s), O-A(m,s) =', &
1574 mean_omb, stdv_omb, mean_oma, stdv_oma
1575
1576 count = 0; count1 = 0
1577 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1578 do n = 1, num_obs
1579 do k = 1, bogus(n) % numlevs
1580 count1 = count1 + 1
1581 call da_increment_stats( bogus(n) % q(k), count, &
1582 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1583 mean_omb, mean_oma, stdv_omb, stdv_oma )
1584 end do
1585 end do
1586 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'bogus q ', count, '/', count1, &
1587 '. O-B(m,s), O-A(m,s) =', &
1588 mean_omb, stdv_omb, mean_oma, stdv_oma
1589 count = 0; count1 = 0
1590 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1591 do n = 1, num_obs
1592 count1 = count1 + 1
1593 call da_increment_stats( bogus(n) % slp, count, &
1594 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1595 mean_omb, mean_oma, stdv_omb, stdv_oma )
1596 end do
1597 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'bogus slp', count, '/', count1, &
1598 '. O-B(m,s), O-A(m,s) =', &
1599 mean_omb, stdv_omb, mean_oma, stdv_oma
1600 write(6,*)
1601
1602 end if
1603
1604 end subroutine da_calc_stats_bogus
1605
1606 subroutine da_calc_stats_airep( num_obs, airep )
1607
1608 implicit none
1609
1610 integer, intent(in) :: num_obs ! Number of obs.
1611 type (airep_type), intent(in) :: airep(1:num_obs) ! Obs data.
1612
1613 integer :: n, k, count, count1
1614 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1615 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1616
1617 if ( num_obs > 0 ) then
1618
1619 count = 0; count1 = 0
1620 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1621 do n = 1, num_obs
1622 do k = 1, airep(n) % numlevs
1623 count1 = count1 + 1
1624 call da_increment_stats( airep(n) % u(k), count, &
1625 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1626 mean_omb, mean_oma, stdv_omb, stdv_oma )
1627 end do
1628 end do
1629 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'airep u ', count, '/', count1, &
1630 '. O-B(m,s), O-A(m,s) =', &
1631 mean_omb, stdv_omb, mean_oma, stdv_oma
1632
1633 count = 0; count1 = 0
1634 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1635 do n = 1, num_obs
1636 do k = 1, airep(n) % numlevs
1637 count1 = count1 + 1
1638 call da_increment_stats( airep(n) % v(k), count, &
1639 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1640 mean_omb, mean_oma, stdv_omb, stdv_oma )
1641 end do
1642 end do
1643 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'airep v ', count, '/', count1, &
1644 '. O-B(m,s), O-A(m,s) =', &
1645 mean_omb, stdv_omb, mean_oma, stdv_oma
1646 count = 0; count1 = 0
1647 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1648 do n = 1, num_obs
1649 do k = 1, airep(n) % numlevs
1650 count1 = count1 + 1
1651 call da_increment_stats( airep(n) % t(k), count, &
1652 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1653 mean_omb, mean_oma, stdv_omb, stdv_oma )
1654 end do
1655 end do
1656 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'airep t ', count, '/', count1, &
1657 '. O-B(m,s), O-A(m,s) =', &
1658 mean_omb, stdv_omb, mean_oma, stdv_oma
1659 write(6,*)
1660
1661 end if
1662
1663 end subroutine da_calc_stats_airep
1664
1665 subroutine da_calc_stats_pilot( num_obs, pilot )
1666
1667 implicit none
1668
1669 integer, intent(in) :: num_obs ! Number of obs.
1670 type (pilot_type), intent(in) :: pilot(1:num_obs) ! Obs data.
1671
1672 integer :: n, k, count, count1
1673 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1674 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1675
1676 if ( num_obs > 0 ) then
1677
1678 count = 0; count1 = 0
1679 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1680 do n = 1, num_obs
1681 do k = 1, pilot(n) % numlevs
1682 count1 = count1 + 1
1683 call da_increment_stats( pilot(n) % u(k), count, &
1684 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1685 mean_omb, mean_oma, stdv_omb, stdv_oma )
1686 end do
1687 end do
1688 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'pilot u ', count, '/', count1, &
1689 '. O-B(m,s), O-A(m,s) =', &
1690 mean_omb, stdv_omb, mean_oma, stdv_oma
1691
1692 count = 0; count1 = 0
1693 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1694 do n = 1, num_obs
1695 do k = 1, pilot(n) % numlevs
1696 count1 = count1 + 1
1697 call da_increment_stats( pilot(n) % v(k), count, &
1698 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1699 mean_omb, mean_oma, stdv_omb, stdv_oma )
1700 end do
1701 end do
1702 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'pilot v ', count, '/', count1, &
1703 '. O-B(m,s), O-A(m,s) =', &
1704 mean_omb, stdv_omb, mean_oma, stdv_oma
1705 write(6,*)
1706
1707 end if
1708 end subroutine da_calc_stats_pilot
1709
1710 subroutine da_calc_stats_profiler( num_obs, profiler )
1711
1712 implicit none
1713
1714 integer, intent(in) :: num_obs ! Number of obs.
1715 type (pilot_type), intent(in) :: profiler(1:num_obs) ! Obs data.
1716
1717 integer :: n, k, count, count1
1718 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1719 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1720
1721 if ( num_obs > 0 ) then
1722
1723 count = 0; count1 = 0
1724 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1725 do n = 1, num_obs
1726 do k = 1, profiler(n) % numlevs
1727 count1 = count1 + 1
1728 call da_increment_stats( profiler(n) % u(k), count, &
1729 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1730 mean_omb, mean_oma, stdv_omb, stdv_oma )
1731 end do
1732 end do
1733 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'profiler u ', count, '/', count1, &
1734 '. O-B(m,s), O-A(m,s) =', &
1735 mean_omb, stdv_omb, mean_oma, stdv_oma
1736
1737 count = 0; count1 = 0
1738 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1739 do n = 1, num_obs
1740 do k = 1, profiler(n) % numlevs
1741 count1 = count1 + 1
1742 call da_increment_stats( profiler(n) % v(k), count, &
1743 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1744 mean_omb, mean_oma, stdv_omb, stdv_oma )
1745 end do
1746 end do
1747 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'profiler v ', count, '/', count1, &
1748 '. O-B(m,s), O-A(m,s) =', &
1749 mean_omb, stdv_omb, mean_oma, stdv_oma
1750 write(6,*)
1751
1752 end if
1753 end subroutine da_calc_stats_profiler
1754
1755 subroutine da_calc_stats_gpsref( num_obs, gpsref )
1756
1757 implicit none
1758
1759 integer, intent(in) :: num_obs ! Number of obs.
1760 type (gpsref_type), intent(in) :: gpsref(1:num_obs) ! Obs data.
1761
1762 integer :: n, k, count, count1
1763 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1764 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1765
1766 if ( num_obs > 0 ) then
1767
1768 count = 0; count1 = 0
1769 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1770 do n = 1, num_obs
1771 do k = 1, gpsref(n) % numlevs
1772 count1 = count1 + 1
1773 call da_increment_stats( gpsref(n) % ref(k), count, &
1774 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1775 mean_omb, mean_oma, stdv_omb, stdv_oma )
1776 end do
1777 end do
1778 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'gpsref re', count, '/', count1, &
1779 '. O-B(m,s), O-A(m,s) =', &
1780 mean_omb, stdv_omb, mean_oma, stdv_oma
1781 write(6,*)
1782
1783 end if
1784 end subroutine da_calc_stats_gpsref
1785
1786 subroutine da_calc_stats_satem( num_obs, satem )
1787
1788 implicit none
1789
1790 integer, intent(in) :: num_obs ! Number of obs.
1791 type (satem_type), intent(in) :: satem(1:num_obs) ! Obs data.
1792
1793 integer :: n, k, count, count1
1794 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1795 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1796
1797 if ( num_obs > 0 ) then
1798
1799 count = 0; count1 = 0
1800 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1801 do n = 1, num_obs
1802 do k = 1, satem(n) % numlevs
1803 count1 = count1 + 1
1804 call da_increment_stats( satem(n) % thickness(k), count, &
1805 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1806 mean_omb, mean_oma, stdv_omb, stdv_oma )
1807 end do
1808 end do
1809 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'satem thk', count, '/', count1, &
1810 '. O-B(m,s), O-A(m,s) =', &
1811 mean_omb, stdv_omb, mean_oma, stdv_oma
1812 write(6,*)
1813
1814 end if
1815
1816 end subroutine da_calc_stats_satem
1817
1818 subroutine da_calc_stats_ssmt1( num_obs, ssmt1 )
1819
1820 implicit none
1821
1822 integer, intent(in) :: num_obs ! Number of obs.
1823 type (ssmt1_type), intent(in) :: ssmt1(1:num_obs) ! Obs data.
1824
1825 integer :: n, k, count, count1
1826 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1827 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1828
1829 if ( num_obs > 0 ) then
1830
1831 count = 0; count1 = 0
1832 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1833 do n = 1, num_obs
1834 do k = 1, ssmt1(n) % numlevs
1835 count1 = count1 + 1
1836 call da_increment_stats( ssmt1(n) % t(k), count, &
1837 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1838 mean_omb, mean_oma, stdv_omb, stdv_oma )
1839 end do
1840 end do
1841 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'ssmt1 t ', count, '/', count1, &
1842 '. O-B(m,s), O-A(m,s) =', &
1843 mean_omb, stdv_omb, mean_oma, stdv_oma
1844 write(6,*)
1845
1846 end if
1847
1848 end subroutine da_calc_stats_ssmt1
1849
1850 subroutine da_calc_stats_ssmt2( num_obs, ssmt2 )
1851
1852 implicit none
1853
1854 integer, intent(in) :: num_obs ! Number of obs.
1855 type (ssmt2_type), intent(in) :: ssmt2(1:num_obs) ! Obs data.
1856
1857 integer :: n, k, count, count1
1858 real :: sum_omb, sum_oma, sum_omb2, sum_oma2
1859 real :: mean_omb, mean_oma, stdv_omb, stdv_oma
1860
1861 if ( num_obs > 0 ) then
1862
1863 count = 0; count1 = 0
1864 sum_omb = 0.0; sum_oma = 0.0; sum_omb2 = 0.0; sum_oma2 = 0.0
1865 do n = 1, num_obs
1866 do k = 1, ssmt2(n) % numlevs
1867 count1 = count1 + 1
1868 call da_increment_stats( ssmt2(n) % rh(k), count, &
1869 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1870 mean_omb, mean_oma, stdv_omb, stdv_oma )
1871 end do
1872 end do
1873 write(6,'(1x,a9,i10,a,i10,a,4f10.4)')'ssmt2 rh ', count, '/', count1, &
1874 '. O-B(m,s), O-A(m,s) =', &
1875 mean_omb, stdv_omb, mean_oma, stdv_oma
1876 write(6,*)
1877
1878 end if
1879
1880 end subroutine da_calc_stats_ssmt2
1881
1882 subroutine da_increment_stats( field, count, &
1883 sum_omb, sum_oma, sum_omb2, sum_oma2, &
1884 mean_omb, mean_oma, stdv_omb, stdv_oma )
1885
1886 implicit none
1887
1888 type (field_type), intent(in) :: field
1889 integer, intent(inout) :: count
1890 real, intent(inout) :: sum_omb, sum_oma, sum_omb2, sum_oma2
1891 real, intent(out) :: mean_omb, mean_oma, stdv_omb, stdv_oma
1892
1893 if ( field % qc >= obs_qc_pointer ) then
1894 count = count + 1
1895 sum_omb = sum_omb + field % omb
1896 sum_oma = sum_oma + field % oma
1897 sum_omb2 = sum_omb2 + field % omb**2
1898 sum_oma2 = sum_oma2 + field % oma**2
1899 end if
1900
1901 if ( count > 0 ) then
1902 mean_omb = sum_omb / real(count)
1903 mean_oma = sum_oma / real(count)
1904 end if
1905
1906 ! Cannot calculate standard deviations for first observation
1907 if ( count > 1 ) then
1908 stdv_omb = sqrt( sum_omb2/ real(count) - mean_omb**2 )
1909 stdv_oma = sqrt( sum_oma2/ real(count) - mean_oma**2 )
1910 end if
1911
1912 end subroutine da_increment_stats
1913
1914 end program da_tune_obs_hollingsworth1