?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
wrf svn trunk commit r4103
[wrffire.git] / wrfv2_fire / var / da / da_verif_obs / da_verif_obs.f90
blob92501474d140470adbc92783577dfb4f6fff2f23
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, stdh
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, height, 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, nht, 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/gts_omb_oma_01'
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, height, & ! Lat/lon, dummy
537 ref_obs, ref_inv, ref_qc, ref_err, ref_inc
538 if (if_write .and. height > 0.0) then
539 call get_std_ht_level(height, nht, stdh, nstdh)
540 if ( ref_qc >= 0) then
541 call update_stats(gpsref%refomb(nht),gpsref%refoma(nht),ref_inv,ref_inc)
542 call update_stats(ggpsref%refomb(nht),ggpsref%refoma(nht),ref_inv,ref_inc)
543 end if
544 end if
545 END DO ! loop over levels
546 END DO ! loop over Obs
547 ENDIF
548 go to 1
549 !---------------------------------------------------------------------
550
551 130 continue ! AIRSRET
552
553 if ( num_obs > 0 ) then
554 do n = 1, num_obs
555 read(diag_unit_in,'(i8)')levels
556 do k = 1, levels
557 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
558 kk,l, stn_id, & ! Station
559 lat, lon, press, & ! Lat/lon, dummy
560 t_obs, t_inv, t_qc, t_error, t_inc, &
561 q_obs, q_inv, q_qc, q_error, q_inc
562 if (if_write .and. press > 0 ) then
563 call get_std_pr_level(press, npr, stdp, nstd)
564 if( t_qc >= 0) then
565 call update_stats(upr%tomb(npr),upr%toma(npr),t_inv,t_inc)
566 call update_stats(gupr%tomb(npr),gupr%toma(npr),t_inv,t_inc)
567 endif
568 if( q_qc >= 0) then
569 call update_stats(upr%qomb(npr),upr%qoma(npr),q_inv,q_inc)
570 call update_stats(gupr%qomb(npr),gupr%qoma(npr),q_inv,q_inc)
571 endif
572 end if
573 end do ! loop over levels
574 end do ! loop over obs
575 end if
576 goto 1
577
578 ! Now process the diagnostics
579 2000 continue
580
581 close (diag_unit_in)
582 ! Write output on outunit
583 out_dir=trim(out_dirs(iexp))
584 if (if_plot_surface ) then
585 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_u',surface%uomb,surface%uoma)
586 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_v',surface%vomb,surface%voma)
587 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_t',surface%tomb,surface%toma)
588 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_p',surface%pomb,surface%poma)
589 call write_diag_single_level(out_dir,diag_unit_out,date,'surface_q',surface%qomb,surface%qoma)
590 end if
591
592 if (if_plot_gpspw ) then
593 call write_diag_single_level(out_dir,diag_unit_out,date,'gpspw_tpw',gpspw%tpwomb,gpspw%tpwoma)
594 end if
595
596 if (if_plot_gpsref ) then
597 call write_diag_multi_level_h(out_dir,diag_unit_out,date,'gps_ref',gpsref%refomb,gpsref%refoma)
598 !rizvi call write_diag_single_level(out_dir,diag_unit_out,date,'avgh_ref',avgh%refomb,avgh%refoma)
599 end if
600
601 if (if_plot_upr ) then
602 call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_u',upr%uomb,upr%uoma)
603 call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_v',upr%vomb,upr%voma)
604 call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_t',upr%tomb,upr%toma)
605 call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_q',upr%qomb,upr%qoma)
606 end if
607 ! Calculate next date:
608 call da_advance_cymdh( date, interval, new_date )
609 date = new_date
610 read(date(1:10), fmt='(i10)')cdate
611 end do ! End loop over date
612 if( if_plot_upr ) then
613 call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_u',gupr%uomb,gupr%uoma)
614 call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_v',gupr%vomb,gupr%voma)
615 call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_t',gupr%tomb,gupr%toma)
616 call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_q',gupr%qomb,gupr%qoma)
617 endif
618 if (if_plot_gpsref ) then
619 call write_diag_multi_level_h(out_dir,diag_unit_out,date,'ggps_ref',ggpsref%refomb,ggpsref%refoma)
620 endif
621
622 end do ! Loop over experiments
623 stop
624
625 1000 print*,' Error while reading unit ',diag_unit_in,' for experiment ',exp_dirs(iexp)
626 stop
627
628 contains
629
630 subroutine get_std_pr_level(prs, npr, stdp, nstd)
631
632 implicit none
633
634 integer, intent(in ) :: nstd
635 real, intent(in) :: stdp(nstd)
636 integer, intent(out) :: npr
637 real, intent(in) :: prs
638
639 real :: pr
640 integer :: k
641
642 pr = prs/100.0
643 if ( pr >= stdp(1) ) then
644 npr = 1
645 return
646 else if ( pr < stdp(nstd-1) ) then
647 npr = nstd
648 return
649 else
650 do k = 2,nstd - 1
651 if (pr >= stdp(k) ) then
652 npr = k
653 return
654 end if
655 end do
656 end if
657
658 end subroutine get_std_pr_level
659
660 subroutine get_std_ht_level(height, nht, stdh, nstdh)
661
662 implicit none
663
664 integer, intent(in ) :: nstdh
665 real, intent(in) :: stdh(nstdh)
666 integer, intent(out) :: nht
667 real, intent(in) :: height
668
669 real :: ht
670 integer :: k
671
672 ht = height*0.001 ! m to km
673 if ( ht <= stdh(1) ) then
674 nht = 1
675 return
676 else if ( ht > stdh(nstdh-1) ) then
677 nht = nstdh
678 return
679 else
680 do k = 2,nstdh - 1
681 if ( ht <= stdh(k) ) then
682 nht = k
683 return
684 end if
685 end do
686 end if
687
688 end subroutine get_std_ht_level
689
690 subroutine update_stats(stats_omb, stats_oma, omb, oma)
691
692 implicit none
693 type(stats_value), intent(inout) :: stats_omb, stats_oma
694 real, intent (in) :: omb, oma
695
696 real :: x1, x2
697
698 stats_omb%num = stats_omb%num + 1
699 stats_oma%num = stats_omb%num
700
701 x1 = 1.0/ stats_omb%num
702 x2 = (stats_omb%num-1)*x1
703
704 stats_omb%bias = x2*stats_omb%bias + omb * x1
705 stats_oma%bias = x2*stats_oma%bias + oma * x1
706
707 stats_omb%abias = x2*stats_omb%abias + abs(omb) * x1
708 stats_oma%abias = x2*stats_oma%abias + abs(oma) * x1
709
710 stats_omb%rmse = x2*stats_omb%rmse + omb*omb * x1
711 stats_oma%rmse = x2*stats_oma%rmse + oma*oma * x1
712
713 end subroutine update_stats
714
715 subroutine write_diag_single_level(out_dir,ounit,ldate,obs_type,omb,oma)
716
717 implicit none
718
719 integer, intent(in) :: ounit
720 character*512,intent(in) :: out_dir
721 character*10,intent(in) :: ldate
722 character*(*),intent(in) :: obs_type
723 type (stats_value),intent(in) :: omb
724 type (stats_value),intent(in) :: oma
725
726 character*512 :: filename
727 integer :: ounit1, ounit2
728 real :: sigt,bar
729
730
731 ounit1 = ounit
732 ounit2 = ounit + 1
733
734 filename = trim(out_dir)//'/'//trim(obs_type)//'_omb.diag'
735 open (ounit1, file = trim(filename), form='formatted',status='unknown',position='append')
736 filename = trim(out_dir)//'/'//trim(obs_type)//'_oma.diag'
737 open (ounit2, file = trim(filename), form='formatted',status='unknown',position='append')
738 if ( omb%num <= 1 ) then
739 sigt=1. ; bar = rmiss
740 write(ounit1,'(1x,a10,1x,6(f6.2,1x))') ldate,rmiss, rmiss, rmiss, rmiss,bar,sigt
741 write(ounit2,'(1x,a10,1x,6(f6.2,1x))') ldate,rmiss, rmiss, rmiss, rmiss,bar,sigt
742 else
743 ! write(ounit1,'(5x,a10,4(2x,a9))') trim(obs_type),' Number','BIAS','ABIAS','RMSE '
744 if (index(obs_type,'_q') > 0 ) then
745 call sig_test(omb%num, omb%bias, omb%rmse, sigt,bar)
746 bar=bar*1000.0
747 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
748 call sig_test(oma%num, oma%bias, oma%rmse, sigt,bar)
749 bar=bar*1000.0
750 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
751 else if( index(obs_type,'_p') > 0 ) then
752 call sig_test(omb%num, omb%bias, omb%rmse, sigt,bar)
753 bar=bar/100.0
754 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
755 call sig_test(oma%num, oma%bias, oma%rmse, sigt,bar)
756 bar=bar/100.0
757 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
758 else
759 call sig_test(omb%num, omb%bias, omb%rmse, sigt,bar)
760 write(ounit1,'(1x,a10,1x,i9,1x,5(f6.2,1x))') ldate,omb%num, omb%bias, omb%abias, sqrt(omb%rmse),bar,sigt
761 call sig_test(oma%num, oma%bias, oma%rmse, sigt,bar)
762 write(ounit2,'(1x,a10,1x,i9,5(1x,f6.2))') ldate,oma%num, oma%bias, oma%abias, sqrt(oma%rmse),bar,sigt
763 endif
764 !
765 end if
766 close(ounit1)
767 close(ounit2)
768
769 end subroutine write_diag_single_level
770
771 subroutine write_diag_multi_level(out_dir,ounit,ldate,obs_type,omb,oma)
772
773 implicit none
774
775 integer, intent(in) :: ounit
776 character*512,intent(in) :: out_dir
777 character*10,intent(in) :: ldate
778 character*(*),intent(in) :: obs_type
779 type (stats_value),intent(in) :: omb(nstd)
780 type (stats_value),intent(in) :: oma(nstd)
781
782 character*512 :: filename
783 integer :: k
784 real :: xnum(nstd)
785 integer :: num(nstd)
786 real, dimension(nstd) :: rmse, bias, abias,sigt,bar
787 integer :: ounit1, ounit2
788
789 ounit1 = ounit
790 ounit2 = ounit + 1
791
792 filename = trim(out_dir)//'/'//trim(obs_type)//'_omb.diag'
793 open (ounit1, file = trim(filename), form='formatted',status='unknown',position='append')
794 filename = trim(out_dir)//'/'//trim(obs_type)//'_oma.diag'
795 open (ounit2, file = trim(filename), form='formatted',status='unknown',position='append')
796
797 do k = 1, nstd
798 num(k) = omb(k)%num
799 if (num(k) <= 1 ) then
800 xnum(k) = rmiss
801 rmse(k) = rmiss
802 bias(k) = rmiss
803 abias(k) = rmiss
804 bar(k) = rmiss
805 sigt(k) = 1.0
806 else
807 if (index(obs_type,'_q') > 0 ) then
808 rmse(k) = sqrt(omb(k)%rmse) * 1000
809 bias(k) = omb(k)%bias * 1000
810 abias(k) = omb(k)%abias * 1000
811 call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
812 bar(k) = bar(k)*1000.
813 else
814 rmse(k) = sqrt(omb(k)%rmse)
815 bias(k) = omb(k)%bias
816 abias(k) = omb(k)%abias
817 call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
818 end if
819 xnum(k) = num(k)
820 end if
821 end do
822
823 write(ounit1,'(1x,a10,1x,16(6(1x,f9.2)))')ldate, (xnum(k), bias(k), abias(k),&
824 rmse(k),bar(k),sigt(k),k=1,nstd)
825
826 do k = 1, nstd
827 num(k) = oma(k)%num
828 if( num(k) <= 1 ) then
829 xnum(k) = rmiss
830 rmse(k) = rmiss
831 bias(k) = rmiss
832 abias(k) = rmiss
833 else
834 if (index(obs_type,'_q') > 0 ) then
835 rmse(k) = sqrt(oma(k)%rmse) * 1000
836 bias(k) = oma(k)%bias * 1000
837 abias(k) = oma(k)%abias * 1000
838 call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
839 bar(k) = bar(k)*1000.
840 else
841 rmse(k) = sqrt(oma(k)%rmse)
842 bias(k) = oma(k)%bias
843 abias(k) = oma(k)%abias
844 call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
845 end if
846 xnum(k) = num(k)
847 end if
848 end do
849 write(ounit2,'(1x,a10,1x,16(6(1x,f9.2)))')ldate, (xnum(k), bias(k), abias(k),&
850 rmse(k),bar(k),sigt(k),k=1,nstd)
851
852 close(ounit1)
853 close(ounit2)
854
855 end subroutine write_diag_multi_level
856 subroutine write_diag_multi_level_h(out_dir,ounit,date,obs_type,omb,oma)
857 implicit none
858 integer, intent(in) :: ounit
859 character*512,intent(in) :: out_dir
860 character*10,intent(in) :: date
861 character*(*),intent(in) :: obs_type
862 type (stats_value),intent(in) :: omb(nstdh)
863 type (stats_value),intent(in) :: oma(nstdh)
864 !
865 character*512 :: filename
866 integer :: k
867 real :: xnum(nstdh)
868 integer :: num(nstdh)
869 real, dimension(nstdh) :: rmse, bias, abias, sigt, bar
870
871 integer :: ounit1, ounit2
872 !
873 ounit1 = ounit
874 ounit2 = ounit + 1
875 !
876
877 filename = trim(out_dir)//'/'//trim(obs_type)//'_omb.diag'
878 open (ounit1, file = trim(filename), form='formatted',status='unknown',position='append')
879 filename = trim(out_dir)//'/'//trim(obs_type)//'_oma.diag'
880 open (ounit2, file = trim(filename), form='formatted',status='unknown',position='append')
881
882 do k = 1, nstdh
883 num(k) = omb(k)%num
884 if( num(k) <= 1 ) then
885 xnum(k) = rmiss
886 rmse(k) = rmiss
887 bias(k) = rmiss
888 abias(k) = rmiss
889 bar(k) = rmiss
890 sigt(k) = 1.0
891 else
892 if( index(obs_type,'_q') > 0 ) then
893
894 rmse(k) = sqrt(omb(k)%rmse) * 1000
895 bias(k) = omb(k)%bias * 1000
896 abias(k) = omb(k)%abias * 1000
897 call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
898 bar(k) = bar(k)*1000.
899 else
900 rmse(k) = sqrt(omb(k)%rmse)
901 bias(k) = omb(k)%bias
902 abias(k) = omb(k)%abias
903 call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
904 endif
905 xnum(k) = num(k)
906 endif
907 enddo
908 write(ounit1,'(1x,a10,1x,125(6(1x,f9.2)))')date, (xnum(k), bias(k), abias(k), &
909 rmse(k),bar(k),sigt(k), k=1,nstdh)
910
911 do k = 1, nstdh
912 num(k) = oma(k)%num
913 if( num(k) <= 1 ) then
914 xnum(k) = rmiss
915 rmse(k) = rmiss
916 bias(k) = rmiss
917 abias(k) = rmiss
918 bar(k) = rmiss
919 sigt(k) = 1.0
920 else
921 if( index(obs_type,'_q') > 0 ) then
922
923 rmse(k) = sqrt(oma(k)%rmse) * 1000
924 bias(k) = oma(k)%bias * 1000
925 abias(k) = oma(k)%abias * 1000
926 call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
927 bar(k) = bar(k)*1000.
928 else
929 rmse(k) = sqrt(oma(k)%rmse)
930 bias(k) = oma(k)%bias
931 abias(k) = oma(k)%abias
932 call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
933 endif
934 xnum(k) = num(k)
935 endif
936 enddo
937 write(ounit2,'(1x,a10,1x,125(6(1x,f9.2)))')date, (xnum(k), bias(k), abias(k), &
938 rmse(k),bar(k),sigt(k), k=1,nstdh)
939
940 !
941 close(ounit1)
942 close(ounit2)
943 !
944 end subroutine write_diag_multi_level_h
945 !
946 subroutine sig_test(num, bias, rmse, sigt,bar)
947 implicit none
948 integer, intent(in) :: num
949 real, intent(in) :: bias, rmse
950 real, intent(out) :: sigt, bar
951
952 real :: t_val, sd, tmp
953 !
954 sigt=0.
955 tmp = num/real(num-1)
956 ! sd = sqrt ( tmp*rmse - bias*bias/real((n-1)) )
957 sd = sqrt ( tmp*( rmse - bias*bias ) )
958 do k=2,34
959 if (real(num-1) < alpha(k,1)) exit
960 end do
961
962 t_val = bias*sqrt( real(num) ) /sd
963 bar = alpha(k-1,2) * sd /sqrt( real(num) )
964
965 if (abs(t_val) >= alpha(k-1,2)) sigt=1.
966
967 end subroutine sig_test
968
969 end program da_verif_obs
A coupled weather-fire forecasting model built on top of Weather Research and Forecasting (WRF).