da_verif_obs.f90
References to this file elsewhere.
1 program da_verif_obs
2 !---------------------------------------------------------------------------
3 ! History:
4 !
5 ! Abstract:
6 ! Program to read diagnostics written in fort.50 by WRFVAR
7 ! and write in proper format to get ploted using PC-XL utility
8 !
9 ! Author: Syed RH Rizvi NCAR/MMM 05/25/2006
10 ! Updates:
11 ! Hui Shao NCAR/RAL/DATC 05/02/2007
12 ! Diagnositics for GPSREF
13 ! Syed RH Rizvi NCAR/MMM 05/08/2007
14 ! Significance test & error bars are added
15 !---------------------------------------------------------------------------
16
17 use da_verif_obs_control, only : surface_type, upr_type, gpspw_type, &
18 gpsref_type, record1, record2, record3, &
19 record4, record5, stats_value, exp_dirs, out_dirs, nstd,nstdh, &
20 rmiss, diag_unit_out, nml_unit, alpha, &
21 diag_unit_in, info_unit, exp_num, end_date, file_path_string, &
22 if_plot_bias, if_plot_airsret, if_plot_airep,if_plot_abias, &
23 if_plot_buoy, if_plot_gpspw, if_plot_gpsref, if_plot_pilot, &
24 if_plot_profiler, if_plot_polaramv, if_plot_qscat, if_plot_rmse, &
25 if_plot_sound, if_plot_sonde_sfc, if_plot_synop, if_plot_surface, &
26 if_plot_upr, if_plot_ships, if_plot_metar, interval, stdp, start_date, &
27 if_plot_geoamv
28 use da_verif_obs_init, only : initialize_surface_type, initialize_upr_type, &
29 initialize_gpspw_type, initialize_gpsref_type, da_advance_cymdh , &
30 initialize_t_tab
31
32
33 implicit none
34 ! Typically 12 hours
35 integer :: num_obs
36 character*20 :: obs_type, dummy_c
37
38 character*5 :: stn_id
39 integer :: n, k, kk, l, levels, dummy_i
40 real :: lat, lon, press, dummy
41 real :: u_obs, u_inv, u_error, u_inc, &
42 v_obs, v_inv, v_error, v_inc, &
43 t_obs, t_inv, t_error, t_inc, &
44 p_obs, p_inv, p_error, p_inc, &
45 q_obs, q_inv, q_error, q_inc, &
46 spd_obs, spd_inv, spd_err, spd_inc
47 real :: tpw_obs, tpw_inv, tpw_err, tpw_inc
48 real :: ref_obs, ref_inv, ref_err, ref_inc
49 integer :: u_qc, v_qc, t_qc, p_qc, q_qc, tpw_qc, spd_qc, ref_qc
50 integer :: npr, ier, iexp
51 character*10 :: date, new_date ! Current date (ccyymmddhh).
52 integer :: sdate, cdate, edate ! Starting, current ending dates.
53 logical :: if_write, is_file
54
55 character(len=512) :: out_dir,filename
56 type (surface_type) :: surface
57 type (upr_type) :: upr, gupr
58 type (gpspw_type) :: gpspw
59 type (gpsref_type) :: gpsref, ggpsref
60
61
62 nml_unit = 10
63 diag_unit_in = 50
64 diag_unit_out = 20
65 info_unit = 30
66
67 exp_num = 0
68 exp_dirs = ''
69 out_dirs = ''
70
71 if_plot_rmse = .false.
72 if_plot_bias = .false.
73 if_plot_abias = .false.
74
75 if_plot_synop = .false.
76 if_plot_sonde_sfc = .false.
77 if_plot_metar = .false.
78 if_plot_ships = .false.
79 if_plot_qscat = .false.
80 if_plot_buoy = .false.
81
82 if_plot_sound = .false.
83 if_plot_geoamv = .false.
84 if_plot_polaramv = .false.
85 if_plot_profiler = .false.
86 if_plot_airep = .false.
87 if_plot_pilot = .false.
88
89 if_plot_gpspw = .false.
90 if_plot_gpsref = .false.
91 if_plot_airsret = .false.
92
93 file_path_string = 'wrfvar/working/gts_omb_oma'
94
95 ! Read in namelist information defined in module define_cons_types
96
97 open ( unit=nml_unit, file='namelist.plot_diag', STATUS='OLD', &
98 form='formatted' )
99 read ( unit=nml_unit, nml=record1, IOSTAT=ier )
100 write ( unit=*, nml = record1 )
101 if ( ier /= 0 ) then
102 write (*,*) 'error in reading namelist record1'
103 stop
104 end if
105
106 read ( unit=nml_unit, nml=record2, iostat=ier )
107 write ( unit=*, nml = record2 )
108 if ( ier /= 0 ) then
109 write (*,*) 'error in reading namelist record2'
110 stop
111 end if
112 read ( unit=nml_unit, nml=record3, iostat=ier )
113 write ( unit=*, nml = record3 )
114 if ( ier /= 0 ) then
115 write (*,*) 'error in reading namelist record3'
116 stop
117 end if
118 read ( unit=nml_unit, nml=record4, iostat=ier )
119 write ( unit=*, nml = record4 )
120 if ( ier /= 0 ) then
121 write (*,*) 'error in reading namelist record4'
122 stop
123 end if
124 read ( unit=nml_unit, nml=record5, iostat=ier )
125 write ( unit=*, nml = record5 )
126 if ( ier /= 0 ) then
127 write (*,*) 'error in reading namelist record5'
128 stop
129 end if
130 close(nml_unit)
131 call initialize_t_tab
132
133 !---------------------------------------------------------------------------
134 ! Loop over experiments
135
136 do iexp =1,exp_num
137
138 if_plot_surface = .false.
139 if (if_plot_synop .or. if_plot_metar .or. if_plot_ships .or. if_plot_buoy .or. &
140 if_plot_sonde_sfc .or. if_plot_qscat ) if_plot_surface = .true.
141
142 if_plot_upr = .false.
143 if (if_plot_sound .or. if_plot_pilot .or. if_plot_profiler .or. &
144 if_plot_geoamv .or. if_plot_polaramv .or. if_plot_airep .or. &
145 if_plot_airsret &
146 ) if_plot_upr= .true.
147
148 ! ! Typically 12 hours
149 read(start_date(1:10), fmt='(i10)')sdate
150 read(end_date(1:10), fmt='(i10)')edate
151 write(6,'(4a)')' Diag Starting date ', start_date, ' Ending date ', end_date
152 write(6,'(a,i8,a)')' Interval between dates = ', interval, ' hours.'
153
154 date = start_date
155 cdate = sdate
156 cdate = sdate
157 call initialize_upr_type(gupr)
158 call initialize_gpsref_type(ggpsref)
159
160
161 do while ( cdate <= edate )
162 ! Initialize various types
163 call initialize_surface_type(surface)
164 call initialize_upr_type(upr)
165 call initialize_gpspw_type(gpspw)
166 call initialize_gpsref_type(gpsref)
167
168 ! construct file name
169
170 filename = TRIM(exp_dirs(iexp))//'/'//date//'/'//trim(file_path_string)
171
172 ! check if the file exists, then open the file
173
174 inquire ( file=trim(filename), exist=is_file)
175 if ( .not. is_file ) then
176 write(0,*) 'can not find the file: ', trim(filename)
177 stop
178 end if
179 open (unit=diag_unit_in, file=trim(filename), form='formatted', &
180 status='old', iostat=ier)
181 1 continue
182
183 if_write = .false.
184 read(diag_unit_in,'(a20,i8)', end=2000, err = 1000)obs_type,num_obs
185 if (index( obs_type,'synop') > 0 ) then
186 if (if_plot_synop ) if_write = .true.
187 goto 10
188
189 elseif (index( obs_type,'metar') > 0 ) then
190 if (if_plot_metar ) if_write = .true.
191 goto 10
192
193 elseif (index( obs_type,'ships') > 0 ) then
194 if (if_plot_ships ) if_write = .true.
195 goto 10
196
197 elseif (index( obs_type,'buoy' ) > 0 ) then
198 if (if_plot_buoy ) if_write = .true.
199 goto 10
200
201 elseif (index( obs_type,'sonde_sfc') > 0 ) then
202 if (if_plot_sonde_sfc ) if_write = .true.
203 goto 10
204
205 elseif (index( obs_type,'polaramv') > 0) then
206 if (if_plot_polaramv ) if_write = .true.
207 goto 20
208
209 elseif (index( obs_type,'geoamv' ) > 0) then
210 if (if_plot_geoamv ) if_write = .true.
211 goto 20
212
213 elseif (index( obs_type,'gpspw') > 0) then
214 if ( if_plot_gpspw ) if_write = .true.
215 goto 30
216
217 elseif (index( obs_type,'sound') > 0) then
218 if (if_plot_sound ) if_write = .true.
219 goto 40
220
221 elseif (index( obs_type,'airep') > 0) then
222 if (if_plot_airep ) if_write = .true.
223 goto 50
224
225 elseif (index( obs_type,'pilot') > 0) then
226 if (if_plot_pilot ) if_write = .true.
227 goto 60
228
229 elseif (index( obs_type,'profiler') > 0) then
230 if (if_plot_profiler ) if_write = .true.
231 goto 60
232
233 elseif (index( obs_type,'ssmir') > 0) then
234 goto 70
235
236 elseif (index( obs_type,'ssmiT') > 0) then
237 goto 80
238
239 elseif (index( obs_type,'satem') > 0) then
240 goto 90
241
242 elseif (index( obs_type,'ssmt1') > 0) then
243 goto 100
244
245 elseif (index( obs_type,'ssmt2') > 0) then
246 goto 100
247
248 elseif (index( obs_type,'qscat') > 0) then
249 if (if_plot_qscat ) if_write = .true.
250 goto 110
251 elseif (index( obs_type,'gpsref' ) > 0) then
252 if (if_plot_gpsref ) if_write = .true.
253 goto 120
254 elseif (index( obs_type,'airsr') > 0) then
255 if (if_plot_airsret ) if_write = .true.
256 goto 130
257 else
258 print*,' Got unknown OBS_TYPE ',obs_type(1:20),' on unit ',diag_unit_in
259 stop
260 end if
261
262 10 continue ! Synop, Metar, Ships, Buoy , Sonde_sfc
263
264 if ( num_obs > 0 ) then
265 do n = 1, num_obs
266 read(diag_unit_in,'(i8)')levels
267 do k = 1, levels
268 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
269 kk, l, stn_id, & ! Station
270 lat, lon, press, & ! Lat/lon, pressure
271 u_obs, u_inv, u_qc, u_error, u_inc, &
272 v_obs, v_inv, v_qc, v_error, v_inc, &
273 t_obs, t_inv, t_qc, t_error, t_inc, &
274 p_obs, p_inv, p_qc, p_error, p_inc, &
275 q_obs, q_inv, q_qc, q_error, q_inc
276 if (if_write) then
277 if (u_qc >= 0) call update_stats(surface%uomb, surface%uoma, u_inv, u_inc)
278 if (v_qc >= 0) call update_stats(surface%vomb, surface%voma, v_inv, v_inc)
279 if (t_qc >= 0) call update_stats(surface%tomb, surface%toma, t_inv, t_inc)
280 if (p_qc >= 0) call update_stats(surface%pomb, surface%poma, p_inv, p_inc)
281 if (q_qc >= 0) call update_stats(surface%qomb, surface%qoma, q_inv, q_inc)
282 end if
283 end do ! loop over levels
284 end do ! loop over Obs
285 end if
286 goto 1
287
288 20 continue ! Polar or Geo AMV's
289
290 if ( num_obs > 0 ) then
291 do n = 1, num_obs
292 read(diag_unit_in,'(i8)')levels
293 do k = 1, levels
294 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
295 kk, l, stn_id, & ! Station
296 lat, lon, press, & ! Lat/lon, pressure
297 u_obs, u_inv, u_qc, u_error, u_inc, &
298 v_obs, v_inv, v_qc, v_error, v_inc
299
300 if (if_write .and. press > 0 ) then
301 call get_std_pr_level(press, npr, stdp, nstd)
302 if( u_qc >= 0) then
303 call update_stats(upr%uomb(npr),upr%uoma(npr),u_inv,u_inc)
304 call update_stats(gupr%uomb(npr),gupr%uoma(npr),u_inv,u_inc)
305 endif
306 if( v_qc >= 0) then
307 call update_stats(upr%vomb(npr),upr%voma(npr),v_inv,v_inc)
308 call update_stats(gupr%vomb(npr),gupr%voma(npr),v_inv,v_inc)
309 endif
310
311 end if
312 end do ! loop over levels
313 end do ! loop over Obs
314 end if
315
316 goto 1
317
318 30 continue ! Gpspw
319
320 if ( num_obs > 0 ) then
321 do n = 1, num_obs
322 read(diag_unit_in,'(i8)')levels
323 do k = 1, levels
324 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
325 kk, l, stn_id, & ! Station
326 lat, lon, dummy, & ! Lat/lon, dummy
327 tpw_obs, tpw_inv, tpw_qc, tpw_err, tpw_inc
328 if (if_write) then
329 if (tpw_qc >= 0) call update_stats(gpspw%tpwomb,gpspw%tpwoma,tpw_inv,tpw_inc)
330 end if
331 end do ! loop over levels
332 end do ! loop over Obs
333 end if
334
335 goto 1
336
337 40 continue ! Sound
338
339 ! [6] Transfer sound obs:
340
341 if ( num_obs > 0 ) then
342 do n = 1, num_obs
343 read(diag_unit_in,'(i8)')levels
344 do k = 1, levels
345 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
346 kk,l, stn_id, & ! Station
347 lat, lon, press, & ! Lat/lon, dummy
348 u_obs, u_inv, u_qc, u_error, u_inc, &
349 v_obs, v_inv, v_qc, v_error, v_inc, &
350 t_obs, t_inv, t_qc, t_error, t_inc, &
351 q_obs, q_inv, q_qc, q_error, q_inc
352 if (if_write .and. press > 0 ) then
353 call get_std_pr_level(press, npr, stdp, nstd)
354
355 if( u_qc >= 0) then
356 call update_stats(upr%uomb(npr),upr%uoma(npr),u_inv,u_inc)
357 call update_stats(gupr%uomb(npr),gupr%uoma(npr),u_inv,u_inc)
358 endif
359 if( v_qc >= 0) then
360 call update_stats(upr%vomb(npr),upr%voma(npr),v_inv,v_inc)
361 call update_stats(gupr%vomb(npr),gupr%voma(npr),v_inv,v_inc)
362 endif
363 if( t_qc >= 0) then
364 call update_stats(upr%tomb(npr),upr%toma(npr),t_inv,t_inc)
365 call update_stats(gupr%tomb(npr),gupr%toma(npr),t_inv,t_inc)
366 endif
367 if( q_qc >= 0) then
368 call update_stats(upr%qomb(npr),upr%qoma(npr),q_inv,q_inc)
369 call update_stats(gupr%qomb(npr),gupr%qoma(npr),q_inv,q_inc)
370 endif
371 !
372 end if
373 end do ! loop over levels
374 end do ! loop over Obs
375 end if
376 goto 1
377
378 50 continue ! Airep
379
380 if ( num_obs > 0 ) then
381 do n = 1, num_obs
382 read(diag_unit_in,'(i8)') levels
383 do k = 1, levels
384 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
385 kk,l, stn_id, & ! Station
386 lat, lon, press, & ! Lat/lon, dummy
387 u_obs, u_inv, u_qc, u_error, u_inc, &
388 v_obs, v_inv, v_qc, v_error, v_inc, &
389 t_obs, t_inv, t_qc, t_error, t_inc
390 if (if_write .and. press > 0 ) then
391 call get_std_pr_level(press, npr, stdp, nstd)
392 if( u_qc >= 0) then
393 call update_stats(upr%uomb(npr),upr%uoma(npr),u_inv,u_inc)
394 call update_stats(gupr%uomb(npr),gupr%uoma(npr),u_inv,u_inc)
395 endif
396 if( v_qc >= 0) then
397 call update_stats(upr%vomb(npr),upr%voma(npr),v_inv,v_inc)
398 call update_stats(gupr%vomb(npr),gupr%voma(npr),v_inv,v_inc)
399 endif
400 if( t_qc >= 0) then
401 call update_stats(upr%tomb(npr),upr%toma(npr),t_inv,t_inc)
402 call update_stats(gupr%tomb(npr),gupr%toma(npr),t_inv,t_inc)
403 endif
404
405 end if
406 end do ! loop over levels
407 end do ! loop over Obs
408 end if
409
410 goto 1
411
412 60 continue ! Pilot & Profiler
413
414 ! [8] Transfer pilot obs:
415 if ( num_obs > 0 ) then
416 do n = 1, num_obs
417 read(diag_unit_in,'(i8)')levels
418 do k = 1, levels
419 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
420 kk,l, stn_id, & ! Station
421 lat, lon, press, & ! Lat/lon, dummy
422 u_obs, u_inv, u_qc, u_error, u_inc, &
423 v_obs, v_inv, v_qc, v_error, v_inc
424 if (if_write .and. press > 0 ) then
425 call get_std_pr_level(press, npr, stdp, nstd)
426 if( u_qc >= 0) then
427 call update_stats(upr%uomb(npr),upr%uoma(npr),u_inv,u_inc)
428 call update_stats(gupr%uomb(npr),gupr%uoma(npr),u_inv,u_inc)
429 endif
430 if( v_qc >= 0) then
431 call update_stats(upr%vomb(npr),upr%voma(npr),v_inv,v_inc)
432 call update_stats(gupr%vomb(npr),gupr%voma(npr),v_inv,v_inc)
433 endif
434
435 end if
436 end do ! loop over levels
437 end do ! loop over Obs
438 end if
439 goto 1
440
441 70 continue ! SSMI retrievals
442
443 if ( num_obs > 0 ) then
444 do n = 1, num_obs
445 read(diag_unit_in,'(i8)')levels
446 do k = 1, levels
447 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
448 kk,l, stn_id, & ! Station
449 lat, lon, dummy, & ! Lat/lon, dummy
450 spd_obs, spd_inv, spd_qc, spd_err, spd_inc
451 end do
452 end do
453 end if
454
455 goto 1
456
457 80 continue ! SSMI radiance
458
459 if ( num_obs > 0 ) then
460 do n = 1, num_obs
461 read(diag_unit_in,*)dummy_c
462 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,7(2f17.7,i8,2f17.7))', err= 1000)&
463 k, l, stn_id, & ! Station
464 lat, lon, dummy, & ! Lat/lon, dummy
465 dummy, dummy, dummy_i, dummy, dummy, &
466 dummy, dummy, dummy_i, dummy, dummy, &
467 dummy, dummy, dummy_i, dummy, dummy, &
468 dummy, dummy, dummy_i, dummy, dummy, &
469 dummy, dummy, dummy_i, dummy, dummy, &
470 dummy, dummy, dummy_i, dummy, dummy, &
471 dummy, dummy, dummy_i, dummy, dummy
472 end do
473 end if
474 goto 1
475
476 90 continue ! SATEM
477
478 if ( num_obs > 0 ) then
479 do n = 1, num_obs
480 read(diag_unit_in,'(i8)') levels
481 do k = 1, levels
482 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
483 kk,l, stn_id, & ! Station
484 lat, lon, dummy, & ! Lat/lon, dummy
485 dummy,dummy, dummy_i, dummy, dummy
486 end do ! loop over levels
487 end do ! loop over Obs
488 end if
489
490 goto 1
491
492 100 continue ! SSMT1 & 2
493
494 if ( num_obs > 0 ) then
495 do n = 1, num_obs
496 read(diag_unit_in,'(i8)') levels
497 do k = 1, levels
498 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
499 kk,l, stn_id, & ! Station
500 lat, lon, dummy, & ! Lat/lon, dummy
501 dummy,dummy, dummy_i, dummy, dummy
502 end do ! loop over levels
503 end do ! loop over obs
504 end if
505
506 goto 1
507
508 110 continue ! Scatrometer winds
509
510 if ( num_obs > 0 ) then
511 do n = 1, num_obs
512 read(diag_unit_in,'(i8)') levels
513 do k = 1, levels
514 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
515 kk, l, stn_id, & ! Station
516 lat, lon, press, & ! Lat/lon, dummy
517 u_obs, u_inv, u_qc, u_error, u_inc, &
518 v_obs, v_inv, v_qc, v_error, v_inc
519 if (if_write) then
520 if (u_qc >= 0) call update_stats(surface%uomb,surface%uoma,u_inv,u_inc)
521 if (v_qc >= 0) call update_stats(surface%vomb,surface%voma,v_inv,v_inc)
522 end if
523 end do ! loop over levels
524 end do ! loop over obs
525 end if
526 goto 1
527
528 120 continue ! Gpsref
529
530 IF ( num_obs > 0 ) THEN
531 DO n = 1, num_obs
532 read(diag_unit_in,'(i8)') levels
533 DO k = 1, levels
534 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
535 kk, l, stn_id, & ! Station
536 lat, lon, press, & ! Lat/lon, dummy
537 ref_obs, ref_inv, ref_qc, ref_err, ref_inc
538 if (if_write) then
539 if( ref_qc >= 0) then
540 call update_stats(gpsref%refomb(k),gpsref%refoma(k),ref_inv,ref_inc)
541 call update_stats(ggpsref%refomb(k),ggpsref%refoma(k),ref_inv,ref_inc)
542 end if
543 endif
544 END DO ! loop over levels
545 END DO ! loop over Obs
546 ENDIF
547 go to 1
548 !---------------------------------------------------------------------
549
550 130 continue ! AIRSRET
551
552 if ( num_obs > 0 ) then
553 do n = 1, num_obs
554 read(diag_unit_in,'(i8)')levels
555 do k = 1, levels
556 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
557 kk,l, stn_id, & ! Station
558 lat, lon, press, & ! Lat/lon, dummy
559 t_obs, t_inv, t_qc, t_error, t_inc, &
560 q_obs, q_inv, q_qc, q_error, q_inc
561 if (if_write .and. press > 0 ) then
562 call get_std_pr_level(press, npr, stdp, nstd)
563 if( t_qc >= 0) then
564 call update_stats(upr%tomb(npr),upr%toma(npr),t_inv,t_inc)
565 call update_stats(gupr%tomb(npr),gupr%toma(npr),t_inv,t_inc)
566 endif
567 if( q_qc >= 0) then
568 call update_stats(upr%qomb(npr),upr%qoma(npr),q_inv,q_inc)
569 call update_stats(gupr%qomb(npr),gupr%qoma(npr),q_inv,q_inc)
570 endif
571 end if
572 end do ! loop over levels
573 end do ! loop over obs
574 end if
575 goto 1
576
577 ! Now process the diagnostics
578 2000 continue
579
580 close (diag_unit_in)
581 ! Write output on outunit
582 out_dir=trim(out_dirs(iexp))
583 if (if_plot_surface ) then
584 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_u',surface%uomb,surface%uoma)
585 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_v',surface%vomb,surface%voma)
586 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_t',surface%tomb,surface%toma)
587 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_p',surface%pomb,surface%poma)
588 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_q',surface%qomb,surface%qoma)
589 end if
590
591 if (if_plot_gpspw ) then
592 call write_diag_single_level(out_dir,diag_unit_out,date,'gpspw_tpw',gpspw%tpwomb,gpspw%tpwoma)
593 end if
594
595 if (if_plot_gpsref ) then
596 call write_diag_multi_level_h(out_dir,diag_unit_out,date,'gps_ref',gpsref%refomb,gpsref%refoma)
597 !rizvi call write_diag_single_level(out_dir,diag_unit_out,date,'avgh_ref',avgh%refomb,avgh%refoma)
598 end if
599
600 if (if_plot_upr ) then
601 call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_u',upr%uomb,upr%uoma)
602 call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_v',upr%vomb,upr%voma)
603 call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_t',upr%tomb,upr%toma)
604 call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_q',upr%qomb,upr%qoma)
605 end if
606 ! Calculate next date:
607 call da_advance_cymdh( date, interval, new_date )
608 date = new_date
609 read(date(1:10), fmt='(i10)')cdate
610 end do ! End loop over date
611 if( if_plot_upr ) then
612 call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_u',gupr%uomb,gupr%uoma)
613 call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_v',gupr%vomb,gupr%voma)
614 call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_t',gupr%tomb,gupr%toma)
615 call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_q',gupr%qomb,gupr%qoma)
616 endif
617 if (if_plot_gpsref ) then
618 call write_diag_multi_level_h(out_dir,diag_unit_out,date,'ggps_ref',ggpsref%refomb,ggpsref%refoma)
619 endif
620
621 end do ! Loop over experiments
622 stop
623
624 1000 print*,' Error while reading unit ',diag_unit_in,' for experiment ',exp_dirs(iexp)
625 stop
626
627 contains
628
629 subroutine get_std_pr_level(prs, npr, stdp, nstd)
630
631 implicit none
632
633 integer, intent(in ) :: nstd
634 real, intent(in) :: stdp(nstd)
635 integer, intent(out) :: npr
636 real, intent(in) :: prs
637
638 real :: pr
639 integer :: k
640
641 pr = prs/100.0
642 if ( pr >= stdp(1) ) then
643 npr = 1
644 return
645 else if ( pr < stdp(nstd-1) ) then
646 npr = nstd
647 return
648 else
649 do k = 2,nstd - 1
650 if (pr >= stdp(k) ) then
651 npr = k
652 return
653 end if
654 end do
655 end if
656
657 end subroutine get_std_pr_level
658
659 subroutine update_stats(stats_omb, stats_oma, omb, oma)
660
661 implicit none
662 type(stats_value), intent(inout) :: stats_omb, stats_oma
663 real, intent (in) :: omb, oma
664
665 real :: x1, x2
666
667 stats_omb%num = stats_omb%num + 1
668 stats_oma%num = stats_omb%num
669
670 x1 = 1.0/ stats_omb%num
671 x2 = (stats_omb%num-1)*x1
672
673 stats_omb%bias = x2*stats_omb%bias + omb * x1
674 stats_oma%bias = x2*stats_oma%bias + oma * x1
675
676 stats_omb%abias = x2*stats_omb%abias + abs(omb) * x1
677 stats_oma%abias = x2*stats_oma%abias + abs(oma) * x1
678
679 stats_omb%rmse = x2*stats_omb%rmse + omb*omb * x1
680 stats_oma%rmse = x2*stats_oma%rmse + oma*oma * x1
681
682 end subroutine update_stats
683
684 subroutine write_diag_single_level(out_dir,ounit,ldate,obs_type,omb,oma)
685
686 implicit none
687
688 integer, intent(in) :: ounit
689 character*512,intent(in) :: out_dir
690 character*10,intent(in) :: ldate
691 character*(*),intent(in) :: obs_type
692 type (stats_value),intent(in) :: omb
693 type (stats_value),intent(in) :: oma
694
695 character*512 :: filename
696 integer :: ounit1, ounit2
697 real :: sigt,bar
698
699
700 ounit1 = ounit
701 ounit2 = ounit + 1
702
703 filename = trim(out_dir)//'/'//trim(obs_type)//'_omb.diag'
704 open (ounit1, file = trim(filename), form='formatted',status='unknown',position='append')
705 filename = trim(out_dir)//'/'//trim(obs_type)//'_oma.diag'
706 open (ounit2, file = trim(filename), form='formatted',status='unknown',position='append')
707 if ( omb%num <= 1 ) then
708 sigt=1. ; bar = rmiss
709 write(ounit1,'(1x,a10,1x,6(f6.2,1x))') ldate,rmiss, rmiss, rmiss, rmiss,bar,sigt
710 write(ounit2,'(1x,a10,1x,6(f6.2,1x))') ldate,rmiss, rmiss, rmiss, rmiss,bar,sigt
711 else
712 ! write(ounit1,'(5x,a10,4(2x,a9))') trim(obs_type),' Number','BIAS','ABIAS','RMSE '
713 if (index(obs_type,'_q') > 0 ) then
714 call sig_test(omb%num, omb%bias, omb%rmse, sigt,bar)
715 bar=bar*1000.0
716 write(ounit1,'(1x,a10,1x,i9,1x,5(f6.2,1x))') ldate,omb%num, 1000.0*omb%bias, 1000.0*omb%abias, 1000.0*sqrt(omb%rmse),bar,sigt
717 call sig_test(oma%num, oma%bias, oma%rmse, sigt,bar)
718 bar=bar*1000.0
719 write(ounit2,'(1x,a10,1x,i9,1x,5(f6.2,1x))') ldate,oma%num, 1000.0*oma%bias, 1000.0*oma%abias, 1000.0*sqrt(oma%rmse),bar,sigt
720 else if( index(obs_type,'_p') > 0 ) then
721 call sig_test(omb%num, omb%bias, omb%rmse, sigt,bar)
722 bar=bar/100.0
723 write(ounit1,'(1x,a10,1x,i9,1x,5(f6.2,1x))')ldate,omb%num, omb%bias/100.0, omb%abias/100.0, sqrt(omb%rmse)/100.0,bar,sigt
724 call sig_test(oma%num, oma%bias, oma%rmse, sigt,bar)
725 bar=bar/100.0
726 write(ounit2,'(1x,a10,1x,i9,5(1x,f6.2))') ldate,oma%num, oma%bias/100.0, oma%abias/100.0, sqrt(oma%rmse)/100.0,bar,sigt
727 else
728 call sig_test(omb%num, omb%bias, omb%rmse, sigt,bar)
729 write(ounit1,'(1x,a10,1x,i9,1x,5(f6.2,1x))') ldate,omb%num, omb%bias, omb%abias, sqrt(omb%rmse),bar,sigt
730 call sig_test(oma%num, oma%bias, oma%rmse, sigt,bar)
731 write(ounit2,'(1x,a10,1x,i9,5(1x,f6.2))') ldate,oma%num, oma%bias, oma%abias, sqrt(oma%rmse),bar,sigt
732 endif
733 !
734 end if
735 close(ounit1)
736 close(ounit2)
737
738 end subroutine write_diag_single_level
739
740 subroutine write_diag_multi_level(out_dir,ounit,ldate,obs_type,omb,oma)
741
742 implicit none
743
744 integer, intent(in) :: ounit
745 character*512,intent(in) :: out_dir
746 character*10,intent(in) :: ldate
747 character*(*),intent(in) :: obs_type
748 type (stats_value),intent(in) :: omb(nstd)
749 type (stats_value),intent(in) :: oma(nstd)
750
751 character*512 :: filename
752 integer :: k
753 real :: xnum(nstd)
754 integer :: num(nstd)
755 real, dimension(nstd) :: rmse, bias, abias,sigt,bar
756 integer :: ounit1, ounit2
757
758 ounit1 = ounit
759 ounit2 = ounit + 1
760
761 filename = trim(out_dir)//'/'//trim(obs_type)//'_omb.diag'
762 open (ounit1, file = trim(filename), form='formatted',status='unknown',position='append')
763 filename = trim(out_dir)//'/'//trim(obs_type)//'_oma.diag'
764 open (ounit2, file = trim(filename), form='formatted',status='unknown',position='append')
765
766 do k = 1, nstd
767 num(k) = omb(k)%num
768 if (num(k) <= 1 ) then
769 xnum(k) = rmiss
770 rmse(k) = rmiss
771 bias(k) = rmiss
772 abias(k) = rmiss
773 bar(k) = rmiss
774 sigt(k) = 1.0
775 else
776 if (index(obs_type,'_q') > 0 ) then
777 rmse(k) = sqrt(omb(k)%rmse) * 1000
778 bias(k) = omb(k)%bias * 1000
779 abias(k) = omb(k)%abias * 1000
780 call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
781 bar(k) = bar(k)*1000.
782 else
783 rmse(k) = sqrt(omb(k)%rmse)
784 bias(k) = omb(k)%bias
785 abias(k) = omb(k)%abias
786 call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
787 end if
788 xnum(k) = num(k)
789 end if
790 end do
791
792 write(ounit1,'(1x,a10,1x,16(6(1x,f9.2)))')ldate, (xnum(k), bias(k), abias(k),&
793 rmse(k),bar(k),sigt(k),k=1,nstd)
794
795 do k = 1, nstd
796 num(k) = oma(k)%num
797 if( num(k) <= 1 ) then
798 xnum(k) = rmiss
799 rmse(k) = rmiss
800 bias(k) = rmiss
801 abias(k) = rmiss
802 else
803 if (index(obs_type,'_q') > 0 ) then
804 rmse(k) = sqrt(oma(k)%rmse) * 1000
805 bias(k) = oma(k)%bias * 1000
806 abias(k) = oma(k)%abias * 1000
807 call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
808 bar(k) = bar(k)*1000.
809 else
810 rmse(k) = sqrt(oma(k)%rmse)
811 bias(k) = oma(k)%bias
812 abias(k) = oma(k)%abias
813 call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
814 end if
815 xnum(k) = num(k)
816 end if
817 end do
818 write(ounit2,'(1x,a10,1x,16(6(1x,f9.2)))')ldate, (xnum(k), bias(k), abias(k),&
819 rmse(k),bar(k),sigt(k),k=1,nstd)
820
821 close(ounit1)
822 close(ounit2)
823
824 end subroutine write_diag_multi_level
825 subroutine write_diag_multi_level_h(out_dir,ounit,date,obs_type,omb,oma)
826 implicit none
827 integer, intent(in) :: ounit
828 character*512,intent(in) :: out_dir
829 character*10,intent(in) :: date
830 character*(*),intent(in) :: obs_type
831 type (stats_value),intent(in) :: omb(nstdh)
832 type (stats_value),intent(in) :: oma(nstdh)
833 !
834 character*512 :: filename
835 integer :: k
836 real :: xnum(nstdh)
837 integer :: num(nstdh)
838 real, dimension(nstdh) :: rmse, bias, abias, sigt, bar
839
840 integer :: ounit1, ounit2
841 !
842 ounit1 = ounit
843 ounit2 = ounit + 1
844 !
845
846 filename = trim(out_dir)//'/'//trim(obs_type)//'_omb.diag'
847 open (ounit1, file = trim(filename), form='formatted',status='unknown',position='append')
848 filename = trim(out_dir)//'/'//trim(obs_type)//'_oma.diag'
849 open (ounit2, file = trim(filename), form='formatted',status='unknown',position='append')
850
851 do k = 1, nstdh
852 num(k) = omb(k)%num
853 if( num(k) <= 1 ) then
854 xnum(k) = rmiss
855 rmse(k) = rmiss
856 bias(k) = rmiss
857 abias(k) = rmiss
858 bar(k) = rmiss
859 sigt(k) = 1.0
860 else
861 if( index(obs_type,'_q') > 0 ) then
862
863 rmse(k) = sqrt(omb(k)%rmse) * 1000
864 bias(k) = omb(k)%bias * 1000
865 abias(k) = omb(k)%abias * 1000
866 call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
867 bar(k) = bar(k)*1000.
868 else
869 rmse(k) = sqrt(omb(k)%rmse)
870 bias(k) = omb(k)%bias
871 abias(k) = omb(k)%abias
872 call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
873 endif
874 xnum(k) = num(k)
875 endif
876 enddo
877 write(ounit1,'(1x,a10,1x,125(6(1x,f9.2)))')date, (xnum(k), bias(k), abias(k), &
878 rmse(k),bar(k),sigt(k), k=1,nstdh)
879
880 do k = 1, nstdh
881 num(k) = oma(k)%num
882 if( num(k) <= 1 ) then
883 xnum(k) = rmiss
884 rmse(k) = rmiss
885 bias(k) = rmiss
886 abias(k) = rmiss
887 bar(k) = rmiss
888 sigt(k) = 1.0
889 else
890 if( index(obs_type,'_q') > 0 ) then
891
892 rmse(k) = sqrt(oma(k)%rmse) * 1000
893 bias(k) = oma(k)%bias * 1000
894 abias(k) = oma(k)%abias * 1000
895 call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
896 bar(k) = bar(k)*1000.
897 else
898 rmse(k) = sqrt(oma(k)%rmse)
899 bias(k) = oma(k)%bias
900 abias(k) = oma(k)%abias
901 call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
902 endif
903 xnum(k) = num(k)
904 endif
905 enddo
906 write(ounit2,'(1x,a10,1x,125(6(1x,f9.2)))')date, (xnum(k), bias(k), abias(k), &
907 rmse(k),bar(k),sigt(k), k=1,nstdh)
908
909 !
910 close(ounit1)
911 close(ounit2)
912 !
913 end subroutine write_diag_multi_level_h
914 !
915 subroutine sig_test(num, bias, rmse, sigt,bar)
916 implicit none
917 integer, intent(in) :: num
918 real, intent(in) :: bias, rmse
919 real, intent(out) :: sigt, bar
920
921 real :: t_val, sd, tmp
922 !
923 sigt=0.
924 tmp = num/real(num-1)
925 ! sd = sqrt ( tmp*rmse - bias*bias/real((n-1)) )
926 sd = sqrt ( tmp*( rmse - bias*bias ) )
927 do k=2,34
928 if (real(num-1) < alpha(k,1)) exit
929 end do
930
931 t_val = bias*sqrt( real(num) ) /sd
932 bar = alpha(k-1,2) * sd /sqrt( real(num) )
933
934 if (abs(t_val) >= alpha(k-1,2)) sigt=1.
935
936 end subroutine sig_test
937
938 end program da_verif_obs