2 !---------------------------------------------------------------------------
6 ! Program to read diagnostics written in fort.50 by WRFVAR
7 ! and write in proper format to get ploted using PC-XL utility
9 ! Author: Syed RH Rizvi NCAR/MMM 05/25/2006
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 !---------------------------------------------------------------------------
17 use da_verif_obs_control
, only
: surface_type
, upr_type
, gpspw_type
, &
18 gpsref_type
, record1
, record2
, record3
, &
19 record4
, record5
, record6
, 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
, if_plot_tamdar
,interval
, stdp
, start_date
, &
27 if_plot_geoamv
, stdh
, num_miss
, &
28 wrf_file
, istart
, iend
, jstart
, jend
29 use da_verif_obs_init
, only
: initialize_surface_type
, initialize_upr_type
, &
30 initialize_gpspw_type
, initialize_gpsref_type
, da_advance_cymdh
, &
32 use da_verif_tools
, only
: map_info
, proj_merc
, proj_ps
,proj_lc
,proj_latlon
, &
33 da_llxy_wrf
,da_xyll
,da_map_set
38 character*20 :: obs_type
, dummy_c
41 integer :: n
, k
, kk
, l
, levels
, dummy_i
42 real :: lat
, lon
, press
, height
, dummy
43 real :: u_obs
, u_inv
, u_error
, u_inc
, &
44 v_obs
, v_inv
, v_error
, v_inc
, &
45 t_obs
, t_inv
, t_error
, t_inc
, &
46 p_obs
, p_inv
, p_error
, p_inc
, &
47 q_obs
, q_inv
, q_error
, q_inc
, &
48 spd_obs
, spd_inv
, spd_err
, spd_inc
49 real :: tpw_obs
, tpw_inv
, tpw_err
, tpw_inc
50 real :: ref_obs
, ref_inv
, ref_err
, ref_inc
51 integer :: u_qc
, v_qc
, t_qc
, p_qc
, q_qc
, tpw_qc
, spd_qc
, ref_qc
52 integer :: npr
, nht
, ier
, iexp
53 character*10 :: date
, new_date
! Current date (ccyymmddhh).
54 integer :: sdate
, cdate
, edate
! Starting, current ending dates.
55 logical :: if_write
, is_file
56 logical, allocatable
:: l_skip(:)
58 character(len
=512) :: out_dir
,filename
59 type (surface_type
) :: surface
60 type (upr_type
) :: upr
, gupr
61 type (gpspw_type
) :: gpspw
62 type (gpsref_type
) :: gpsref
, ggpsref
64 integer :: nx
, ny
, nz
, num_date
, idate
65 real :: dx
, cen_lat
, cen_lon
, truelat1
, truelat2
, stand_lon
66 integer :: map_proj_wrf
67 logical :: l_got_info
, inside
80 if_plot_rmse
= .false
.
81 if_plot_bias
= .false
.
82 if_plot_abias
= .false
.
84 if_plot_synop
= .false
.
85 if_plot_sonde_sfc
= .false
.
86 if_plot_metar
= .false
.
87 if_plot_ships
= .false
.
88 if_plot_qscat
= .false
.
89 if_plot_buoy
= .false
.
91 if_plot_sound
= .false
.
92 if_plot_geoamv
= .false
.
93 if_plot_polaramv
= .false
.
94 if_plot_profiler
= .false
.
95 if_plot_airep
= .false
.
96 if_plot_pilot
= .false
.
98 if_plot_gpspw
= .false
.
99 if_plot_gpsref
= .false
.
100 if_plot_airsret
= .false
.
102 if_plot_tamdar
= .false
.
104 file_path_string
= 'wrfvar/gts_omb_oma_01'
112 ! Read in namelist information defined in module define_cons_types
114 open ( unit
=nml_unit
, file
='namelist.plot_diag', STATUS
='OLD', &
116 read ( unit
=nml_unit
, nml
=record1
, IOSTAT
=ier
)
117 write ( unit
=*, nml
= record1
)
119 write (*,*) 'error in reading namelist record1'
123 read ( unit
=nml_unit
, nml
=record2
, iostat
=ier
)
124 write ( unit
=*, nml
= record2
)
126 write (*,*) 'error in reading namelist record2'
129 read ( unit
=nml_unit
, nml
=record3
, iostat
=ier
)
130 write ( unit
=*, nml
= record3
)
132 write (*,*) 'error in reading namelist record3'
135 read ( unit
=nml_unit
, nml
=record4
, iostat
=ier
)
136 write ( unit
=*, nml
= record4
)
138 write (*,*) 'error in reading namelist record4'
141 read ( unit
=nml_unit
, nml
=record5
, iostat
=ier
)
142 write ( unit
=*, nml
= record5
)
144 write (*,*) 'error in reading namelist record5'
147 read ( unit
=nml_unit
, nml
=record6
, iostat
=ier
)
149 write (*,*) 'error in reading namelist record6'
152 write ( unit
=*, nml
= record6
)
155 call initialize_t_tab
158 if ( l_got_info
) then
160 istart
= max(1, istart
)
162 jstart
= max(1, jstart
)
166 if_plot_surface
= .false
.
167 if (if_plot_synop
.or
. if_plot_metar
.or
. if_plot_ships
.or
. if_plot_buoy
.or
. &
168 if_plot_sonde_sfc
.or
. if_plot_qscat
) if_plot_surface
= .true
.
170 if_plot_upr
= .false
.
171 if (if_plot_sound
.or
. if_plot_pilot
.or
. if_plot_profiler
.or
. &
172 if_plot_geoamv
.or
. if_plot_polaramv
.or
. if_plot_airep
.or
. &
173 if_plot_airsret
.or
. if_plot_tamdar
) if_plot_upr
= .true
.
175 read(start_date(1:10), fmt
='(i10)')sdate
176 read(end_date(1:10), fmt
='(i10)')edate
177 write(6,'(4a)')' Diag Starting date ', start_date
, ' Ending date ', end_date
178 write(6,'(a,i8,a)')' Interval between dates = ', interval
, ' hours.'
182 do while ( date
<= end_date
)
183 num_date
= num_date
+ 1
184 call da_advance_cymdh(date
, interval
, date
)
186 allocate(l_skip(num_date
))
189 ! check for missing dates
190 idate
= 0 ! index of date
192 do while ( date
<= end_date
)
195 filename
= TRIM(exp_dirs(iexp
))//'/'//date
//'/'//trim(file_path_string
)
196 inquire ( file
=trim(filename
), exist
=is_file
)
197 if ( .not
. is_file
) then
198 l_skip(idate
) = .true
.
200 if ( l_skip(idate
) ) exit
202 call da_advance_cymdh(date
, interval
, date
)
205 !---------------------------------------------------------------------------
206 ! Loop over experiments
213 call initialize_upr_type(gupr
)
214 call initialize_gpsref_type(ggpsref
)
216 do while ( cdate
<= edate
)
217 ! Initialize various types
218 call initialize_surface_type(surface
)
219 call initialize_upr_type(upr
)
220 call initialize_gpspw_type(gpspw
)
221 call initialize_gpsref_type(gpsref
)
225 ! construct file name
226 filename
= TRIM(exp_dirs(iexp
))//'/'//date
//'/'//trim(file_path_string
)
228 inquire ( file
=trim(filename
), exist
=is_file
)
229 if ( l_skip(idate
) .or
. .not
. is_file
) then
230 print*, 'skipping file ', trim(filename
)
232 ! Write output on outunit
233 out_dir
=trim(out_dirs(iexp
))
234 if (if_plot_surface
) then
235 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_u',surface
%uomb
,surface
%uoma
)
236 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_v',surface
%vomb
,surface
%voma
)
237 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_t',surface
%tomb
,surface
%toma
)
238 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_p',surface
%pomb
,surface
%poma
)
239 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_q',surface
%qomb
,surface
%qoma
)
242 if (if_plot_gpspw
) then
243 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'gpspw_tpw',gpspw
%tpwomb
,gpspw
%tpwoma
)
246 if (if_plot_gpsref
) then
247 call write_diag_multi_level_h(out_dir
,diag_unit_out
,date
,'gps_ref',gpsref
%refomb
,gpsref
%refoma
)
250 if (if_plot_upr
) then
251 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'upr_u',upr
%uomb
,upr
%uoma
)
252 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'upr_v',upr
%vomb
,upr
%voma
)
253 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'upr_t',upr
%tomb
,upr
%toma
)
254 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'upr_q',upr
%qomb
,upr
%qoma
)
256 ! Calculate next date:
257 call da_advance_cymdh( date
, interval
, new_date
)
259 read(date(1:10), fmt
='(i10)')cdate
263 open (unit
=diag_unit_in
, file
=trim(filename
), form
='formatted', &
264 status
='old', iostat
=ier
)
268 read(diag_unit_in
,'(a20,i8)', end=2000, err
= 1000)obs_type
,num_obs
269 if (index( obs_type
,'synop') > 0 ) then
270 if (if_plot_synop
) if_write
= .true
.
273 elseif (index( obs_type
,'metar') > 0 ) then
274 if (if_plot_metar
) if_write
= .true
.
277 elseif (index( obs_type
,'ships') > 0 ) then
278 if (if_plot_ships
) if_write
= .true
.
281 elseif (index( obs_type
,'buoy' ) > 0 ) then
282 if (if_plot_buoy
) if_write
= .true
.
285 elseif (index( obs_type
,'sonde_sfc') > 0 ) then
286 if (if_plot_sonde_sfc
) if_write
= .true
.
289 elseif (index( obs_type
,'polaramv') > 0) then
290 if (if_plot_polaramv
) if_write
= .true
.
293 elseif (index( obs_type
,'geoamv' ) > 0) then
294 if (if_plot_geoamv
) if_write
= .true
.
297 elseif (index( obs_type
,'gpspw') > 0) then
298 if ( if_plot_gpspw
) if_write
= .true
.
301 elseif (index( obs_type
,'sound') > 0) then
302 if (if_plot_sound
) if_write
= .true
.
305 elseif (index( obs_type
,'airep') > 0) then
306 if (if_plot_airep
) if_write
= .true
.
309 elseif (index( obs_type
,'pilot') > 0) then
310 if (if_plot_pilot
) if_write
= .true
.
313 elseif (index( obs_type
,'profiler') > 0) then
314 if (if_plot_profiler
) if_write
= .true
.
317 elseif (index( obs_type
,'ssmir') > 0) then
320 elseif (index( obs_type
,'ssmiT') > 0) then
323 elseif (index( obs_type
,'satem') > 0) then
326 elseif (index( obs_type
,'ssmt1') > 0) then
329 elseif (index( obs_type
,'ssmt2') > 0) then
332 elseif (index( obs_type
,'qscat') > 0) then
333 if (if_plot_qscat
) if_write
= .true
.
335 elseif (index( obs_type
,'gpsref' ) > 0) then
336 if (if_plot_gpsref
) if_write
= .true
.
338 elseif (index( obs_type
,'airsr') > 0) then
339 if (if_plot_airsret
) if_write
= .true
.
342 elseif (index( obs_type
,'tamdar') > 0) then
343 if (if_plot_tamdar
) if_write
= .true
.
347 print*,' Got unknown OBS_TYPE ',obs_type(1:20),' on unit ',diag_unit_in
351 10 continue ! Synop, Metar, Ships, Buoy , Sonde_sfc
353 if ( num_obs
> 0 ) then
355 read(diag_unit_in
,'(i8)')levels
357 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
358 kk
, l
, stn_id
, & ! Station
359 lat
, lon
, press
, & ! Lat/lon, pressure
360 u_obs
, u_inv
, u_qc
, u_error
, u_inc
, &
361 v_obs
, v_inv
, v_qc
, v_error
, v_inc
, &
362 t_obs
, t_inv
, t_qc
, t_error
, t_inc
, &
363 p_obs
, p_inv
, p_qc
, p_error
, p_inc
, &
364 q_obs
, q_inv
, q_qc
, q_error
, q_inc
366 if ( l_got_info
) call check_domain(lat
, lon
, inside
)
368 if (u_qc
>= 0) call update_stats(surface
%uomb
, surface
%uoma
, u_inv
, u_inc
)
369 if (v_qc
>= 0) call update_stats(surface
%vomb
, surface
%voma
, v_inv
, v_inc
)
370 if (t_qc
>= 0) call update_stats(surface
%tomb
, surface
%toma
, t_inv
, t_inc
)
371 if (p_qc
>= 0) call update_stats(surface
%pomb
, surface
%poma
, p_inv
, p_inc
)
372 if (q_qc
>= 0) call update_stats(surface
%qomb
, surface
%qoma
, q_inv
, q_inc
)
375 end do ! loop over levels
376 end do ! loop over Obs
380 20 continue ! Polar or Geo AMV's
382 if ( num_obs
> 0 ) then
384 read(diag_unit_in
,'(i8)')levels
386 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
387 kk
, l
, stn_id
, & ! Station
388 lat
, lon
, press
, & ! Lat/lon, pressure
389 u_obs
, u_inv
, u_qc
, u_error
, u_inc
, &
390 v_obs
, v_inv
, v_qc
, v_error
, v_inc
392 if (if_write
.and
. press
> 0 ) then
393 call get_std_pr_level(press
, npr
, stdp
, nstd
)
394 if ( l_got_info
) call check_domain(lat
, lon
, inside
)
396 if( u_qc
>= 0 .and
. npr
> 0 ) then
397 call update_stats(upr
%uomb(npr
),upr
%uoma(npr
),u_inv
,u_inc
)
398 call update_stats(gupr
%uomb(npr
),gupr
%uoma(npr
),u_inv
,u_inc
)
400 if( v_qc
>= 0 .and
. npr
> 0 ) then
401 call update_stats(upr
%vomb(npr
),upr
%voma(npr
),v_inv
,v_inc
)
402 call update_stats(gupr
%vomb(npr
),gupr
%voma(npr
),v_inv
,v_inc
)
406 end do ! loop over levels
407 end do ! loop over Obs
414 if ( num_obs
> 0 ) then
416 read(diag_unit_in
,'(i8)')levels
418 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
419 kk
, l
, stn_id
, & ! Station
420 lat
, lon
, dummy
, & ! Lat/lon, dummy
421 tpw_obs
, tpw_inv
, tpw_qc
, tpw_err
, tpw_inc
423 if ( l_got_info
) call check_domain(lat
, lon
, inside
)
425 if (tpw_qc
>= 0) call update_stats(gpspw
%tpwomb
,gpspw
%tpwoma
,tpw_inv
,tpw_inc
)
428 end do ! loop over levels
429 end do ! loop over Obs
436 ! [6] Transfer sound obs:
438 if ( num_obs
> 0 ) then
440 read(diag_unit_in
,'(i8)')levels
442 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
443 kk
,l
, stn_id
, & ! Station
444 lat
, lon
, press
, & ! Lat/lon, dummy
445 u_obs
, u_inv
, u_qc
, u_error
, u_inc
, &
446 v_obs
, v_inv
, v_qc
, v_error
, v_inc
, &
447 t_obs
, t_inv
, t_qc
, t_error
, t_inc
, &
448 q_obs
, q_inv
, q_qc
, q_error
, q_inc
449 if (if_write
.and
. press
> 0 ) then
450 call get_std_pr_level(press
, npr
, stdp
, nstd
)
451 if ( l_got_info
) call check_domain(lat
, lon
, inside
)
453 if( u_qc
>= 0 .and
. npr
> 0 ) then
454 call update_stats(upr
%uomb(npr
),upr
%uoma(npr
),u_inv
,u_inc
)
455 call update_stats(gupr
%uomb(npr
),gupr
%uoma(npr
),u_inv
,u_inc
)
457 if( v_qc
>= 0 .and
. npr
> 0 ) then
458 call update_stats(upr
%vomb(npr
),upr
%voma(npr
),v_inv
,v_inc
)
459 call update_stats(gupr
%vomb(npr
),gupr
%voma(npr
),v_inv
,v_inc
)
461 if( t_qc
>= 0 .and
. npr
> 0 ) then
462 call update_stats(upr
%tomb(npr
),upr
%toma(npr
),t_inv
,t_inc
)
463 call update_stats(gupr
%tomb(npr
),gupr
%toma(npr
),t_inv
,t_inc
)
465 if( q_qc
>= 0 .and
. npr
> 0 ) then
466 call update_stats(upr
%qomb(npr
),upr
%qoma(npr
),q_inv
,q_inc
)
467 call update_stats(gupr
%qomb(npr
),gupr
%qoma(npr
),q_inv
,q_inc
)
472 end do ! loop over levels
473 end do ! loop over Obs
479 if ( num_obs
> 0 ) then
481 read(diag_unit_in
,'(i8)') levels
483 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
484 kk
,l
, stn_id
, & ! Station
485 lat
, lon
, press
, & ! Lat/lon, dummy
486 u_obs
, u_inv
, u_qc
, u_error
, u_inc
, &
487 v_obs
, v_inv
, v_qc
, v_error
, v_inc
, &
488 t_obs
, t_inv
, t_qc
, t_error
, t_inc
489 if (if_write
.and
. press
> 0 ) then
490 call get_std_pr_level(press
, npr
, stdp
, nstd
)
491 if ( l_got_info
) call check_domain(lat
, lon
, inside
)
493 if( u_qc
>= 0 .and
. npr
> 0 ) then
494 call update_stats(upr
%uomb(npr
),upr
%uoma(npr
),u_inv
,u_inc
)
495 call update_stats(gupr
%uomb(npr
),gupr
%uoma(npr
),u_inv
,u_inc
)
497 if( v_qc
>= 0 .and
. npr
> 0 ) then
498 call update_stats(upr
%vomb(npr
),upr
%voma(npr
),v_inv
,v_inc
)
499 call update_stats(gupr
%vomb(npr
),gupr
%voma(npr
),v_inv
,v_inc
)
501 if( t_qc
>= 0 .and
. npr
> 0 ) then
502 call update_stats(upr
%tomb(npr
),upr
%toma(npr
),t_inv
,t_inc
)
503 call update_stats(gupr
%tomb(npr
),gupr
%toma(npr
),t_inv
,t_inc
)
508 end do ! loop over levels
509 end do ! loop over Obs
514 60 continue ! Pilot & Profiler
516 ! [8] Transfer pilot obs:
517 if ( num_obs
> 0 ) then
519 read(diag_unit_in
,'(i8)')levels
521 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
522 kk
,l
, stn_id
, & ! Station
523 lat
, lon
, press
, & ! Lat/lon, dummy
524 u_obs
, u_inv
, u_qc
, u_error
, u_inc
, &
525 v_obs
, v_inv
, v_qc
, v_error
, v_inc
526 if (if_write
.and
. press
> 0 ) then
527 call get_std_pr_level(press
, npr
, stdp
, nstd
)
528 if ( l_got_info
) call check_domain(lat
, lon
, inside
)
530 if( u_qc
>= 0 .and
. npr
> 0 ) then
531 call update_stats(upr
%uomb(npr
),upr
%uoma(npr
),u_inv
,u_inc
)
532 call update_stats(gupr
%uomb(npr
),gupr
%uoma(npr
),u_inv
,u_inc
)
534 if( v_qc
>= 0 .and
. npr
> 0 ) then
535 call update_stats(upr
%vomb(npr
),upr
%voma(npr
),v_inv
,v_inc
)
536 call update_stats(gupr
%vomb(npr
),gupr
%voma(npr
),v_inv
,v_inc
)
541 end do ! loop over levels
542 end do ! loop over Obs
546 70 continue ! SSMI retrievals
548 if ( num_obs
> 0 ) then
550 read(diag_unit_in
,'(i8)')levels
552 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
553 kk
,l
, stn_id
, & ! Station
554 lat
, lon
, dummy
, & ! Lat/lon, dummy
555 spd_obs
, spd_inv
, spd_qc
, spd_err
, spd_inc
562 80 continue ! SSMI radiance
564 if ( num_obs
> 0 ) then
566 read(diag_unit_in
,*)dummy_c
567 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,7(2f17.7,i8,2f17.7))', err
= 1000)&
568 k
, l
, stn_id
, & ! Station
569 lat
, lon
, dummy
, & ! Lat/lon, dummy
570 dummy
, dummy
, dummy_i
, dummy
, dummy
, &
571 dummy
, dummy
, dummy_i
, dummy
, dummy
, &
572 dummy
, dummy
, dummy_i
, dummy
, dummy
, &
573 dummy
, dummy
, dummy_i
, dummy
, dummy
, &
574 dummy
, dummy
, dummy_i
, dummy
, dummy
, &
575 dummy
, dummy
, dummy_i
, dummy
, dummy
, &
576 dummy
, dummy
, dummy_i
, dummy
, dummy
583 if ( num_obs
> 0 ) then
585 read(diag_unit_in
,'(i8)') levels
587 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
588 kk
,l
, stn_id
, & ! Station
589 lat
, lon
, dummy
, & ! Lat/lon, dummy
590 dummy
,dummy
, dummy_i
, dummy
, dummy
591 end do ! loop over levels
592 end do ! loop over Obs
597 100 continue ! SSMT1 & 2
599 if ( num_obs
> 0 ) then
601 read(diag_unit_in
,'(i8)') levels
603 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
604 kk
,l
, stn_id
, & ! Station
605 lat
, lon
, dummy
, & ! Lat/lon, dummy
606 dummy
,dummy
, dummy_i
, dummy
, dummy
607 end do ! loop over levels
608 end do ! loop over obs
613 110 continue ! Scatrometer winds
615 if ( num_obs
> 0 ) then
617 read(diag_unit_in
,'(i8)') levels
619 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
620 kk
, l
, stn_id
, & ! Station
621 lat
, lon
, press
, & ! Lat/lon, dummy
622 u_obs
, u_inv
, u_qc
, u_error
, u_inc
, &
623 v_obs
, v_inv
, v_qc
, v_error
, v_inc
625 if ( l_got_info
) call check_domain(lat
, lon
, inside
)
627 if (u_qc
>= 0) call update_stats(surface
%uomb
,surface
%uoma
,u_inv
,u_inc
)
628 if (v_qc
>= 0) call update_stats(surface
%vomb
,surface
%voma
,v_inv
,v_inc
)
631 end do ! loop over levels
632 end do ! loop over obs
636 120 continue ! Gpsref
638 IF ( num_obs
> 0 ) THEN
640 read(diag_unit_in
,'(i8)') levels
642 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
643 kk
, l
, stn_id
, & ! Station
644 lat
, lon
, height
, & ! Lat/lon, dummy
645 ref_obs
, ref_inv
, ref_qc
, ref_err
, ref_inc
646 if (if_write
.and
. height
> 0.0) then
647 call get_std_ht_level(height
, nht
, stdh
, nstdh
)
648 if ( l_got_info
) call check_domain(lat
, lon
, inside
)
650 if ( ref_qc
>= 0) then
651 call update_stats(gpsref
%refomb(nht
),gpsref
%refoma(nht
),ref_inv
,ref_inc
)
652 call update_stats(ggpsref
%refomb(nht
),ggpsref
%refoma(nht
),ref_inv
,ref_inc
)
656 END DO ! loop over levels
657 END DO ! loop over Obs
660 !---------------------------------------------------------------------
662 130 continue ! AIRSRET
664 if ( num_obs
> 0 ) then
666 read(diag_unit_in
,'(i8)')levels
668 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
669 kk
,l
, stn_id
, & ! Station
670 lat
, lon
, press
, & ! Lat/lon, dummy
671 t_obs
, t_inv
, t_qc
, t_error
, t_inc
, &
672 q_obs
, q_inv
, q_qc
, q_error
, q_inc
673 if (if_write
.and
. press
> 0 ) then
674 call get_std_pr_level(press
, npr
, stdp
, nstd
)
675 if ( l_got_info
) call check_domain(lat
, lon
, inside
)
677 if( t_qc
>= 0 .and
. npr
> 0 ) then
678 call update_stats(upr
%tomb(npr
),upr
%toma(npr
),t_inv
,t_inc
)
679 call update_stats(gupr
%tomb(npr
),gupr
%toma(npr
),t_inv
,t_inc
)
681 if( q_qc
>= 0 .and
. npr
> 0 ) then
682 call update_stats(upr
%qomb(npr
),upr
%qoma(npr
),q_inv
,q_inc
)
683 call update_stats(gupr
%qomb(npr
),gupr
%qoma(npr
),q_inv
,q_inc
)
687 end do ! loop over levels
688 end do ! loop over obs
694 if ( num_obs
> 0 ) then
696 read(diag_unit_in
,'(i8)')levels
698 read(diag_unit_in
,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err
= 1000)&
701 u_obs
, u_inv
, u_qc
, u_error
, u_inc
, &
702 v_obs
, v_inv
, v_qc
, v_error
, v_inc
, &
703 t_obs
, t_inv
, t_qc
, t_error
, t_inc
, &
704 q_obs
, q_inv
, q_qc
, q_error
, q_inc
705 if (if_write
.and
. press
> 0 ) then
706 call get_std_pr_level(press
, npr
, stdp
, nstd
)
709 call update_stats(upr
%uomb(npr
),upr
%uoma(npr
),u_inv
,u_inc
)
710 call update_stats(gupr
%uomb(npr
),gupr
%uoma(npr
),u_inv
,u_inc
)
713 call update_stats(upr
%vomb(npr
),upr
%voma(npr
),v_inv
,v_inc
)
714 call update_stats(gupr
%vomb(npr
),gupr
%voma(npr
),v_inv
,v_inc
)
717 call update_stats(upr
%tomb(npr
),upr
%toma(npr
),t_inv
,t_inc
)
718 call update_stats(gupr
%tomb(npr
),gupr
%toma(npr
),t_inv
,t_inc
)
721 call update_stats(upr
%qomb(npr
),upr
%qoma(npr
),q_inv
,q_inc
)
722 call update_stats(gupr
%qomb(npr
),gupr
%qoma(npr
),q_inv
,q_inc
)
732 ! Now process the diagnostics
736 ! Write output on outunit
737 out_dir
=trim(out_dirs(iexp
))
738 if (if_plot_surface
) then
739 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_u',surface
%uomb
,surface
%uoma
)
740 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_v',surface
%vomb
,surface
%voma
)
741 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_t',surface
%tomb
,surface
%toma
)
742 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_p',surface
%pomb
,surface
%poma
)
743 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'surface_q',surface
%qomb
,surface
%qoma
)
746 if (if_plot_gpspw
) then
747 call write_diag_single_level(out_dir
,diag_unit_out
,date
,'gpspw_tpw',gpspw
%tpwomb
,gpspw
%tpwoma
)
750 if (if_plot_gpsref
) then
751 call write_diag_multi_level_h(out_dir
,diag_unit_out
,date
,'gps_ref',gpsref
%refomb
,gpsref
%refoma
)
752 !rizvi call write_diag_single_level(out_dir,diag_unit_out,date,'avgh_ref',avgh%refomb,avgh%refoma)
755 if (if_plot_upr
) then
756 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'upr_u',upr
%uomb
,upr
%uoma
)
757 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'upr_v',upr
%vomb
,upr
%voma
)
758 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'upr_t',upr
%tomb
,upr
%toma
)
759 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'upr_q',upr
%qomb
,upr
%qoma
)
761 ! Calculate next date:
762 call da_advance_cymdh( date
, interval
, new_date
)
764 read(date(1:10), fmt
='(i10)')cdate
765 end do ! End loop over date
766 if( if_plot_upr
) then
767 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'gupr_u',gupr
%uomb
,gupr
%uoma
)
768 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'gupr_v',gupr
%vomb
,gupr
%voma
)
769 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'gupr_t',gupr
%tomb
,gupr
%toma
)
770 call write_diag_multi_level(out_dir
,diag_unit_out
,date
,'gupr_q',gupr
%qomb
,gupr
%qoma
)
772 if (if_plot_gpsref
) then
773 call write_diag_multi_level_h(out_dir
,diag_unit_out
,date
,'ggps_ref',ggpsref
%refomb
,ggpsref
%refoma
)
776 end do ! Loop over experiments
779 1000 print*,' Error while reading unit ',diag_unit_in
,' for experiment ',exp_dirs(iexp
)
784 subroutine get_std_pr_level(prs
, npr
, stdp
, nstd
)
788 integer, intent(in
) :: nstd
789 real, intent(in
) :: stdp(nstd
)
790 integer, intent(out
) :: npr
791 real, intent(in
) :: prs
796 npr
= num_miss
! initialize as a missing value
798 if ( pr
>= stdp(1) ) then
801 else if ( pr
< stdp(nstd
-1) .and
. pr
>= stdp(nstd
) ) then
806 if (pr
>= stdp(k
) ) then
813 end subroutine get_std_pr_level
815 subroutine get_std_ht_level(height
, nht
, stdh
, nstdh
)
819 integer, intent(in
) :: nstdh
820 real, intent(in
) :: stdh(nstdh
)
821 integer, intent(out
) :: nht
822 real, intent(in
) :: height
827 ht
= height
*0.001 ! m to km
828 if ( ht
<= stdh(1) ) then
831 else if ( ht
> stdh(nstdh
-1) ) then
836 if ( ht
<= stdh(k
) ) then
843 end subroutine get_std_ht_level
845 subroutine update_stats(stats_omb
, stats_oma
, omb
, oma
)
848 type(stats_value
), intent(inout
) :: stats_omb
, stats_oma
849 real, intent (in
) :: omb
, oma
853 stats_omb
%num
= stats_omb
%num
+ 1
854 stats_oma
%num
= stats_omb
%num
856 x1
= 1.0/ stats_omb
%num
857 x2
= (stats_omb
%num
-1)*x1
859 stats_omb
%bias
= x2
*stats_omb
%bias
+ omb
* x1
860 stats_oma
%bias
= x2
*stats_oma
%bias
+ oma
* x1
862 stats_omb
%abias
= x2
*stats_omb
%abias
+ abs(omb
) * x1
863 stats_oma
%abias
= x2
*stats_oma
%abias
+ abs(oma
) * x1
865 stats_omb
%rmse
= x2
*stats_omb
%rmse
+ omb
*omb
* x1
866 stats_oma
%rmse
= x2
*stats_oma
%rmse
+ oma
*oma
* x1
868 end subroutine update_stats
870 subroutine write_diag_single_level(out_dir
,ounit
,ldate
,obs_type
,omb
,oma
)
874 integer, intent(in
) :: ounit
875 character*512,intent(in
) :: out_dir
876 character*10,intent(in
) :: ldate
877 character*(*),intent(in
) :: obs_type
878 type (stats_value
),intent(in
) :: omb
879 type (stats_value
),intent(in
) :: oma
881 character*512 :: filename
882 integer :: ounit1
, ounit2
889 filename
= trim(out_dir
)//'/'//trim(obs_type
)//'_omb.diag'
890 open (ounit1
, file
= trim(filename
), form
='formatted',status
='unknown',position
='append')
891 filename
= trim(out_dir
)//'/'//trim(obs_type
)//'_oma.diag'
892 open (ounit2
, file
= trim(filename
), form
='formatted',status
='unknown',position
='append')
893 if ( omb
%num
<= 1 ) then
894 sigt
=1. ; bar
= rmiss
895 write(ounit1
,'(1x,a10,1x,6(f6.2,1x))') ldate
,rmiss
, rmiss
, rmiss
, rmiss
,bar
,sigt
896 write(ounit2
,'(1x,a10,1x,6(f6.2,1x))') ldate
,rmiss
, rmiss
, rmiss
, rmiss
,bar
,sigt
898 ! write(ounit1,'(5x,a10,4(2x,a9))') trim(obs_type),' Number','BIAS','ABIAS','RMSE '
899 if (index(obs_type
,'_q') > 0 ) then
900 call sig_test(omb
%num
, omb
%bias
, omb
%rmse
, sigt
,bar
)
902 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
903 call sig_test(oma
%num
, oma
%bias
, oma
%rmse
, sigt
,bar
)
905 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
906 else if( index(obs_type
,'_p') > 0 ) then
907 call sig_test(omb
%num
, omb
%bias
, omb
%rmse
, sigt
,bar
)
909 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
910 call sig_test(oma
%num
, oma
%bias
, oma
%rmse
, sigt
,bar
)
912 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
914 call sig_test(omb
%num
, omb
%bias
, omb
%rmse
, sigt
,bar
)
915 write(ounit1
,'(1x,a10,1x,i9,1x,5(f6.2,1x))') ldate
,omb
%num
, omb
%bias
, omb
%abias
, sqrt(omb
%rmse
),bar
,sigt
916 call sig_test(oma
%num
, oma
%bias
, oma
%rmse
, sigt
,bar
)
917 write(ounit2
,'(1x,a10,1x,i9,5(1x,f6.2))') ldate
,oma
%num
, oma
%bias
, oma
%abias
, sqrt(oma
%rmse
),bar
,sigt
924 end subroutine write_diag_single_level
926 subroutine write_diag_multi_level(out_dir
,ounit
,ldate
,obs_type
,omb
,oma
)
930 integer, intent(in
) :: ounit
931 character*512,intent(in
) :: out_dir
932 character*10,intent(in
) :: ldate
933 character*(*),intent(in
) :: obs_type
934 type (stats_value
),intent(in
) :: omb(nstd
)
935 type (stats_value
),intent(in
) :: oma(nstd
)
937 character*512 :: filename
941 real, dimension(nstd
) :: rmse
, bias
, abias
,sigt
,bar
942 integer :: ounit1
, ounit2
947 filename
= trim(out_dir
)//'/'//trim(obs_type
)//'_omb.diag'
948 open (ounit1
, file
= trim(filename
), form
='formatted',status
='unknown',position
='append')
949 filename
= trim(out_dir
)//'/'//trim(obs_type
)//'_oma.diag'
950 open (ounit2
, file
= trim(filename
), form
='formatted',status
='unknown',position
='append')
954 if (num(k
) <= 1 ) then
962 if (index(obs_type
,'_q') > 0 ) then
963 rmse(k
) = sqrt(omb(k
)%rmse
) * 1000
964 bias(k
) = omb(k
)%bias
* 1000
965 abias(k
) = omb(k
)%abias
* 1000
966 call sig_test(num(k
), omb(k
)%bias
, omb(k
)%rmse
, sigt(k
),bar(k
))
967 bar(k
) = bar(k
)*1000.
969 rmse(k
) = sqrt(omb(k
)%rmse
)
970 bias(k
) = omb(k
)%bias
971 abias(k
) = omb(k
)%abias
972 call sig_test(num(k
), omb(k
)%bias
, omb(k
)%rmse
, sigt(k
),bar(k
))
978 write(ounit1
,'(1x,a10,1x,16(6(1x,f12.2)))')ldate
, (xnum(k
), bias(k
), abias(k
),&
979 rmse(k
),bar(k
),sigt(k
),k
=1,nstd
)
983 if( num(k
) <= 1 ) then
989 if (index(obs_type
,'_q') > 0 ) then
990 rmse(k
) = sqrt(oma(k
)%rmse
) * 1000
991 bias(k
) = oma(k
)%bias
* 1000
992 abias(k
) = oma(k
)%abias
* 1000
993 call sig_test(num(k
), oma(k
)%bias
, oma(k
)%rmse
, sigt(k
),bar(k
))
994 bar(k
) = bar(k
)*1000.
996 rmse(k
) = sqrt(oma(k
)%rmse
)
997 bias(k
) = oma(k
)%bias
998 abias(k
) = oma(k
)%abias
999 call sig_test(num(k
), oma(k
)%bias
, oma(k
)%rmse
, sigt(k
),bar(k
))
1004 write(ounit2
,'(1x,a10,1x,16(6(1x,f12.2)))')ldate
, (xnum(k
), bias(k
), abias(k
),&
1005 rmse(k
),bar(k
),sigt(k
),k
=1,nstd
)
1010 end subroutine write_diag_multi_level
1011 subroutine write_diag_multi_level_h(out_dir
,ounit
,date
,obs_type
,omb
,oma
)
1013 integer, intent(in
) :: ounit
1014 character*512,intent(in
) :: out_dir
1015 character*10,intent(in
) :: date
1016 character*(*),intent(in
) :: obs_type
1017 type (stats_value
),intent(in
) :: omb(nstdh
)
1018 type (stats_value
),intent(in
) :: oma(nstdh
)
1020 character*512 :: filename
1023 integer :: num(nstdh
)
1024 real, dimension(nstdh
) :: rmse
, bias
, abias
, sigt
, bar
1026 integer :: ounit1
, ounit2
1032 filename
= trim(out_dir
)//'/'//trim(obs_type
)//'_omb.diag'
1033 open (ounit1
, file
= trim(filename
), form
='formatted',status
='unknown',position
='append')
1034 filename
= trim(out_dir
)//'/'//trim(obs_type
)//'_oma.diag'
1035 open (ounit2
, file
= trim(filename
), form
='formatted',status
='unknown',position
='append')
1039 if( num(k
) <= 1 ) then
1047 if( index(obs_type
,'_q') > 0 ) then
1049 rmse(k
) = sqrt(omb(k
)%rmse
) * 1000
1050 bias(k
) = omb(k
)%bias
* 1000
1051 abias(k
) = omb(k
)%abias
* 1000
1052 call sig_test(num(k
), omb(k
)%bias
, omb(k
)%rmse
, sigt(k
),bar(k
))
1053 bar(k
) = bar(k
)*1000.
1055 rmse(k
) = sqrt(omb(k
)%rmse
)
1056 bias(k
) = omb(k
)%bias
1057 abias(k
) = omb(k
)%abias
1058 call sig_test(num(k
), omb(k
)%bias
, omb(k
)%rmse
, sigt(k
),bar(k
))
1063 write(ounit1
,'(1x,a10,1x,150(6(1x,f12.2)))')date
, (xnum(k
), bias(k
), abias(k
), &
1064 rmse(k
),bar(k
),sigt(k
), k
=1,nstdh
)
1068 if( num(k
) <= 1 ) then
1076 if( index(obs_type
,'_q') > 0 ) then
1078 rmse(k
) = sqrt(oma(k
)%rmse
) * 1000
1079 bias(k
) = oma(k
)%bias
* 1000
1080 abias(k
) = oma(k
)%abias
* 1000
1081 call sig_test(num(k
), oma(k
)%bias
, oma(k
)%rmse
, sigt(k
),bar(k
))
1082 bar(k
) = bar(k
)*1000.
1084 rmse(k
) = sqrt(oma(k
)%rmse
)
1085 bias(k
) = oma(k
)%bias
1086 abias(k
) = oma(k
)%abias
1087 call sig_test(num(k
), oma(k
)%bias
, oma(k
)%rmse
, sigt(k
),bar(k
))
1092 write(ounit2
,'(1x,a10,1x,150(6(1x,f12.2)))')date
, (xnum(k
), bias(k
), abias(k
), &
1093 rmse(k
),bar(k
),sigt(k
), k
=1,nstdh
)
1099 end subroutine write_diag_multi_level_h
1101 subroutine sig_test(num
, bias
, rmse
, sigt
,bar
)
1103 integer, intent(in
) :: num
1104 real, intent(in
) :: bias
, rmse
1105 real, intent(out
) :: sigt
, bar
1107 real :: t_val
, sd
, tmp
1110 tmp
= num
/real(num
-1)
1111 ! sd = sqrt ( tmp*rmse - bias*bias/real((n-1)) )
1112 sd
= sqrt ( tmp
*( rmse
- bias
*bias
) )
1114 if (real(num
-1) < alpha(k
,1)) exit
1117 t_val
= bias
*sqrt( real(num
) ) /sd
1118 bar
= alpha(k
-1,2) * sd
/sqrt( real(num
) )
1120 if (abs(t_val
) >= alpha(k
-1,2)) sigt
=1.
1122 end subroutine sig_test
1125 subroutine get_fileinfo
1127 #
include "netcdf.inc"
1130 l_got_info
= .false
.
1131 iost(1) = nf_open(trim(wrf_file
), NF_NOWRITE
, ncid
)
1132 if ( iost(1) /= NF_NOERR
) then
1133 print*, 'INFO: wrf_file: '//trim(wrf_file
)//' does not exist for retrieving mapping info'
1136 print*, 'Retrieving mapping info from wrf_file: ',trim(wrf_file
)
1138 iost(2) = nf_get_att_int(ncid
, NF_GLOBAL
, 'WEST-EAST_GRID_DIMENSION', nx
)
1139 iost(3) = nf_get_att_int(ncid
, NF_GLOBAL
, 'SOUTH-NORTH_GRID_DIMENSION', ny
)
1140 iost(4) = nf_get_att_int(ncid
, NF_GLOBAL
, 'BOTTOM-TOP_GRID_DIMENSION', nz
)
1141 iost(5) = nf_get_att_double(ncid
, NF_GLOBAL
, 'DX', dx
)
1142 iost(6) = nf_get_att_double(ncid
, NF_GLOBAL
, 'CEN_LAT', cen_lat
)
1143 iost(7) = nf_get_att_double(ncid
, NF_GLOBAL
, 'CEN_LON', cen_lon
)
1144 iost(8) = nf_get_att_double(ncid
, NF_GLOBAL
, 'TRUELAT1', truelat1
)
1145 iost(9) = nf_get_att_double(ncid
, NF_GLOBAL
, 'TRUELAT2', truelat2
)
1146 iost(10) = nf_get_att_double(ncid
, NF_GLOBAL
, 'STAND_LON', stand_lon
)
1147 iost(11) = nf_get_att_int(ncid
, NF_GLOBAL
, 'MAP_PROJ', map_proj_wrf
)
1148 iost(12) = nf_close(ncid
)
1149 if ( .not
. any(iost
/=NF_NOERR
) ) then
1152 print*, 'nx, ny, nz, dx, map_proj, cen_lat, cen_lon, truelat1, truelat2, stand_lon = ', &
1153 nx
, ny
, nz
, dx
, map_proj_wrf
, cen_lat
, cen_lon
, truelat1
, truelat2
, stand_lon
1154 end subroutine get_fileinfo
1156 subroutine set_mapinfo
1158 integer :: map_proj_util
1160 real :: start_x
, start_y
1163 if ( map_proj_wrf
== 0 .or
. map_proj_wrf
== 6 ) then
1164 map_proj_util
= proj_latlon
1165 else if ( map_proj_wrf
== 3 ) then
1166 map_proj_util
= proj_merc
1167 else if ( map_proj_wrf
== 1 ) then
1168 map_proj_util
= proj_lc
1169 else if ( map_proj_wrf
== 2 ) then
1170 map_proj_util
= proj_ps
1172 call da_map_set(map_proj_util
, cen_lat
,cen_lon
, xref
, yref
, dx
, &
1173 stand_lon
, truelat1
, truelat2
, truelat1
, stand_lon
, map_info
)
1174 !call da_llxy_wrf(map_info, cen_lat, cen_lon, start_x, start_y)
1175 end subroutine set_mapinfo
1177 subroutine check_domain(lat
, lon
, inside
)
1179 real, intent(in
) :: lat
, lon
1180 logical, intent(out
) :: inside
1182 call da_llxy_wrf(map_info
, lat
, lon
, xx
, yy
)
1184 if ( xx
>= istart
.and
. xx
<= iend
.and
. &
1185 yy
>= jstart
.and
. yy
<= jend
) then
1188 end subroutine check_domain
1190 end program da_verif_obs