?  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_anal / da_verif_anal.f90
blob4a71f9a39ac64faad8a1fda7cb8ff6522ce23e8c
1 program da_verif_anal !---------------------------------------------------------------------------
2 ! History:
3 !
4 ! Abstract: Program to calculate statistics for various experiment
5 ! for verification againsr its own analysis
6 !
7 ! Author: Syed RH Rizvi NCAR/MMM 05/25/2006
8 !---------------------------------------------------------------------------
9
10 use da_verif_anal_control, only : control_main, control_times, control_vars, &
11 max_3d_variables, max_2d_variables,num_vert_levels,verification_file_string,&
12 missing,namelist_unit,time_series_unit,time_average_unit,&
13 ncl_info_unit, grads_ctl_unit, out_dat_unit, time_series_2d, profile_time_series_3d,&
14 time_average_2d, profile_time_average_3d, filename, stime, etime, &
15 hstart, hend, hdate, date, pdate, desc3d, desc2d, var_to_get, var_to_plot,&
16 length_var, length_plot, output_input_grid, use_lowest_heights, vert_args, &
17 nx, ny, nz, number_of_levels, io_status, debug1, debug2, verify_its_own_analysis, &
18 num_verifying_experiments, verify_forecast_hour, domain, control_exp_dir, verif_dirs, &
19 out_dirs,start_year, end_year, start_month, end_month, start_day, end_day, &
20 start_hour, end_hour,start_minutes, end_minutes, start_seconds, end_seconds,interval_hour, &
21 num3dvar, num2dvar, var3d, var2d, num_scores, score_names, vertical_type
22
23
24 use da_netcdf_interface, only : da_get_dims_cdf, da_get_gl_att_int_cdf, da_get_gl_att_real_cdf, &
25 da_get_var_3d_real_cdf, da_get_var_2d_real_cdf, da_get_var_2d_int_cdf
26
27 implicit none
28
29 character (len=512) :: control_file, verif_file
30 character (len=512) :: out_dir
31 character (len=512) :: namelist_file, grads_file
32
33 integer :: time_loop_count
34 integer :: time(6), ptime(6)
35 integer :: nx1, ny1, nz1
36 integer :: nx2, ny2, nz2
37 integer :: i,k
38 integer :: ivar, iexp, iscore
39 integer, allocatable,dimension(:) :: count_recs
40 integer :: irec, dat_unit
41 character (len=10) :: sdate
42 character (len=20) :: file_string, domain_string, out_hr
43 logical, allocatable,dimension(:) :: first_score
44
45 real, allocatable, dimension(:,:,:) :: data_out1, data_out2
46 real, allocatable, dimension(:,:,:) :: data_out1_z, data_out2_z
47
48 real, allocatable, dimension(:,:,:) :: sum3d, asum3d, sqr3d, diff, absdiff, sqdiff
49 real, allocatable, dimension(:,:) :: score_avg_prof
50 real, allocatable, dimension(:) :: avg_prof
51 integer :: num_grads_levels
52 real, dimension( 100) :: grads_levels
53 integer, allocatable, dimension(:) :: num_counter
54
55 data grads_levels /1000.0, 925.0, 850.0, 700.0, 500.0, 400.0, 300.0, 250.0, &
56 200.0, 150.0, 100.0, 70.0, 50.0, 30.0, 20.0, 10.0, &
57 84*9999.9/
58
59 num_grads_levels=16
60
61
62
63 !---------------------------------------------------------------------
64 verify_forecast_hour = 0
65
66 namelist_file = 'namelist.in'
67 grads_file = 'statistics'
68
69 !----------------------------------------------------------------------------
70 debug1 = .false.
71 debug2 = .false.
72 vertical_type = 'n'
73
74 !--3D need update
75 num3dvar=max_3d_variables
76 var3d(1)='U'
77 var3d(2)='V'
78 var3d(3)='W'
79 var3d(4)='TK'
80 var3d(5)='PH'
81 var3d(6)='RH'
82
83 desc3d(1)='U component of wind '
84 desc3d(2)='V component of wind '
85 desc3d(3)='W component of wind '
86 desc3d(4)='Temperature of air '
87 desc3d(5)='Geopotential Temperature'
88 desc3d(6)='Specific Humidity '
89
90 !--2D need update
91 num2dvar=max_2d_variables
92 var2d(1)='SLP'
93 var2d(2)='PSFC'
94 var2d(3)='U10M'
95 var2d(4)='V10M'
96 var2d(5)='T2M'
97 var2d(6)='Q2M'
98
99 desc2d(1)='Sea Level Pressure '
100 desc2d(2)='Surface Pressure '
101 desc2d(3)='10 meter Wind U compoment'
102 desc2d(4)='10 meter Wind V compoment'
103 desc2d(5)='2 meter Temperature '
104 desc2d(6)='2 meter Specific Humidity'
105
106 !--Score names
107 num_scores = 3
108 score_names(1) = 'BIAS'
109 score_names(2) = 'RMSE'
110 score_names(3) = 'ABIAS'
111
112
113 domain = 1
114 verification_file_string = 'wrfout'
115 out_hr ='_00'
116 !---------------------------------------------------------------------
117 ! Read namelist
118 !----------------------------------------------------------------------------
119 io_status = 0
120
121
122 open(unit = namelist_unit, file = trim(namelist_file), &
123 status = 'old' , access = 'sequential', &
124 form = 'formatted', action = 'read', &
125 iostat = io_status )
126
127
128 if(io_status /= 0) then
129 print *, 'Error to open namelist file: ', trim(namelist_file)
130 else
131 read(unit=namelist_unit, nml = control_main , iostat = io_status)
132
133 if(io_status /= 0) then
134 print *, 'Error to read control_main. Stopped.'
135 stop
136 endif
137
138
139
140
141 read(unit=namelist_unit, nml = control_times , iostat = io_status )
142
143 if(io_status /= 0) then
144 print *, 'Error to read control_times Stopped.'
145 stop
146 endif
147
148 read(unit=namelist_unit, nml = control_vars , iostat = io_status )
149
150 if(io_status /= 0) then
151 print *, 'Error to read control_vars Stopped.'
152 stop
153 endif
154
155 close(unit=namelist_unit)
156 endif
157 !----------------------------------------------------------------------------
158 !---------------------------------------------------------------------
159 if( num_grads_levels == 0) then ! output cartesian grid
160 write(6,*) ' generic cartesian output grid '
161 use_lowest_heights = .false.
162 output_input_grid = .true.
163 num_grads_levels = 1
164 else if ( num_grads_levels < 0) then
165 write(6,*) ' interp to z profile at lowest terrain column '
166 use_lowest_heights = .true. ! use z distribution from
167 output_input_grid = .true. ! lowest point of terrain
168 num_grads_levels = 1
169 else
170 output_input_grid = .false.
171 use_lowest_heights = .false.
172 endif
173 !
174 write(domain_string, fmt ='("_d",i2.2,"_")') domain
175 file_string = trim(verification_file_string)//trim(domain_string)
176 write(out_hr, fmt ='("_",i2.2)') verify_forecast_hour
177 allocate(count_recs(num_scores))
178 allocate(first_score(num_scores))
179 !---------------------------------------------------------------------
180
181 stime(1) = start_year
182 stime(2) = start_month
183 stime(3) = start_day
184 stime(4) = start_hour
185 stime(5) = start_minutes
186 stime(6) = start_seconds
187
188
189 call build_hdate(hstart, stime )
190
191 etime(1) = end_year
192 etime(2) = end_month
193 etime(3) = end_day
194 etime(4) = end_hour
195 etime(5) = end_minutes
196 etime(6) = end_seconds
197
198 call build_hdate(hend, etime )
199
200
201 if ( hend < hstart ) then
202 print*, '****************************************************************'
203 print*, 'End time is before the start time'
204 print*, ' Start time = ', hstart ,' & End time = ', hend
205 print*, '****************************************************************'
206 stop
207 endif
208
209 hdate = hstart
210 call build_hdate(sdate, stime )
211
212 filename = trim(file_string)//hdate
213
214 loop_verif : do iexp = 1, num_verifying_experiments
215 verif_file = trim(verif_dirs(iexp))//'/'//sdate//'/'//trim(filename)
216 out_dir = trim(out_dirs(iexp))//'/'
217 do iscore = 1, num_scores
218 grads_file = trim(out_dir)//trim(score_names(iscore))//trim(out_hr)
219 call create_grads_ctl (grads_file, verif_file, 1 , hdate, 1, &
220 var3d, num3dvar, var2d, num2dvar, desc3d, desc2d, &
221 output_input_grid, use_lowest_heights, &
222 grads_levels, num_grads_levels, number_of_levels, &
223 vert_args, vertical_type, debug1, debug2, &
224 ncl_info_unit, grads_ctl_unit, missing )
225 enddo
226 close(grads_ctl_unit)
227 close(ncl_info_unit)
228 enddo loop_verif
229 number_of_levels=num_grads_levels
230 !---------------------------------------------------------------------
231
232 !---------------------------------------------------------------------
233 loop_verif_exp : do iexp = 1, num_verifying_experiments
234
235 count_recs = 1
236 first_score = .true.
237 out_dir = trim(out_dirs(iexp))//'/'
238 !---------------------------------------------------------------------
239 !--For 3D variables
240 !---------------------------------------------------------------------
241 loop_3d : do ivar=1,num3dvar
242
243 ! Open profile units
244
245 ! profile_time_average_3d = trim(out_dir)//trim(var3d(ivar))//'_average'//trim(out_hr)
246 profile_time_average_3d = trim(out_dir)//trim(var3d(ivar))//'_profile'//trim(out_hr)
247 open(time_average_unit, file=trim(profile_time_average_3d), form='formatted', &
248 status='unknown')
249 profile_time_series_3d = trim(out_dir)//trim(var3d(ivar))//'_time_series'//trim(out_hr)
250 open(time_series_unit, file=profile_time_series_3d,form='formatted', &
251 status='unknown')
252
253 !----------------------------------------------------------------------------
254 var_to_get = var3d(ivar)
255 var_to_plot = var_to_get
256 call check_special_variable( var_to_get )
257
258 length_var = len_trim(var_to_get)
259 length_plot = len_trim(var_to_plot)
260 !
261 !----------------------------------------------------------------------------
262 time_loop_count = 0
263 hdate = hstart
264 time = stime
265
266 time_loop_3d : do
267
268 print*,' processing exp: ',iexp,' 3d var: ',trim(var_to_get),' for time: ',hdate
269 !----------------------------------------------------------------------------
270 call build_hdate(hdate, time)
271
272 if ( hdate > hend ) exit time_loop_3d
273
274 call build_hdate(date, time )
275 ptime = time
276 call advance_date(ptime,-verify_forecast_hour)
277 call build_hdate(pdate, ptime)
278
279
280 filename = trim(file_string)//hdate
281
282 if( verify_its_own_analysis) then
283 control_file = trim(verif_dirs(iexp))//'/'//date//'/'//trim(filename)
284 else
285 control_file = trim(control_exp_dir)//'/'//date//'/'//trim(filename)
286 endif
287
288 verif_file = trim(verif_dirs(iexp))//'/'//pdate//'/'//trim(filename)
289
290
291 ! Get the dimensions of the both files and check
292 call get_dimensions(control_file,nx1,ny1,nz1)
293 call get_dimensions(verif_file,nx2,ny2,nz2)
294
295 if ( nx1 /= nx2 .or. ny1 /= ny2 .or. nz1 /= nz2 ) then
296 print*, '********************************************************'
297 print*, 'Dimension mismatch between files of the experiments ....'
298 print*, '********************************************************'
299 stop
300 else
301 nx = nx1
302 ny = ny1
303 nz = nz1
304 if (time_loop_count == 0 ) then
305 allocate( sum3d(nx,ny,number_of_levels))
306 allocate( asum3d(nx,ny,number_of_levels))
307 allocate( sqr3d(nx,ny,number_of_levels))
308 sum3d = 0.0
309 asum3d = 0.0
310 sqr3d = 0.0
311 endif
312 endif
313 ! first, get control data
314 allocate ( data_out1_z (nx, ny, number_of_levels) )
315 call compute_data_3d( control_file, var_to_plot, length_plot, &
316 nx, ny, nz, number_of_levels, 1, vert_args, &
317 vertical_type, missing, data_out1_z, debug1 )
318 ! second, get verifying data
319 allocate ( data_out2_z (nx, ny, number_of_levels) )
320 call compute_data_3d( verif_file, var_to_plot, length_plot, &
321 nx, ny, nz, number_of_levels, 1, vert_args, &
322 vertical_type, missing, data_out2_z, debug2 )
323
324 allocate(diff(nx,ny,number_of_levels))
325 allocate(absdiff(nx,ny,number_of_levels))
326 allocate(sqdiff(nx,ny,number_of_levels))
327 call get_diffs(data_out1_z, data_out2_z, diff, absdiff, sqdiff, nx, ny, &
328 number_of_levels, missing)
329 deallocate(data_out1_z)
330 deallocate(data_out2_z)
331
332 allocate( avg_prof(number_of_levels) )
333 allocate( num_counter(number_of_levels) )
334 allocate( score_avg_prof(num_scores, number_of_levels) )
335 do iscore = 1, num_scores
336 if ( trim(score_names(iscore)) == 'BIAS' ) then
337 call domain_average( diff, avg_prof, num_counter, nx, ny, number_of_levels, missing,0)
338 elseif ( trim(score_names(iscore)) == 'RMSE' ) then
339 call domain_average( sqdiff, avg_prof, num_counter, nx, ny, number_of_levels, missing,1)
340 elseif ( trim(score_names(iscore)) == 'ABIAS' ) then
341 call domain_average( absdiff, avg_prof, num_counter, nx, ny, number_of_levels, missing,0)
342 endif
343 score_avg_prof(iscore,:) = avg_prof(:)
344 enddo
345 call write_profile(date, score_avg_prof, num_counter, number_of_levels, num_scores, &
346 time_series_unit)
347 deallocate( avg_prof )
348 deallocate( num_counter )
349 deallocate( score_avg_prof )
350
351 call get_sum(sum3d,diff,nx,ny,number_of_levels,missing)
352 call get_sum(asum3d,absdiff,nx,ny,number_of_levels,missing)
353 call get_sum(sqr3d,sqdiff,nx,ny,number_of_levels,missing)
354
355 deallocate(diff)
356 deallocate(absdiff)
357 deallocate(sqdiff)
358
359 time_loop_count = time_loop_count + 1
360
361 call advance_date(time,interval_hour)
362
363 enddo time_loop_3d ! time loop over
364
365
366
367 close(time_series_unit)
368 call time_average(sum3d,nx,ny,number_of_levels,time_loop_count,missing,0)
369 call time_average(asum3d,nx,ny,number_of_levels,time_loop_count,missing,0)
370 call time_average(sqr3d,nx,ny,number_of_levels,time_loop_count,missing,1)
371
372 allocate( avg_prof(number_of_levels) )
373 allocate( num_counter(number_of_levels) )
374 allocate( score_avg_prof(num_scores, number_of_levels) )
375
376 do iscore = 1, num_scores
377
378 grads_file = trim(out_dir)//trim(score_names(iscore))//out_hr
379 if ( trim(score_names(iscore)) == 'BIAS' ) then
380 dat_unit = out_dat_unit + (iscore -1)
381 irec = count_recs(iscore)
382 ! call write_out_data( grads_file, irec, sum3d, nx, ny, number_of_levels, &
383 ! first_score(iscore), dat_unit )
384 call domain_average( sum3d, avg_prof, num_counter,nx, ny, number_of_levels, missing,0)
385 if (first_score(iscore) ) first_score(iscore) = .false.
386 count_recs(iscore) = irec
387 elseif ( trim(score_names(iscore)) == 'RMSE' ) then
388 dat_unit = out_dat_unit + (iscore -1)
389 irec = count_recs(iscore)
390 ! call write_out_data( grads_file, irec, sqr3d, nx, ny, number_of_levels, &
391 ! first_score(iscore), dat_unit )
392 call domain_average( sqr3d, avg_prof, num_counter,nx, ny, number_of_levels, missing,1)
393 if (first_score(iscore) ) first_score(iscore) = .false.
394 count_recs(iscore) = irec
395 elseif ( trim(score_names(iscore)) == 'ABIAS' ) then
396 dat_unit = out_dat_unit + (iscore -1)
397 irec = count_recs(iscore)
398 ! call write_out_data( grads_file, irec, asum3d, nx, ny, number_of_levels, &
399 ! first_score(iscore), dat_unit )
400 call domain_average( asum3d, avg_prof, num_counter,nx, ny, number_of_levels, missing,0)
401 if (first_score(iscore) ) first_score(iscore) = .false.
402 count_recs(iscore) = irec
403 endif
404 score_avg_prof(iscore,:) = avg_prof(:)
405
406 enddo
407 deallocate(sum3d)
408 deallocate(asum3d)
409 deallocate(sqr3d)
410
411 call write_profile(sdate, score_avg_prof, num_counter, number_of_levels, num_scores, &
412 time_average_unit)
413 deallocate(score_avg_prof)
414 deallocate(avg_prof)
415 deallocate( num_counter )
416
417 enddo loop_3d
418 print*, ' successful completion of loop_3d '
419
420 !--------------------------------------------------------------------------------
421 !--Loop For 2D variables
422 !--------------------------------------------------------------------------------
423 loop_2d : do ivar = 1, num2dvar
424
425 ! Open profile units
426
427 time_average_2d = trim(out_dir)//trim(var2d(ivar))//'_average'//trim(out_hr)
428 open(time_average_unit, file=trim(time_average_2d), form='formatted', &
429 status='unknown')
430 time_series_2d = trim(out_dir)//trim(var2d(ivar))//'_time_series'//trim(out_hr)
431 open(time_series_unit, file=time_series_2d,form='formatted', &
432 status='unknown')
433
434 !----------------------------------------------------------------------------
435 var_to_get = var2d(ivar)
436 var_to_plot = var_to_get
437
438 call check_special_variable( var_to_get )
439
440 length_var = len_trim(var_to_get)
441 length_plot = len_trim(var_to_plot)
442 !
443 !----------------------------------------------------------------------------
444
445 time_loop_count = 0
446 hdate = hstart
447 time = stime
448
449 time_loop_2d : do
450
451 print*,' processing exp: ',iexp,' 2d var: ',trim(var_to_get),' for time: ',hdate
452 !----------------------------------------------------------------------------
453 call build_hdate(hdate, time )
454
455 if ( hdate > hend ) exit time_loop_2d
456
457 call build_hdate(date, time )
458 ptime = time
459 call advance_date(ptime,-verify_forecast_hour)
460 call build_hdate(pdate, ptime)
461
462 filename = trim(file_string)//hdate
463 if( verify_its_own_analysis) then
464 control_file = trim(verif_dirs(iexp))//'/'//date//'/'//trim(filename)
465 else
466 control_file = trim(control_exp_dir)//'/'//date//'/'//trim(filename)
467 endif
468
469 verif_file = trim(verif_dirs(iexp))//'/'//pdate//'/'//trim(filename)
470 !
471 ! Get the dimensions of the both files and check
472 call get_dimensions(control_file,nx1,ny1,nz1)
473 call get_dimensions(verif_file,nx2,ny2,nz2)
474
475 if ( nx1 /= nx2 .or. ny1 /= ny2 .or. nz1 /= nz2 ) then
476 print*, '********************************************************'
477 print*, 'Dimension mismatch between files of the experiments ....'
478 print*, '********************************************************'
479 stop
480 else
481 nx = nx1
482 ny = ny1
483 nz = nz1
484 if (time_loop_count == 0 ) then
485 allocate( sum3d(nx,ny,1))
486 allocate( asum3d(nx,ny,1))
487 allocate( sqr3d(nx,ny,1))
488 sum3d = 0.0
489 asum3d = 0.0
490 sqr3d = 0.0
491 endif
492 endif
493
494 allocate(data_out1(nx, ny, 1))
495 allocate(data_out2(nx, ny, 1))
496
497 call g_output_2d (control_file, 1, var_to_plot, length_plot, &
498 nx, ny, nz, data_out1, debug1)
499
500 call g_output_2d (verif_file, 1, var_to_plot, length_plot, &
501 nx, ny, nz, data_out2, debug2)
502
503 allocate(diff(nx,ny,1))
504 allocate(absdiff(nx,ny,1))
505 allocate(sqdiff(nx,ny,1))
506 call get_diffs(data_out1, data_out2, diff, absdiff, sqdiff, nx, ny, 1, missing)
507 deallocate(data_out1)
508 deallocate(data_out2)
509
510 allocate( avg_prof(1) )
511 allocate( num_counter(1) )
512 allocate( score_avg_prof(num_scores, 1) )
513 do iscore = 1, num_scores
514 if ( trim(score_names(iscore)) == 'BIAS' ) then
515 call domain_average( diff, avg_prof, num_counter, nx, ny, 1, missing,0)
516 elseif ( trim(score_names(iscore)) == 'RMSE' ) then
517 call domain_average( sqdiff, avg_prof, num_counter, nx, ny, 1, missing,1)
518 elseif ( trim(score_names(iscore)) == 'ABIAS' ) then
519 call domain_average( absdiff, avg_prof, num_counter, nx, ny, 1, missing,0)
520 endif
521 score_avg_prof(iscore,:) = avg_prof(:)
522 enddo
523 call write_profile(date, score_avg_prof, num_counter, 1, num_scores, &
524 time_series_unit)
525 deallocate( avg_prof )
526 deallocate( num_counter )
527 deallocate( score_avg_prof )
528
529 call get_sum(sum3d,diff,nx,ny,1,missing)
530 call get_sum(asum3d,absdiff,nx,ny,1,missing)
531 call get_sum(sqr3d,sqdiff,nx,ny,1,missing)
532 ! write(80+time_loop_count,*) sum3d
533 ! write(90+time_loop_count,*) asum3d
534
535 deallocate(diff)
536 deallocate(absdiff)
537 deallocate(sqdiff)
538 time_loop_count = time_loop_count + 1
539
540 call advance_date(time,interval_hour)
541
542
543 enddo time_loop_2d
544
545 close(time_series_unit)
546
547 !---------------------------------------------------------------------
548 ! calculate bias and RMSE
549 !---------------------------------------------------------------------
550 call time_average(sum3d,nx,ny,1,time_loop_count,missing,0)
551 call time_average(asum3d,nx,ny,1,time_loop_count,missing,0)
552 call time_average(sqr3d,nx,ny,1,time_loop_count,missing,1)
553
554 allocate( avg_prof(1) )
555 allocate( num_counter(1) )
556 allocate( score_avg_prof(num_scores, 1) )
557
558 !---------------------------------------------------------------------
559 ! Writting the results in grads file
560 !---------------------------------------------------------------------
561 do iscore = 1, num_scores
562 ! grads_file = trim(out_dir)//trim(score_names(iscore))//out_hr
563 if ( trim(score_names(iscore)) == 'BIAS' ) then
564 ! dat_unit = out_dat_unit + (iscore -1)
565 ! irec = count_recs(iscore)
566 ! call write_out_data( grads_file, irec, sum3d, nx, ny, 1, &
567 ! first_score(iscore), dat_unit )
568 ! count_recs(iscore) = irec
569 call domain_average( sum3d, avg_prof, num_counter,nx, ny, 1, missing,0)
570
571 elseif ( trim(score_names(iscore)) == 'RMSE' ) then
572 ! dat_unit = out_dat_unit + (iscore -1)
573 ! irec = count_recs(iscore)
574 ! call write_out_data( grads_file, irec, sqr3d, nx, ny, 1, &
575 ! first_score(iscore), dat_unit )
576 ! count_recs(iscore) = irec
577 call domain_average( sqr3d, avg_prof, num_counter,nx, ny, 1, missing,1)
578 else
579 ! dat_unit = out_dat_unit + (iscore -1)
580 ! irec = count_recs(iscore)
581 ! call write_out_data( grads_file, irec, asum3d, nx, ny, 1, &
582 ! first_score(iscore), dat_unit )
583 ! count_recs(iscore) = irec
584 call domain_average( asum3d, avg_prof, num_counter,nx, ny, 1, missing,0)
585 endif
586 score_avg_prof(iscore,:) = avg_prof(:)
587 enddo
588 !---------------------------------------------------------------------
589 deallocate(sum3d)
590 deallocate(asum3d)
591 deallocate(sqr3d)
592
593 call write_profile(sdate, score_avg_prof, num_counter, 1, num_scores, &
594 time_average_unit)
595
596 deallocate(score_avg_prof)
597 deallocate(avg_prof)
598 deallocate( num_counter )
599
600
601 enddo loop_2d
602 print*, ' successful completion of loop_2d '
603 !
604 ! Writting Latitude and Longitude at the end of the file
605 !
606
607 filename = trim(file_string)//hstart
608
609 if( verify_its_own_analysis) then
610 control_file = trim(verif_dirs(iexp))//'/'//sdate//'/'//trim(filename)
611 else
612 control_file = trim(control_exp_dir)//'/'//sdate//'/'//trim(filename)
613 endif
614
615 call get_dimensions(control_file,nx,ny,nz)
616 allocate(data_out1(nx, ny, 1))
617 allocate(data_out2(nx, ny, 1))
618
619 var_to_plot = 'XLAT'
620 length_plot = len_trim(var_to_plot)
621 call g_output_2d (control_file, 1, var_to_plot, length_plot, &
622 nx, ny, nz, data_out1, debug1)
623 var_to_plot = 'XLONG'
624 length_plot = len_trim(var_to_plot)
625 call g_output_2d (control_file, 1, var_to_plot, length_plot, &
626 nx, ny, nz, data_out2, debug1)
627 ! do iscore = 1, num_scores
628 ! grads_file = trim(out_dir)//trim(score_names(iscore))//out_hr
629 ! if ( trim(score_names(iscore)) == 'BIAS' ) then
630 ! dat_unit = out_dat_unit + (iscore -1)
631 ! irec = count_recs(iscore)
632 ! call write_out_data( grads_file, irec, data_out1, nx, ny, 1, &
633 ! first_score(iscore), dat_unit )
634 ! call write_out_data( grads_file, irec, data_out2, nx, ny, 1, &
635 ! first_score(iscore), dat_unit )
636 ! count_recs(iscore) = irec
637 ! elseif ( trim(score_names(iscore)) == 'RMSE' ) then
638 ! dat_unit = out_dat_unit + (iscore -1)
639 ! irec = count_recs(iscore)
640 ! call write_out_data( grads_file, irec, data_out1, nx, ny, 1, &
641 ! first_score(iscore), dat_unit )
642 ! call write_out_data( grads_file, irec, data_out2, nx, ny, 1, &
643 ! first_score(iscore), dat_unit )
644 ! count_recs(iscore) = irec
645 ! else
646 ! dat_unit = out_dat_unit + (iscore -1)
647 ! irec = count_recs(iscore)
648 ! call write_out_data( grads_file, irec, data_out1, nx, ny, 1, &
649 ! first_score(iscore), dat_unit )
650 ! call write_out_data( grads_file, irec, data_out2, nx, ny, 1, &
651 ! first_score(iscore), dat_unit )
652 ! count_recs(iscore) = irec
653 ! endif
654 ! enddo
655 deallocate(data_out1)
656 deallocate(data_out2)
657
658 close(time_average_unit)
659 do iscore = 1, num_scores
660 close(out_dat_unit + (iscore-1))
661 enddo
662 print*,' Finished Experiment : ', trim(verif_dirs(iexp))
663 enddo loop_verif_exp
664
665 !-----------------------------------------------------
666
667 contains
668 !-----------------------------------------------------
669 subroutine advance_date( time, delta )
670
671 implicit none
672
673 integer, intent(inout) :: time(6)
674 integer, intent(in) :: delta
675
676 integer :: ccyy, mm, dd, hh
677 integer, dimension(12) :: mmday
678 ! character(len=10) :: ccyymmddhh
679
680 !-----------------------------------------------------
681 mmday = (/31,28,31,30,31,30,31,31,30,31,30,31/)
682 mmday(2) = 28
683 !-----------------------------------------------------
684 ccyy = time(1)
685 mm = time(2)
686 dd = time(3)
687 hh = time(4)
688 !-----------------------------------------------------
689
690 hh = hh + delta
691
692 do while (hh < 0)
693 hh = hh + 24
694 call change_date ( ccyy, mm, dd, -1 )
695 end do
696
697 do while (hh > 23)
698 hh = hh - 24
699 call change_date ( ccyy, mm, dd, 1 )
700 end do
701
702 ! write(ccyymmddhh(1:10), fmt='(i4, 3i2.2)') ccyy, mm, dd, hh
703 time(1) = ccyy
704 time(2) = mm
705 time(3) = dd
706 time(4) = hh
707
708 end subroutine advance_date
709 !---------------------------------------------------------------------------
710 subroutine change_date ( ccyy, mm, dd, delta)
711 integer, intent(inout) :: ccyy, mm, dd
712 integer, intent(in) :: delta
713 !
714 integer, dimension(12) :: mmday
715 mmday = (/31,28,31,30,31,30,31,31,30,31,30,31/)
716
717 mmday(2) = 28
718
719 if (mod(ccyy,4) == 0) then
720 mmday(2) = 29
721
722 if ( mod(ccyy,100) == 0) then
723 mmday(2) = 28
724 endif
725
726 if(mod(ccyy,400) == 0) then
727 mmday(2) = 29
728 end if
729 endif
730
731 dd = dd + delta
732
733 if(dd == 0) then
734 mm = mm - 1
735
736 if(mm == 0) then
737 mm = 12
738 ccyy = ccyy - 1
739 endif
740
741 dd = mmday(mm)
742 elseif ( dd .gt. mmday(mm) ) then
743 dd = 1
744 mm = mm + 1
745 if(mm > 12 ) then
746 mm = 1
747 ccyy = ccyy + 1
748 end if
749 end if
750 end subroutine change_date
751
752 subroutine build_hdate(hdate, time)
753
754 ! PURPOSE:
755 ! From the Year, Month, Day, Hour, Minute, and Second values,
756 ! creates a 19-character string representing the date, in the
757 ! format: "YYYY-MM-DD hh:mm:ss"
758
759 ! INPUT:
760 integer, intent(in) :: time(6) ! all time specs in it
761 ! OUTPUT:
762 character*(*), intent(out) :: hdate ! 'YYYY-MM-DD hh:mm:ss'
763
764 ! LOCAL:
765 integer iyr ! year (e.g., 1997, 2001)
766 integer imo ! month (01 - 12)
767 integer idy ! day of the month (01 - 31)
768 integer ihr ! hour (00-23)
769 integer imi ! minute (00-59)
770 integer isc ! second (00-59)
771 !
772 ! integer i ! Loop counter.
773 integer hlen ! Length of hdate string
774
775 hlen = len(hdate)
776 iyr = time(1)
777 imo = time(2)
778 idy = time(3)
779 ihr = time(4)
780 imi = time(5)
781 isc = time(6)
782
783 if (hlen.eq.19) then
784 write(hdate,19) iyr, imo, idy, ihr, imi, isc
785 19 format(i4,'-',i2.2,'-',i2.2,'_',i2.2,':',i2.2,':',i2.2)
786
787 elseif (hlen.eq.16) then
788 write(hdate,16) iyr, imo, idy, ihr, imi
789 16 format(i4,'-',i2.2,'-',i2.2,'_',i2.2,':',i2.2)
790
791 elseif (hlen.eq.13) then
792 write(hdate,13) iyr, imo, idy, ihr
793 13 format(i4,'-',i2.2,'-',i2.2,'_',i2.2)
794
795 elseif (hlen.eq.10) then
796 write(hdate,10) iyr, imo, idy, ihr
797 10 format(i4,i2.2,i2.2,i2.2)
798 endif
799
800 return
801 end subroutine build_hdate
802
803
804 !---------------------------------------------------------------------------------
805
806 subroutine create_grads_ctl( grads_file, file_for_time, file_time_index, times, output_times, &
807 variables3d, number_of_3dvar, variables2d, number_of_2dvar, &
808 desc3d, desc2d, output_input_grid, use_lowest_heights, &
809 grads_levels, num_grads_levels, number_of_levels, &
810 vert_args, vertical_type, debug1, debug2, &
811 ncl_info_unit, grads_ctl_unit, missing)
812
813
814 implicit none
815 #include "netcdf.inc"
816
817 integer, intent(in) :: output_times
818 integer, intent(in) :: ncl_info_unit
819 integer, intent(in) :: grads_ctl_unit
820 character (len=512), intent(in) :: grads_file
821 character (len=*), intent(in) :: file_for_time
822 character (len=19), intent(in) :: times
823 integer, intent(in) :: file_time_index
824 integer, intent(inout) :: number_of_3dvar, number_of_2dvar
825 character (len=20), intent(inout), dimension(number_of_3dvar) :: variables3d
826 character (len=20), intent(inout), dimension(number_of_2dvar) :: variables2d
827 character (len=50), intent(inout), dimension(number_of_3dvar) :: desc3d
828 character (len=50), intent(inout), dimension(number_of_2dvar) :: desc2d
829 character (len=50) :: descdumm
830
831 logical, intent(in) :: output_input_grid, use_lowest_heights
832 integer, intent(in) :: num_grads_levels
833 integer, intent (inout) :: number_of_levels
834 real, dimension( num_grads_levels ), intent(in) :: grads_levels
835 real, intent(inout) :: vert_args(100)
836 character (len=1), intent(inout) :: vertical_type
837 logical, intent(in) :: debug1, debug2
838
839
840 real, allocatable, dimension(:,:,:) :: z, ph, phb
841 ! real, allocatable, dimension(:,:,:) :: p, pb
842 ! real, allocatable, dimension(:,:,:) :: data_out, data_out_z
843 character (len=30) :: var_to_get, var_to_plot
844 integer :: length_var, length_plot
845 integer :: ivar
846 integer :: num_vert
847 integer, dimension(2) :: loc_of_min_z
848 real , intent(in) :: missing
849
850
851 integer :: count_var
852 integer :: ncid, dimid, nf_status
853 integer :: rcode, trcode
854 real :: value_real
855 integer :: nx, ny, nz
856 integer :: nlgen
857 integer :: ndims, dims(4)
858 integer :: i, k
859 integer :: ilon
860
861
862 character (len=180) :: nclfile,ctlfile, datfile
863 character (len=35) :: tdef
864 integer :: timestamp, datestamp
865 ! integer :: file_recl
866
867
868 real, allocatable, dimension(:,:) :: xlat, xlon
869 real :: xlat_a(4), xlon_a(4)
870 ! real :: xlat_n_max, xlat_s_max
871 ! real :: xlon_w
872 ! real :: xlon_e
873 real :: abslatmin, abslatmax
874 real :: abslonmin, abslonmax
875 real :: truelat1, truelat2, temp
876 real :: cen_lat, cen_lon
877 ! real :: centeri, centerj
878 integer :: ipoints, jpoints
879 integer :: ncref, nrref
880 real :: dx, dy
881 real :: dxll
882 integer :: map_proj
883 logical :: debug
884
885
886 !==================================================================================
887 ! need to pull out some data to set up dimensions, etc.
888 !
889 call get_dimensions (file_for_time, nx, ny, nz )
890 nlgen = nz
891 !
892 !==================================================================================
893 ! open output files
894 ctlfile = trim(grads_file)//".ctl"
895 datfile = trim(grads_file)//".dat"
896
897 open (grads_ctl_unit, file=trim(ctlfile),status='unknown')
898 write (grads_ctl_unit, '("dset ^",a50)') datfile
899 #ifdef bytesw
900 write (grads_ctl_unit, '("options byteswapped")')
901 #endif
902 write (grads_ctl_unit, '("undef ",e7.1)') missing
903
904 !==================================================================================
905 ! How will the vertical coordinate look like
906
907 IF ( (.not. output_input_grid) .and. (.not. use_lowest_heights)) THEN
908 ! we have user supplied vertical levels - CAN WE DO IT?
909
910 nf_status = nf_open (file_for_time, nf_nowrite, ncid)
911
912 call handle_err('Error opening file',nf_status)
913 if ( vertical_type == 'p' ) then
914 rcode = nf_inq_varid ( ncid, "P", dimid )
915 trcode = rcode
916 rcode = nf_inq_varid ( ncid, "PB", dimid )
917 if ( nf_status == 0 ) rcode = trcode
918 else if ( vertical_type == 'z' ) then
919 rcode = nf_inq_varid ( ncid, "PH", dimid )
920 trcode = rcode
921 rcode = nf_inq_varid ( ncid, "PHB", dimid )
922 if ( nf_status == 0 ) rcode = trcode
923 endif
924 nf_status = nf_close (ncid)
925 call handle_err('Error closing file',nf_status)
926
927 if ( rcode == 0 ) then
928 ! we can do it
929 write(6,*) ' '
930 write(6,*) " Asking to interpolate to ",vertical_type," levels - we can do that"
931 write(6,*) ' '
932 number_of_levels = num_grads_levels
933 vert_args(1:number_of_levels) = grads_levels(1:number_of_levels)
934 else
935 ! no interp, just put out computational grid
936 write(6,*) ' '
937 write(6,*) ' FIELDS MISSING TO INTERPOLATE TO USER SPECIFIED VERTICAL LEVELS'
938 write(6,*) ' WILL OUTPUT ON MODEL GRID'
939 write(6,*) ' '
940 number_of_levels = nz
941 vertical_type = 'n'
942 endif
943
944 END IF
945
946 IF ( (output_input_grid) .and. (use_lowest_heights)) THEN
947 ! use lowest column for z heights of grids - CAN WE DO IT?
948
949 nf_status = nf_open (file_for_time, nf_nowrite, ncid)
950 call handle_err('Error opening file',nf_status)
951 rcode = nf_inq_varid ( ncid, "P", dimid )
952 trcode = rcode
953 rcode = nf_inq_varid ( ncid, "PB", dimid )
954 if ( nf_status == 0 ) rcode = trcode
955 nf_status = nf_close (ncid)
956 call handle_err('Error closing file',nf_status)
957
958 if ( rcode == 0 ) then
959 ! we can do it
960 write(6,*) ' '
961 write(6,*) " Asking to interpolate to lowerst h level - we can do that"
962 write(6,*) ' '
963 allocate( z(nx,ny,nz) )
964 allocate( ph(nx,ny,nz+1) )
965 allocate( phb(nx,ny,nz+1) )
966
967 ! get base and perturbation height at full eta levels:
968
969 call da_get_var_3d_real_cdf( file_for_time,'PH',ph, &
970 nx,ny,nz+1,1,debug2 )
971 call da_get_var_3d_real_cdf( file_for_time,'PHB',phb, &
972 nx,ny,nz+1,1,debug2 )
973
974 ! compute Z at half levels:
975
976 ph = (ph+phb)/9.81
977 z = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1))
978 z = z/1000. ! convert to kilometers
979
980 number_of_levels = nz
981 vertical_type = 'z'
982 loc_of_min_z = minloc(z(:,:,1))
983 vert_args(1:number_of_levels) = &
984 z(loc_of_min_z(1),loc_of_min_z(2),1:number_of_levels)
985 vert_args(1) = vert_args(1) + 0.002
986 vert_args(nz) = vert_args(nz) - 0.002
987
988 deallocate( z )
989 deallocate( ph )
990 deallocate( phb )
991 else
992 ! no interp, just put out computational grid
993 write(6,*) ' '
994 write(6,*) ' FIELDS MISSING TO INTERPOLATE TO HEIGHT LEVELS'
995 write(6,*) ' WILL OUTPUT ON MODEL GRID'
996 write(6,*) ' '
997 number_of_levels = nz
998 vertical_type = 'n'
999 endif
1000
1001 END IF
1002
1003
1004 IF ( output_input_grid .and. (.not. use_lowest_heights)) THEN
1005 ! no interp, just put out computational grid
1006
1007 write(6,*) " Will use model levels for output"
1008 number_of_levels = nz
1009
1010 ENDIF
1011
1012 !==================================================================================
1013
1014 if(debug1) then
1015 write(6,*) ' number of levels = ',number_of_levels
1016 do k=1, number_of_levels
1017 write(6,*) ' k, vert_args(k) ',k,vert_args(k)
1018 enddo
1019 end if
1020
1021 !==================================================================================
1022 ! work out times and time differences
1023
1024 tdef = ' 11 linear 00z01jan2000 1hr'
1025 call time_calc(times, timestamp, datestamp, debug2, tdef, 1 )
1026 write (tdef(9:10),'(i2)') output_times
1027
1028 !==================================================================================
1029 ! try to get map information
1030
1031 call da_get_gl_att_int_cdf( file_for_time, 'MAP_PROJ', map_proj, debug2 )
1032
1033
1034 if ( map_proj /= 0 ) then
1035 ! get more map parameters first
1036 call da_get_gl_att_real_cdf( file_for_time, 'DX', dx, debug2 )
1037 call da_get_gl_att_real_cdf( file_for_time, 'DY', dy, debug2 )
1038 call da_get_gl_att_real_cdf( file_for_time, 'CEN_LAT', cen_lat, debug2 )
1039 call da_get_gl_att_real_cdf( file_for_time, 'TRUELAT1', truelat1, debug2 )
1040 call da_get_gl_att_real_cdf( file_for_time, 'TRUELAT2', truelat2, debug2 )
1041
1042 nf_status = nf_open (file_for_time, nf_nowrite, ncid)
1043 call handle_err('Error opening file',nf_status)
1044 rcode = NF_GET_ATT_REAL(ncid, nf_global, 'STAND_LON', value_real )
1045 nf_status = nf_close (ncid)
1046 call handle_err('Error closing file',nf_status)
1047 if ( rcode == 0) then
1048 call da_get_gl_att_real_cdf( file_for_time, 'STAND_LON', cen_lon, debug2 )
1049 else
1050 write(6,*) ' ##### #####'
1051 write(6,*) ' ##### NOTE probably dealing with version 1 data #####'
1052 write(6,*) ' ##### Using CEN_LON in calculations #####'
1053 write(6,*) ' ##### Please check project of GrADS data #####'
1054 write(6,*) ' ##### #####'
1055 write(6,*) ' '
1056 call da_get_gl_att_real_cdf( file_for_time, 'CEN_LON', cen_lon, debug2 )
1057 endif
1058
1059 allocate( xlat(nx,ny) )
1060 allocate( xlon(nx,ny) )
1061 ! debug = .true.
1062 debug = debug2
1063 call da_get_var_2d_real_cdf( file_for_time, 'XLAT', xlat, nx,ny, 1, debug )
1064 call da_get_var_2d_real_cdf( file_for_time, 'XLONG',xlon, nx,ny, 1, debug )
1065 ! debug = .false.
1066
1067 end if
1068
1069 if (map_proj == 0 .OR. map_proj == 3) then
1070 ! NO or MERCATOR
1071
1072 ! check for dateline
1073 ilon = 0
1074 if ( abs(xlon(1,1) - xlon(nx,ny)) .GT. 180. ) ilon = 1
1075
1076 IF ( ilon == 1 ) THEN
1077 write(grads_ctl_unit,'("xdef ",i4," linear ",f9.4," ",f8.4)') &
1078 nx,xlon(1,1),(abs(xlon(1,1)-(360.+xlon(nx,ny)))/(nx-1))
1079 ELSE
1080 write(grads_ctl_unit,'("xdef ",i4," linear ",f9.4," ",f8.4)') &
1081 nx,xlon(1,1),(abs(xlon(1,1)-xlon(nx,ny))/(nx-1))
1082 ENDIF
1083 write(grads_ctl_unit,'("ydef ",i4," linear ",f9.4," ",f8.4)') &
1084 ny,xlat(1,1),(abs(xlat(1,1)-xlat(nx,ny))/(ny-1))
1085 if (vertical_type == 'n' ) then
1086 write (grads_ctl_unit, '("zdef ",i3, " linear 1 1")') number_of_levels
1087 else
1088 write(grads_ctl_unit,'("zdef ",i3, " levels ")') number_of_levels
1089 do k = 1,number_of_levels
1090 write(grads_ctl_unit,'(" ",f10.5)') vert_args(k)
1091 enddo
1092 endif
1093
1094
1095 else if (map_proj == 1) then
1096 ! LAMBERT-CONFORMAL
1097
1098
1099 ! make sure truelat1 is the larger number
1100 if (truelat1 < truelat2) then
1101 if (debug2) write (6,*) ' switching true lat values'
1102 temp = truelat1
1103 truelat1 = truelat2
1104 truelat2 = temp
1105 endif
1106
1107 xlat_a(1) = xlat(1,1)
1108 xlat_a(2) = xlat(1,ny)
1109 xlat_a(3) = xlat(nx,1)
1110 xlat_a(4) = xlat(nx,ny)
1111 xlon_a(1) = xlon(1,1)
1112 xlon_a(2) = xlon(1,ny)
1113 xlon_a(3) = xlon(nx,1)
1114 xlon_a(4) = xlon(nx,ny)
1115 abslatmin = 99999.
1116 abslatmax = -99999.
1117 abslonmin = 99999.
1118 abslonmax = -99999.
1119
1120 ! check for dateline
1121 ilon = 0
1122 if ( abs(xlon_a(1) - xlon_a(3)) .GT. 180. ) ilon = 1
1123 if ( abs(xlon_a(2) - xlon_a(4)) .GT. 180. ) ilon = 1
1124
1125 do i=1,4
1126 abslatmin=min(abslatmin,xlat_a(i))
1127 abslatmax=max(abslatmax,xlat_a(i))
1128 IF ( xlon_a(i) .lt. 0.0 .AND. ilon .eq. 1 ) THEN
1129 abslonmin=min(abslonmin,360.+xlon_a(i))
1130 abslonmax=max(abslonmax,360.+xlon_a(i))
1131 ELSE
1132 abslonmin=min(abslonmin,xlon_a(i))
1133 abslonmax=max(abslonmax,xlon_a(i))
1134 ENDIF
1135 enddo
1136 !
1137 ! xlat_s_max = -90.
1138 ! xlat_n_max = -90.
1139 !
1140 ! do i = 1, nx
1141 ! xlat_s_max = max (xlat_s_max,xlat(i,1))
1142 ! xlat_n_max = max (xlat_n_max,xlat(i,ny))
1143 ! enddo
1144
1145 ! xlon_w = xlon(1, ny)
1146 ! xlon_e = xlon(nx, ny)
1147 ! centeri = int((cen_lon-xlon_w)*((nx-1)/(xlon_e-xlon_w))+1)
1148 ! centerj = ((cen_lat-xlat_s_max)*((ny)/(xlat_n_max-xlat_s_max)))
1149
1150 dxll = (dx/1000.)/111./2.
1151 ipoints = int((abslatmax-abslatmin+2)/dxll)
1152 jpoints = int((abslonmax-abslonmin+2)/dxll)
1153
1154 write(grads_ctl_unit,'("pdef ",i3," ",i3," lcc ",f7.3," ",f8.3," ",&
1155 & f8.3," ",f8.3," ",f4.0," ",f4.0," ",f8.3," ",&
1156 & f7.0," ",f7.0)')&
1157 ! nx,ny,cen_lat,cen_lon,centeri,centerj,&
1158 nx,ny,xlat(1,1),xlon(1,1),1.0,1.0,&
1159 truelat1,truelat2,cen_lon,dx,dy
1160 write(grads_ctl_unit,'("xdef ",i4," linear ",f6.1," ",f12.8)') jpoints, &
1161 (abslonmin-1.),dxll
1162 write(grads_ctl_unit,'("ydef ",i4," linear ",f6.1," ",f12.8)') ipoints, &
1163 (abslatmin-1.),dxll
1164 if (vertical_type == 'n' ) then
1165 write (grads_ctl_unit, '("zdef ",i3, " linear 1 1")') number_of_levels
1166 else
1167 write(grads_ctl_unit,'("zdef ",i3, " levels ")') number_of_levels
1168 do k = 1,number_of_levels
1169 write(grads_ctl_unit,'(" ",f10.5)') vert_args(k)
1170 enddo
1171 endif
1172
1173
1174 elseif (map_proj == 2) then
1175 ! POLAR STEREO
1176
1177
1178 xlat_a(1) = xlat(1,1)
1179 xlat_a(2) = xlat(1,ny)
1180 xlat_a(3) = xlat(nx,1)
1181 xlat_a(4) = xlat(nx,ny)
1182 xlon_a(1) = xlon(1,1)
1183 xlon_a(2) = xlon(1,ny)
1184 xlon_a(3) = xlon(nx,1)
1185 xlon_a(4) = xlon(nx,ny)
1186 abslatmin = 99999.
1187 abslatmax = -99999.
1188 abslonmin = 99999.
1189 abslonmax = -99999.
1190
1191 do i=1,4
1192 abslatmin=min(abslatmin,xlat_a(i))
1193 abslonmin=min(abslonmin,xlon_a(i))
1194 abslatmax=max(abslatmax,xlat_a(i))
1195 abslonmax=max(abslonmax,xlon_a(i))
1196 enddo
1197
1198 dxll = (dx/1000.)/111./2.
1199 ipoints = int((abslatmax-abslatmin+2)/dxll) + 20
1200 jpoints = int((abslonmax-abslonmin+2)/dxll) + 20
1201
1202 ncref = nx/2
1203 nrref = ny/2
1204
1205 write(grads_ctl_unit,'("pdef ",i3," ",i3," ops ",f8.3," ",f8.3," ",f12.4," ", &
1206 & f12.4," ",i4.1," ",i4.1," ",f12.2," ",f12.2)') &
1207 nx,ny,xlat(ncref,nrref), xlon(ncref,nrref),dx*0.1,dy*0.1, &
1208 ncref,nrref,dx,dy
1209 write(grads_ctl_unit,'("xdef ",i4," linear ",f6.1," ",f12.8)') jpoints, &
1210 (abslonmin-1.),dxll
1211 write(grads_ctl_unit,'("ydef ",i4," linear ",f6.1," ",f12.8)') ipoints, &
1212 (abslatmin-1.),dxll
1213 if (vertical_type == 'n' ) then
1214 write (grads_ctl_unit, '("zdef ",i3, " linear 1 1")') number_of_levels
1215 else
1216 write(grads_ctl_unit,'("zdef ",i3, " levels ")') number_of_levels
1217 do k = 1,number_of_levels
1218 write(grads_ctl_unit,'(" ",f10.5)') vert_args(k)
1219 enddo
1220 endif
1221
1222
1223 endif ! END of map projections
1224
1225
1226 write(grads_ctl_unit, '("tdef",a35)') tdef
1227
1228 !==================================================================================
1229 !
1230 ! Write all required information for NCL plot
1231 !
1232 nclfile = trim(grads_file)//".info"
1233 nclfile = trim(grads_file)//".info"
1234 open (ncl_info_unit, file=trim(nclfile),status='unknown')
1235 nclfile = trim(grads_file)//".info"
1236
1237 open (ncl_info_unit, file=trim(nclfile),status='unknown')
1238
1239 write(ncl_info_unit, '(6(i4,1x))') number_of_3dvar, number_of_2dvar, map_proj, nx, ny, number_of_levels
1240 open (ncl_info_unit, file=trim(nclfile),status='unknown')
1241 write(ncl_info_unit, '(3(f7.2,1x))') truelat1,truelat2,cen_lon
1242 open (ncl_info_unit, file=trim(nclfile),status='unknown')
1243 ! do k = 1,number_of_levels
1244 write(ncl_info_unit,'(100f10.5)') (vert_args(k),k=1,number_of_levels)
1245 ! enddo
1246
1247 !
1248 ! Information writting for NCL is over
1249 !
1250 !==================================================================================
1251
1252 ! First check to see if we have all the variables
1253
1254 ! write(6,*) ' CHECK to make sure we have all the variables in the input file'
1255 call check_all_variables ( file_for_time, &
1256 variables3d, desc3d, number_of_3dvar, &
1257 variables2d, desc2d, number_of_2dvar, &
1258 debug1 )
1259 count_var = number_of_3dvar+number_of_2dvar+2
1260
1261 write(grads_ctl_unit, '("vars ",i3)') count_var
1262
1263 !==================================================================================
1264 do ivar = 1, number_of_3dvar
1265
1266 var_to_get = variables3d(ivar)
1267 var_to_plot = var_to_get
1268
1269 call check_special_variable( var_to_get )
1270 length_var = len_trim(var_to_get)
1271 length_plot = len_trim(var_to_plot)
1272
1273 call da_get_dims_cdf( file_for_time, var_to_get(1:length_var), &
1274 dims, ndims, debug2 )
1275
1276 if (dims(3) < number_of_levels ) then
1277 Write(6,*) 'No of vertical level is less here for: ', var_to_get
1278 else
1279 num_vert = number_of_levels
1280 endif
1281
1282 write(ncl_info_unit, '(a)' ) var_to_plot
1283 write(grads_ctl_unit,'(a15,i3," 0 ",a50)') var_to_plot, num_vert, desc3d(ivar)
1284
1285 enddo
1286
1287 do ivar = 1, number_of_2dvar
1288
1289 var_to_get = variables2d(ivar)
1290 var_to_plot = var_to_get
1291
1292 length_var = len_trim(var_to_get)
1293 length_plot = len_trim(var_to_plot)
1294
1295 write(ncl_info_unit, '(a)' ) var_to_plot
1296 write(grads_ctl_unit,'(a15," 0 0 ",a50)') var_to_plot, desc2d(ivar)
1297
1298 enddo
1299 var_to_plot = 'XLAT'
1300 descdumm = 'Latitude array'
1301 write(ncl_info_unit, '(a)' ) var_to_plot
1302 write(grads_ctl_unit,'(a15," 0 0 ",a50)') var_to_plot, descdumm
1303
1304 var_to_plot = 'XLONG'
1305 descdumm = 'Longitude array'
1306 write(ncl_info_unit, '(a)' ) var_to_plot
1307 write(grads_ctl_unit,'(a15," 0 0 ",a50)') var_to_plot, descdumm
1308
1309 write(grads_ctl_unit, '("endvars")' )
1310
1311 deallocate( xlat )
1312 deallocate( xlon )
1313
1314 end subroutine create_grads_ctl
1315 !==================================================================================
1316 #if 0
1317 subroutine write_out_data( grads_file, irec, data_in, nx, ny, number_of_levels, &
1318 first, grads_dat_unit )
1319
1320 implicit none
1321 #include "netcdf.inc"
1322
1323
1324 character (len=512), intent(in) :: grads_file
1325 integer, intent(inout) :: irec
1326 integer, intent(in) :: grads_dat_unit
1327 integer, intent(in) :: nx, ny, number_of_levels
1328 logical, intent(in) :: first
1329
1330 real, intent(in), dimension(:,:,:) :: data_in
1331
1332 integer :: file_recl, rec_length
1333 integer :: ii, jj, kk
1334
1335 character (len=512) :: datfile
1336
1337 ! open output files
1338 datfile = trim(grads_file)//".dat"
1339
1340 #ifdef RECL4
1341 file_recl = 4
1342 #endif
1343 #ifdef RECL1
1344 file_recl = 1
1345 #endif
1346
1347 if ( first ) then
1348
1349 if ( nx == 2 .and. ny /= 2 ) then
1350 rec_length = ny*file_recl
1351 elseif ( nx /= 2 .and. ny == 2 ) then
1352 rec_length = nx*file_recl
1353 elseif ( nx == 2 .and. ny == 2 ) then
1354 rec_length = file_recl
1355 else
1356 rec_length = nx*ny*file_recl
1357 endif
1358 open (grads_dat_unit, file=trim(datfile), form="unformatted",access="direct", &
1359 recl=rec_length,status='unknown')
1360
1361 endif
1362
1363 do kk=1,number_of_levels
1364 write(grads_dat_unit,rec=irec) ((data_in(ii,jj,kk),ii=1,nx),jj=1,ny)
1365 irec = irec + 1
1366 enddo
1367
1368 end subroutine write_out_data
1369 #endif
1370
1371 !----------------------------------------------------------------------------------
1372 subroutine write_profile(date, profile, counter, nlevel, nscore, out_unit)
1373
1374 integer, intent(in) :: nlevel, nscore, out_unit
1375 real, intent(in), dimension(:,:) :: profile
1376 integer, intent(in), dimension(:) :: counter
1377 character (len=10), intent(in) :: date
1378 write(out_unit,fmt='(a10,1x,100(i8,1x,3(f14.8,1x)))') date, &
1379 (counter(k), (profile(i,k),i=1,nscore),k=1,nlevel)
1380
1381 end subroutine write_profile
1382
1383 !---------------------------------------------------------------------------------
1384 subroutine time_calc( time, timestamp, datestamp, debug , tdef,it)
1385
1386 implicit none
1387
1388 character (len=19), intent(in) :: time
1389 character (len=35), intent(inout) :: tdef
1390 integer, intent(out) :: timestamp, datestamp
1391 logical, intent(in) :: debug
1392
1393 integer, intent(in) :: it
1394 integer :: hours, minutes, seconds, year, month, day,hour1,hourint
1395 integer :: mins1,minsint
1396
1397 save hourint
1398 save minsint
1399
1400 read(time(18:19),*) seconds
1401 read(time(15:16),*) minutes
1402 read(time(12:13),*) hours
1403 read(time(1:4),*) year
1404 read(time(6:7),*) month
1405 read(time(9:10),*) day
1406
1407 if(debug) write(6,*) ' day, month, year, hours, minutes, seconds '
1408 if(debug) write(6,*) day, month, year, hours, minutes, seconds
1409
1410 if ( it == 1) then
1411 write (tdef(19:20),'(i2)') hours
1412 if ( day < 10 ) then
1413 write (tdef(23:23),'(i1)') day
1414 else
1415 write (tdef(22:23),'(i2)') day
1416 endif
1417 write (tdef(27:30),'(i4)') year
1418 if (month == 1) write (tdef(24:26),'(a3)') 'jan'
1419 if (month == 2) write (tdef(24:26),'(a3)') 'feb'
1420 if (month == 3) write (tdef(24:26),'(a3)') 'mar'
1421 if (month == 4) write (tdef(24:26),'(a3)') 'apr'
1422 if (month == 5) write (tdef(24:26),'(a3)') 'may'
1423 if (month == 6) write (tdef(24:26),'(a3)') 'jun'
1424 if (month == 7) write (tdef(24:26),'(a3)') 'jul'
1425 if (month == 8) write (tdef(24:26),'(a3)') 'aug'
1426 if (month == 9) write (tdef(24:26),'(a3)') 'sep'
1427 if (month ==10) write (tdef(24:26),'(a3)') 'oct'
1428 if (month ==11) write (tdef(24:26),'(a3)') 'nov'
1429 if (month ==12) write (tdef(24:26),'(a3)') 'dec'
1430 hour1=hours
1431 mins1=minutes
1432 elseif ( it == 2) then
1433 hourint = abs(hours-hour1)
1434 minsint = abs(minutes-mins1)
1435 if (hourint == 0 ) then
1436 if (minsint == 0 ) minsint = 1
1437 if(debug) write(6,*) "interval is",minsint
1438 write (tdef(34:35),'(a2)') "mn"
1439 write (tdef(32:33),'(i2)') minsint
1440 if(debug) write(6,*) "TDEF is",tdef
1441 else
1442 if(debug) write(6,*) "Interval is",hourint
1443 write (tdef(32:33),'(i2)') hourint
1444 if(debug) write(6,*) "TDEF is",tdef
1445 endif
1446 endif
1447
1448 timestamp = seconds+100*minutes+10000*hours
1449
1450 if((year > 1800) .and. (year < 2000)) year = year-1900
1451 if((year >= 2000)) year = year-2000
1452
1453 if(month >= 2) day = day+31 ! add january
1454 if(month >= 3) day = day+28 ! add february
1455 if(month >= 4) day = day+31 ! add march
1456 if(month >= 5) day = day+30 ! add april
1457 if(month >= 6) day = day+31 ! add may
1458 if(month >= 7) day = day+30 ! add june
1459 if(month >= 8) day = day+31 ! add july
1460 if(month >= 9) day = day+31 ! add august
1461 if(month >= 10) day = day+30 ! add september
1462 if(month >= 11) day = day+31 ! add october
1463 if(month >= 12) day = day+30 ! add november
1464 if((month > 2) .and. (mod(year,4) == 0)) day = day+1 ! get leap year day
1465
1466 datestamp = (year)*1000 + day
1467 ! datestamp = (year+2100)*1000 + day
1468
1469 if(debug) then
1470 write(6,*) ' time, timestamp, datestamp ',time(1:19),timestamp,datestamp
1471 endif
1472
1473 end subroutine time_calc
1474
1475 !-------------------------------------------------------------------------
1476
1477 subroutine g_output_3d (file, file_time_index, var, length_var, &
1478 nx, ny, nz, data_out, debug)
1479 implicit none
1480
1481 character (len=512), intent(in) :: file
1482 integer, intent(in) :: file_time_index
1483 character (len=30), intent(in) :: var
1484 integer, intent(in) :: length_var
1485 integer , intent(in) :: nx, ny, nz
1486 real, intent(out), dimension(:,:,:) :: data_out
1487 logical, intent(in) :: debug
1488 real, allocatable, dimension(:,:,:) :: data_tmp, data_tmp2
1489 real, allocatable, dimension(:,:,:) :: u, v
1490 real, allocatable, dimension(:,:) :: xlat, xlon
1491 ! real, allocatable, dimension(:,:,:) :: z
1492 real, allocatable, dimension(:,:,:) :: ph, phb
1493 real, allocatable, dimension(:,:,:) :: p, pb
1494 real, allocatable, dimension(:,:,:) :: t, qv
1495 integer :: map_proj
1496 real :: cen_lon, truelat1, truelat2
1497
1498
1499 REAL , PARAMETER :: g = 9.81 ! acceleration due to gravity (m {s}^-2)
1500 REAL , PARAMETER :: r_d = 287.
1501 REAL , PARAMETER :: r_v = 461.6
1502 REAL , PARAMETER :: cp = 7.*r_d/2.
1503 REAL , PARAMETER :: cv = cp-r_d
1504 REAL , PARAMETER :: cliq = 4190.
1505 REAL , PARAMETER :: cice = 2106.
1506 REAL , PARAMETER :: psat = 610.78
1507 REAL , PARAMETER :: rcv = r_d/cv
1508 REAL , PARAMETER :: rcp = r_d/cp
1509 REAL , PARAMETER :: c2 = cp * rcv
1510 REAL , PARAMETER :: T0 = 273.16
1511
1512 REAL , PARAMETER :: p1000mb = 100000.
1513 REAL , PARAMETER :: cpovcv = cp/(cp-r_d)
1514 REAL , PARAMETER :: cvovcp = 1./cpovcv
1515 ! REAL :: pp
1516
1517 if(debug) then
1518 write(6,*) ' calculations for variable ',var
1519 end if
1520
1521 if(var == 'U' ) then
1522
1523 allocate ( data_tmp(nx+1,ny,nz) )
1524 call da_get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, &
1525 file_time_index, debug )
1526 data_out = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:))
1527
1528 deallocate ( data_tmp )
1529
1530 else if(var == 'V' ) then
1531
1532 allocate ( data_tmp(nx,ny+1,nz) )
1533 call da_get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, &
1534 file_time_index, debug )
1535 data_out = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:))
1536 deallocate ( data_tmp )
1537
1538 else if(var == 'UMET' ) then
1539
1540 call da_get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug )
1541
1542 IF ( map_proj == 1 .OR. map_proj == 2 ) THEN
1543
1544 allocate ( u(nx,ny,nz) )
1545 allocate ( v(nx,ny,nz) )
1546 allocate ( xlat(nx,ny) )
1547 allocate ( xlon(nx,ny) )
1548
1549 allocate ( data_tmp(nx+1,ny,nz) )
1550 call da_get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, &
1551 file_time_index, debug )
1552 u = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:))
1553 deallocate ( data_tmp )
1554
1555 allocate ( data_tmp(nx,ny+1,nz) )
1556 call da_get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, &
1557 file_time_index, debug )
1558 v = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:))
1559 deallocate ( data_tmp )
1560
1561 call da_get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug )
1562 call da_get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug )
1563 call da_get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug )
1564 call da_get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug )
1565 call da_get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug )
1566
1567 call rotate_wind (u,v,nx,ny,nz,var, &
1568 map_proj,cen_lon,xlat,xlon, &
1569 truelat1,truelat2,data_out)
1570
1571 deallocate ( xlat )
1572 deallocate ( xlon )
1573 deallocate ( u )
1574 deallocate ( v )
1575
1576 ELSE
1577
1578 allocate ( data_tmp(nx+1,ny,nz) )
1579 call da_get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, &
1580 file_time_index, debug )
1581 data_out = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:))
1582 deallocate ( data_tmp )
1583
1584 ENDIF
1585
1586 else if(var == 'VMET' ) then
1587
1588 call da_get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug )
1589
1590 IF ( map_proj == 1 .OR. map_proj == 2 ) THEN
1591
1592 allocate ( u(nx,ny,nz) )
1593 allocate ( v(nx,ny,nz) )
1594 allocate ( xlat(nx,ny) )
1595 allocate ( xlon(nx,ny) )
1596
1597 allocate ( data_tmp(nx+1,ny,nz) )
1598 call da_get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, &
1599 file_time_index, debug )
1600 u = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:))
1601 deallocate ( data_tmp )
1602
1603 allocate ( data_tmp(nx,ny+1,nz) )
1604 call da_get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, &
1605 file_time_index, debug )
1606 v = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:))
1607 deallocate ( data_tmp )
1608
1609 call da_get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug )
1610 call da_get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug )
1611 call da_get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug )
1612 call da_get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug )
1613 call da_get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug )
1614
1615 call rotate_wind (u,v,nx,ny,nz,var, &
1616 map_proj,cen_lon,xlat,xlon, &
1617 truelat1,truelat2,data_out)
1618
1619 deallocate ( xlat )
1620 deallocate ( xlon )
1621 deallocate ( u )
1622 deallocate ( v )
1623
1624 ELSE
1625
1626 allocate ( data_tmp(nx,ny+1,nz) )
1627 call da_get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, &
1628 file_time_index, debug )
1629 data_out = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:))
1630 deallocate ( data_tmp )
1631
1632 ENDIF
1633
1634 else if(var == 'W' ) then
1635
1636 allocate ( data_tmp(nx,ny,nz+1) )
1637 call da_get_var_3d_real_cdf( file,"W", data_tmp, nx, ny, nz+1, &
1638 file_time_index, debug )
1639 data_out = 0.5*(data_tmp(:,:,1:nz)+data_tmp(:,:,2:nz+1))
1640 deallocate ( data_tmp )
1641
1642 else if(var == 'P' ) then
1643
1644 allocate ( p(nx,ny,nz) )
1645 allocate ( pb(nx,ny,nz) )
1646
1647 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
1648 file_time_index, debug )
1649 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
1650 file_time_index, debug )
1651 data_out = (p+pb)*.01
1652
1653 deallocate ( p )
1654 deallocate ( pb )
1655
1656 else if(var == 'Z' ) then
1657
1658 allocate ( ph(nx,ny,nz+1) )
1659 allocate ( phb(nx,ny,nz+1) )
1660
1661 call da_get_var_3d_real_cdf( file,"PH", ph, nx, ny, nz+1, &
1662 file_time_index, debug )
1663 call da_get_var_3d_real_cdf( file,"PHB", phb, nx, ny, nz+1, &
1664 file_time_index, debug )
1665 ph = (ph+phb)/9.81
1666 data_out = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1))
1667
1668 deallocate ( ph )
1669 deallocate ( phb )
1670
1671 else if(var == 'THETA' ) then
1672
1673 call da_get_var_3d_real_cdf( file,"T", data_out, nx, ny, nz, &
1674 file_time_index, debug )
1675 data_out = data_out + 300.
1676
1677 else if(var == 'TK' ) then
1678
1679 allocate ( p(nx,ny,nz) )
1680 allocate ( pb(nx,ny,nz) )
1681 allocate ( data_tmp(nx,ny,nz) )
1682
1683 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
1684 file_time_index, debug )
1685 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
1686 file_time_index, debug )
1687 p = p+pb
1688
1689 call da_get_var_3d_real_cdf( file,"T", data_tmp, nx, ny, nz, &
1690 file_time_index, debug )
1691 data_out = (data_tmp+300.)*(p/p1000mb)**rcp
1692
1693 deallocate ( p )
1694 deallocate ( pb )
1695 deallocate ( data_tmp )
1696
1697 else if(var == 'TC' ) then
1698
1699 allocate ( p(nx,ny,nz) )
1700 allocate ( pb(nx,ny,nz) )
1701 allocate ( data_tmp(nx,ny,nz) )
1702
1703 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
1704 file_time_index, debug )
1705 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
1706 file_time_index, debug )
1707 p = p+pb
1708
1709 call da_get_var_3d_real_cdf( file,"T", data_tmp, nx, ny, nz, &
1710 file_time_index, debug )
1711 data_out = (data_tmp+300.)*(p/p1000mb)**rcp -T0
1712
1713 deallocate ( p )
1714 deallocate ( pb )
1715 deallocate ( data_tmp )
1716
1717 else if(var == 'TD' ) then
1718
1719 allocate ( p(nx,ny,nz) )
1720 allocate ( pb(nx,ny,nz) )
1721 allocate ( qv(nx,ny,nz) )
1722 allocate ( data_tmp(nx,ny,nz) )
1723
1724 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
1725 file_time_index, debug )
1726 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
1727 file_time_index, debug )
1728 p = p+pb
1729
1730 call da_get_var_3d_real_cdf( file,"QVAPOR", qv, nx, ny, nz, &
1731 file_time_index, debug )
1732
1733 data_tmp = qv*(p/100.)/(0.622+qv)
1734 data_tmp = AMAX1(data_tmp,0.001)
1735 data_out = (243.5*log(data_tmp)-440.8)/(19.48-log(data_tmp))
1736
1737 deallocate ( p )
1738 deallocate ( pb )
1739 deallocate ( qv )
1740 deallocate ( data_tmp )
1741
1742 else if(var == 'RH' ) then
1743
1744 allocate ( p(nx,ny,nz) )
1745 allocate ( pb(nx,ny,nz) )
1746 allocate ( qv(nx,ny,nz) )
1747 allocate ( t(nx,ny,nz) )
1748 allocate ( data_tmp(nx,ny,nz) )
1749 allocate ( data_tmp2(nx,ny,nz) )
1750
1751 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
1752 file_time_index, debug )
1753 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
1754 file_time_index, debug )
1755 p = p+pb
1756
1757 call da_get_var_3d_real_cdf( file,"T", t, nx, ny, nz, &
1758 file_time_index, debug )
1759 call da_get_var_3d_real_cdf( file,"QVAPOR", qv, nx, ny, nz, &
1760 file_time_index, debug )
1761
1762 t = (t+300.)*(p/p1000mb)**rcp
1763 data_tmp2 = 10.*0.6112*exp(17.67*(t-T0)/(t-29.65))
1764 data_tmp = 0.622*data_tmp2/(0.01 * p - (1.-0.622)*data_tmp2)
1765 data_out = 100.*AMAX1(AMIN1(qv/data_tmp,1.0),0.0)
1766
1767
1768 deallocate ( p )
1769 deallocate ( pb )
1770 deallocate ( qv )
1771 deallocate ( t )
1772 deallocate ( data_tmp )
1773 deallocate ( data_tmp2 )
1774
1775 else
1776 call da_get_var_3d_real_cdf( file,var(1:length_var), &
1777 data_out, nx,ny,nz, &
1778 file_time_index, debug )
1779 endif
1780
1781
1782 end subroutine g_output_3d
1783
1784 !-------------------------------------------------------------------------
1785
1786 subroutine g_output_2d (file, file_time_index, var, length_var, &
1787 nx, ny, nz, data_out, debug)
1788 implicit none
1789
1790 character (len=512), intent(in) :: file
1791 integer, intent(in) :: file_time_index
1792 character (len=30), intent(in) :: var
1793 integer, intent(in) :: length_var
1794 integer, intent(in) :: nx, ny, nz
1795 real, intent(out), dimension(:,:,:) :: data_out
1796 logical, intent(in) :: debug
1797 integer, allocatable, dimension(:,:,:) :: data_int
1798 real, allocatable, dimension(:,:,:) :: u10, v10
1799 real, allocatable, dimension(:,:) :: psfc,t2m,q2m
1800 real, allocatable, dimension(:,:) :: xlat, xlon
1801 real, allocatable, dimension(:,:,:) :: z,ph,phb
1802 real, allocatable, dimension(:,:,:) :: p,pb
1803 real, allocatable, dimension(:,:,:) :: ts,qv
1804 integer :: map_proj
1805 real :: cen_lon, truelat1, truelat2
1806
1807 if(debug) then
1808 write(6,*) ' calculations for variable ',var
1809 end if
1810
1811 if(var == 'SLP') then
1812
1813 allocate ( z(nx,ny,nz) )
1814 allocate ( ph(nx,ny,nz+1) )
1815 allocate ( phb(nx,ny,nz+1) )
1816 allocate ( p(nx,ny,nz) )
1817 allocate ( pb(nx,ny,nz) )
1818 allocate ( ts(nx,ny,nz) )
1819 allocate ( qv(nx,ny,nz) )
1820
1821 call da_get_var_3d_real_cdf( file,"PH", ph, nx, ny,nz+1, &
1822 file_time_index, debug )
1823 call da_get_var_3d_real_cdf( file,"PHB", phb, nx, ny,nz+1, &
1824 file_time_index, debug )
1825 ph = (ph+phb)/9.81
1826 z = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1))
1827
1828 call da_get_var_3d_real_cdf( file,"P", p, nx, ny,nz, &
1829 file_time_index, debug )
1830 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny,nz, &
1831 file_time_index, debug )
1832 p = p+pb
1833
1834 call da_get_var_3d_real_cdf( file,"T", ts, nx, ny,nz, &
1835 file_time_index, debug )
1836 call da_get_var_3d_real_cdf( file,"QVAPOR", qv, nx, ny,nz, &
1837 file_time_index, debug )
1838
1839 call compute_seaprs (nx, ny, nz, z, ts, p, qv, data_out, debug)
1840
1841
1842 deallocate ( z )
1843 deallocate ( ph )
1844 deallocate ( phb )
1845 deallocate ( p )
1846 deallocate ( pb )
1847 deallocate ( ts )
1848 deallocate ( qv )
1849
1850 else if(var == 'PSFC' ) then
1851 allocate ( psfc(nx,ny) )
1852 call da_get_var_2d_real_cdf( file,"PSFC", psfc, nx, ny, &
1853 file_time_index, debug )
1854 data_out(:,:,1) = psfc(:,:)
1855 deallocate ( psfc )
1856
1857 else if(var == 'T2M' ) then
1858 allocate ( t2m(nx,ny) )
1859 call da_get_var_2d_real_cdf( file,"T2", t2m, nx, ny, &
1860 file_time_index, debug )
1861 data_out(:,:,1) = t2m(:,:)
1862 deallocate ( t2m )
1863
1864 else if(var == 'Q2M' ) then
1865 allocate ( q2m(nx,ny) )
1866 call da_get_var_2d_real_cdf( file,"Q2", q2m, nx, ny, &
1867 file_time_index, debug )
1868 data_out(:,:,1) = q2m(:,:)
1869 deallocate ( q2m )
1870
1871 else if(var == 'U10M' ) then
1872
1873 call da_get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug )
1874
1875 IF ( map_proj == 1 .OR. map_proj == 2 ) THEN
1876
1877 allocate ( u10(nx,ny,1) )
1878 allocate ( v10(nx,ny,1) )
1879 allocate ( xlat(nx, ny) )
1880 allocate ( xlon(nx, ny) )
1881 call da_get_var_2d_real_cdf( file,"U10", u10, nx, ny, &
1882 file_time_index, debug )
1883 call da_get_var_2d_real_cdf( file,"V10", v10, nx, ny, &
1884 file_time_index, debug )
1885
1886 call da_get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug )
1887 call da_get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug )
1888 call da_get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug )
1889 call da_get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug )
1890 call da_get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug )
1891
1892 call rotate_wind (u10,v10,nx,ny,1,var, &
1893 map_proj,cen_lon,xlat,xlon, &
1894 truelat1,truelat2,data_out)
1895
1896 deallocate ( xlat )
1897 deallocate ( xlon )
1898 deallocate ( u10 )
1899 deallocate ( v10 )
1900
1901 ELSE
1902
1903 call da_get_var_2d_real_cdf( file,"U10", data_out, nx, ny, &
1904 file_time_index, debug )
1905
1906 ENDIF
1907
1908 else if(var == 'V10M' ) then
1909
1910 call da_get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug )
1911
1912 IF ( map_proj == 1 .OR. map_proj == 2 ) THEN
1913
1914 allocate ( u10(nx,ny,1) )
1915 allocate ( v10(nx,ny,1) )
1916 allocate ( xlat(nx, ny) )
1917 allocate ( xlon(nx, ny) )
1918 call da_get_var_2d_real_cdf( file,"U10", u10, nx, ny, &
1919 file_time_index, debug )
1920 call da_get_var_2d_real_cdf( file,"V10", v10, nx, ny, &
1921 file_time_index, debug )
1922
1923 call da_get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug )
1924 call da_get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug )
1925 call da_get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug )
1926 call da_get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug )
1927 call da_get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug )
1928
1929 call rotate_wind (u10,v10,nx,ny,1,var, &
1930 map_proj,cen_lon,xlat,xlon, &
1931 truelat1,truelat2,data_out)
1932
1933 deallocate ( xlat )
1934 deallocate ( xlon )
1935 deallocate ( u10 )
1936 deallocate ( v10 )
1937
1938 ELSE
1939
1940 call da_get_var_2d_real_cdf( file,"V10", data_out, nx, ny, &
1941 file_time_index, debug )
1942
1943 ENDIF
1944
1945 else if(var == 'XLONG' ) then
1946 call da_get_var_2d_real_cdf( file,var(1:length_var), &
1947 data_out, nx,ny, &
1948 file_time_index, debug )
1949 WHERE ( data_out < 0.0 )
1950 data_out = data_out + 360.0
1951 ENDWHERE
1952
1953 else if(var == 'IVGTYP' .or. var == 'ISLTYP') then
1954
1955 allocate (data_int(nx,ny,1))
1956 call da_get_var_2d_int_cdf( file,var(1:length_var), &
1957 data_int, nx,ny, &
1958 file_time_index, debug )
1959 data_out = data_int
1960 deallocate (data_int)
1961
1962 else
1963 call da_get_var_2d_real_cdf( file,var(1:length_var), &
1964 data_out, nx,ny, &
1965 file_time_index, debug )
1966 endif
1967
1968
1969 end subroutine g_output_2d
1970
1971 !------------------------------------------------------------------
1972
1973 subroutine check_special_variable( var_to_get )
1974
1975 implicit none
1976 character (len=20), intent(inout) :: var_to_get
1977
1978 if(var_to_get(1:6) == 'THETA ' .or. var_to_get(1:6) == 'TC ' &
1979 .or. var_to_get(1:6) == 'TK ') then
1980 var_to_get(1:6) = 'T '
1981 else if(var_to_get(1:6) == 'TD ' .or. var_to_get(1:6) == 'RH ' ) then
1982 var_to_get(1:6) = 'QVAPOR'
1983 else if(var_to_get(1:2) == 'Z ') then
1984 var_to_get(1:6) = 'PH '
1985 else if(var_to_get(1:4) == 'UMET') then
1986 var_to_get(1:6) = 'U '
1987 else if(var_to_get(1:4) == 'VMET') then
1988 var_to_get(1:6) = 'V '
1989 end if
1990
1991 end subroutine check_special_variable
1992
1993 !--------------------------------------------------------
1994
1995 subroutine interp_to_z( data_in , nx_in , ny_in , nz_in , &
1996 data_out, nx_out, ny_out, nz_out, &
1997 z_in, z_out, missing_value, &
1998 vertical_type, debug )
1999 implicit none
2000 integer, intent(in) :: nx_in , ny_in , nz_in
2001 integer, intent(in) :: nx_out, ny_out, nz_out
2002 real, intent(in) :: missing_value
2003 real, dimension(nx_in , ny_in , nz_in ), intent(in ) :: data_in, z_in
2004 real, dimension(nx_out, ny_out, nz_out), intent(out) :: data_out
2005 real, dimension(nz_out), intent(in) :: z_out
2006 logical, intent(in) :: debug
2007 character (len=1) , intent(in) :: vertical_type
2008
2009 real, dimension(nz_in) :: data_in_z, zz_in
2010 real, dimension(nz_out) :: data_out_z
2011
2012 integer :: i,j,k
2013
2014 do i=1,nx_in
2015 do j=1,ny_in
2016
2017 do k=1,nz_in
2018 data_in_z(k) = data_in(i,j,k)
2019 zz_in(k) = z_in(i,j,k)
2020 enddo
2021
2022 !Hui do k=1,nz_out
2023 !Hui data_out_z(k) = data_out(i,j,k)
2024 !Hui enddo
2025
2026 call interp_1d( data_in_z, zz_in, nz_in, &
2027 data_out_z, z_out, nz_out, &
2028 vertical_type, missing_value )
2029
2030 do k=1,nz_out
2031 data_out(i,j,k) = data_out_z(k)
2032 enddo
2033
2034
2035 enddo
2036 enddo
2037
2038 end subroutine interp_to_z
2039
2040 !----------------------------------------------
2041
2042 subroutine interp_1d( a, xa, na, &
2043 b, xb, nb, vertical_type, missing_value )
2044 implicit none
2045 integer, intent(in) :: na, nb
2046 real, intent(in), dimension(na) :: a, xa
2047 real, intent(in), dimension(nb) :: xb
2048 real, intent(out), dimension(nb) :: b
2049 real, intent(in) :: missing_value
2050
2051 integer :: n_in, n_out
2052 logical :: interp
2053 real :: w1, w2
2054 character (len=1) ,intent(in) :: vertical_type
2055
2056
2057 if ( vertical_type == 'p' ) then
2058
2059 do n_out = 1, nb
2060
2061 b(n_out) = missing_value
2062 interp = .false.
2063 n_in = 1
2064
2065 do while ( (.not.interp) .and. (n_in < na) )
2066
2067 if( (xa(n_in) >= xb(n_out)) .and. &
2068 (xa(n_in+1) <= xb(n_out)) ) then
2069 interp = .true.
2070 w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
2071 w2 = 1. - w1
2072 b(n_out) = w1*a(n_in) + w2*a(n_in+1)
2073 end if
2074 n_in = n_in +1
2075
2076 enddo
2077
2078 enddo
2079
2080 else
2081
2082 do n_out = 1, nb
2083
2084 b(n_out) = missing_value
2085 interp = .false.
2086 n_in = 1
2087
2088 do while ( (.not.interp) .and. (n_in < na) )
2089
2090 if( (xa(n_in) <= xb(n_out)) .and. &
2091 (xa(n_in+1) >= xb(n_out)) ) then
2092 interp = .true.
2093 w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
2094 w2 = 1. - w1
2095 b(n_out) = w1*a(n_in) + w2*a(n_in+1)
2096 end if
2097 n_in = n_in +1
2098
2099 enddo
2100
2101 enddo
2102
2103 endif
2104
2105 end subroutine interp_1d
2106
2107 !-------------------------------------------------------------------------
2108 !
2109 ! This routines has been taken "as is" from wrf_user_fortran_util_0.f
2110 !
2111 ! This routine assumes
2112 ! index order is (i,j,k)
2113 ! wrf staggering
2114 ! units: pressure (Pa), temperature(K), height (m), mixing ratio (kg kg{-1})
2115 ! availability of 3d p, t, and qv; 2d terrain; 1d half-level zeta string
2116 ! output units of SLP are Pa, but you should divide that by 100 for the
2117 ! weather weenies.
2118 ! virtual effects are included
2119 !
2120 ! Dave
2121
2122 subroutine compute_seaprs ( nx , ny , nz , &
2123 z, t , p , q , &
2124 sea_level_pressure,debug)
2125 ! & t_sea_level, t_surf, level )
2126 IMPLICIT NONE
2127 ! Estimate sea level pressure.
2128 INTEGER, intent(in) :: nx , ny , nz
2129 REAL, intent(in) :: z(nx,ny,nz)
2130 REAL, intent(in) :: p(nx,ny,nz) , q(nx,ny,nz)
2131 REAL, intent(inout) :: t(nx,ny,nz)
2132 ! The output is the 2d sea level pressure.
2133 REAL, intent(out) :: sea_level_pressure(nx,ny)
2134 INTEGER level(nx,ny)
2135 REAL t_surf(nx,ny) , t_sea_level(nx,ny)
2136 LOGICAL, intent(in) :: debug
2137
2138 ! Some required physical constants:
2139
2140 REAL R, G, GAMMA
2141 PARAMETER (R=287.04, G=9.81, GAMMA=0.0065)
2142
2143 ! Specific constants for assumptions made in this routine:
2144
2145 REAL TC, PCONST
2146 PARAMETER (TC=273.16+17.5, PCONST = 10000)
2147 LOGICAL ridiculous_mm5_test
2148 PARAMETER (ridiculous_mm5_test = .TRUE.)
2149 ! PARAMETER (ridiculous_mm5_test = .false.)
2150
2151 ! Local variables:
2152
2153 INTEGER i , j , k
2154 INTEGER klo , khi
2155
2156
2157 REAL plo , phi , tlo, thi , zlo , zhi
2158 REAL p_at_pconst , t_at_pconst , z_at_pconst
2159 REAL z_half_lowest
2160
2161 REAL , PARAMETER :: cp = 7.*R/2.
2162 REAL , PARAMETER :: rcp = R/cp
2163 REAL , PARAMETER :: p1000mb = 100000.
2164
2165 LOGICAL l1 , l2 , l3, found
2166
2167 ! Find least zeta level that is PCONST Pa above the surface. We later use this
2168 ! level to extrapolate a surface pressure and temperature, which is supposed
2169 ! to reduce the effect of the diurnal heating cycle in the pressure field.
2170
2171 t(:,:,:) = (t(:,:,:)+300.)*(p(:,:,:)/p1000mb)**rcp
2172
2173 DO j = 1 , ny
2174 DO i = 1 , nx
2175 level(i,j) = -1
2176
2177 k = 1
2178 found = .false.
2179 do while( (.not. found) .and. (k.le.nz))
2180 IF ( p(i,j,k) .LT. p(i,j,1)-PCONST ) THEN
2181 level(i,j) = k
2182 found = .true.
2183 END IF
2184 k = k+1
2185 END DO
2186
2187 IF ( level(i,j) .EQ. -1 ) THEN
2188 PRINT '(A,I4,A)','Troubles finding level ', &
2189 NINT(PCONST)/100,' above ground.'
2190 PRINT '(A,I4,A,I4,A)', &
2191 'Problems first occur at (',i,',',j,')'
2192 PRINT '(A,F6.1,A)', &
2193 'Surface pressure = ',p(i,j,1)/100,' hPa.'
2194 STOP 'Error_in_finding_100_hPa_up'
2195 END IF
2196
2197
2198 END DO
2199 END DO
2200
2201 ! Get temperature PCONST Pa above surface. Use this to extrapolate
2202 ! the temperature at the surface and down to sea level.
2203
2204 DO j = 1 , ny
2205 DO i = 1 , nx
2206
2207 klo = MAX ( level(i,j) - 1 , 1 )
2208 khi = MIN ( klo + 1 , nz - 1 )
2209
2210 IF ( klo .EQ. khi ) THEN
2211 PRINT '(A)','Trapping levels are weird.'
2212 PRINT '(A,I3,A,I3,A)','klo = ',klo,', khi = ',khi, &
2213 ': and they should not be equal.'
2214 STOP 'Error_trapping_levels'
2215 END IF
2216
2217 plo = p(i,j,klo)
2218 phi = p(i,j,khi)
2219 tlo = t(i,j,klo)*(1. + 0.608 * q(i,j,klo) )
2220 thi = t(i,j,khi)*(1. + 0.608 * q(i,j,khi) )
2221 ! zlo = zetahalf(klo)/ztop*(ztop-terrain(i,j))+terrain(i,j)
2222 ! zhi = zetahalf(khi)/ztop*(ztop-terrain(i,j))+terrain(i,j)
2223 zlo = z(i,j,klo)
2224 zhi = z(i,j,khi)
2225
2226 p_at_pconst = p(i,j,1) - pconst
2227 t_at_pconst = thi-(thi-tlo)*LOG(p_at_pconst/phi)*LOG(plo/phi)
2228 z_at_pconst = zhi-(zhi-zlo)*LOG(p_at_pconst/phi)*LOG(plo/phi)
2229
2230 t_surf(i,j) = t_at_pconst*(p(i,j,1)/p_at_pconst)**(gamma*R/g)
2231 t_sea_level(i,j) = t_at_pconst+gamma*z_at_pconst
2232
2233 END DO
2234 END DO
2235
2236 ! If we follow a traditional computation, there is a correction to the sea level
2237 ! temperature if both the surface and sea level temnperatures are *too* hot.
2238
2239 IF ( ridiculous_mm5_test ) THEN
2240 DO j = 1 , ny
2241 DO i = 1 , nx
2242 l1 = t_sea_level(i,j) .LT. TC
2243 l2 = t_surf (i,j) .LE. TC
2244 l3 = .NOT. l1
2245 IF ( l2 .AND. l3 ) THEN
2246 t_sea_level(i,j) = TC
2247 ELSE
2248 t_sea_level(i,j) = TC - 0.005*(t_surf(i,j)-TC)**2
2249 END IF
2250 END DO
2251 END DO
2252 END IF
2253
2254 ! The grand finale: ta da!
2255
2256 DO j = 1 , ny
2257 DO i = 1 , nx
2258 ! z_half_lowest=zetahalf(1)/ztop*(ztop-terrain(i,j))+terrain(i,j)
2259 z_half_lowest=z(i,j,1)
2260 sea_level_pressure(i,j) = p(i,j,1) * &
2261 EXP((2.*g*z_half_lowest)/ &
2262 (R*(t_sea_level(i,j)+t_surf(i,j))))
2263 END DO
2264 END DO
2265
2266 if (debug) then
2267 print *,'sea pres input at weird location i=20,j=1,k=1'
2268 print *,'t=',t(20,1,1),t(20,2,1),t(20,3,1)
2269 print *,'z=',z(20,1,1),z(20,2,1),z(20,3,1)
2270 print *,'p=',p(20,1,1),p(20,2,1),p(20,3,1)
2271 print *,'slp=',sea_level_pressure(20,1), &
2272 sea_level_pressure(20,2),sea_level_pressure(20,3)
2273 endif
2274 ! print *,'t=',t(10:15,10:15,1),t(10:15,2,1),t(10:15,3,1)
2275 ! print *,'z=',z(10:15,1,1),z(10:15,2,1),z(10:15,3,1)
2276 ! print *,'p=',p(10:15,1,1),p(10:15,2,1),p(10:15,3,1)
2277 ! print *,'slp=',sea_level_pressure(10:15,10:15), &
2278 ! sea_level_pressure(10:15,10:15),sea_level_pressure(20,10:15)
2279
2280 end subroutine compute_seaprs
2281
2282 !------------------------------------------------------------------
2283
2284
2285 subroutine rotate_wind (u,v,d1,d2,d3,var, &
2286 map_proj,cen_lon,xlat,xlon, &
2287 truelat1,truelat2,data_out)
2288
2289 implicit none
2290
2291 integer, intent(in) :: d1, d2, d3
2292
2293 real, dimension(d1,d2,d3), intent(out) :: data_out
2294 integer, intent(in) :: map_proj
2295 integer ::i,j,k
2296 real, intent(in) :: cen_lon, truelat1, truelat2
2297 real :: cone
2298 real, dimension(d1,d2,d3), intent(in) :: u,v
2299 real, dimension(d1,d2), intent(in) :: xlat, xlon
2300 real, dimension(d1,d2) :: diff, alpha
2301
2302 character (len=10), intent(in) :: var
2303
2304 REAL , PARAMETER :: pii = 3.14159265
2305 REAL , PARAMETER :: radians_per_degree = pii/180.
2306
2307
2308
2309
2310 cone = 1. ! PS
2311 if( map_proj .eq. 1) then ! Lambert Conformal mapping
2312 IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
2313 cone=(ALOG(COS(truelat1*radians_per_degree))- &
2314 ALOG(COS(truelat2*radians_per_degree))) / &
2315 (ALOG(TAN((90.-ABS(truelat1))*radians_per_degree*0.5 ))- &
2316 ALOG(TAN((90.-ABS(truelat2))*radians_per_degree*0.5 )) )
2317 ELSE
2318 cone = SIN(ABS(truelat1)*radians_per_degree )
2319 ENDIF
2320 end if
2321
2322
2323 diff = xlon - cen_lon
2324
2325 do i = 1, d1
2326 do j = 1, d2
2327 if(diff(i,j) .gt. 180.) then
2328 diff(i,j) = diff(i,j) - 360.
2329 end if
2330 if(diff(i,j) .lt. -180.) then
2331 diff(i,j) = diff(i,j) + 360.
2332 end if
2333 end do
2334 end do
2335
2336 do i = 1, d1
2337 do j = 1, d2
2338 if(xlat(i,j) .lt. 0.) then
2339 alpha(i,j) = - diff(i,j) * cone * radians_per_degree
2340 else
2341 alpha(i,j) = diff(i,j) * cone * radians_per_degree
2342 end if
2343 end do
2344 end do
2345
2346
2347 if(var(1:1) .eq. "U") then
2348 do k=1,d3
2349 data_out(:,:,k) = v(:,:,k)*sin(alpha) + u(:,:,k)*cos(alpha)
2350 end do
2351 else if(var(1:1) .eq. "V") then
2352 do k=1,d3
2353 data_out(:,:,k) = v(:,:,k)*cos(alpha) - u(:,:,k)*sin(alpha)
2354 end do
2355 end if
2356
2357
2358 end subroutine rotate_wind
2359
2360
2361 !------------------------------------------------------------------
2362 subroutine handle_err(rmarker,nf_status)
2363
2364 #include "netcdf.inc"
2365 integer, intent(in) :: nf_status
2366 character*(*), intent(in) :: rmarker
2367 if (nf_status .ne. nf_noerr) then
2368 write(*,*) 'NetCDF error : ',rmarker
2369 write(*,*) ' ',nf_strerror(nf_status)
2370 stop
2371 endif
2372
2373 end subroutine handle_err
2374
2375 !------------------------------------------------------------------
2376
2377 subroutine check_all_variables ( infile, &
2378 variables3d, desc3d, number_of_3dvar, &
2379 variables2d, desc2d, number_of_2dvar, &
2380 debug )
2381
2382 #include "netcdf.inc"
2383
2384 character (len=512), intent(in) :: infile
2385 integer, intent(inout) :: number_of_3dvar ,number_of_2dvar
2386 character (len=20), dimension(number_of_3dvar), intent(inout) :: variables3d
2387 character (len=20), dimension(number_of_2dvar), intent(inout) :: variables2d
2388 character (len=50), dimension(number_of_3dvar), intent(inout) :: desc3d
2389 character (len=50), dimension(number_of_2dvar), intent(inout) :: desc2d
2390 logical, intent(in) :: debug
2391 integer :: nf_status, ncid, rcode, id_var, trcode
2392 integer :: missing3d, missing2d
2393 integer :: newi
2394
2395 nf_status = nf_open (infile, nf_nowrite, ncid)
2396 call handle_err('Error opening file',nf_status)
2397
2398
2399 missing3d = 0
2400 do i = 1,number_of_3dvar
2401 if ( variables3d(i) == 'UMET' ) then
2402 rcode = nf_inq_varid ( ncid, "U", id_var )
2403 trcode = rcode
2404 rcode = nf_inq_varid ( ncid, "V", id_var )
2405 if ( rcode == 0 ) rcode = trcode
2406 else if ( variables3d(i) == 'VMET' ) then
2407 rcode = nf_inq_varid ( ncid, "U", id_var )
2408 trcode = rcode
2409 rcode = nf_inq_varid ( ncid, "V", id_var )
2410 if ( rcode == 0 ) rcode = trcode
2411 else if ( variables3d(i) == 'Z' ) then
2412 rcode = nf_inq_varid ( ncid, "PH", id_var )
2413 trcode = rcode
2414 rcode = nf_inq_varid ( ncid, "PHB", id_var )
2415 if ( rcode == 0 ) rcode = trcode
2416 else if ( variables3d(i) == 'P' ) then
2417 rcode = nf_inq_varid ( ncid, "P", id_var )
2418 trcode = rcode
2419 rcode = nf_inq_varid ( ncid, "PB", id_var )
2420 if ( rcode == 0 ) rcode = trcode
2421 else if ( variables3d(i) == 'THETA' ) then
2422 rcode = nf_inq_varid ( ncid, "T", id_var )
2423 else if ( variables3d(i) == 'TK' ) then
2424 rcode = nf_inq_varid ( ncid, "P", id_var )
2425 trcode = rcode
2426 rcode = nf_inq_varid ( ncid, "PB", id_var )
2427 if ( rcode == 0 ) rcode = trcode
2428 trcode = rcode
2429 rcode = nf_inq_varid ( ncid, "T", id_var )
2430 if ( rcode == 0 ) rcode = trcode
2431 else if ( variables3d(i) == 'TC' ) then
2432 rcode = nf_inq_varid ( ncid, "P", id_var )
2433 trcode = rcode
2434 rcode = nf_inq_varid ( ncid, "PB", id_var )
2435 if ( rcode == 0 ) rcode = trcode
2436 trcode = rcode
2437 rcode = nf_inq_varid ( ncid, "T", id_var )
2438 if ( rcode == 0 ) rcode = trcode
2439 else if ( variables3d(i) == 'TD' ) then
2440 rcode = nf_inq_varid ( ncid, "P", id_var )
2441 trcode = rcode
2442 rcode = nf_inq_varid ( ncid, "PB", id_var )
2443 if ( rcode == 0 ) rcode = trcode
2444 trcode = rcode
2445 rcode = nf_inq_varid ( ncid, "QVAPOR", id_var )
2446 if ( rcode == 0 ) rcode = trcode
2447 else if ( variables3d(i) == 'RH' ) then
2448 rcode = nf_inq_varid ( ncid, "P", id_var )
2449 trcode = rcode
2450 rcode = nf_inq_varid ( ncid, "PB", id_var )
2451 if ( rcode == 0 ) rcode = trcode
2452 trcode = rcode
2453 rcode = nf_inq_varid ( ncid, "T", id_var )
2454 if ( rcode == 0 ) rcode = trcode
2455 trcode = rcode
2456 rcode = nf_inq_varid ( ncid, "QVAPOR", id_var )
2457 if ( rcode == 0 ) rcode = trcode
2458 else
2459 rcode = nf_inq_varid ( ncid, variables3d(i), id_var )
2460 endif
2461 if (rcode .ne. 0) then
2462 write(6,*) ' Not in file, remove from list: ',trim(variables3d(i))
2463 variables3d(i) = ' '
2464 desc3d(i) = ' '
2465 missing3d = missing3d+1
2466 endif
2467 enddo
2468
2469
2470 missing2d = 0
2471 do i = 1,number_of_2dvar
2472 if ( variables2d(i) == 'U10M' ) then
2473 rcode = nf_inq_varid ( ncid, "U10", id_var )
2474 trcode = rcode
2475 rcode = nf_inq_varid ( ncid, "V10", id_var )
2476 if ( rcode == 0 ) rcode = trcode
2477 else if ( variables2d(i) == 'V10M' ) then
2478 rcode = nf_inq_varid ( ncid, "U10", id_var )
2479 trcode = rcode
2480 rcode = nf_inq_varid ( ncid, "V10", id_var )
2481 if ( rcode == 0 ) rcode = trcode
2482 else if ( variables2d(i) == 'SLP' ) then
2483 rcode = nf_inq_varid ( ncid, "P", id_var )
2484 trcode = rcode
2485 rcode = nf_inq_varid ( ncid, "PB", id_var )
2486 if ( rcode == 0 ) rcode = trcode
2487 trcode = rcode
2488 rcode = nf_inq_varid ( ncid, "PH", id_var )
2489 if ( rcode == 0 ) rcode = trcode
2490 trcode = rcode
2491 rcode = nf_inq_varid ( ncid, "PHB", id_var )
2492 if ( rcode == 0 ) rcode = trcode
2493 trcode = rcode
2494 rcode = nf_inq_varid ( ncid, "T", id_var )
2495 if ( rcode == 0 ) rcode = trcode
2496 trcode = rcode
2497 rcode = nf_inq_varid ( ncid, "QVAPOR", id_var )
2498 if ( rcode == 0 ) rcode = trcode
2499 else if ( variables2d(i) == 'T2M' ) then
2500 rcode = nf_inq_varid ( ncid, "T2", id_var )
2501 if ( rcode == 0 ) rcode = trcode
2502 else if ( variables2d(i) == 'Q2M' ) then
2503 rcode = nf_inq_varid ( ncid, "Q2", id_var )
2504 if ( rcode == 0 ) rcode = trcode
2505 else
2506 rcode = nf_inq_varid ( ncid, variables2d(i), id_var )
2507 endif
2508 if (rcode .ne. 0) then
2509 write(6,*) ' Not in file, remove from list: ',trim(variables2d(i))
2510 variables2d(i) = ' '
2511 desc2d(i) = ' '
2512 missing2d = missing2d+1
2513 endif
2514 enddo
2515
2516
2517 newi = 0
2518 do i = 1,number_of_3dvar
2519 if ( variables3d(i) /= ' ' ) then
2520 newi = newi+1
2521 variables3d(newi) = variables3d(i)
2522 desc3d(newi) = desc3d(i)
2523 endif
2524 enddo
2525 number_of_3dvar = number_of_3dvar - missing3d
2526 newi = 0
2527 do i = 1,number_of_2dvar
2528 if ( variables2d(i) /= ' ' ) then
2529 newi = newi+1
2530 variables2d(newi) = variables2d(i)
2531 desc2d(newi) = desc2d(i)
2532 endif
2533 enddo
2534 number_of_2dvar = number_of_2dvar - missing2d
2535
2536
2537 nf_status = nf_close (ncid)
2538 call handle_err('Error closing file',nf_status)
2539
2540 end subroutine check_all_variables
2541
2542 !------------------------------------------------------------------
2543 subroutine get_dimensions(infile,nx,ny,nz)
2544
2545 #include "netcdf.inc"
2546 character (len=512), intent(in) :: infile
2547 integer :: ncid, dimid, nf_status
2548 integer, intent(inout) :: nx, ny, nz
2549 integer :: nlgen
2550
2551 ! need to pull out some data to set up dimensions, etc.
2552 nf_status = nf_open (infile, nf_nowrite, ncid)
2553 call handle_err('Error opening file',nf_status)
2554 !
2555 nf_status = nf_inq_dimid (ncid, 'west_east_stag', dimid)
2556 call handle_err('west_east_stag',nf_status)
2557 nf_status = nf_inq_dimlen (ncid, dimid, nx)
2558 call handle_err('Get NX',nf_status)
2559 nx = nx-1
2560 !
2561 nf_status = nf_inq_dimid (ncid, 'south_north_stag', dimid)
2562 call handle_err('south_north_stag',nf_status)
2563 nf_status = nf_inq_dimlen (ncid, dimid, ny)
2564 call handle_err('Get NY',nf_status)
2565 ny = ny-1
2566 !
2567 nf_status = nf_inq_dimid (ncid, 'bottom_top', dimid)
2568 call handle_err('bottom_top',nf_status)
2569 nf_status = nf_inq_dimlen (ncid, dimid, nz)
2570 call handle_err('Get NZ',nf_status)
2571 nlgen = nz
2572 !
2573 nf_status = nf_close (ncid)
2574 call handle_err('Error closing file',nf_status)
2575
2576 end subroutine get_dimensions
2577 !------------------------------------------------------------------
2578
2579 subroutine get_diffs(var1, var2, diff, absdiff, sqdiff, nx, ny, nz, missing)
2580
2581 real, intent(in), dimension(:,:,:) :: var1, var2
2582 real, intent(out), dimension(:,:,:) :: diff, absdiff, sqdiff
2583 integer, intent(in) :: nx, ny, nz
2584 integer :: i,j,k
2585 real, intent(in) :: missing
2586
2587 do k = 1, nz
2588 do j = 1, ny
2589 do i = 1, nx
2590 if ( var1(i,j,k) /= missing .and. var2(i,j,k) /= missing ) then
2591 diff(i,j,k) = var2(i,j,k) - var1(i,j,k)
2592 absdiff(i,j,k) = abs(var2(i,j,k) - var1(i,j,k))
2593 sqdiff(i,j,k) = (var2(i,j,k) - var1(i,j,k) )*(var2(i,j,k) - var1(i,j,k) )
2594 else
2595 diff(i,j,k) = missing
2596 absdiff(i,j,k) = missing
2597 sqdiff(i,j,k) = missing
2598 endif
2599 enddo
2600 enddo
2601 enddo
2602
2603 end subroutine get_diffs
2604 !------------------------------------------------------------------
2605 subroutine domain_average(var, avg_prof, counter, nx, ny, nz, missing,isq)
2606
2607 integer, intent(in) :: nx,ny,nz,isq
2608 real, intent(in), dimension(:,:,:) :: var
2609 integer, intent(out), dimension(:) :: counter
2610 real, intent(out), dimension(:) :: avg_prof
2611 real, intent(in) :: missing
2612
2613 integer :: i,j,k,icount
2614 real :: dsum, dmiss
2615 integer :: imiss
2616
2617 !9999
2618 !Hui: set dmiss value consistent with plot script
2619 ! dmiss = -9999.999
2620 dmiss = -99.99
2621 imiss = -99
2622 do k = 1, nz
2623 icount = 0
2624 dsum = 0
2625 do j = 1, ny
2626 do i = 1, nx
2627 if ( var(i,j,k) /= missing ) then
2628 icount = icount + 1
2629 dsum = dsum + var(i,j,k)
2630 endif
2631 enddo
2632 enddo
2633 avg_prof(k) = dmiss
2634 !Hui counter(k) = 0
2635 counter(k) = imiss
2636 if (icount /= 0 ) then
2637 counter(k) = icount
2638 if ( isq .eq. 0 ) then
2639 avg_prof(k) = dsum /float(icount)
2640 else
2641 avg_prof(k) = sqrt(dsum /float(icount))
2642 endif
2643 endif
2644 enddo
2645
2646 end subroutine domain_average
2647 !------------------------------------------------------------------
2648
2649 subroutine get_sum(dsum, dvar, nx, ny, nz, missing)
2650
2651 integer, intent(in) :: nx, ny, nz
2652 real, intent(in) :: missing
2653 real, intent(in),dimension(:,:,:) :: dvar
2654 real, intent(inout),dimension(:,:,:) :: dsum
2655
2656 integer :: i,j,k
2657
2658 do k = 1, nz
2659 do j = 1, ny
2660 do i = 1, nx
2661 if ( dvar(i,j,k) /= missing .and. dsum(i,j,k) /= missing ) then
2662 dsum(i,j,k) = dsum(i,j,k) + dvar(i,j,k)
2663 else
2664 dsum(i,j,k) = missing
2665 endif
2666 enddo
2667 enddo
2668 enddo
2669
2670 end subroutine get_sum
2671 !------------------------------------------------------------------
2672 subroutine time_average(dvar, nx, ny, nz, time_count,missing, isqr)
2673
2674 integer, intent(in) :: nx, ny, nz,time_count,isqr
2675 real, intent(in) :: missing
2676 real, intent(inout), dimension(:,:,:) :: dvar
2677
2678 integer :: i,j,k
2679
2680 do k = 1, nz
2681 do j = 1, ny
2682 do i = 1, nx
2683 if ( dvar(i,j,k) /= missing ) then
2684 if (isqr .eq. 1 ) then
2685 dvar(i,j,k) = sqrt(dvar(i,j,k)/float(time_count))
2686 else
2687 dvar(i,j,k) = dvar(i,j,k)/float(time_count)
2688 endif
2689 else
2690 dvar(i,j,k) = missing
2691 endif
2692 enddo
2693 enddo
2694 enddo
2695
2696 end subroutine time_average
2697 !------------------------------------------------------------------
2698
2699 subroutine compute_data_3d( infile, var, length, nx, ny, nz, levels, time_idx, &
2700 vert_args, vertical_type, missing, data_out_z, debug )
2701
2702 integer, intent(in) :: time_idx
2703 integer, intent(in) :: nx, ny, nz, levels
2704 integer, intent(in) :: length
2705 real, intent(in) :: missing
2706 real, intent(in) :: vert_args(100)
2707 character (len=1), intent(in) :: vertical_type
2708 character (len=30), intent(in) :: var
2709 character (len=512), intent(in) :: infile
2710 logical, intent(in) :: debug
2711
2712 real, intent(out), dimension(:,:,:) :: data_out_z
2713 real, allocatable, dimension(:,:,:) :: data_out
2714 real, allocatable, dimension(:,:,:) :: z, ph, phb
2715 real, allocatable, dimension(:,:,:) :: p, pb
2716
2717
2718 ! first, get some base-state-stuff
2719
2720 if ( vertical_type == 'p' ) then
2721 allocate( p (nx, ny, nz) )
2722 allocate( pb(nx, ny, nz) )
2723 call da_get_var_3d_real_cdf( infile,'PB',pb,nx,ny,nz,time_idx,debug )
2724 endif
2725 if ( vertical_type == 'z' ) then
2726 allocate( z (nx, ny, nz) )
2727 allocate( ph (nx, ny, nz+1) )
2728 allocate( phb(nx, ny, nz+1) )
2729 call da_get_var_3d_real_cdf( infile,'PHB',phb,nx,ny,nz+1,time_idx,debug )
2730 endif
2731 ! first, get p/z if needed
2732 if ( vertical_type == 'p' ) then
2733 call da_get_var_3d_real_cdf( infile,'P',p, nx, ny, nz, time_idx,debug )
2734 p = p+pb
2735 endif
2736 if ( vertical_type == 'z' ) then
2737 call da_get_var_3d_real_cdf( infile,'PH',ph, nx, ny, nz+1, time_idx,debug )
2738 ph = (ph+phb)/9.81
2739 z = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1))
2740 ! need to convert to kilometers for coordinate
2741 z = z/1000.
2742 endif
2743
2744 allocate ( data_out (nx, ny, nz) )
2745
2746 call g_output_3d (infile, time_idx, var, length, nx, ny, nz, data_out, debug)
2747
2748 if ( vertical_type == 'p' ) then
2749 call interp_to_z( data_out, nx, ny, nz, data_out_z, nx, ny, levels, &
2750 p/100., vert_args, missing, vertical_type, debug )
2751
2752 else if ( vertical_type == 'z' ) then
2753 call interp_to_z( data_out, nx, ny, nz, data_out_z, nx, ny, levels, &
2754 z, vert_args, missing, vertical_type, debug )
2755 else
2756 data_out_z = data_out
2757 endif
2758 deallocate ( data_out )
2759 if ( vertical_type == 'p' ) then
2760 deallocate( p )
2761 deallocate( pb )
2762 endif
2763 if ( vertical_type == 'z' ) then
2764 deallocate( z )
2765 deallocate( ph )
2766 deallocate( phb )
2767 endif
2768
2769 end subroutine compute_data_3d
2770
2771 !---------------------------------------------------------------------
2772 end program da_verif_anal
Mirror of the official repository at https://github.com/wrf-model/WRF