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