?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge remote-tracking branch 'origin/release-v4.5.2'
[WRF.git] / var / da / da_verif_grid / da_verif_grid.f90
bloba7bf2752a99d662b9913f34a24bbdb0003ed38e2
1 program da_verif_grid
2 !---------------------------------------------------------------------------
3 ! History:
4 !
5 ! Abstract: Program to calculate statistics for multiple experiments for
6 ! verification analayis/forecasts against any desired analysis
7 !
8 ! Author: Syed RH Rizvi NCAR/MMM 05/25/2006
9 ! Udates:
10 ! Syed RH Rizvi NCAR/MMM 06/05/2009
11 ! a) Added verification for wind vector & geopotentials
12 ! b) For vertical profiles, added the choice to display the desired
13 ! number of pressure levels
14 !---------------------------------------------------------------------------
15
16 use da_verif_grid_control, only : control_main, control_times, control_vars, sub_domain, &
17 max_3d_variables, max_2d_variables,num_vert_levels,verification_file_string,&
18 missing,namelist_unit,time_series_unit, time_series_2d, profile_time_series_3d,&
19 filename, stime, etime, hstart, hend, hdate, date, pdate, vert_levels, &
20 nx, ny, nz, io_status, debug1, debug2, verify_its_own_analysis, &
21 num_verifying_experiments, verify_forecast_hour, domain, control_exp_dir, verif_dirs, &
22 out_dirs,start_year, end_year, start_month, end_month, start_day, end_day, &
23 start_hour, end_hour,start_minutes, end_minutes, start_seconds, end_seconds,interval_hour, &
24 num3dvar, num2dvar, var3d, var2d, num_scores, score_names, vertical_type, &
25 istart, iend, jstart, jend, unit_all, unit_land, unit_water
26
27 use da_netcdf_interface, only : da_get_dims_cdf, da_get_gl_att_int_cdf, da_get_gl_att_real_cdf, &
28 da_get_var_3d_real_cdf, da_get_var_2d_real_cdf, da_get_var_2d_int_cdf
29
30 implicit none
31
32 character (len=512) :: control_file, verif_file
33 character (len=512) :: out_dir, outname_all, outname_land, outname_water
34
35 integer, parameter :: imiss = -99
36 real, parameter :: rmiss = -99.99
37 integer :: time_loop_count
38 integer :: time(6), ptime(6)
39 integer :: nx1, ny1, nz1
40 integer :: nx2, ny2, nz2
41 integer :: i,k
42 integer :: ivar, iexp, iscore
43 character (len=10) :: sdate
44 character (len=20) :: file_string, domain_string, out_hr
45 logical, allocatable,dimension(:) :: first_score
46 logical :: exist_file1, exist_file2
47
48 real, allocatable, dimension(:,:,:) :: data_out1, data_out2
49 real, allocatable, dimension(:,:,:) :: u1, u2, v1, v2
50
51 real, allocatable, dimension(:,:,:) :: sum3d, asum3d, sqr3d, diff, absdiff, sqdiff
52 real, allocatable, dimension(:,:) :: score_avg_prof
53 real, allocatable, dimension(:) :: avg_prof
54 real, allocatable, dimension(:,:) :: landmask ! 1 for land, 0 for water
55 integer, parameter :: island = 1, iswater = 0
56 integer, allocatable, dimension(:) :: num_counter
57 integer, allocatable, dimension(:,:,:) :: mask_count
58 real, allocatable, dimension(:,:,:) :: mask_stats ! 1st dimension: 1:bias, 2:rmse, 3:abias
59 ! 2nd dimension: 1:all, 2:land, 3:water
60 ! 3rd dimension: number of vertical levels
61
62 !----------------------------------------------------------------------------
63 verify_forecast_hour = 0
64 !----------------------------------------------------------------------------
65 debug1 = .false.
66 debug2 = .false.
67 vertical_type = 'p'
68
69 !--3D need update
70 num3dvar=max_3d_variables
71 var3d(1)='U'
72 var3d(2)='V'
73 var3d(3)='T'
74 var3d(4)='QVAPOR'
75 var3d(5)='Z'
76 var3d(6)='WV'
77
78 !--2D need update
79 num2dvar=max_2d_variables
80 var2d(1)='SLP'
81 var2d(2)='PSFC'
82 var2d(3)='U10M'
83 var2d(4)='V10M'
84 var2d(5)='T2M'
85 var2d(6)='Q2M'
86 var2d(7)='MU'
87
88 !--Score names
89 num_scores = 3
90 score_names(1) = 'BIAS'
91 score_names(2) = 'RMSE'
92 score_names(3) = 'ABIAS'
93
94 domain = 1
95 verification_file_string = 'wrfout'
96 out_hr ='_00'
97 !---------------------------------------------------------------------
98 ! Read namelist
99 !----------------------------------------------------------------------------
100 io_status = 0
101
102
103 open(unit = namelist_unit, file = 'namelist.in', &
104 status = 'old' , access = 'sequential', &
105 form = 'formatted', action = 'read', &
106 iostat = io_status )
107
108
109 if(io_status /= 0) then
110 print *, 'Error to open namelist.in file: '
111 else
112 read(unit=namelist_unit, nml = control_main , iostat = io_status)
113
114 if(io_status /= 0) then
115 print *, 'Error to read control_main. Stopped.'
116 stop
117 endif
118
119
120 read(unit=namelist_unit, nml = control_times , iostat = io_status )
121
122 if(io_status /= 0) then
123 print *, 'Error to read control_times Stopped.'
124 stop
125 endif
126
127 read(unit=namelist_unit, nml = control_vars , iostat = io_status )
128
129 if(io_status /= 0) then
130 print *, 'Error to read control_vars Stopped.'
131 stop
132 endif
133
134 read(unit=namelist_unit, nml = sub_domain , iostat = io_status )
135
136 if(io_status /= 0) then
137 print *, 'Error reading sub_domain'
138 istart = 1
139 iend = 10000
140 jstart = 1
141 jend = 10000
142 !stop
143 endif
144
145 close(unit=namelist_unit)
146 endif
147 !----------------------------------------------------------------------------
148 !---------------------------------------------------------------------
149 !
150 write(domain_string, fmt ='("_d",i2.2,"_")') domain
151 file_string = trim(verification_file_string)//trim(domain_string)
152 write(out_hr, fmt ='("_",i2.2)') verify_forecast_hour
153 allocate(first_score(num_scores))
154 !---------------------------------------------------------------------
155
156 stime(1) = start_year
157 stime(2) = start_month
158 stime(3) = start_day
159 stime(4) = start_hour
160 stime(5) = start_minutes
161 stime(6) = start_seconds
162
163
164 call build_hdate(hstart, stime )
165
166 etime(1) = end_year
167 etime(2) = end_month
168 etime(3) = end_day
169 etime(4) = end_hour
170 etime(5) = end_minutes
171 etime(6) = end_seconds
172
173 call build_hdate(hend, etime )
174
175
176 if ( hend < hstart ) then
177 print*, '****************************************************************'
178 print*, 'End time is before the start time'
179 print*, ' Start time = ', hstart ,' & End time = ', hend
180 print*, '****************************************************************'
181 stop
182 endif
183
184 hdate = hstart
185 call build_hdate(sdate, stime )
186
187 filename = trim(file_string)//hdate
188
189 !---------------------------------------------------------------------
190 loop_verif_exp : do iexp = 1, num_verifying_experiments
191
192 first_score = .true.
193 out_dir = trim(out_dirs(iexp))//'/'
194 !---------------------------------------------------------------------
195 !--For 3D variables
196 !---------------------------------------------------------------------
197 allocate( avg_prof(num_vert_levels) )
198 allocate( num_counter(num_vert_levels) )
199 allocate( score_avg_prof(num_scores, num_vert_levels) )
200
201 loop_3d : do ivar=1,num3dvar
202
203 ! Open profile units
204
205 profile_time_series_3d = trim(out_dir)//trim(var3d(ivar))//'_time_series'//trim(out_hr)
206 open(time_series_unit, file=profile_time_series_3d,form='formatted', &
207 status='unknown')
208
209 outname_all = trim(profile_time_series_3d)//'_sub'
210 outname_land = trim(profile_time_series_3d)//'_sub_land'
211 outname_water = trim(profile_time_series_3d)//'_sub_water'
212 open(unit_all, file=trim(outname_all), form='formatted', status='unknown')
213 open(unit_land, file=trim(outname_land), form='formatted', status='unknown')
214 open(unit_water, file=trim(outname_water),form='formatted', status='unknown')
215 !----------------------------------------------------------------------------
216 !
217 !----------------------------------------------------------------------------
218 time_loop_count = 0
219 hdate = hstart
220 time = stime
221
222 time_loop_3d : do
223
224 !----------------------------------------------------------------------------
225 call build_hdate(hdate, time)
226 print*,' processing exp: ',iexp,' 3d var: ',trim(var3d(ivar)),' for time: ',hdate
227
228 if ( hdate > hend ) exit time_loop_3d
229
230 call build_hdate(date, time )
231 ptime = time
232 call advance_date(ptime,-verify_forecast_hour)
233 call build_hdate(pdate, ptime)
234
235
236 filename = trim(file_string)//hdate
237
238 if( verify_its_own_analysis) then
239 control_file = trim(verif_dirs(iexp))//'/'//date//'/'//trim(filename)
240 else
241 control_file = trim(control_exp_dir)//'/'//date//'/'//trim(filename)
242 endif
243
244 verif_file = trim(verif_dirs(iexp))//'/'//pdate//'/'//trim(filename)
245
246 inquire(file=trim(control_file), exist=exist_file1)
247 inquire(file=trim(verif_file), exist=exist_file2)
248 if ( (.not. exist_file1) .or. (.not. exist_file2) ) then
249 avg_prof = rmiss
250 num_counter = imiss
251 score_avg_prof = rmiss
252 call write_profile(date, score_avg_prof, num_counter, num_vert_levels, num_scores, &
253 time_series_unit)
254 call write_profile(date, score_avg_prof, num_counter, num_vert_levels, num_scores, unit_all)
255 call write_profile(date, score_avg_prof, num_counter, num_vert_levels, num_scores, unit_land)
256 call write_profile(date, score_avg_prof, num_counter, num_vert_levels, num_scores, unit_water)
257 call advance_date(time,interval_hour)
258 if ( .not. exist_file1 ) print*, trim(control_file), ' file not found'
259 if ( .not. exist_file2 ) print*, trim(verif_file), ' file not found'
260 print*, 'skipping date ', date
261 cycle time_loop_3d
262 end if
263
264 ! print*,' control_file : ', trim(control_file)
265 ! print*,' verif_file : ', trim(verif_file)
266
267 ! Get the dimensions of the both files and check
268 call get_dimensions(control_file,nx1,ny1,nz1)
269 call get_dimensions(verif_file,nx2,ny2,nz2)
270
271 if ( nx1 /= nx2 .or. ny1 /= ny2 .or. nz1 /= nz2 ) then
272 print*, '********************************************************'
273 print*, 'Dimension mismatch between files of the experiments ....'
274 print*, '********************************************************'
275 stop
276 else
277 nx = nx1
278 ny = ny1
279 nz = nz1
280 if (time_loop_count == 0 ) then
281 allocate( sum3d(nx,ny,num_vert_levels))
282 allocate( asum3d(nx,ny,num_vert_levels))
283 allocate( sqr3d(nx,ny,num_vert_levels))
284 sum3d = 0.0
285 asum3d = 0.0
286 sqr3d = 0.0
287
288 if ( ivar == 1 ) then
289 allocate(landmask(nx,ny))
290 call da_get_var_2d_real_cdf(verif_file,"LANDMASK", landmask, nx, ny, 1, .false.)
291 istart = max(1, istart)
292 iend = min(nx, iend)
293 jstart = max(1, jstart)
294 jend = min(ny, jend)
295 end if
296
297 endif
298 endif
299 allocate(diff(nx,ny,num_vert_levels))
300 allocate(absdiff(nx,ny,num_vert_levels))
301 allocate(sqdiff(nx,ny,num_vert_levels))
302 ! first, get control data
303 if( trim(var3d(ivar)) .eq. "WV" ) then
304 allocate ( u1 (nx, ny, num_vert_levels) )
305 allocate ( u2 (nx, ny, num_vert_levels) )
306 allocate ( v1 (nx, ny, num_vert_levels) )
307 allocate ( v2 (nx, ny, num_vert_levels) )
308 call compute_wind_3d( control_file, nx, ny, nz, num_vert_levels, 1, &
309 vert_levels, vertical_type, missing, u1, v1, debug1 )
310 call compute_wind_3d( verif_file, nx, ny, nz, num_vert_levels, 1, &
311 vert_levels, vertical_type, missing, u2, v2, debug1 )
312 call get_diffs_wind(u1, v1, u2, v2, diff, absdiff, sqdiff, nx, ny, &
313 num_vert_levels, missing)
314 deallocate(u1, v1, u2, v2)
315 else
316 allocate ( data_out1 (nx, ny, num_vert_levels) )
317 allocate ( data_out2 (nx, ny, num_vert_levels) )
318 call compute_data_3d( control_file, trim(var3d(ivar)), len_trim(var3d(ivar)), &
319 nx, ny, nz, num_vert_levels, 1, vert_levels, &
320 vertical_type, missing, data_out1, debug1 )
321 ! second, get verifying data
322 call compute_data_3d( verif_file, trim(var3d(ivar)), len_trim(var3d(ivar)), &
323 nx, ny, nz, num_vert_levels, 1, vert_levels, &
324 vertical_type, missing, data_out2, debug2 )
325
326 call get_diffs(data_out1, data_out2, diff, absdiff, sqdiff, nx, ny, &
327 num_vert_levels, missing)
328 deallocate(data_out1)
329 deallocate(data_out2)
330 end if
331
332 do iscore = 1, num_scores
333 if ( trim(score_names(iscore)) == 'BIAS' ) then
334 call domain_average( diff, avg_prof, num_counter, nx, ny, num_vert_levels, missing,0)
335 elseif ( trim(score_names(iscore)) == 'RMSE' ) then
336 call domain_average( sqdiff, avg_prof, num_counter, nx, ny, num_vert_levels, missing,1)
337 elseif ( trim(score_names(iscore)) == 'ABIAS' ) then
338 call domain_average( absdiff, avg_prof, num_counter, nx, ny, num_vert_levels, missing,0)
339 endif
340 score_avg_prof(iscore,:) = avg_prof(:)
341 enddo
342 call write_profile(date, score_avg_prof, num_counter, num_vert_levels, num_scores, &
343 time_series_unit)
344
345 allocate(mask_stats(3,3,nz))
346 allocate(mask_count(3,3,nz))
347 call mask_domain_average( diff, mask_stats(1,:,:), mask_count(1,:,:), nx, ny, num_vert_levels, missing,0)
348 call mask_domain_average( sqdiff, mask_stats(2,:,:), mask_count(2,:,:), nx, ny, num_vert_levels, missing,1)
349 call mask_domain_average( absdiff, mask_stats(3,:,:), mask_count(3,:,:), nx, ny, num_vert_levels, missing,0)
350 call write_profile(date, mask_stats(:,1,:), mask_count(1,1,:), num_vert_levels, 3, unit_all)
351 call write_profile(date, mask_stats(:,2,:), mask_count(1,2,:), num_vert_levels, 3, unit_land)
352 call write_profile(date, mask_stats(:,3,:), mask_count(1,3,:), num_vert_levels, 3, unit_water)
353 deallocate(mask_stats)
354 deallocate(mask_count)
355
356 call get_sum(sum3d,diff,nx,ny,num_vert_levels,missing)
357 call get_sum(asum3d,absdiff,nx,ny,num_vert_levels,missing)
358 call get_sum(sqr3d,sqdiff,nx,ny,num_vert_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 close(time_series_unit)
371 close(unit_all)
372 close(unit_land)
373 close(unit_water)
374 deallocate(sum3d)
375 deallocate(asum3d)
376 deallocate(sqr3d)
377
378 enddo loop_3d
379 print*, ' successful completion of loop_3d '
380
381 deallocate( avg_prof )
382 deallocate( num_counter )
383 deallocate( score_avg_prof )
384
385 !--------------------------------------------------------------------------------
386 !--Loop For 2D variables
387 !--------------------------------------------------------------------------------
388 allocate( avg_prof(1) )
389 allocate( num_counter(1) )
390 allocate( score_avg_prof(num_scores, 1) )
391
392 loop_2d : do ivar = 1, num2dvar
393
394 ! Open profile units
395
396 time_series_2d = trim(out_dir)//trim(var2d(ivar))//'_time_series'//trim(out_hr)
397 open(time_series_unit, file=time_series_2d,form='formatted', &
398 status='unknown')
399
400 outname_all = trim(time_series_2d)//'_sub'
401 outname_land = trim(time_series_2d)//'_sub_land'
402 outname_water = trim(time_series_2d)//'_sub_water'
403 open(unit_all, file=trim(outname_all), form='formatted', status='unknown')
404 open(unit_land, file=trim(outname_land), form='formatted', status='unknown')
405 open(unit_water, file=trim(outname_water),form='formatted', status='unknown')
406 !----------------------------------------------------------------------------
407 !
408 !----------------------------------------------------------------------------
409
410 time_loop_count = 0
411 hdate = hstart
412 time = stime
413
414 time_loop_2d : do
415
416 !----------------------------------------------------------------------------
417 call build_hdate(hdate, time )
418 print*,' processing exp: ',iexp,' 2d var: ',trim(var2d(ivar)),' for time: ',hdate
419
420 if ( hdate > hend ) exit time_loop_2d
421
422 call build_hdate(date, time )
423 ptime = time
424 call advance_date(ptime,-verify_forecast_hour)
425 call build_hdate(pdate, ptime)
426
427 filename = trim(file_string)//hdate
428 if( verify_its_own_analysis) then
429 control_file = trim(verif_dirs(iexp))//'/'//date//'/'//trim(filename)
430 else
431 control_file = trim(control_exp_dir)//'/'//date//'/'//trim(filename)
432 endif
433
434 verif_file = trim(verif_dirs(iexp))//'/'//pdate//'/'//trim(filename)
435
436 inquire(file=trim(control_file), exist=exist_file1)
437 inquire(file=trim(verif_file), exist=exist_file2)
438 if ( (.not. exist_file1) .or. (.not. exist_file2) ) then
439 avg_prof = rmiss
440 num_counter = imiss
441 score_avg_prof = rmiss
442 call write_profile(date, score_avg_prof, num_counter, 1, num_scores, &
443 time_series_unit)
444 call write_profile(date, score_avg_prof, num_counter, 1, num_scores, unit_all)
445 call write_profile(date, score_avg_prof, num_counter, 1, num_scores, unit_land)
446 call write_profile(date, score_avg_prof, num_counter, 1, num_scores, unit_water)
447 call advance_date(time,interval_hour)
448 if ( .not. exist_file1 ) print*, trim(control_file), ' file not found'
449 if ( .not. exist_file2 ) print*, trim(verif_file), ' file not found'
450 print*, 'skipping date ', date
451 cycle time_loop_2d
452 end if
453
454 !
455 ! Get the dimensions of the both files and check
456 call get_dimensions(control_file,nx1,ny1,nz1)
457 call get_dimensions(verif_file,nx2,ny2,nz2)
458
459 if ( nx1 /= nx2 .or. ny1 /= ny2 .or. nz1 /= nz2 ) then
460 print*, '********************************************************'
461 print*, 'Dimension mismatch between files of the experiments ....'
462 print*, '********************************************************'
463 stop
464 else
465 nx = nx1
466 ny = ny1
467 nz = nz1
468 if (time_loop_count == 0 ) then
469 allocate( sum3d(nx,ny,1))
470 allocate( asum3d(nx,ny,1))
471 allocate( sqr3d(nx,ny,1))
472 sum3d = 0.0
473 asum3d = 0.0
474 sqr3d = 0.0
475 endif
476 endif
477
478 allocate(data_out1(nx, ny, 1))
479 allocate(data_out2(nx, ny, 1))
480
481 call g_output_2d (control_file, 1, trim(var2d(ivar)), len_trim(var2d(ivar)), &
482 nx, ny, nz, data_out1, debug1)
483
484 call g_output_2d (verif_file, 1, trim(var2d(ivar)), len_trim(var2d(ivar)), &
485 nx, ny, nz, data_out2, debug2)
486
487 allocate(diff(nx,ny,1))
488 allocate(absdiff(nx,ny,1))
489 allocate(sqdiff(nx,ny,1))
490 call get_diffs(data_out1, data_out2, diff, absdiff, sqdiff, nx, ny, 1, missing)
491 deallocate(data_out1)
492 deallocate(data_out2)
493
494 do iscore = 1, num_scores
495 if ( trim(score_names(iscore)) == 'BIAS' ) then
496 call domain_average( diff, avg_prof, num_counter, nx, ny, 1, missing,0)
497 elseif ( trim(score_names(iscore)) == 'RMSE' ) then
498 call domain_average( sqdiff, avg_prof, num_counter, nx, ny, 1, missing,1)
499 elseif ( trim(score_names(iscore)) == 'ABIAS' ) then
500 call domain_average( absdiff, avg_prof, num_counter, nx, ny, 1, missing,0)
501 endif
502 score_avg_prof(iscore,:) = avg_prof(:)
503 enddo
504 call write_profile(date, score_avg_prof, num_counter, 1, num_scores, &
505 time_series_unit)
506
507 allocate(mask_stats(3,3,1))
508 allocate(mask_count(3,3,1))
509 call mask_domain_average( diff, mask_stats(1,:,:), mask_count(1,:,:), nx, ny, 1, missing,0)
510 call mask_domain_average( sqdiff, mask_stats(2,:,:), mask_count(2,:,:), nx, ny, 1, missing,1)
511 call mask_domain_average( absdiff, mask_stats(3,:,:), mask_count(3,:,:), nx, ny, 1, missing,0)
512 call write_profile(date, mask_stats(:,1,:), mask_count(1,1,:), 1, 3, unit_all)
513 call write_profile(date, mask_stats(:,2,:), mask_count(1,2,:), 1, 3, unit_land)
514 call write_profile(date, mask_stats(:,3,:), mask_count(1,3,:), 1, 3, unit_water)
515 deallocate(mask_stats)
516 deallocate(mask_count)
517
518 call get_sum(sum3d,diff,nx,ny,1,missing)
519 call get_sum(asum3d,absdiff,nx,ny,1,missing)
520 call get_sum(sqr3d,sqdiff,nx,ny,1,missing)
521
522 deallocate(diff)
523 deallocate(absdiff)
524 deallocate(sqdiff)
525 time_loop_count = time_loop_count + 1
526
527 call advance_date(time,interval_hour)
528
529
530 enddo time_loop_2d
531
532 close(time_series_unit)
533 close(unit_all)
534 close(unit_land)
535 close(unit_water)
536 deallocate(sum3d)
537 deallocate(asum3d)
538 deallocate(sqr3d)
539
540 enddo loop_2d
541
542 deallocate( avg_prof )
543 deallocate( num_counter )
544 deallocate( score_avg_prof )
545 deallocate( landmask )
546
547 print*, ' successful completion of loop_2d '
548 print*,' Finished Experiment : ', trim(verif_dirs(iexp))
549 enddo loop_verif_exp
550
551 !-----------------------------------------------------
552
553 contains
554 !-----------------------------------------------------
555 subroutine advance_date( time, delta )
556
557 implicit none
558
559 integer, intent(inout) :: time(6)
560 integer, intent(in) :: delta
561
562 integer :: ccyy, mm, dd, hh
563 integer, dimension(12) :: mmday
564 ! character(len=10) :: ccyymmddhh
565
566 !-----------------------------------------------------
567 mmday = (/31,28,31,30,31,30,31,31,30,31,30,31/)
568 mmday(2) = 28
569 !-----------------------------------------------------
570 ccyy = time(1)
571 mm = time(2)
572 dd = time(3)
573 hh = time(4)
574 !-----------------------------------------------------
575
576 hh = hh + delta
577
578 do while (hh < 0)
579 hh = hh + 24
580 call change_date ( ccyy, mm, dd, -1 )
581 end do
582
583 do while (hh > 23)
584 hh = hh - 24
585 call change_date ( ccyy, mm, dd, 1 )
586 end do
587
588 ! write(ccyymmddhh(1:10), fmt='(i4, 3i2.2)') ccyy, mm, dd, hh
589 time(1) = ccyy
590 time(2) = mm
591 time(3) = dd
592 time(4) = hh
593
594 end subroutine advance_date
595 !---------------------------------------------------------------------------
596 subroutine change_date ( ccyy, mm, dd, delta)
597 integer, intent(inout) :: ccyy, mm, dd
598 integer, intent(in) :: delta
599 !
600 integer, dimension(12) :: mmday
601 mmday = (/31,28,31,30,31,30,31,31,30,31,30,31/)
602
603 mmday(2) = 28
604
605 if (mod(ccyy,4) == 0) then
606 mmday(2) = 29
607
608 if ( mod(ccyy,100) == 0) then
609 mmday(2) = 28
610 endif
611
612 if(mod(ccyy,400) == 0) then
613 mmday(2) = 29
614 end if
615 endif
616
617 dd = dd + delta
618
619 if(dd == 0) then
620 mm = mm - 1
621
622 if(mm == 0) then
623 mm = 12
624 ccyy = ccyy - 1
625 endif
626
627 dd = mmday(mm)
628 elseif ( dd .gt. mmday(mm) ) then
629 dd = 1
630 mm = mm + 1
631 if(mm > 12 ) then
632 mm = 1
633 ccyy = ccyy + 1
634 end if
635 end if
636 end subroutine change_date
637
638 subroutine build_hdate(hdate, time)
639
640 ! PURPOSE:
641 ! From the Year, Month, Day, Hour, Minute, and Second values,
642 ! creates a 19-character string representing the date, in the
643 ! format: "YYYY-MM-DD hh:mm:ss"
644
645 ! INPUT:
646 integer, intent(in) :: time(6) ! all time specs in it
647 ! OUTPUT:
648 character*(*), intent(out) :: hdate ! 'YYYY-MM-DD hh:mm:ss'
649
650 ! LOCAL:
651 integer iyr ! year (e.g., 1997, 2001)
652 integer imo ! month (01 - 12)
653 integer idy ! day of the month (01 - 31)
654 integer ihr ! hour (00-23)
655 integer imi ! minute (00-59)
656 integer isc ! second (00-59)
657 !
658 ! integer i ! Loop counter.
659 integer hlen ! Length of hdate string
660
661 hlen = len(hdate)
662 iyr = time(1)
663 imo = time(2)
664 idy = time(3)
665 ihr = time(4)
666 imi = time(5)
667 isc = time(6)
668
669 if (hlen.eq.19) then
670 write(hdate,19) iyr, imo, idy, ihr, imi, isc
671 19 format(i4,'-',i2.2,'-',i2.2,'_',i2.2,':',i2.2,':',i2.2)
672
673 elseif (hlen.eq.16) then
674 write(hdate,16) iyr, imo, idy, ihr, imi
675 16 format(i4,'-',i2.2,'-',i2.2,'_',i2.2,':',i2.2)
676
677 elseif (hlen.eq.13) then
678 write(hdate,13) iyr, imo, idy, ihr
679 13 format(i4,'-',i2.2,'-',i2.2,'_',i2.2)
680
681 elseif (hlen.eq.10) then
682 write(hdate,10) iyr, imo, idy, ihr
683 10 format(i4,i2.2,i2.2,i2.2)
684 endif
685
686 return
687 end subroutine build_hdate
688
689 !----------------------------------------------------------------------------------
690 subroutine write_profile(date, profile, counter, nlevel, nscore, out_unit)
691
692 integer, intent(in) :: nlevel, nscore, out_unit
693 real, intent(in), dimension(:,:) :: profile
694 integer, intent(in), dimension(:) :: counter
695 character (len=10), intent(in) :: date
696 write(out_unit,fmt='(a10,1x,100(i8,1x,3(f14.8,1x)))') date, &
697 (counter(k), (profile(i,k),i=1,nscore),k=1,nlevel)
698
699 end subroutine write_profile
700
701 !---------------------------------------------------------------------------------
702 subroutine time_calc( time, timestamp, datestamp, debug , tdef,it)
703
704 implicit none
705
706 character (len=19), intent(in) :: time
707 character (len=35), intent(inout) :: tdef
708 integer, intent(out) :: timestamp, datestamp
709 logical, intent(in) :: debug
710
711 integer, intent(in) :: it
712 integer :: hours, minutes, seconds, year, month, day,hour1,hourint
713 integer :: mins1,minsint
714
715 save hourint
716 save minsint
717
718 read(time(18:19),*) seconds
719 read(time(15:16),*) minutes
720 read(time(12:13),*) hours
721 read(time(1:4),*) year
722 read(time(6:7),*) month
723 read(time(9:10),*) day
724
725 if(debug) write(6,*) ' day, month, year, hours, minutes, seconds '
726 if(debug) write(6,*) day, month, year, hours, minutes, seconds
727
728 if ( it == 1) then
729 write (tdef(19:20),'(i2)') hours
730 if ( day < 10 ) then
731 write (tdef(23:23),'(i1)') day
732 else
733 write (tdef(22:23),'(i2)') day
734 endif
735 write (tdef(27:30),'(i4)') year
736 if (month == 1) write (tdef(24:26),'(a3)') 'jan'
737 if (month == 2) write (tdef(24:26),'(a3)') 'feb'
738 if (month == 3) write (tdef(24:26),'(a3)') 'mar'
739 if (month == 4) write (tdef(24:26),'(a3)') 'apr'
740 if (month == 5) write (tdef(24:26),'(a3)') 'may'
741 if (month == 6) write (tdef(24:26),'(a3)') 'jun'
742 if (month == 7) write (tdef(24:26),'(a3)') 'jul'
743 if (month == 8) write (tdef(24:26),'(a3)') 'aug'
744 if (month == 9) write (tdef(24:26),'(a3)') 'sep'
745 if (month ==10) write (tdef(24:26),'(a3)') 'oct'
746 if (month ==11) write (tdef(24:26),'(a3)') 'nov'
747 if (month ==12) write (tdef(24:26),'(a3)') 'dec'
748 hour1=hours
749 mins1=minutes
750 elseif ( it == 2) then
751 hourint = abs(hours-hour1)
752 minsint = abs(minutes-mins1)
753 if (hourint == 0 ) then
754 if (minsint == 0 ) minsint = 1
755 if(debug) write(6,*) "interval is",minsint
756 write (tdef(34:35),'(a2)') "mn"
757 write (tdef(32:33),'(i2)') minsint
758 if(debug) write(6,*) "TDEF is",tdef
759 else
760 if(debug) write(6,*) "Interval is",hourint
761 write (tdef(32:33),'(i2)') hourint
762 if(debug) write(6,*) "TDEF is",tdef
763 endif
764 endif
765
766 timestamp = seconds+100*minutes+10000*hours
767
768 if((year > 1800) .and. (year < 2000)) year = year-1900
769 if((year >= 2000)) year = year-2000
770
771 if(month >= 2) day = day+31 ! add january
772 if(month >= 3) day = day+28 ! add february
773 if(month >= 4) day = day+31 ! add march
774 if(month >= 5) day = day+30 ! add april
775 if(month >= 6) day = day+31 ! add may
776 if(month >= 7) day = day+30 ! add june
777 if(month >= 8) day = day+31 ! add july
778 if(month >= 9) day = day+31 ! add august
779 if(month >= 10) day = day+30 ! add september
780 if(month >= 11) day = day+31 ! add october
781 if(month >= 12) day = day+30 ! add november
782 if((month > 2) .and. (mod(year,4) == 0)) day = day+1 ! get leap year day
783
784 datestamp = (year)*1000 + day
785 ! datestamp = (year+2100)*1000 + day
786
787 if(debug) then
788 write(6,*) ' time, timestamp, datestamp ',time(1:19),timestamp,datestamp
789 endif
790
791 end subroutine time_calc
792
793 !-------------------------------------------------------------------------
794
795 subroutine g_output_3d (file, file_time_index, var, length_var, &
796 nx, ny, nz, data_out, debug)
797 implicit none
798
799 character (len=*), intent(in) :: file
800 integer, intent(in) :: file_time_index
801 character (len=*), intent(in) :: var
802 integer, intent(in) :: length_var
803 integer , intent(in) :: nx, ny, nz
804 real, intent(out), dimension(:,:,:) :: data_out
805 logical, intent(in) :: debug
806 real, allocatable, dimension(:,:,:) :: data_tmp, data_tmp2
807 real, allocatable, dimension(:,:,:) :: u, v
808 real, allocatable, dimension(:,:) :: xlat, xlon
809 ! real, allocatable, dimension(:,:,:) :: z
810 real, allocatable, dimension(:,:,:) :: ph, phb
811 real, allocatable, dimension(:,:,:) :: p, pb
812 real, allocatable, dimension(:,:,:) :: t, qv
813 integer :: map_proj
814 real :: cen_lon, truelat1, truelat2
815
816
817 REAL , PARAMETER :: g = 9.81 ! acceleration due to gravity (m {s}^-2)
818 REAL , PARAMETER :: r_d = 287.
819 REAL , PARAMETER :: r_v = 461.6
820 REAL , PARAMETER :: cp = 7.*r_d/2.
821 REAL , PARAMETER :: cv = cp-r_d
822 REAL , PARAMETER :: cliq = 4190.
823 REAL , PARAMETER :: cice = 2106.
824 REAL , PARAMETER :: psat = 610.78
825 REAL , PARAMETER :: rcv = r_d/cv
826 REAL , PARAMETER :: rcp = r_d/cp
827 REAL , PARAMETER :: c2 = cp * rcv
828 REAL , PARAMETER :: T0 = 273.16
829
830 REAL , PARAMETER :: p1000mb = 100000.
831 REAL , PARAMETER :: cpovcv = cp/(cp-r_d)
832 REAL , PARAMETER :: cvovcp = 1./cpovcv
833 ! REAL :: pp
834
835 if(debug) then
836 write(6,*) ' calculations for variable ',var
837 end if
838
839 if(var == 'U' ) then
840
841 allocate ( data_tmp(nx+1,ny,nz) )
842 call da_get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, &
843 file_time_index, debug )
844 data_out = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:))
845
846 deallocate ( data_tmp )
847
848 else if(var == 'V' ) then
849
850 allocate ( data_tmp(nx,ny+1,nz) )
851 call da_get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, &
852 file_time_index, debug )
853 data_out = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:))
854 deallocate ( data_tmp )
855
856 else if(var == 'UMET' ) then
857
858 call da_get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug )
859
860 IF ( map_proj == 1 .OR. map_proj == 2 ) THEN
861
862 allocate ( u(nx,ny,nz) )
863 allocate ( v(nx,ny,nz) )
864 allocate ( xlat(nx,ny) )
865 allocate ( xlon(nx,ny) )
866
867 allocate ( data_tmp(nx+1,ny,nz) )
868 call da_get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, &
869 file_time_index, debug )
870 u = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:))
871 deallocate ( data_tmp )
872
873 allocate ( data_tmp(nx,ny+1,nz) )
874 call da_get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, &
875 file_time_index, debug )
876 v = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:))
877 deallocate ( data_tmp )
878
879 call da_get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug )
880 call da_get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug )
881 call da_get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug )
882 call da_get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug )
883 call da_get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug )
884
885 call rotate_wind (u,v,nx,ny,nz,var, &
886 map_proj,cen_lon,xlat,xlon, &
887 truelat1,truelat2,data_out)
888
889 deallocate ( xlat )
890 deallocate ( xlon )
891 deallocate ( u )
892 deallocate ( v )
893
894 ELSE
895
896 allocate ( data_tmp(nx+1,ny,nz) )
897 call da_get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, &
898 file_time_index, debug )
899 data_out = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:))
900 deallocate ( data_tmp )
901
902 ENDIF
903
904 else if(var == 'VMET' ) then
905
906 call da_get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug )
907
908 IF ( map_proj == 1 .OR. map_proj == 2 ) THEN
909
910 allocate ( u(nx,ny,nz) )
911 allocate ( v(nx,ny,nz) )
912 allocate ( xlat(nx,ny) )
913 allocate ( xlon(nx,ny) )
914
915 allocate ( data_tmp(nx+1,ny,nz) )
916 call da_get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, &
917 file_time_index, debug )
918 u = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:))
919 deallocate ( data_tmp )
920
921 allocate ( data_tmp(nx,ny+1,nz) )
922 call da_get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, &
923 file_time_index, debug )
924 v = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:))
925 deallocate ( data_tmp )
926
927 call da_get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug )
928 call da_get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug )
929 call da_get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug )
930 call da_get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug )
931 call da_get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug )
932
933 call rotate_wind (u,v,nx,ny,nz,var, &
934 map_proj,cen_lon,xlat,xlon, &
935 truelat1,truelat2,data_out)
936
937 deallocate ( xlat )
938 deallocate ( xlon )
939 deallocate ( u )
940 deallocate ( v )
941
942 ELSE
943
944 allocate ( data_tmp(nx,ny+1,nz) )
945 call da_get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, &
946 file_time_index, debug )
947 data_out = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:))
948 deallocate ( data_tmp )
949
950 ENDIF
951
952 else if(var == 'W' ) then
953
954 allocate ( data_tmp(nx,ny,nz+1) )
955 call da_get_var_3d_real_cdf( file,"W", data_tmp, nx, ny, nz+1, &
956 file_time_index, debug )
957 data_out = 0.5*(data_tmp(:,:,1:nz)+data_tmp(:,:,2:nz+1))
958 deallocate ( data_tmp )
959
960 else if(var == 'P' ) then
961
962 allocate ( p(nx,ny,nz) )
963 allocate ( pb(nx,ny,nz) )
964
965 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
966 file_time_index, debug )
967 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
968 file_time_index, debug )
969 data_out = (p+pb)*.01
970
971 deallocate ( p )
972 deallocate ( pb )
973
974 else if(var == 'Z' ) then
975
976 allocate ( ph(nx,ny,nz+1) )
977 allocate ( phb(nx,ny,nz+1) )
978
979 call da_get_var_3d_real_cdf( file,"PH", ph, nx, ny, nz+1, &
980 file_time_index, debug )
981 call da_get_var_3d_real_cdf( file,"PHB", phb, nx, ny, nz+1, &
982 file_time_index, debug )
983 ph = (ph+phb)/9.81
984 data_out = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1))
985
986 deallocate ( ph )
987 deallocate ( phb )
988
989 else if(var == 'THETA' ) then
990
991 call da_get_var_3d_real_cdf( file,"T", data_out, nx, ny, nz, &
992 file_time_index, debug )
993 data_out = data_out + 300.
994
995 else if(var == 'T' ) then
996
997 allocate ( p(nx,ny,nz) )
998 allocate ( pb(nx,ny,nz) )
999 allocate ( data_tmp(nx,ny,nz) )
1000
1001 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
1002 file_time_index, debug )
1003 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
1004 file_time_index, debug )
1005 p = p+pb
1006
1007 call da_get_var_3d_real_cdf( file,"T", data_tmp, nx, ny, nz, &
1008 file_time_index, debug )
1009 data_out = (data_tmp+300.)*(p/p1000mb)**rcp
1010
1011 deallocate ( p )
1012 deallocate ( pb )
1013 deallocate ( data_tmp )
1014
1015 else if(var == 'TC' ) then
1016
1017 allocate ( p(nx,ny,nz) )
1018 allocate ( pb(nx,ny,nz) )
1019 allocate ( data_tmp(nx,ny,nz) )
1020
1021 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
1022 file_time_index, debug )
1023 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
1024 file_time_index, debug )
1025 p = p+pb
1026
1027 call da_get_var_3d_real_cdf( file,"T", data_tmp, nx, ny, nz, &
1028 file_time_index, debug )
1029 data_out = (data_tmp+300.)*(p/p1000mb)**rcp -T0
1030
1031 deallocate ( p )
1032 deallocate ( pb )
1033 deallocate ( data_tmp )
1034
1035 else if(var == 'TD' ) then
1036
1037 allocate ( p(nx,ny,nz) )
1038 allocate ( pb(nx,ny,nz) )
1039 allocate ( qv(nx,ny,nz) )
1040 allocate ( data_tmp(nx,ny,nz) )
1041
1042 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
1043 file_time_index, debug )
1044 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
1045 file_time_index, debug )
1046 p = p+pb
1047
1048 call da_get_var_3d_real_cdf( file,"QVAPOR", qv, nx, ny, nz, &
1049 file_time_index, debug )
1050
1051 data_tmp = qv*(p/100.)/(0.622+qv)
1052 data_tmp = AMAX1(data_tmp,0.001)
1053 data_out = (243.5*log(data_tmp)-440.8)/(19.48-log(data_tmp))
1054
1055 deallocate ( p )
1056 deallocate ( pb )
1057 deallocate ( qv )
1058 deallocate ( data_tmp )
1059
1060 else if(var == 'RH' ) then
1061
1062 allocate ( p(nx,ny,nz) )
1063 allocate ( pb(nx,ny,nz) )
1064 allocate ( qv(nx,ny,nz) )
1065 allocate ( t(nx,ny,nz) )
1066 allocate ( data_tmp(nx,ny,nz) )
1067 allocate ( data_tmp2(nx,ny,nz) )
1068
1069 call da_get_var_3d_real_cdf( file,"P", p, nx, ny, nz, &
1070 file_time_index, debug )
1071 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, &
1072 file_time_index, debug )
1073 p = p+pb
1074
1075 call da_get_var_3d_real_cdf( file,"T", t, nx, ny, nz, &
1076 file_time_index, debug )
1077 call da_get_var_3d_real_cdf( file,"QVAPOR", qv, nx, ny, nz, &
1078 file_time_index, debug )
1079
1080 t = (t+300.)*(p/p1000mb)**rcp
1081 data_tmp2 = 10.*0.6112*exp(17.67*(t-T0)/(t-29.65))
1082 data_tmp = 0.622*data_tmp2/(0.01 * p - (1.-0.622)*data_tmp2)
1083 data_out = 100.*AMAX1(AMIN1(qv/data_tmp,1.0),0.0)
1084
1085
1086 deallocate ( p )
1087 deallocate ( pb )
1088 deallocate ( qv )
1089 deallocate ( t )
1090 deallocate ( data_tmp )
1091 deallocate ( data_tmp2 )
1092
1093 else
1094 call da_get_var_3d_real_cdf( file,var(1:length_var), &
1095 data_out, nx,ny,nz, &
1096 file_time_index, debug )
1097 endif
1098
1099
1100 end subroutine g_output_3d
1101
1102 !-------------------------------------------------------------------------
1103
1104 subroutine g_output_2d (file, file_time_index, var, length_var, &
1105 nx, ny, nz, data_out, debug)
1106 implicit none
1107
1108 character (len=*), intent(in) :: file
1109 integer, intent(in) :: file_time_index
1110 character (len=*), intent(in) :: var
1111 integer, intent(in) :: length_var
1112 integer, intent(in) :: nx, ny, nz
1113 real, intent(out), dimension(:,:,:) :: data_out
1114 logical, intent(in) :: debug
1115 integer, allocatable, dimension(:,:,:) :: data_int
1116 real, allocatable, dimension(:,:,:) :: u10, v10
1117 real, allocatable, dimension(:,:) :: psfc,t2m,q2m,mu
1118 real, allocatable, dimension(:,:) :: xlat, xlon
1119 real, allocatable, dimension(:,:,:) :: z,ph,phb
1120 real, allocatable, dimension(:,:,:) :: p,pb
1121 real, allocatable, dimension(:,:,:) :: ts,qv
1122 integer :: map_proj
1123 real :: cen_lon, truelat1, truelat2
1124
1125 if(debug) then
1126 write(6,*) ' calculations for variable ',var
1127 end if
1128
1129 if(var == 'SLP') then
1130
1131 allocate ( z(nx,ny,nz) )
1132 allocate ( ph(nx,ny,nz+1) )
1133 allocate ( phb(nx,ny,nz+1) )
1134 allocate ( p(nx,ny,nz) )
1135 allocate ( pb(nx,ny,nz) )
1136 allocate ( ts(nx,ny,nz) )
1137 allocate ( qv(nx,ny,nz) )
1138
1139 call da_get_var_3d_real_cdf( file,"PH", ph, nx, ny,nz+1, &
1140 file_time_index, debug )
1141 call da_get_var_3d_real_cdf( file,"PHB", phb, nx, ny,nz+1, &
1142 file_time_index, debug )
1143 ph = (ph+phb)/9.81
1144 z = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1))
1145
1146 call da_get_var_3d_real_cdf( file,"P", p, nx, ny,nz, &
1147 file_time_index, debug )
1148 call da_get_var_3d_real_cdf( file,"PB", pb, nx, ny,nz, &
1149 file_time_index, debug )
1150 p = p+pb
1151
1152 call da_get_var_3d_real_cdf( file,"T", ts, nx, ny,nz, &
1153 file_time_index, debug )
1154 call da_get_var_3d_real_cdf( file,"QVAPOR", qv, nx, ny,nz, &
1155 file_time_index, debug )
1156
1157 call compute_seaprs (nx, ny, nz, z, ts, p, qv, data_out, debug)
1158
1159
1160 deallocate ( z )
1161 deallocate ( ph )
1162 deallocate ( phb )
1163 deallocate ( p )
1164 deallocate ( pb )
1165 deallocate ( ts )
1166 deallocate ( qv )
1167
1168 else if(var == 'MU' ) then
1169 allocate ( mu(nx,ny) )
1170 call da_get_var_2d_real_cdf( file,"MU", mu, nx, ny, &
1171 file_time_index, debug )
1172 data_out(:,:,1) = mu(:,:)
1173 deallocate ( mu )
1174
1175 else if(var == 'PSFC' ) then
1176 allocate ( psfc(nx,ny) )
1177 call da_get_var_2d_real_cdf( file,"PSFC", psfc, nx, ny, &
1178 file_time_index, debug )
1179 data_out(:,:,1) = psfc(:,:)
1180 deallocate ( psfc )
1181
1182 else if(var == 'T2M' ) then
1183 allocate ( t2m(nx,ny) )
1184 call da_get_var_2d_real_cdf( file,"T2", t2m, nx, ny, &
1185 file_time_index, debug )
1186 data_out(:,:,1) = t2m(:,:)
1187 deallocate ( t2m )
1188
1189 else if(var == 'Q2M' ) then
1190 allocate ( q2m(nx,ny) )
1191 call da_get_var_2d_real_cdf( file,"Q2", q2m, nx, ny, &
1192 file_time_index, debug )
1193 data_out(:,:,1) = q2m(:,:)
1194 deallocate ( q2m )
1195
1196 else if(var == 'U10M' ) then
1197
1198 call da_get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug )
1199
1200 IF ( map_proj == 1 .OR. map_proj == 2 ) THEN
1201
1202 allocate ( u10(nx,ny,1) )
1203 allocate ( v10(nx,ny,1) )
1204 allocate ( xlat(nx, ny) )
1205 allocate ( xlon(nx, ny) )
1206 call da_get_var_2d_real_cdf( file,"U10", u10, nx, ny, &
1207 file_time_index, debug )
1208 call da_get_var_2d_real_cdf( file,"V10", v10, nx, ny, &
1209 file_time_index, debug )
1210
1211 call da_get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug )
1212 call da_get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug )
1213 call da_get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug )
1214 call da_get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug )
1215 call da_get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug )
1216
1217 call rotate_wind (u10,v10,nx,ny,1,var, &
1218 map_proj,cen_lon,xlat,xlon, &
1219 truelat1,truelat2,data_out)
1220
1221 deallocate ( xlat )
1222 deallocate ( xlon )
1223 deallocate ( u10 )
1224 deallocate ( v10 )
1225
1226 ELSE
1227
1228 call da_get_var_2d_real_cdf( file,"U10", data_out, nx, ny, &
1229 file_time_index, debug )
1230
1231 ENDIF
1232
1233 else if(var == 'V10M' ) then
1234
1235 call da_get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug )
1236
1237 IF ( map_proj == 1 .OR. map_proj == 2 ) THEN
1238
1239 allocate ( u10(nx,ny,1) )
1240 allocate ( v10(nx,ny,1) )
1241 allocate ( xlat(nx, ny) )
1242 allocate ( xlon(nx, ny) )
1243 call da_get_var_2d_real_cdf( file,"U10", u10, nx, ny, &
1244 file_time_index, debug )
1245 call da_get_var_2d_real_cdf( file,"V10", v10, nx, ny, &
1246 file_time_index, debug )
1247
1248 call da_get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug )
1249 call da_get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug )
1250 call da_get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug )
1251 call da_get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug )
1252 call da_get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug )
1253
1254 call rotate_wind (u10,v10,nx,ny,1,var, &
1255 map_proj,cen_lon,xlat,xlon, &
1256 truelat1,truelat2,data_out)
1257
1258 deallocate ( xlat )
1259 deallocate ( xlon )
1260 deallocate ( u10 )
1261 deallocate ( v10 )
1262
1263 ELSE
1264
1265 call da_get_var_2d_real_cdf( file,"V10", data_out, nx, ny, &
1266 file_time_index, debug )
1267
1268 ENDIF
1269
1270 else if(var == 'XLONG' ) then
1271 call da_get_var_2d_real_cdf( file,var(1:length_var), &
1272 data_out, nx,ny, &
1273 file_time_index, debug )
1274 WHERE ( data_out < 0.0 )
1275 data_out = data_out + 360.0
1276 ENDWHERE
1277
1278 else if(var == 'IVGTYP' .or. var == 'ISLTYP') then
1279
1280 allocate (data_int(nx,ny,1))
1281 call da_get_var_2d_int_cdf( file,var(1:length_var), &
1282 data_int, nx,ny, &
1283 file_time_index, debug )
1284 data_out = data_int
1285 deallocate (data_int)
1286
1287 else
1288 call da_get_var_2d_real_cdf( file,var(1:length_var), &
1289 data_out, nx,ny, &
1290 file_time_index, debug )
1291 endif
1292
1293
1294 end subroutine g_output_2d
1295
1296 !------------------------------------------------------------------
1297
1298 subroutine interp_to_z( data_in , nx_in , ny_in , nz_in , &
1299 data_out, nx_out, ny_out, nz_out, &
1300 z_in, z_out, missing_value, &
1301 vertical_type, debug )
1302 implicit none
1303 integer, intent(in) :: nx_in , ny_in , nz_in
1304 integer, intent(in) :: nx_out, ny_out, nz_out
1305 real, intent(in) :: missing_value
1306 real, dimension(nx_in , ny_in , nz_in ), intent(in ) :: data_in, z_in
1307 real, dimension(nx_out, ny_out, nz_out), intent(out) :: data_out
1308 real, dimension(nz_out), intent(in) :: z_out
1309 logical, intent(in) :: debug
1310 character (len=1) , intent(in) :: vertical_type
1311
1312 real, dimension(nz_in) :: data_in_z, zz_in
1313 real, dimension(nz_out) :: data_out_z
1314
1315 integer :: i,j,k
1316
1317 do i=1,nx_in
1318 do j=1,ny_in
1319
1320 do k=1,nz_in
1321 data_in_z(k) = data_in(i,j,k)
1322 zz_in(k) = z_in(i,j,k)
1323 enddo
1324
1325 !Hui do k=1,nz_out
1326 !Hui data_out_z(k) = data_out(i,j,k)
1327 !Hui enddo
1328
1329 call interp_1d( data_in_z, zz_in, nz_in, &
1330 data_out_z, z_out, nz_out, &
1331 vertical_type, missing_value )
1332
1333 do k=1,nz_out
1334 data_out(i,j,k) = data_out_z(k)
1335 enddo
1336
1337
1338 enddo
1339 enddo
1340
1341 end subroutine interp_to_z
1342
1343 !----------------------------------------------
1344
1345 subroutine interp_1d( a, xa, na, &
1346 b, xb, nb, vertical_type, missing_value )
1347 implicit none
1348 integer, intent(in) :: na, nb
1349 real, intent(in), dimension(na) :: a, xa
1350 real, intent(in), dimension(nb) :: xb
1351 real, intent(out), dimension(nb) :: b
1352 real, intent(in) :: missing_value
1353
1354 integer :: n_in, n_out
1355 logical :: interp
1356 real :: w1, w2
1357 character (len=1) ,intent(in) :: vertical_type
1358
1359
1360 if ( vertical_type == 'p' ) then
1361
1362 do n_out = 1, nb
1363
1364 b(n_out) = missing_value
1365 interp = .false.
1366 n_in = 1
1367
1368 do while ( (.not.interp) .and. (n_in < na) )
1369
1370 if( (xa(n_in) >= xb(n_out)) .and. &
1371 (xa(n_in+1) <= xb(n_out)) ) then
1372 interp = .true.
1373 w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
1374 w2 = 1. - w1
1375 b(n_out) = w1*a(n_in) + w2*a(n_in+1)
1376 end if
1377 n_in = n_in +1
1378
1379 enddo
1380
1381 enddo
1382
1383 else
1384
1385 do n_out = 1, nb
1386
1387 b(n_out) = missing_value
1388 interp = .false.
1389 n_in = 1
1390
1391 do while ( (.not.interp) .and. (n_in < na) )
1392
1393 if( (xa(n_in) <= xb(n_out)) .and. &
1394 (xa(n_in+1) >= xb(n_out)) ) then
1395 interp = .true.
1396 w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
1397 w2 = 1. - w1
1398 b(n_out) = w1*a(n_in) + w2*a(n_in+1)
1399 end if
1400 n_in = n_in +1
1401
1402 enddo
1403
1404 enddo
1405
1406 endif
1407
1408 end subroutine interp_1d
1409
1410 !-------------------------------------------------------------------------
1411 !
1412 ! This routines has been taken "as is" from wrf_user_fortran_util_0.f
1413 !
1414 ! This routine assumes
1415 ! index order is (i,j,k)
1416 ! wrf staggering
1417 ! units: pressure (Pa), temperature(K), height (m), mixing ratio (kg kg{-1})
1418 ! availability of 3d p, t, and qv; 2d terrain; 1d half-level zeta string
1419 ! output units of SLP are Pa, but you should divide that by 100 for the
1420 ! weather weenies.
1421 ! virtual effects are included
1422 !
1423 ! Dave
1424
1425 subroutine compute_seaprs ( nx , ny , nz , &
1426 z, t , p , q , &
1427 sea_level_pressure,debug)
1428 ! & t_sea_level, t_surf, level )
1429 IMPLICIT NONE
1430 ! Estimate sea level pressure.
1431 INTEGER, intent(in) :: nx , ny , nz
1432 REAL, intent(in) :: z(nx,ny,nz)
1433 REAL, intent(in) :: p(nx,ny,nz) , q(nx,ny,nz)
1434 REAL, intent(inout) :: t(nx,ny,nz)
1435 ! The output is the 2d sea level pressure.
1436 REAL, intent(out) :: sea_level_pressure(nx,ny)
1437 INTEGER level(nx,ny)
1438 REAL t_surf(nx,ny) , t_sea_level(nx,ny)
1439 LOGICAL, intent(in) :: debug
1440
1441 ! Some required physical constants:
1442
1443 REAL R, G, GAMMA
1444 PARAMETER (R=287.04, G=9.81, GAMMA=0.0065)
1445
1446 ! Specific constants for assumptions made in this routine:
1447
1448 REAL TC, PCONST
1449 PARAMETER (TC=273.16+17.5, PCONST = 10000)
1450 LOGICAL ridiculous_mm5_test
1451 PARAMETER (ridiculous_mm5_test = .TRUE.)
1452 ! PARAMETER (ridiculous_mm5_test = .false.)
1453
1454 ! Local variables:
1455
1456 INTEGER i , j , k
1457 INTEGER klo , khi
1458
1459
1460 REAL plo , phi , tlo, thi , zlo , zhi
1461 REAL p_at_pconst , t_at_pconst , z_at_pconst
1462 REAL z_half_lowest
1463
1464 REAL , PARAMETER :: cp = 7.*R/2.
1465 REAL , PARAMETER :: rcp = R/cp
1466 REAL , PARAMETER :: p1000mb = 100000.
1467
1468 LOGICAL l1 , l2 , l3, found
1469
1470 ! Find least zeta level that is PCONST Pa above the surface. We later use this
1471 ! level to extrapolate a surface pressure and temperature, which is supposed
1472 ! to reduce the effect of the diurnal heating cycle in the pressure field.
1473
1474 t(:,:,:) = (t(:,:,:)+300.)*(p(:,:,:)/p1000mb)**rcp
1475
1476 DO j = 1 , ny
1477 DO i = 1 , nx
1478 level(i,j) = -1
1479
1480 k = 1
1481 found = .false.
1482 do while( (.not. found) .and. (k.le.nz))
1483 IF ( p(i,j,k) .LT. p(i,j,1)-PCONST ) THEN
1484 level(i,j) = k
1485 found = .true.
1486 END IF
1487 k = k+1
1488 END DO
1489
1490 IF ( level(i,j) .EQ. -1 ) THEN
1491 PRINT '(A,I4,A)','Troubles finding level ', &
1492 NINT(PCONST)/100,' above ground.'
1493 PRINT '(A,I4,A,I4,A)', &
1494 'Problems first occur at (',i,',',j,')'
1495 PRINT '(A,F6.1,A)', &
1496 'Surface pressure = ',p(i,j,1)/100,' hPa.'
1497 STOP 'Error_in_finding_100_hPa_up'
1498 END IF
1499
1500
1501 END DO
1502 END DO
1503
1504 ! Get temperature PCONST Pa above surface. Use this to extrapolate
1505 ! the temperature at the surface and down to sea level.
1506
1507 DO j = 1 , ny
1508 DO i = 1 , nx
1509
1510 klo = MAX ( level(i,j) - 1 , 1 )
1511 khi = MIN ( klo + 1 , nz - 1 )
1512
1513 IF ( klo .EQ. khi ) THEN
1514 PRINT '(A)','Trapping levels are weird.'
1515 PRINT '(A,I3,A,I3,A)','klo = ',klo,', khi = ',khi, &
1516 ': and they should not be equal.'
1517 STOP 'Error_trapping_levels'
1518 END IF
1519
1520 plo = p(i,j,klo)
1521 phi = p(i,j,khi)
1522 tlo = t(i,j,klo)*(1. + 0.608 * q(i,j,klo) )
1523 thi = t(i,j,khi)*(1. + 0.608 * q(i,j,khi) )
1524 ! zlo = zetahalf(klo)/ztop*(ztop-terrain(i,j))+terrain(i,j)
1525 ! zhi = zetahalf(khi)/ztop*(ztop-terrain(i,j))+terrain(i,j)
1526 zlo = z(i,j,klo)
1527 zhi = z(i,j,khi)
1528
1529 p_at_pconst = p(i,j,1) - pconst
1530 t_at_pconst = thi-(thi-tlo)*LOG(p_at_pconst/phi)*LOG(plo/phi)
1531 z_at_pconst = zhi-(zhi-zlo)*LOG(p_at_pconst/phi)*LOG(plo/phi)
1532
1533 t_surf(i,j) = t_at_pconst*(p(i,j,1)/p_at_pconst)**(gamma*R/g)
1534 t_sea_level(i,j) = t_at_pconst+gamma*z_at_pconst
1535
1536 END DO
1537 END DO
1538
1539 ! If we follow a traditional computation, there is a correction to the sea level
1540 ! temperature if both the surface and sea level temnperatures are *too* hot.
1541
1542 IF ( ridiculous_mm5_test ) THEN
1543 DO j = 1 , ny
1544 DO i = 1 , nx
1545 l1 = t_sea_level(i,j) .LT. TC
1546 l2 = t_surf (i,j) .LE. TC
1547 l3 = .NOT. l1
1548 IF ( l2 .AND. l3 ) THEN
1549 t_sea_level(i,j) = TC
1550 ELSE
1551 t_sea_level(i,j) = TC - 0.005*(t_surf(i,j)-TC)**2
1552 END IF
1553 END DO
1554 END DO
1555 END IF
1556
1557 ! The grand finale: ta da!
1558
1559 DO j = 1 , ny
1560 DO i = 1 , nx
1561 ! z_half_lowest=zetahalf(1)/ztop*(ztop-terrain(i,j))+terrain(i,j)
1562 z_half_lowest=z(i,j,1)
1563 sea_level_pressure(i,j) = p(i,j,1) * &
1564 EXP((2.*g*z_half_lowest)/ &
1565 (R*(t_sea_level(i,j)+t_surf(i,j))))
1566 END DO
1567 END DO
1568
1569 if (debug) then
1570 print *,'sea pres input at weird location i=20,j=1,k=1'
1571 print *,'t=',t(20,1,1),t(20,2,1),t(20,3,1)
1572 print *,'z=',z(20,1,1),z(20,2,1),z(20,3,1)
1573 print *,'p=',p(20,1,1),p(20,2,1),p(20,3,1)
1574 print *,'slp=',sea_level_pressure(20,1), &
1575 sea_level_pressure(20,2),sea_level_pressure(20,3)
1576 endif
1577 ! print *,'t=',t(10:15,10:15,1),t(10:15,2,1),t(10:15,3,1)
1578 ! print *,'z=',z(10:15,1,1),z(10:15,2,1),z(10:15,3,1)
1579 ! print *,'p=',p(10:15,1,1),p(10:15,2,1),p(10:15,3,1)
1580 ! print *,'slp=',sea_level_pressure(10:15,10:15), &
1581 ! sea_level_pressure(10:15,10:15),sea_level_pressure(20,10:15)
1582
1583 end subroutine compute_seaprs
1584
1585 !------------------------------------------------------------------
1586
1587
1588 subroutine rotate_wind (u,v,d1,d2,d3,var, &
1589 map_proj,cen_lon,xlat,xlon, &
1590 truelat1,truelat2,data_out)
1591
1592 implicit none
1593
1594 integer, intent(in) :: d1, d2, d3
1595
1596 real, dimension(d1,d2,d3), intent(out) :: data_out
1597 integer, intent(in) :: map_proj
1598 integer ::i,j,k
1599 real, intent(in) :: cen_lon, truelat1, truelat2
1600 real :: cone
1601 real, dimension(d1,d2,d3), intent(in) :: u,v
1602 real, dimension(d1,d2), intent(in) :: xlat, xlon
1603 real, dimension(d1,d2) :: diff, alpha
1604
1605 character (len=10), intent(in) :: var
1606
1607 REAL , PARAMETER :: pii = 3.14159265
1608 REAL , PARAMETER :: radians_per_degree = pii/180.
1609
1610
1611
1612
1613 cone = 1. ! PS
1614 if( map_proj .eq. 1) then ! Lambert Conformal mapping
1615 IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
1616 cone=(ALOG(COS(truelat1*radians_per_degree))- &
1617 ALOG(COS(truelat2*radians_per_degree))) / &
1618 (ALOG(TAN((90.-ABS(truelat1))*radians_per_degree*0.5 ))- &
1619 ALOG(TAN((90.-ABS(truelat2))*radians_per_degree*0.5 )) )
1620 ELSE
1621 cone = SIN(ABS(truelat1)*radians_per_degree )
1622 ENDIF
1623 end if
1624
1625
1626 diff = xlon - cen_lon
1627
1628 do i = 1, d1
1629 do j = 1, d2
1630 if(diff(i,j) .gt. 180.) then
1631 diff(i,j) = diff(i,j) - 360.
1632 end if
1633 if(diff(i,j) .lt. -180.) then
1634 diff(i,j) = diff(i,j) + 360.
1635 end if
1636 end do
1637 end do
1638
1639 do i = 1, d1
1640 do j = 1, d2
1641 if(xlat(i,j) .lt. 0.) then
1642 alpha(i,j) = - diff(i,j) * cone * radians_per_degree
1643 else
1644 alpha(i,j) = diff(i,j) * cone * radians_per_degree
1645 end if
1646 end do
1647 end do
1648
1649
1650 if(var(1:1) .eq. "U") then
1651 do k=1,d3
1652 data_out(:,:,k) = v(:,:,k)*sin(alpha) + u(:,:,k)*cos(alpha)
1653 end do
1654 else if(var(1:1) .eq. "V") then
1655 do k=1,d3
1656 data_out(:,:,k) = v(:,:,k)*cos(alpha) - u(:,:,k)*sin(alpha)
1657 end do
1658 end if
1659
1660
1661 end subroutine rotate_wind
1662
1663
1664 !------------------------------------------------------------------
1665 subroutine handle_err(rmarker,nf_status)
1666
1667 #include "netcdf.inc"
1668 integer, intent(in) :: nf_status
1669 character*(*), intent(in) :: rmarker
1670 if (nf_status .ne. nf_noerr) then
1671 write(*,*) 'NetCDF error : ',rmarker
1672 write(*,*) ' ',nf_strerror(nf_status)
1673 stop
1674 endif
1675
1676 end subroutine handle_err
1677
1678 !------------------------------------------------------------------
1679 subroutine get_dimensions(infile,nx,ny,nz)
1680
1681 #include "netcdf.inc"
1682 character (len=512), intent(in) :: infile
1683 integer :: ncid, dimid, nf_status
1684 integer, intent(inout) :: nx, ny, nz
1685 integer :: nlgen
1686
1687 ! need to pull out some data to set up dimensions, etc.
1688 nf_status = nf_open (infile, nf_nowrite, ncid)
1689 call handle_err('Error opening file: '//trim(infile),nf_status)
1690 !
1691 nf_status = nf_inq_dimid (ncid, 'west_east_stag', dimid)
1692 call handle_err('west_east_stag',nf_status)
1693 nf_status = nf_inq_dimlen (ncid, dimid, nx)
1694 call handle_err('Get NX',nf_status)
1695 nx = nx-1
1696 !
1697 nf_status = nf_inq_dimid (ncid, 'south_north_stag', dimid)
1698 call handle_err('south_north_stag',nf_status)
1699 nf_status = nf_inq_dimlen (ncid, dimid, ny)
1700 call handle_err('Get NY',nf_status)
1701 ny = ny-1
1702 !
1703 nf_status = nf_inq_dimid (ncid, 'bottom_top', dimid)
1704 call handle_err('bottom_top',nf_status)
1705 nf_status = nf_inq_dimlen (ncid, dimid, nz)
1706 call handle_err('Get NZ',nf_status)
1707 nlgen = nz
1708 !
1709 nf_status = nf_close (ncid)
1710 call handle_err('Error closing file',nf_status)
1711
1712 end subroutine get_dimensions
1713 !------------------------------------------------------------------
1714
1715 subroutine get_diffs(var1, var2, diff, absdiff, sqdiff, nx, ny, nz, missing)
1716
1717 real, intent(in), dimension(:,:,:) :: var1, var2
1718 real, intent(out), dimension(:,:,:) :: diff, absdiff, sqdiff
1719 integer, intent(in) :: nx, ny, nz
1720 integer :: i,j,k
1721 real, intent(in) :: missing
1722
1723 do k = 1, nz
1724 do j = 1, ny
1725 do i = 1, nx
1726 if ( var1(i,j,k) /= missing .and. var2(i,j,k) /= missing ) then
1727 diff(i,j,k) = var2(i,j,k) - var1(i,j,k)
1728 absdiff(i,j,k) = abs(var2(i,j,k) - var1(i,j,k))
1729 sqdiff(i,j,k) = (var2(i,j,k) - var1(i,j,k) )*(var2(i,j,k) - var1(i,j,k) )
1730 else
1731 diff(i,j,k) = missing
1732 absdiff(i,j,k) = missing
1733 sqdiff(i,j,k) = missing
1734 endif
1735 enddo
1736 enddo
1737 enddo
1738
1739 end subroutine get_diffs
1740
1741 !------------------------------------------------------------------
1742 subroutine get_diffs_wind(u1, v1, u2,v2,diff, absdiff, sqdiff, nx, ny, nz, missing)
1743
1744 real, intent(in), dimension(:,:,:) :: u1, u2, v1, v2
1745 real, intent(out), dimension(:,:,:) :: diff, absdiff, sqdiff
1746 integer, intent(in) :: nx, ny, nz
1747 integer :: i,j,k
1748 real, intent(in) :: missing
1749 real :: u, v
1750
1751 do k = 1, nz
1752 do j = 1, ny
1753 do i = 1, nx
1754 if ( u1(i,j,k) /= missing .and. v1(i,j,k) /= missing .and. &
1755 u2(i,j,k) /= missing .and. v2(i,j,k) /= missing ) then
1756 u = u1(i,j,k) - u2(i,j,k)
1757 v = v1(i,j,k) - v2(i,j,k)
1758 diff(i,j,k) = u + v
1759 absdiff(i,j,k) = abs(u) + abs (v)
1760 sqdiff(i,j,k) = u*u + v*v
1761 else
1762 diff(i,j,k) = missing
1763 absdiff(i,j,k) = missing
1764 sqdiff(i,j,k) = missing
1765 endif
1766 enddo
1767 enddo
1768 enddo
1769
1770 end subroutine get_diffs_wind
1771
1772 !------------------------------------------------------------------
1773 subroutine domain_average(var, avg_prof, counter, nx, ny, nz, missing,isq)
1774
1775 integer, intent(in) :: nx,ny,nz,isq
1776 real, intent(in), dimension(:,:,:) :: var
1777 integer, intent(out), dimension(:) :: counter
1778 real, intent(out), dimension(:) :: avg_prof
1779 real, intent(in) :: missing
1780
1781 integer :: i,j,k,icount
1782 real :: dsum, dmiss
1783 integer :: imiss
1784
1785 !9999
1786 !Hui: set dmiss value consistent with plot script
1787 ! dmiss = -9999.999
1788 dmiss = -99.99
1789 imiss = -99
1790 do k = 1, nz
1791 icount = 0
1792 dsum = 0
1793 do j = 1, ny
1794 do i = 1, nx
1795 if ( var(i,j,k) /= missing ) then
1796 icount = icount + 1
1797 dsum = dsum + var(i,j,k)
1798 endif
1799 enddo
1800 enddo
1801 avg_prof(k) = dmiss
1802 !Hui counter(k) = 0
1803 counter(k) = imiss
1804 if (icount /= 0 ) then
1805 counter(k) = icount
1806 if ( isq .eq. 0 ) then
1807 avg_prof(k) = dsum /float(icount)
1808 else
1809 avg_prof(k) = sqrt(dsum /float(icount))
1810 endif
1811 endif
1812 enddo
1813
1814 end subroutine domain_average
1815 !------------------------------------------------------------------
1816
1817 subroutine mask_domain_average(var, avg_prof, counter, nx, ny, nz, missing,isq)
1818
1819 integer, intent(in) :: nx,ny,nz,isq
1820 real, intent(in), dimension(nx,ny,nz) :: var
1821 integer, intent(out), dimension(3,nz) :: counter
1822 real, intent(out), dimension(3,nz) :: avg_prof
1823 real, intent(in) :: missing
1824
1825 integer :: i,j,k,imask
1826 integer :: icount(3)
1827 real :: dsum(3)
1828 real :: dmiss
1829 integer :: imiss
1830
1831 ! 1: all, 2:land, 3:water
1832
1833 dmiss = -99.99
1834 imiss = -99
1835 do k = 1, nz
1836 icount = 0
1837 dsum = 0.0
1838 do j = 1, ny
1839 do i = 1, nx
1840 if ( var(i,j,k) /= missing ) then
1841 if ( i>=istart .and. i<=iend .and. &
1842 j>=jstart .and. j<=jend ) then
1843 icount(1) = icount(1) + 1
1844 dsum(1) = dsum(1) + var(i,j,k)
1845 if ( int(landmask(i,j)) == island ) then
1846 icount(2) = icount(2) + 1
1847 dsum(2) = dsum(2) + var(i,j,k)
1848 else if ( int(landmask(i,j)) == iswater ) then
1849 icount(3) = icount(3) + 1
1850 dsum(3) = dsum(3) + var(i,j,k)
1851 end if
1852 end if
1853 end if
1854 end do
1855 end do
1856 avg_prof(:,k) = dmiss
1857 counter(:,k) = imiss
1858 do imask = 1, 3
1859 if ( icount(imask) /= 0 ) then
1860 counter(imask,k) = icount(imask)
1861 if ( isq .eq. 0 ) then
1862 avg_prof(imask,k) = dsum(imask)/float(icount(imask))
1863 else
1864 avg_prof(imask,k) = sqrt(dsum(imask)/float(icount(imask)))
1865 end if
1866 end if
1867 end do
1868 enddo
1869
1870 end subroutine mask_domain_average
1871 !------------------------------------------------------------------
1872
1873 subroutine get_sum(dsum, dvar, nx, ny, nz, missing)
1874
1875 integer, intent(in) :: nx, ny, nz
1876 real, intent(in) :: missing
1877 real, intent(in),dimension(:,:,:) :: dvar
1878 real, intent(inout),dimension(:,:,:) :: dsum
1879
1880 integer :: i,j,k
1881
1882 do k = 1, nz
1883 do j = 1, ny
1884 do i = 1, nx
1885 if ( dvar(i,j,k) /= missing .and. dsum(i,j,k) /= missing ) then
1886 dsum(i,j,k) = dsum(i,j,k) + dvar(i,j,k)
1887 else
1888 dsum(i,j,k) = missing
1889 endif
1890 enddo
1891 enddo
1892 enddo
1893
1894 end subroutine get_sum
1895 !------------------------------------------------------------------
1896 subroutine time_average(dvar, nx, ny, nz, time_count,missing, isqr)
1897
1898 integer, intent(in) :: nx, ny, nz,time_count,isqr
1899 real, intent(in) :: missing
1900 real, intent(inout), dimension(:,:,:) :: dvar
1901
1902 integer :: i,j,k
1903
1904 do k = 1, nz
1905 do j = 1, ny
1906 do i = 1, nx
1907 if ( dvar(i,j,k) /= missing ) then
1908 if (isqr .eq. 1 ) then
1909 dvar(i,j,k) = sqrt(dvar(i,j,k)/float(time_count))
1910 else
1911 dvar(i,j,k) = dvar(i,j,k)/float(time_count)
1912 endif
1913 else
1914 dvar(i,j,k) = missing
1915 endif
1916 enddo
1917 enddo
1918 enddo
1919
1920 end subroutine time_average
1921 !------------------------------------------------------------------
1922 subroutine compute_data_3d( infile, var, length, nx, ny, nz, levels, time_idx, &
1923 vert_levels, vertical_type, missing, data_out_z, debug )
1924
1925 integer, intent(in) :: time_idx
1926 integer, intent(in) :: nx, ny, nz, levels
1927 integer, intent(in) :: length
1928 real, intent(in) :: missing
1929 real, intent(in) :: vert_levels(:)
1930 character (len=1), intent(in) :: vertical_type
1931 character (len=*), intent(in) :: var
1932 character (len=*), intent(in) :: infile
1933 logical, intent(in) :: debug
1934
1935 real, intent(out), dimension(:,:,:) :: data_out_z
1936 real, allocatable, dimension(:,:,:) :: data_out
1937 real, allocatable, dimension(:,:,:) :: z, ph, phb
1938 real, allocatable, dimension(:,:,:) :: p, pb
1939
1940
1941 ! first, get some base-state-stuff
1942
1943 if ( vertical_type == 'p' ) then
1944 allocate( p (nx, ny, nz) )
1945 allocate( pb(nx, ny, nz) )
1946 call da_get_var_3d_real_cdf( infile,'PB',pb,nx,ny,nz,time_idx,debug )
1947 call da_get_var_3d_real_cdf( infile,'P',p, nx, ny, nz, time_idx,debug )
1948 p = (p+pb)/100.0 ! Convert to hPa
1949 endif
1950 if ( vertical_type == 'z' ) then
1951 allocate( z (nx, ny, nz) )
1952 allocate( ph (nx, ny, nz+1) )
1953 allocate( phb(nx, ny, nz+1) )
1954 call da_get_var_3d_real_cdf( infile,'PHB',phb,nx,ny,nz+1,time_idx,debug )
1955 call da_get_var_3d_real_cdf( infile,'PH',ph, nx, ny, nz+1, time_idx,debug )
1956 ph = (ph+phb)/9.81
1957 z = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1))
1958 z = z/1000. ! Convert to kilometers
1959 endif
1960
1961 allocate ( data_out (nx, ny, nz) )
1962
1963 call g_output_3d (infile, time_idx, var, length, nx, ny, nz, data_out, debug)
1964
1965 if ( vertical_type == 'p' ) then
1966 call interp_to_z( data_out, nx, ny, nz, data_out_z, nx, ny, levels, &
1967 p, vert_levels, missing, vertical_type, debug )
1968
1969 else if ( vertical_type == 'z' ) then
1970 call interp_to_z( data_out, nx, ny, nz, data_out_z, nx, ny, levels, &
1971 z, vert_levels, missing, vertical_type, debug )
1972 else
1973 data_out_z = data_out
1974 endif
1975 deallocate ( data_out )
1976 if ( vertical_type == 'p' ) then
1977 deallocate( p )
1978 deallocate( pb )
1979 endif
1980 if ( vertical_type == 'z' ) then
1981 deallocate( z )
1982 deallocate( ph )
1983 deallocate( phb )
1984 endif
1985
1986 end subroutine compute_data_3d
1987
1988 subroutine compute_wind_3d( infile, nx, ny, nz, levels, time_idx, &
1989 vert_levels, vertical_type, missing, data_out_z1, data_out_z2, debug )
1990
1991 integer, intent(in) :: time_idx
1992 integer, intent(in) :: nx, ny, nz, levels
1993 real, intent(in) :: missing
1994 real, intent(in) :: vert_levels(:)
1995 character (len=1), intent(in) :: vertical_type
1996 character (len=*), intent(in) :: infile
1997 logical, intent(in) :: debug
1998
1999 real, intent(out), dimension(:,:,:) :: data_out_z1
2000 real, intent(out), dimension(:,:,:) :: data_out_z2
2001 real, allocatable, dimension(:,:,:) :: data_out
2002 real, allocatable, dimension(:,:,:) :: z, ph, phb
2003 real, allocatable, dimension(:,:,:) :: p, pb
2004 character (len=1) :: var
2005
2006
2007 ! first, get some base-state-stuff
2008
2009 if ( vertical_type == 'p' ) then
2010 allocate( p (nx, ny, nz) )
2011 allocate( pb(nx, ny, nz) )
2012 call da_get_var_3d_real_cdf( infile,'PB',pb,nx,ny,nz,time_idx,debug )
2013 call da_get_var_3d_real_cdf( infile,'P',p, nx, ny, nz, time_idx,debug )
2014 p = (p+pb)/100.0 ! Convert to hPa
2015 endif
2016 if ( vertical_type == 'z' ) then
2017 allocate( z (nx, ny, nz) )
2018 allocate( ph (nx, ny, nz+1) )
2019 allocate( phb(nx, ny, nz+1) )
2020 call da_get_var_3d_real_cdf( infile,'PHB',phb,nx,ny,nz+1,time_idx,debug )
2021 call da_get_var_3d_real_cdf( infile,'PH',ph, nx, ny, nz+1, time_idx,debug )
2022 ph = (ph+phb)/9.81
2023 z = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1))
2024 z = z/1000. ! Convert to kilometers
2025 endif
2026
2027 allocate ( data_out2 (nx, ny, nz) )
2028 allocate ( data_out1 (nx, ny, nz) )
2029 var='U'
2030 call g_output_3d (infile, time_idx, var, 1, nx, ny, nz, data_out1, debug)
2031 var='V'
2032 call g_output_3d (infile, time_idx, var, 1, nx, ny, nz, data_out2, debug)
2033
2034 if ( vertical_type == 'p' ) then
2035 call interp_to_z( data_out1, nx, ny, nz, data_out_z1, nx, ny, levels, &
2036 p, vert_levels, missing, vertical_type, debug )
2037 call interp_to_z( data_out2, nx, ny, nz, data_out_z2, nx, ny, levels, &
2038 p, vert_levels, missing, vertical_type, debug )
2039
2040 else if ( vertical_type == 'z' ) then
2041 call interp_to_z( data_out1, nx, ny, nz, data_out_z1, nx, ny, levels, &
2042 z, vert_levels, missing, vertical_type, debug )
2043 call interp_to_z( data_out2, nx, ny, nz, data_out_z2, nx, ny, levels, &
2044 z, vert_levels, missing, vertical_type, debug )
2045 else
2046 data_out_z1 = data_out1
2047 data_out_z2 = data_out2
2048 endif
2049 deallocate ( data_out1 )
2050 deallocate ( data_out2 )
2051 if ( vertical_type == 'p' ) then
2052 deallocate( p )
2053 deallocate( pb )
2054 endif
2055 if ( vertical_type == 'z' ) then
2056 deallocate( z )
2057 deallocate( ph )
2058 deallocate( phb )
2059 endif
2060
2061 end subroutine compute_wind_3d
2062
2063 !---------------------------------------------------------------------
2064 end program da_verif_grid
Mirror of the official repository at https://github.com/wrf-model/WRF