?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge remote-tracking branch 'origin/release-v4.6.1'
[WRF.git] / var / da / da_verif_obs / da_verif_obs.f90
blob1cdf27b81b1a919b443e3c574f027827a9c57c66
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, 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 , &
31 initialize_t_tab
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
34
35 implicit none
36
37 integer :: num_obs
38 character*20 :: obs_type, dummy_c
39
40 character*5 :: stn_id
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(:)
57
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
63
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
68
69 inside = .true.
70
71 nml_unit = 10
72 diag_unit_in = 50
73 diag_unit_out = 20
74 info_unit = 30
75
76 exp_num = 0
77 exp_dirs = ''
78 out_dirs = ''
79
80 if_plot_rmse = .false.
81 if_plot_bias = .false.
82 if_plot_abias = .false.
83
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.
90
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.
97
98 if_plot_gpspw = .false.
99 if_plot_gpsref = .false.
100 if_plot_airsret = .false.
101
102 if_plot_tamdar = .false.
103
104 file_path_string = 'wrfvar/gts_omb_oma_01'
105 wrf_file = 'foo.nc'
106
107 istart = 1
108 iend = 10000
109 jstart = 1
110 jend = 10000
111
112 ! Read in namelist information defined in module define_cons_types
113
114 open ( unit=nml_unit, file='namelist.plot_diag', STATUS='OLD', &
115 form='formatted' )
116 read ( unit=nml_unit, nml=record1, IOSTAT=ier )
117 write ( unit=*, nml = record1 )
118 if ( ier /= 0 ) then
119 write (*,*) 'error in reading namelist record1'
120 stop
121 end if
122
123 read ( unit=nml_unit, nml=record2, iostat=ier )
124 write ( unit=*, nml = record2 )
125 if ( ier /= 0 ) then
126 write (*,*) 'error in reading namelist record2'
127 stop
128 end if
129 read ( unit=nml_unit, nml=record3, iostat=ier )
130 write ( unit=*, nml = record3 )
131 if ( ier /= 0 ) then
132 write (*,*) 'error in reading namelist record3'
133 stop
134 end if
135 read ( unit=nml_unit, nml=record4, iostat=ier )
136 write ( unit=*, nml = record4 )
137 if ( ier /= 0 ) then
138 write (*,*) 'error in reading namelist record4'
139 stop
140 end if
141 read ( unit=nml_unit, nml=record5, iostat=ier )
142 write ( unit=*, nml = record5 )
143 if ( ier /= 0 ) then
144 write (*,*) 'error in reading namelist record5'
145 stop
146 end if
147 read ( unit=nml_unit, nml=record6, iostat=ier )
148 if ( ier /= 0 ) then
149 write (*,*) 'error in reading namelist record6'
150 !stop
151 else
152 write ( unit=*, nml = record6 )
153 end if
154 close(nml_unit)
155 call initialize_t_tab
156
157 call get_fileinfo
158 if ( l_got_info ) then
159 call set_mapinfo
160 istart = max(1, istart)
161 iend = min(nx, iend)
162 jstart = max(1, jstart)
163 jend = min(ny, jend)
164 end if
165
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.
169
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.
174
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.'
179
180 num_date = 0
181 date = start_date
182 do while ( date <= end_date )
183 num_date = num_date + 1
184 call da_advance_cymdh(date, interval, date)
185 end do
186 allocate(l_skip(num_date))
187 l_skip(:) = .false.
188
189 ! check for missing dates
190 idate = 0 ! index of date
191 date = start_date
192 do while ( date <= end_date )
193 idate = idate + 1
194 do iexp = 1, exp_num
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.
199 end if
200 if ( l_skip(idate) ) exit
201 end do
202 call da_advance_cymdh(date, interval, date)
203 end do
204
205 !---------------------------------------------------------------------------
206 ! Loop over experiments
207
208 do iexp =1,exp_num
209
210 idate = 0
211 date = start_date
212 cdate = sdate
213 call initialize_upr_type(gupr)
214 call initialize_gpsref_type(ggpsref)
215
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)
222
223 idate = idate + 1
224
225 ! construct file name
226 filename = TRIM(exp_dirs(iexp))//'/'//date//'/'//trim(file_path_string)
227
228 inquire ( file=trim(filename), exist=is_file)
229 if ( l_skip(idate) .or. .not. is_file ) then
230 print*, 'skipping file ', trim(filename)
231 !stop
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)
240 end if
241
242 if (if_plot_gpspw ) then
243 call write_diag_single_level(out_dir,diag_unit_out,date,'gpspw_tpw',gpspw%tpwomb,gpspw%tpwoma)
244 end if
245
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)
248 end if
249
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)
255 end if
256 ! Calculate next date:
257 call da_advance_cymdh( date, interval, new_date )
258 date = new_date
259 read(date(1:10), fmt='(i10)')cdate
260 cycle
261 end if
262
263 open (unit=diag_unit_in, file=trim(filename), form='formatted', &
264 status='old', iostat=ier)
265 1 continue
266
267 if_write = .false.
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.
271 goto 10
272
273 elseif (index( obs_type,'metar') > 0 ) then
274 if (if_plot_metar ) if_write = .true.
275 goto 10
276
277 elseif (index( obs_type,'ships') > 0 ) then
278 if (if_plot_ships ) if_write = .true.
279 goto 10
280
281 elseif (index( obs_type,'buoy' ) > 0 ) then
282 if (if_plot_buoy ) if_write = .true.
283 goto 10
284
285 elseif (index( obs_type,'sonde_sfc') > 0 ) then
286 if (if_plot_sonde_sfc ) if_write = .true.
287 goto 10
288
289 elseif (index( obs_type,'polaramv') > 0) then
290 if (if_plot_polaramv ) if_write = .true.
291 goto 20
292
293 elseif (index( obs_type,'geoamv' ) > 0) then
294 if (if_plot_geoamv ) if_write = .true.
295 goto 20
296
297 elseif (index( obs_type,'gpspw') > 0) then
298 if ( if_plot_gpspw ) if_write = .true.
299 goto 30
300
301 elseif (index( obs_type,'sound') > 0) then
302 if (if_plot_sound ) if_write = .true.
303 goto 40
304
305 elseif (index( obs_type,'airep') > 0) then
306 if (if_plot_airep ) if_write = .true.
307 goto 50
308
309 elseif (index( obs_type,'pilot') > 0) then
310 if (if_plot_pilot ) if_write = .true.
311 goto 60
312
313 elseif (index( obs_type,'profiler') > 0) then
314 if (if_plot_profiler ) if_write = .true.
315 goto 60
316
317 elseif (index( obs_type,'ssmir') > 0) then
318 goto 70
319
320 elseif (index( obs_type,'ssmiT') > 0) then
321 goto 80
322
323 elseif (index( obs_type,'satem') > 0) then
324 goto 90
325
326 elseif (index( obs_type,'ssmt1') > 0) then
327 goto 100
328
329 elseif (index( obs_type,'ssmt2') > 0) then
330 goto 100
331
332 elseif (index( obs_type,'qscat') > 0) then
333 if (if_plot_qscat ) if_write = .true.
334 goto 110
335 elseif (index( obs_type,'gpsref' ) > 0) then
336 if (if_plot_gpsref ) if_write = .true.
337 goto 120
338 elseif (index( obs_type,'airsr') > 0) then
339 if (if_plot_airsret ) if_write = .true.
340 goto 130
341
342 elseif (index( obs_type,'tamdar') > 0) then
343 if (if_plot_tamdar ) if_write = .true.
344 goto 140
345
346 else
347 print*,' Got unknown OBS_TYPE ',obs_type(1:20),' on unit ',diag_unit_in
348 stop
349 end if
350
351 10 continue ! Synop, Metar, Ships, Buoy , Sonde_sfc
352
353 if ( num_obs > 0 ) then
354 do n = 1, num_obs
355 read(diag_unit_in,'(i8)')levels
356 do k = 1, 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
365 if (if_write) then
366 if ( l_got_info ) call check_domain(lat, lon, inside)
367 if ( inside ) then
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)
373 end if
374 end if
375 end do ! loop over levels
376 end do ! loop over Obs
377 end if
378 goto 1
379
380 20 continue ! Polar or Geo AMV's
381
382 if ( num_obs > 0 ) then
383 do n = 1, num_obs
384 read(diag_unit_in,'(i8)')levels
385 do k = 1, 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
391
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)
395 if ( inside ) then
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)
399 endif
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)
403 endif
404 end if
405 end if
406 end do ! loop over levels
407 end do ! loop over Obs
408 end if
409
410 goto 1
411
412 30 continue ! Gpspw
413
414 if ( num_obs > 0 ) then
415 do n = 1, num_obs
416 read(diag_unit_in,'(i8)')levels
417 do k = 1, 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
422 if (if_write) then
423 if ( l_got_info ) call check_domain(lat, lon, inside)
424 if ( inside ) then
425 if (tpw_qc >= 0) call update_stats(gpspw%tpwomb,gpspw%tpwoma,tpw_inv,tpw_inc)
426 end if
427 end if
428 end do ! loop over levels
429 end do ! loop over Obs
430 end if
431
432 goto 1
433
434 40 continue ! Sound
435
436 ! [6] Transfer sound obs:
437
438 if ( num_obs > 0 ) then
439 do n = 1, num_obs
440 read(diag_unit_in,'(i8)')levels
441 do k = 1, 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)
452 if ( inside ) then
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)
456 endif
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)
460 endif
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)
464 endif
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)
468 endif
469 end if
470 !
471 end if
472 end do ! loop over levels
473 end do ! loop over Obs
474 end if
475 goto 1
476
477 50 continue ! Airep
478
479 if ( num_obs > 0 ) then
480 do n = 1, num_obs
481 read(diag_unit_in,'(i8)') levels
482 do k = 1, 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)
492 if ( inside ) then
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)
496 endif
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)
500 endif
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)
504 endif
505 end if
506
507 end if
508 end do ! loop over levels
509 end do ! loop over Obs
510 end if
511
512 goto 1
513
514 60 continue ! Pilot & Profiler
515
516 ! [8] Transfer pilot obs:
517 if ( num_obs > 0 ) then
518 do n = 1, num_obs
519 read(diag_unit_in,'(i8)')levels
520 do k = 1, 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)
529 if ( inside ) then
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)
533 endif
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)
537 endif
538 end if
539
540 end if
541 end do ! loop over levels
542 end do ! loop over Obs
543 end if
544 goto 1
545
546 70 continue ! SSMI retrievals
547
548 if ( num_obs > 0 ) then
549 do n = 1, num_obs
550 read(diag_unit_in,'(i8)')levels
551 do k = 1, 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
556 end do
557 end do
558 end if
559
560 goto 1
561
562 80 continue ! SSMI radiance
563
564 if ( num_obs > 0 ) then
565 do n = 1, num_obs
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
577 end do
578 end if
579 goto 1
580
581 90 continue ! SATEM
582
583 if ( num_obs > 0 ) then
584 do n = 1, num_obs
585 read(diag_unit_in,'(i8)') levels
586 do k = 1, 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
593 end if
594
595 goto 1
596
597 100 continue ! SSMT1 & 2
598
599 if ( num_obs > 0 ) then
600 do n = 1, num_obs
601 read(diag_unit_in,'(i8)') levels
602 do k = 1, 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
609 end if
610
611 goto 1
612
613 110 continue ! Scatrometer winds
614
615 if ( num_obs > 0 ) then
616 do n = 1, num_obs
617 read(diag_unit_in,'(i8)') levels
618 do k = 1, 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
624 if (if_write) then
625 if ( l_got_info ) call check_domain(lat, lon, inside)
626 if ( inside ) then
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)
629 end if
630 end if
631 end do ! loop over levels
632 end do ! loop over obs
633 end if
634 goto 1
635
636 120 continue ! Gpsref
637
638 IF ( num_obs > 0 ) THEN
639 DO n = 1, num_obs
640 read(diag_unit_in,'(i8)') levels
641 DO k = 1, 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)
649 if ( inside ) then
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)
653 end if
654 end if
655 end if
656 END DO ! loop over levels
657 END DO ! loop over Obs
658 ENDIF
659 go to 1
660 !---------------------------------------------------------------------
661
662 130 continue ! AIRSRET
663
664 if ( num_obs > 0 ) then
665 do n = 1, num_obs
666 read(diag_unit_in,'(i8)')levels
667 do k = 1, 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)
676 if ( inside ) then
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)
680 endif
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)
684 endif
685 end if
686 end if
687 end do ! loop over levels
688 end do ! loop over obs
689 end if
690 goto 1
691
692 140 continue
693
694 if ( num_obs > 0 ) then
695 do n = 1, num_obs
696 read(diag_unit_in,'(i8)')levels
697 do k = 1, levels
698 read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
699 kk,l, stn_id, &
700 lat, lon, press, &
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)
707
708 if( u_qc >= 0) then
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)
711 endif
712 if( v_qc >= 0) then
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)
715 endif
716 if( t_qc >= 0) then
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)
719 endif
720 if( q_qc >= 0) then
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)
723 endif
724
725 end if
726 end do
727 end do
728 end if
729 goto 1
730
731
732 ! Now process the diagnostics
733 2000 continue
734
735 close (diag_unit_in)
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)
744 end if
745
746 if (if_plot_gpspw ) then
747 call write_diag_single_level(out_dir,diag_unit_out,date,'gpspw_tpw',gpspw%tpwomb,gpspw%tpwoma)
748 end if
749
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)
753 end if
754
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)
760 end if
761 ! Calculate next date:
762 call da_advance_cymdh( date, interval, new_date )
763 date = 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)
771 endif
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)
774 endif
775
776 end do ! Loop over experiments
777 stop
778
779 1000 print*,' Error while reading unit ',diag_unit_in,' for experiment ',exp_dirs(iexp)
780 stop
781
782 contains
783
784 subroutine get_std_pr_level(prs, npr, stdp, nstd)
785
786 implicit none
787
788 integer, intent(in ) :: nstd
789 real, intent(in) :: stdp(nstd)
790 integer, intent(out) :: npr
791 real, intent(in) :: prs
792
793 real :: pr
794 integer :: k
795
796 npr = num_miss ! initialize as a missing value
797 pr = prs/100.0
798 if ( pr >= stdp(1) ) then
799 npr = 1
800 return
801 else if ( pr < stdp(nstd-1) .and. pr >= stdp(nstd) ) then
802 npr = nstd
803 return
804 else
805 do k = 2,nstd - 1
806 if (pr >= stdp(k) ) then
807 npr = k
808 return
809 end if
810 end do
811 end if
812
813 end subroutine get_std_pr_level
814
815 subroutine get_std_ht_level(height, nht, stdh, nstdh)
816
817 implicit none
818
819 integer, intent(in ) :: nstdh
820 real, intent(in) :: stdh(nstdh)
821 integer, intent(out) :: nht
822 real, intent(in) :: height
823
824 real :: ht
825 integer :: k
826
827 ht = height*0.001 ! m to km
828 if ( ht <= stdh(1) ) then
829 nht = 1
830 return
831 else if ( ht > stdh(nstdh-1) ) then
832 nht = nstdh
833 return
834 else
835 do k = 2,nstdh - 1
836 if ( ht <= stdh(k) ) then
837 nht = k
838 return
839 end if
840 end do
841 end if
842
843 end subroutine get_std_ht_level
844
845 subroutine update_stats(stats_omb, stats_oma, omb, oma)
846
847 implicit none
848 type(stats_value), intent(inout) :: stats_omb, stats_oma
849 real, intent (in) :: omb, oma
850
851 real :: x1, x2
852
853 stats_omb%num = stats_omb%num + 1
854 stats_oma%num = stats_omb%num
855
856 x1 = 1.0/ stats_omb%num
857 x2 = (stats_omb%num-1)*x1
858
859 stats_omb%bias = x2*stats_omb%bias + omb * x1
860 stats_oma%bias = x2*stats_oma%bias + oma * x1
861
862 stats_omb%abias = x2*stats_omb%abias + abs(omb) * x1
863 stats_oma%abias = x2*stats_oma%abias + abs(oma) * x1
864
865 stats_omb%rmse = x2*stats_omb%rmse + omb*omb * x1
866 stats_oma%rmse = x2*stats_oma%rmse + oma*oma * x1
867
868 end subroutine update_stats
869
870 subroutine write_diag_single_level(out_dir,ounit,ldate,obs_type,omb,oma)
871
872 implicit none
873
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
880
881 character*512 :: filename
882 integer :: ounit1, ounit2
883 real :: sigt,bar
884
885
886 ounit1 = ounit
887 ounit2 = ounit + 1
888
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
897 else
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)
901 bar=bar*1000.0
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)
904 bar=bar*1000.0
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)
908 bar=bar/100.0
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)
911 bar=bar/100.0
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
913 else
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
918 endif
919 !
920 end if
921 close(ounit1)
922 close(ounit2)
923
924 end subroutine write_diag_single_level
925
926 subroutine write_diag_multi_level(out_dir,ounit,ldate,obs_type,omb,oma)
927
928 implicit none
929
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)
936
937 character*512 :: filename
938 integer :: k
939 real :: xnum(nstd)
940 integer :: num(nstd)
941 real, dimension(nstd) :: rmse, bias, abias,sigt,bar
942 integer :: ounit1, ounit2
943
944 ounit1 = ounit
945 ounit2 = ounit + 1
946
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')
951
952 do k = 1, nstd
953 num(k) = omb(k)%num
954 if (num(k) <= 1 ) then
955 xnum(k) = rmiss
956 rmse(k) = rmiss
957 bias(k) = rmiss
958 abias(k) = rmiss
959 bar(k) = rmiss
960 sigt(k) = 1.0
961 else
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.
968 else
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))
973 end if
974 xnum(k) = num(k)
975 end if
976 end do
977
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)
980
981 do k = 1, nstd
982 num(k) = oma(k)%num
983 if( num(k) <= 1 ) then
984 xnum(k) = rmiss
985 rmse(k) = rmiss
986 bias(k) = rmiss
987 abias(k) = rmiss
988 else
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.
995 else
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))
1000 end if
1001 xnum(k) = num(k)
1002 end if
1003 end do
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)
1006
1007 close(ounit1)
1008 close(ounit2)
1009
1010 end subroutine write_diag_multi_level
1011 subroutine write_diag_multi_level_h(out_dir,ounit,date,obs_type,omb,oma)
1012 implicit none
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)
1019 !
1020 character*512 :: filename
1021 integer :: k
1022 real :: xnum(nstdh)
1023 integer :: num(nstdh)
1024 real, dimension(nstdh) :: rmse, bias, abias, sigt, bar
1025
1026 integer :: ounit1, ounit2
1027 !
1028 ounit1 = ounit
1029 ounit2 = ounit + 1
1030 !
1031
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')
1036
1037 do k = 1, nstdh
1038 num(k) = omb(k)%num
1039 if( num(k) <= 1 ) then
1040 xnum(k) = rmiss
1041 rmse(k) = rmiss
1042 bias(k) = rmiss
1043 abias(k) = rmiss
1044 bar(k) = rmiss
1045 sigt(k) = 1.0
1046 else
1047 if( index(obs_type,'_q') > 0 ) then
1048
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.
1054 else
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))
1059 endif
1060 xnum(k) = num(k)
1061 endif
1062 enddo
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)
1065
1066 do k = 1, nstdh
1067 num(k) = oma(k)%num
1068 if( num(k) <= 1 ) then
1069 xnum(k) = rmiss
1070 rmse(k) = rmiss
1071 bias(k) = rmiss
1072 abias(k) = rmiss
1073 bar(k) = rmiss
1074 sigt(k) = 1.0
1075 else
1076 if( index(obs_type,'_q') > 0 ) then
1077
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.
1083 else
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))
1088 endif
1089 xnum(k) = num(k)
1090 endif
1091 enddo
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)
1094
1095 !
1096 close(ounit1)
1097 close(ounit2)
1098 !
1099 end subroutine write_diag_multi_level_h
1100 !
1101 subroutine sig_test(num, bias, rmse, sigt,bar)
1102 implicit none
1103 integer, intent(in) :: num
1104 real, intent(in) :: bias, rmse
1105 real, intent(out) :: sigt, bar
1106
1107 real :: t_val, sd, tmp
1108 !
1109 sigt=0.
1110 tmp = num/real(num-1)
1111 ! sd = sqrt ( tmp*rmse - bias*bias/real((n-1)) )
1112 sd = sqrt ( tmp*( rmse - bias*bias ) )
1113 do k=2,34
1114 if (real(num-1) < alpha(k,1)) exit
1115 end do
1116
1117 t_val = bias*sqrt( real(num) ) /sd
1118 bar = alpha(k-1,2) * sd /sqrt( real(num) )
1119
1120 if (abs(t_val) >= alpha(k-1,2)) sigt=1.
1121
1122 end subroutine sig_test
1123
1124
1125 subroutine get_fileinfo
1126 implicit none
1127 #include "netcdf.inc"
1128 integer :: iost(12)
1129 integer :: ncid
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'
1134 return
1135 else
1136 print*, 'Retrieving mapping info from wrf_file: ',trim(wrf_file)
1137 end if
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
1150 l_got_info = .true.
1151 end if
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
1155
1156 subroutine set_mapinfo
1157 implicit none
1158 integer :: map_proj_util
1159 real :: xref, yref
1160 real :: start_x, start_y
1161 xref = nx/2.0
1162 yref = ny/2.0
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
1171 end if
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
1176
1177 subroutine check_domain(lat, lon, inside)
1178 implicit none
1179 real, intent(in) :: lat, lon
1180 logical, intent(out) :: inside
1181 real :: xx, yy
1182 call da_llxy_wrf(map_info, lat, lon, xx, yy)
1183 inside = .false.
1184 if ( xx >= istart .and. xx <= iend .and. &
1185 yy >= jstart .and. yy <= jend ) then
1186 inside = .true.
1187 end if
1188 end subroutine check_domain
1189
1190 end program da_verif_obs
Mirror of the official repository at https://github.com/wrf-model/WRF