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