da_plot_gen_be.f90
References to this file elsewhere.
1 program da_plot_gen_be
2 !--------------------------------------------------------------------------!
3 ! Purpose: 1) Compute ratio of chi_b, temp_b and ps_b explained by psi to !
4 ! its full fields, and plot the global and local figures of !
5 ! the ratio. !
6 ! 2) Plot the global eigenvalues, the first 5 global eigenvectors !
7 ! and the first 5 local eigenvalues for psi, chi_u, t_u, and !
8 ! rh. !
9 ! 3) Plot the scale-lengths for psi, chi_u, t_u, rh, and ps_u. !
10 ! !
11 ! Note: this program is working for gen_be.NMC.dat and gen_be.ENS.dat for !
12 ! for cv_options=5. !
13 ! !
14 ! History: !
15 ! 06/09/2005 Original developer. Yong-Run Guo !
16 ! 06/20/2005 Modified for plotting the BES file !
17 ! including bin information and ps_u !
18 ! covariance. Yong-Run Guo !
19 ! !
20 !--------------------------------------------------------------------------!
21
22 implicit none
23
24 type eigen_type
25 real(kind=8), pointer :: g_val(:)
26 real(kind=8), pointer :: g_vec(:,:)
27 real(kind=8), pointer :: l_val(:,:)
28 real(kind=8), pointer :: l_vec(:,:,:)
29 real(kind=8), pointer :: sl(:)
30 real(kind=8), pointer :: power(:,:)
31 end type eigen_type
32
33 type (eigen_type) :: cv
34
35 integer, parameter :: iunit = 10, input_unit = 11
36
37 !
38 integer :: ni, nj, nk, nkdum, num_bins, num_bins2d
39 integer :: bin_type ! Type of bin to average over. !!!DALE ADD.
40 real(kind=8) :: lat_min, lat_max ! Used if bin_type = 2 (degrees).
41 real(kind=8) :: binwidth_lat ! Used if bin_type = 2 (degrees). !!!DALE ADD..
42 real(kind=8) :: hgt_min, hgt_max ! Used if bin_type = 2 (m).
43 real(kind=8) :: binwidth_hgt ! Used if bin_type = 2 (m). !!!DALE ADD..
44 integer :: num_bins_lat ! Used if bin_type = 2.
45 integer :: num_bins_hgt ! Used if bin_type = 2.
46
47 logical :: ldum1, ldum2 ! Dummy logicals.
48
49 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
50 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
51 real(kind=8), allocatable :: regcoeff_chi(:) ! psi/chi regression cooefficient.
52 real(kind=8), allocatable :: regcoeff_ps (:,:) ! psi/ps regression cooefficient.
53 real(kind=8), allocatable :: regcoeff_t (:,:,:)! psi/T regression cooefficient.
54 real(kind=8), allocatable :: reg_chi(:,:)
55 real(kind=8), allocatable :: reg_ps (:,:)
56 real(kind=8), allocatable :: reg_t(:,:,:)
57
58 ! Namelist variables:
59
60 character(len= 10) :: start_date, end_date ! Starting and ending dates.
61 character(len=120) :: main_title1
62 character(len=120) :: main_title2
63 character(len= 3) :: be_method ! Be method ('NMC', or 'ENS')
64 character(len= 8) :: uh_method ! Uh_method (power, scale)
65 character(len= 80) :: gen_be_dir
66 character(len= 8) :: code_version
67 integer :: stride ! in s-N direction
68 integer :: interval ! file interval
69 integer :: ne ! number of ENS members
70 real :: Resolution_km
71 logical :: domain_averaged
72
73 namelist /plot_title/start_date, end_date, interval, &
74 main_title1, main_title2, &
75 Resolution_km, stride, ne, &
76 be_method, Uh_method, gen_be_dir, &
77 code_version, domain_averaged
78
79 ! working variables
80 integer :: i, j, k, n, m, b, member, n_times, ier
81 integer :: sdate, cdate, edate ! Starting, current ending dates.
82 character(len= 10) :: variable ! Variable name
83 character(len= 80) :: filename
84 character(len= 10) :: date, new_date ! Current date (ccyymmddhh).
85 character(len= 3) :: ce ! Member index -> characte
86 character(len= 5) :: cvv ! control variable name
87
88 ! Variable fields
89 real(kind=8), allocatable :: psi (:,:,:) ! psi.
90 real(kind=8), allocatable :: chi (:,:,:) ! chi
91 real(kind=8), allocatable :: chi_u (:,:,:) ! chi_u
92
93 real(kind=8), allocatable :: temp(:,:,:) ! Temperature.
94 real(kind=8), allocatable :: t_u (:,:,:) ! Temperature.
95
96 real(kind=8), allocatable :: ps (:,:) ! Surface pressure.
97 real(kind=8), allocatable :: ps_u(:,:) ! Surface pressure.
98
99 ! Arrays for Ratio computation:
100 real(kind=4), allocatable :: chi_loc (:,:), chi_global(:) ! chi_b covariance
101 real(kind=4), allocatable :: temp_loc(:,:), temp_global(:)! Temp_b covariance
102 real(kind=4), allocatable :: ps_loc (:) ! ps_b .covariance
103 real(kind=4), allocatable :: chi_var (:,:), chi_gvar(:) ! chi variance
104 real(kind=4), allocatable :: temp_var(:,:), temp_gvar(:) ! Temp variance
105 real(kind=4), allocatable :: ps_var (:) ! ps variance
106 real(kind=8) :: balance, unbalan, avg, avg2, avg3
107
108 integer :: max_wavenumber
109 real(kind=8), allocatable :: total_power(:,:) ! Total Power spectrum.
110 !
111 ! For contour plot:
112
113 integer :: IOFFP, ndot, jj, interval_j, len_title
114 real :: flo, hi, cint, spval0
115
116 character (len=120) :: Title
117 common /conre1/ IOFFP, SPVAL0
118 IOFFP = 1
119 SPVAL0 = -99999.9
120 !
121 ! 1.0 Initialize GKS
122 ! ------------------
123
124 call opngks
125 call da_setup_color_table
126
127 ! 2.0 Read namelist
128 ! -----------------
129
130 start_date = '2005061000'
131 end_date = '2005061000'
132 interval = 12
133 main_title1 = ' '
134 main_title2 = ' '
135 stride = 1
136 Resolution_km = 100.0
137 be_method = 'NMC'
138 NE = 1
139 Uh_method = 'scale'
140 gen_be_dir = '.'
141 code_version = 'wrfvar'
142
143 open(unit=5, file='namelist.title', status='old')
144 read(unit=5, nml=plot_title,iostat=ier)
145 if (ier/=0) then
146 print '(a,i2,a)','iostat=',ier,' Error in Namelist file...'
147 stop
148 end if
149 print plot_title
150 close (5)
151 read(start_date(1:10), fmt='(i10)')sdate
152 read(end_date (1:10), fmt='(i10)')edate
153
154 ! 3.0 Open the gen_be file and read in the information
155 ! ----------------------------------------------------
156
157 filename = trim(gen_be_dir)//'/gen_be.'//trim(be_method)//'.dat'
158 print '("*** Unit=",i3,3X,"filename=",a60)',input_unit, filename
159 open (input_unit, file = filename, form='unformatted')
160
161 ! 3.1 Read the domain size
162
163 read(input_unit) ni, nj, nk
164
165 ! 3.2. Setup the bin information
166
167 if (code_version == 'wrf3dvar') then
168
169 ! 3.2.1 the bin information from wrf3dvar/gen_be
170
171 read(input_unit)bin_type, num_bins_hgt, binwidth_hgt, binwidth_lat
172 read(input_unit)num_bins, num_bins2d
173
174 call da_create_bins( ni, nj, nk, bin_type, num_bins, num_bins2d, bin, bin2d )
175
176 else if (code_version == 'wrfvar ') then
177
178 ! 3.2.2 Read the bin information from wrfvar/gen_be
179
180 allocate( bin(1:ni,1:nj,1:nk) )
181 allocate( bin2d(1:ni,1:nj) )
182
183 read(input_unit)bin_type
184 read(input_unit)lat_min, lat_max, binwidth_lat
185 read(input_unit)hgt_min, hgt_max, binwidth_hgt
186 read(input_unit)num_bins, num_bins2d
187 read(input_unit)bin(1:ni,1:nj,1:nk)
188 read(input_unit)bin2d(1:ni,1:nj)
189
190 else
191
192 print '(a,a,a)', "Not support this version of ",code_version," code."
193 stop 22222
194
195 end if
196
197 ! ----------------------------------------------------------------------------------
198 ! PART I. Regression coefficient
199 ! ---------------------------------------------------------------------------------
200
201 ! 4.0 Read in the regression coefficients
202 ! ---------------------------------------
203
204 allocate( regcoeff_chi(1:num_bins) )
205 allocate( regcoeff_ps (1:nk,1:num_bins2d) )
206 allocate( regcoeff_t (1:nk,1:nk,1:num_bins2d) )
207
208 read(input_unit) regcoeff_chi
209 read(input_unit) regcoeff_ps
210 read(input_unit) regcoeff_t
211
212 print '(a,5I8)', &
213 "1) Finished Ceofficients reading:: ni, nj, nk, num_bins, num_bins2d:", &
214 ni, nj, nk, num_bins, num_bins2d
215
216 ! .. Skip over processing the regression coefficient ......................
217 ! goto 1000
218 ! ..........................................................................
219
220 ! 4.1 Re-assign the regression coefficients
221
222 allocate(reg_chi(nj, nk))
223 allocate(reg_ps (nj, nk))
224 allocate(reg_t (nj, ni, nk))
225
226 do k=1,nk
227 do j =1, nj
228 b = bin(1,j,k)
229 reg_chi(j,k) = regcoeff_chi(b)
230 end do
231 end do
232
233 !
234 do j=1,nj
235 b = bin2d(1,j)
236 do k=1,nk
237 reg_ps(j,k) = regcoeff_ps(k,b)
238 end do
239 end do
240 !
241 do j=1,nj
242 b = bin2d(1,j)
243 do i=1,nk
244 do k=1,nk
245 reg_t(j,i,k) = regcoeff_t(i,k,b)
246 end do
247 end do
248 end do
249
250 ! 4.1.1 domain averaged regression coefficients
251
252 if (domain_averaged) then
253
254 print '(/5x, a/)', &
255 '*** Using the averaged regression coefficients for balanced part ***'
256
257 ! .. ps and chi:
258 do k=1,nk
259 avg= 0.0
260 avg2=0.0
261 do j=1,nj
262 avg = avg + reg_ps (j,k)/float(nj)
263 avg2= avg2 + reg_chi (j,k)/float(nj)
264 end do
265 !
266 do j=1,nj
267 reg_ps (j,k)=avg
268 reg_chi(j,k)=avg2
269 end do
270 end do
271
272 ! .. temperature:
273 do i=1,nk
274 do k=1,nk
275 avg3= 0.0
276
277 do j=1,nj
278 avg3= avg3 + reg_t (j,k,i)/float(nj)
279 end do
280
281 do j=1,nj
282 reg_t(j,k,i)=avg3
283 end do
284
285 end do
286 end do
287
288 end if
289
290 print '(a)', "2) Re-assign the regression coefficients."
291
292 ! 4.2 Compute the balanced part of chi, ps and t
293
294 print '("ni,nj,nk:",3I8)', ni,nj,nk
295
296 allocate(psi (ni, nj, nk))
297 allocate(chi (ni, nj, nk))
298 allocate(chi_u (ni, nj, nk))
299 allocate(temp(ni, nj, nk))
300 allocate(t_u (ni, nj, nk))
301 allocate(ps (ni, nj))
302 allocate(ps_u(ni, nj))
303
304 allocate( chi_loc(nj, nk)); allocate( chi_global(nk))
305 allocate(temp_loc(nj, nk)); allocate(temp_global(nk))
306 allocate( ps_loc(nj))
307 allocate( chi_var(nj, nk)); allocate( chi_gvar(nk))
308 allocate(temp_var(nj, nk)); allocate(temp_gvar(nk))
309 allocate( ps_var(nj))
310
311 chi_loc = 0.0
312 temp_loc = 0.0
313 ps_loc = 0.0
314
315 date = start_date
316 cdate = sdate
317 n_times = 0
318
319 do while ( cdate <= edate )
320 n_times = n_times + 1
321
322 write(6,'(i4.4,a,a)') n_times,' Calculating unbalanced fields for date ', date
323
324 do member = 1, ne
325
326 write(ce,'(i3)')member
327 if ( member < 10 ) ce = '00'//ce(3:3)
328 if ( member >= 10 .and. member < 100 ) ce = '0'//ce(2:3)
329
330 ! 4.2.1 Read psi (mean-removed):
331 variable = 'psi'
332 filename = trim(gen_be_dir)//'/'//trim(variable)//'/'//date(1:10)
333 filename = trim(filename)//'.'//trim(variable)//'.'//trim(be_method)//'.e'//ce
334 ! print '("Read file:",a60)', filename
335 open (iunit, file = filename, form='unformatted')
336 read(iunit)ni, nj, nk
337 read(iunit)psi
338 close(iunit)
339
340 ! 4.2.2 Read chi (mean-removed):
341 variable = 'chi'
342 filename = trim(gen_be_dir)//'/'//trim(variable)//'/'//date(1:10)
343 filename = trim(filename)//'.'//trim(variable)//'.'//trim(be_method)//'.e'//ce
344 ! print '("Read file:",a60)', filename
345 open (iunit, file = filename, form='unformatted')
346 read(iunit)ni, nj, nk
347 read(iunit)chi
348 close(iunit)
349
350 ! 4.2.2.1 Read chi_u (mean-removed):
351 ! variable = 'chi_u'
352 ! filename = trim(gen_be_dir)//'/'//trim(variable)//'/'//date(1:10)
353 ! filename = trim(filename)//'.'//trim(variable)//'.'//trim(be_method)//'.e'//ce
354 ! print '("Read file:",a60)', filename
355 ! open (iunit, file = filename, form='unformatted')
356 ! read(iunit)ni, nj, nk
357 ! read(iunit)chi_u
358 ! close(iunit)
359
360 do k=1,nk
361 do j =1, nj
362 do i = 1,ni
363 balance = reg_chi(j,k)*psi(i,j,k)
364 chi_loc(j,k) = chi_loc(j,k) + balance * chi(i,j,k)
365 chi_var(j,k) = chi_var(j,k) + chi(i,j,k) * chi(i,j,k)
366
367 ! unbalan = chi(i,j,k) - balance
368 ! if (abs(chi_u(i,j,k)-unbalan) > 1.E-25) &
369 ! print '("n,i,j,k,chi_u,unbalan:",4i5,2e20.12)', n_times,i,j,k, chi_u(i,j,k), unbalan
370 end do
371 end do
372 end do
373
374 ! 4.2.3 Read ps (mean-removed):
375 variable = 'ps'
376 filename = trim(gen_be_dir)//'/'//trim(variable)//'/'//date(1:10)
377 filename = trim(filename)//'.'//trim(variable)//'.'//trim(be_method)//'.e'//ce//'.01'
378 ! print '("Read file:",a60)', filename
379 open (iunit, file = filename, form='unformatted')
380 read(iunit)ni, nj, nkdum
381 read(iunit)ldum1, ldum2 ! Dummy logicals.
382 read(iunit)ps
383 close(iunit)
384
385 ! 4.2.3.1 Read ps_u (mean-removed):
386 ! variable = 'ps_u'
387 ! filename = trim(gen_be_dir)//'/'//trim(variable)//'/'//date(1:10)
388 ! filename = trim(filename)//'.'//trim(variable)//'.'//trim(be_method)//'.e'//ce//'.01'
389 ! print '("Read file:",a60)', filename
390 ! open (iunit, file = filename, form='unformatted')
391 ! read(iunit)ni, nj, nkdum
392 ! read(iunit)ldum1, ldum2 ! Dummy logicals.
393 ! read(iunit)ps_u
394 ! close(iunit)
395
396 do j =1, nj
397 do i = 1,ni
398 balance = sum(reg_ps(j,1:nk) * psi(i,j,1:nk))
399 ps_loc(j) = ps_loc(j) + balance * ps(i,j)
400 ps_var(j) = ps_var(j) + ps(i,j) * ps(i,j)
401
402 ! unbalan = ps(i,j) - balance
403 ! if (abs(ps_u(i,j)-unbalan) > 1.E-12) &
404 ! print '("n,i,j,ps_u,unbalan:",3i5,2e20.12)', n_times,i,j,k, ps_u(i,j), unbalan
405 end do
406 end do
407
408 ! 4.2.4 Read t (mean-removed):
409 variable = 't'
410 filename = trim(gen_be_dir)//'/'//trim(variable)//'/'//date(1:10)
411 filename = trim(filename)//'.'//trim(variable)//'.'//trim(be_method)//'.e'//ce
412 ! print '("Read file:",a60)', filename
413 open (iunit, file = filename, form='unformatted')
414 read(iunit)ni, nj, nk
415 read(iunit)temp
416 close(iunit)
417
418 ! 4.2.4.1 Read t_u (mean-removed):
419 ! variable = 't_u'
420 ! filename = trim(gen_be_dir)//'/'//trim(variable)//'/'//date(1:10)
421 ! filename = trim(filename)//'.'//trim(variable)//'.'//trim(be_method)//'.e'//ce
422 ! print '("Read file:",a60)', filename
423 ! open (iunit, file = filename, form='unformatted')
424 ! read(iunit)ni, nj, nk
425 ! read(iunit)t_u
426 ! close(iunit)
427
428 do k=1,nk
429 do j =1, nj
430 do i = 1,ni
431 balance = sum(reg_t(j,k,1:nk) * psi(i,j,1:nk))
432 temp_loc(j,k) = temp_loc(j,k) + balance * temp(i,j,k)
433 temp_var(j,k) = temp_var(j,k) + temp(i,j,k) * temp(i,j,k)
434
435 ! unbalan = temp(i,j,k) - balance
436 ! if (abs(t_u(i,j,k)-unbalan) > 1.E-12) &
437 ! print '("n,i,j,k,temp_u,unbalan:",4i5,2e20.12)', n_times,i,j,k, t_u(i,j,k), unbalan
438 end do
439 end do
440 end do
441
442 if (member == ne) print '(10x,i5,a)', member, ' members has been processed.'
443
444 end do ! End loop over ensemble members.
445
446 ! Calculate next date:
447 call da_advance_cymdh( date, interval, new_date )
448 date = new_date
449 read(date(1:10), fmt='(i10)')cdate
450
451 end do ! End loop over times
452
453 ! 5.0 print the ratio of balance / full_field::
454
455 print '(/10x,a,3I8/)', 'Balanced part:', ni, nj, n_times
456
457 ! 5.1 Chi and temperature
458
459 chi_global = 0.0
460 temp_global = 0.0
461 do k=1,nk
462 do j =1, nj
463 chi_global(k) = chi_global(k) + chi_loc(j,k)
464 chi_gvar (k) = chi_gvar (k) + chi_var(j,k)
465 chi_loc(j,k) = chi_loc(j,k) / chi_var(j,k)
466
467 temp_global(k) = temp_global(k) + temp_loc(j,k)
468 temp_gvar (k) = temp_gvar (k) + temp_var(j,k)
469 temp_loc(j,k) = temp_loc(j,k) / temp_var(j,k)
470
471 end do
472 chi_global (k) = chi_global (k) / chi_gvar (k)
473 temp_global(k) = temp_global(k) / temp_gvar (k)
474 print '(i3,a,f12.5,a,f12.5)', &
475 k, " Ratio of chi, temp explained by psi: <chi_b*chi>/<chi*chi>=", &
476 chi_global(k), " and <t_b*t>/<t*t>=", temp_global(k)
477 end do
478
479 ! 5.2 Surface pressure:
480
481 print '(/a)', "Surface pressure explained by psi: <ps_b*ps>/<ps*ps>"
482 do j = 1,nj
483 ps_loc(j) = ps_loc(j) / ps_var(j)
484 print '(10x,i6,f12.5)', j,ps_loc(j)
485 end do
486
487 ! 6.0 Plot the <Xb*X> / <X*X>
488
489 ! 6.1 Surface pressure and chi, Temp globally explained by psi
490
491 call da_balance_ratio_plot (ps_loc, chi_global, temp_global, nj, nk, main_title1)
492
493 ! 6.2 chi and Temp expalined by psi locally
494
495 NDOT = -682 ! dashed-line pattern
496 flo = 0.0
497 hi = 0.0
498 call setusv('LW',2000)
499
500 ! 6.2.1 Balanced chi
501
502 Print '(a)', 'PLOT the local <chi_b*chi>/<chi*chi>...'
503
504 call set(0.20,1.0,0.2,0.9, 1.0,float(nj),1.0,float(nk),1)
505
506 ! Contour:
507 cint= 0.025
508 call CONREC(chi_loc,nj,nk,nk,flo,hi, cint,-1, 0, NDOT)
509 call PERIML( nj/10, 5, nk/5, 5)
510
511 write(Title,'("CHI_B LOcallY EXPLAinED BY PSI")')
512 print '(a)', Title
513 call set(0.1,0.97,0.1,0.97, 0.0,100.,0.0,100.,1)
514 call PWRITX(50.0, 95.0,Title,23,20,0,0)
515 call PWRITX(40.0, 2.0,"GRID in Y-DIRECTION",19,20,0,0)
516 call PWRITX(1.5,30.0,"<CHI_B*CHI> / <CHI*CHI>",23,20,90,-1)
517 call FRAME
518
519 ! 6.2.1 Balanced temperature
520
521 Print '(a)', 'PLOT the local <temp_b*temp>/<temp*temp>...'
522
523 call set(0.20,1.0,0.2,0.9, 1.0,float(nj),1.0,float(nk),1)
524
525 ! Contour:
526 cint= 0.050
527 call CONREC(temp_loc,nj,nk,nk,flo,hi, cint,-1, 0, NDOT)
528 call PERIML( nj/10, 5, nk/5, 5)
529
530 write(Title,'("TEMP_B LOcallY EXPLAinED BY PSI")')
531 print '(a)', Title
532 call set(0.1,0.97,0.1,0.97, 0.0,100.,0.0,100.,1)
533 call PWRITX(50.0, 95.0,Title,24,20,0,0)
534 call PWRITX(40.0, 2.0,"GRID in Y-DIRECTION",19,20,0,0)
535 call PWRITX(1.5,30.0,"<TEMP_B*TEMP> / <TEMP*TEMP>",26,20,90,-1)
536 call FRAME
537
538 ! --------------------------------------------------------------------------
539 ! PART II. Eigenvector and eigenvalue
540 ! --------------------------------------------------------------------------
541 1000 continue
542
543 call da_allocate_gen_be(cv, ni, nj, nk)
544
545 do n = 1,5
546
547 read(input_unit) variable
548 read(input_unit) nk, num_bins2d
549 print*,'n=',n,' Eigenvector/values for variable ',variable,nk,num_bins2d
550 cvv = trim(variable)
551 !
552 if ( n <=4) then
553
554 read(input_unit) cv%g_vec
555 read(input_unit) cv%g_val
556 read(input_unit) cv%l_vec
557 read(input_unit) cv%l_val
558
559 ! .. Skip over plotting eigenvectors ...................................
560 ! goto 2000
561 ! ......................................................................
562
563 if (cvv == 'psi') CVV = 'PSI'
564 if (cvv == 'chi_u') CVV = 'CHI_U'
565 if (cvv == 't_u') CVV = 'T_U'
566 if (cvv == 'rh') CVV = 'RH'
567
568 ! .. Global eigenvectors:
569 call da_eigen_plot(cvv, cv, main_title1)
570
571 ! .. local eigenvalues:
572 call da_eigen_plot2(cvv, cv, main_title1, 99)
573
574 else
575
576 read(input_unit) cv%g_vec(1,1)
577 read(input_unit) cv%g_val(1)
578 read(input_unit) cv%l_vec(1,1,1:nj)
579 read(input_unit) cv%l_val(1,1:nj)
580
581 if (cvv == 'ps_u') CVV = 'PS_U'
582 if (cvv == 'ps') CVV = 'PS'
583
584 do j = 1,nj
585 ps_var(j) = cv%l_val(1,j)
586 end do
587
588 call da_covariance_plot (cvv, ps_var, nj, main_title1)
589
590 end if
591
592 2000 continue
593
594 end do
595
596
597 ! -------------------------------------------------------------------------
598
599 do n = 1,5
600 read(input_unit) variable
601 print*,'n=',n,' scale length for variable ',variable,nk,num_bins2d
602 if (n==5) then
603 cv%sl = 0.0
604 read(input_unit) cv%sl(1:1)
605 else
606 read(input_unit) cv%sl
607 end if
608
609 cvv = trim(variable)
610
611 if (cvv == 'psi') CVV = 'PSI'
612 if (cvv == 'chi_u') CVV = 'CHI_U'
613 if (cvv == 't_u') CVV = 'T_U'
614 if (cvv == 'rh') CVV = 'RH'
615
616 call da_scale_plot(cvv, cv, main_title1, resolution_km)
617
618
619 end do
620
621 ! -------------------------------------------------------------------------
622
623 call clsgks
624
625 stop
626
627 contains
628
629 subroutine da_setup_color_table
630
631 implicit none
632
633 call gscr(1, 0, 1.00, 1.00, 1.00) ! White
634 call gscr(1, 1, 0.00, 0.00, 0.00) ! BLACK
635 call gscr(1, 2, 1.00, 0.00, 0.00) ! Red
636 call gscr(1, 3, 0.00, 1.00, 0.00) ! Green
637 call gscr(1, 4, 0.00, 0.00, 1.00) ! Blue
638 call gscr(1, 5, 0.60, 0.00, 0.80) ! Dark violet
639 call gscr(1, 6, 0.00, 1.00, 1.00) ! Cyran
640 call gscr(1, 7, 1.00, 0.00, 1.00) ! Magenta
641 call gscr(1, 8, 0.90, 0.25, 0.00) ! FreshRed
642 call gscr(1, 9, 0.40, 0.30, 0.20) ! Tan
643 call gscr(1, 10, 1.00, 1.00, 0.00) ! Yellow
644 call gscr(1, 11, 0.60, 0.60, 0.60) ! Gray
645
646 end subroutine da_setup_color_table
647
648 subroutine da_allocate_gen_be(be, ix, jy, kz)
649
650 implicit none
651
652 type (eigen_type), intent(inout) :: be
653 integer, intent(in) :: ix, jy, kz
654 integer :: max_waveno
655
656 max_waveno = ix/2 - 1
657
658 allocate(be%g_val(kz))
659 allocate(be%g_vec(kz, kz))
660 allocate(be%l_vec(kz, kz,jy))
661 allocate(be%l_val(kz,jy))
662 allocate(be%sl(kz))
663 allocate(be%power(max_waveno,kz))
664
665 be%g_val = 0.0
666 be%g_vec = 0.0
667 be%l_val = 0.0
668 be%l_vec = 0.0
669 be%sl = 0.0
670
671 end subroutine da_allocate_gen_be
672
673 subroutine da_advance_cymdh( start_date, dh, end_date )
674
675 implicit none
676
677 character (len=10), intent(in) :: start_date ! In date (ccyymmddhh).
678 integer, intent(in) :: dh ! Period to advance (-ve for past).
679 character (len=10), intent(out) :: end_date ! Out date (ccyymmddhh).
680
681 integer :: ccyy, mm, dd, hh
682
683 read(start_date(1:10), fmt='(i4, 3i2)') ccyy, mm, dd, hh
684
685 hh = hh + dh
686
687 do while (hh < 0)
688 hh = hh + 24
689 call da_change_date ( ccyy, mm, dd, -1 )
690 end do
691
692 do while (hh > 23)
693 hh = hh - 24
694 call da_change_date ( ccyy, mm, dd, 1 )
695 end do
696
697 write(end_date(1:10), fmt='(i4, 3i2.2)') ccyy, mm, dd, hh
698
699 end subroutine da_advance_cymdh
700
701 subroutine da_change_date( ccyy, mm, dd, delta )
702
703 implicit none
704
705 integer, intent(inout) :: ccyy, mm, dd
706 integer, intent(in) :: delta
707
708 integer, dimension(12) :: mmday
709
710 mmday = (/31,28,31,30,31,30,31,31,30,31,30,31/)
711
712 mmday(2) = 28
713
714 if (mod(ccyy,4) == 0) then
715 mmday(2) = 29
716
717 if ( mod(ccyy,100) == 0) then
718 mmday(2) = 28
719 end if
720
721 if(mod(ccyy,400) == 0) then
722 mmday(2) = 29
723 end if
724 end if
725
726 dd = dd + delta
727
728 if(dd == 0) then
729 mm = mm - 1
730
731 if(mm == 0) then
732 mm = 12
733 ccyy = ccyy - 1
734 end if
735
736 dd = mmday(mm)
737 elseif ( dd .gt. mmday(mm) ) then
738 dd = 1
739 mm = mm + 1
740 if(mm > 12 ) then
741 mm = 1
742 ccyy = ccyy + 1
743 end if
744 end if
745 end subroutine da_change_date
746
747 subroutine da_create_bins( ni, nj, nk, bin_type, num_bins, num_bins2d, &
748 bin, bin2d, &
749 lat_min, lat_max, binwidth_lat, &
750 hgt_min, hgt_max, binwidth_hgt, latitude, height )
751 !----------------------------------------------------------------------------------------
752 !
753 ! Purpose: To create the bins for calculation of statistics.
754 !
755 ! Input : ni, nj, nk ----- Dimensions
756 ! bin_type ----- 0: No binning;
757 ! 1: bins for X-direction mean;
758 ! 2: bins for each of binwidth_lat/binwidth_hgt.
759 ! 3: bins for each of binwidth_lat/nk.
760 ! 4: num_hor_bins horizontal bins /nk.
761 ! 5: Average over all horizontal points (nk bins for 3D fields).
762 ! 6: Average over all points (only 1 bin).
763 ! Optional for bin_type = 2:
764 ! lat_min, lat_max-- Minimum/maximum latitude ranges for bin_type = 2
765 ! binwidth_lat ----- interval between bins for latitude in degree for bin_type = 2
766 ! binwidth_hgt ----- interval between bins for height in meter for bin_type = 2
767 ! num_bins_hgt ----- the number of height bins for bin_type = 2
768 ! latitude ----- 3d field latitude in degree for bin_type = 2
769 ! height ----- 3d field height in meter for bin_type = 2
770 !
771 ! Output: num_bins,num_bins2d ---- the number of bins for 3d and 2d fields
772 ! bin , bin2d ---- Assigned bin to a gridpoints for 3d and 2d fields
773 !
774 !-----------------------------------------------------------------------------------------
775 implicit none
776
777 integer, intent(in) :: ni, nj, nk ! Dimensions read in.
778 integer, intent(in) :: bin_type ! Type of bin to average over.
779 integer, intent(out) :: num_bins ! Number of bins.
780 integer, intent(out) :: num_bins2d ! Number of bins for 2D fields.
781 integer, intent(out) :: bin(1:ni,1:nj,1:nk) ! Bin assigned to each 3D point.
782 integer, intent(out) :: bin2d(1:ni,1:nj) ! Bin assigned to each 2D point.
783
784 real(kind=8), intent(in),optional:: lat_min, lat_max ! Used if bin_type = 2 (degrees).
785 real(kind=8), intent(in),optional:: binwidth_lat ! Used if bin_type = 2 (degrees).
786 real(kind=8), intent(in),optional:: hgt_min, hgt_max ! Used if bin_type = 2 (degrees).
787 real(kind=8), intent(in),optional:: binwidth_hgt ! Used if bin_type = 2 (m).
788 real(kind=8), intent(in),optional:: latitude(1:ni,1:nj) ! Latitude (degrees).
789 real(kind=8), intent(in),optional:: height(1:ni,1:nj,1:nk) ! Height (m).
790
791 integer :: b, i, j, k ! Loop counters.
792 integer :: count ! Counter
793 integer :: num_bins_lat ! Used if bin_type = 2.
794 integer :: num_bins_hgt ! Used if bin_type = 2.
795 integer :: bin_lat ! Latitude bin.
796 integer :: bin_hgt ! Height bin.
797 integer :: num_bins_i, num_bins_j ! Used if bin_type = 4.
798 integer :: nii, njj ! Used if bin_type = 4.
799 integer :: bin_i(1:ni), bin_j(1:nj) ! Used if bin_type = 4.
800 real(kind=8), allocatable :: binstart_lat(:) ! Used if bin_type = 2 (degrees).
801 real(kind=8), allocatable :: binstart_hgt(:) ! Used if bin_type = 2 (degrees).
802
803 if ( bin_type == 0 ) then ! No averaging in space
804
805 num_bins = nk * nj * ni
806 num_bins2d = nj * ni ! Equals number of horizontal points.
807
808 count = 1
809 do k = 1, nk
810 do j = 1, nj
811 do i = 1, ni
812 bin(i,j,k) = count
813 count = count + 1
814 end do
815 end do
816 end do
817 bin2d(:,:) = bin(:,:,1)
818
819 else if ( bin_type == 1 ) then ! Average over x-direction.
820
821 num_bins = nj * nk
822 num_bins2d = nj
823
824 count = 1
825 do k = 1, nk
826 do j = 1, nj
827 bin(1:ni,j,k) = count
828 count = count + 1
829 end do
830 end do
831 bin2d(:,:) = bin(:,:,1)
832
833 else if ( bin_type == 2 ) then ! Global latitude/height bins:
834
835 ! Setup latitude bins:
836 write(6,'(/a,f12.5)')' Minimum latitude = ', lat_min
837 write(6,'(a,f12.5)')' Maximum latitude = ', lat_max
838 write(6,'(a,f12.5)')' Latitude bin width = ', binwidth_lat
839 num_bins_lat = nint( ( lat_max - lat_min ) / binwidth_lat)
840 write(6,'(a,i8)') ' Number of latitude bins = ', num_bins_lat
841
842 allocate( binstart_lat(1:num_bins_lat) )
843 do b = 1, num_bins_lat ! Assume south to north (as in WRF).
844 binstart_lat(b) = lat_min + real(b-1) * binwidth_lat
845 ! write(6,'(i4,f15.5)')b, binstart_lat(b)
846 end do
847
848 ! Setup height bins:
849 write(6,'(/a,f12.5)')' Minimum height = ', hgt_min
850 write(6,'(a,f12.5)')' Maximum height = ', hgt_max
851 write(6,'(a,f12.5)')' Height bin width = ', binwidth_hgt
852 num_bins_hgt = nint( ( hgt_max - hgt_min ) / binwidth_hgt)
853 write(6,'(a,i8)')' Number of height bins = ', num_bins_hgt
854
855 allocate( binstart_hgt(1:num_bins_hgt) )
856 do b = 1, num_bins_hgt
857 binstart_hgt(b) = hgt_min + real(b-1) * binwidth_hgt
858 ! write(6,'(i4,f15.5)')b, binstart_hgt(b)
859 end do
860
861 num_bins = num_bins_lat * num_bins_hgt
862 num_bins2d = num_bins_lat
863
864 ! Select height bins:
865 do j = 1, nj
866 do i = 1, ni
867 do k = 1, nk
868 if ( height(i,j,k) < binstart_hgt(1) ) then
869 bin_hgt = 1 ! In first bin.
870 else if ( height(i,j,k) >= binstart_hgt(num_bins_hgt) ) then
871 bin_hgt = num_bins_hgt ! In final bin.
872 else
873 do b = 1, num_bins_hgt-1
874 if ( height(i,j,k) >= binstart_hgt(b) .and. &
875 height(i,j,k) < binstart_hgt(b+1) ) then
876 bin_hgt = b
877 exit
878 end if
879 end do
880 end if
881
882 ! Select latitude bin that point falls in:
883 if ( k == 1 ) then
884 do b = 1, num_bins_lat-1
885 if ( latitude(i,j) >= binstart_lat(b) .and. &
886 latitude(i,j) < binstart_lat(b+1) ) then
887 bin_lat = b
888 exit
889 end if
890 end do
891 if ( latitude(i,j) >= binstart_lat(num_bins_lat) ) then ! In final bin.
892 bin_lat = num_bins_lat
893 end if
894 bin2d(i,j) = bin_lat
895 end if
896 bin(i,j,k) = bin_lat + num_bins_lat * ( bin_hgt - 1 )
897 end do
898 end do
899 end do
900
901 deallocate( binstart_lat )
902 deallocate( binstart_hgt )
903
904 else if ( bin_type == 3 ) then ! Latitude/nk bins:
905
906 ! Setup latitude bins:
907 write(6,'(/a,f12.5)')' Minimum latitude = ', lat_min
908 write(6,'(a,f12.5)')' Maximum latitude = ', lat_max
909 write(6,'(a,f12.5)')' Latitude bin width = ', binwidth_lat
910 num_bins_lat = nint( ( lat_max - lat_min ) / binwidth_lat)
911 write(6,'(a,i8)') ' Number of latitude bins = ', num_bins_lat
912
913 allocate( binstart_lat(1:num_bins_lat) )
914 do b = 1, num_bins_lat ! Assume south to north (as in WRF).
915 binstart_lat(b) = lat_min + real(b-1) * binwidth_lat
916 ! write(6,'(i4,f15.5)')b, binstart_lat(b)
917 end do
918
919 num_bins = num_bins_lat * nk
920 num_bins2d = num_bins_lat
921
922 ! Select bins:
923 do j = 1, nj
924 do i = 1, ni
925 do k = 1, nk
926
927 ! Select latitude bin that point falls in:
928 if ( k == 1 ) then
929 do b = 1, num_bins_lat-1
930 if ( latitude(i,j) >= binstart_lat(b) .and. &
931 latitude(i,j) < binstart_lat(b+1) ) then
932 bin_lat = b
933 exit
934 end if
935 end do
936 if ( latitude(i,j) >= binstart_lat(num_bins_lat) ) then ! In final bin.
937 bin_lat = num_bins_lat
938 end if
939 bin2d(i,j) = bin_lat
940 end if
941 bin(i,j,k) = bin_lat + num_bins_lat * ( k - 1 )
942 end do
943 end do
944 end do
945
946 deallocate( binstart_lat )
947
948 else if ( bin_type == 4 ) then ! binwidth_lat/nk bins:
949
950 ! Setup horizontal bins:
951 write(6,'(/a,f12.5)')' Number of grid-cells to average over = ', binwidth_lat
952 ! use binwidth_lat, but actually an integer number of points.
953
954 num_bins_j = int( real(nj) / real(binwidth_lat) )
955 njj = int(binwidth_lat) * num_bins_j
956 do j = 1, njj
957 bin_j(j) = 1 + int( real(j-1) / binwidth_lat)
958 end do
959 if ( nj > njj ) bin_j(njj+1:nj) = bin_j(njj)
960
961 num_bins_i = int( real(ni) / real(binwidth_lat) )
962 nii = int(binwidth_lat) * num_bins_i
963 do i = 1, nii
964 bin_i(i) = 1 + int( real(i-1) / binwidth_lat )
965 end do
966 if ( ni > nii ) bin_i(nii+1:ni) = bin_i(nii)
967
968 num_bins2d = num_bins_i * num_bins_j
969 num_bins = num_bins2d * nk
970
971 do j = 1, nj
972 do i = 1, ni
973 bin2d(i,j) = bin_i(i) + ( bin_j(j) - 1 ) * num_bins_i
974 do k = 1, nk
975 bin(i,j,k) = bin2d(i,j) + (k - 1) * num_bins2d
976 end do
977 end do
978 end do
979
980 else if ( bin_type == 5 ) then ! Average over all horizontal points.
981
982 num_bins = nk
983 num_bins2d = 1
984
985 do k = 1, nk
986 bin(:,:,k) = k
987 end do
988 bin2d(:,:) = 1
989
990 else if ( bin_type == 6 ) then ! Average over all points.
991
992 num_bins = 1
993 num_bins2d = 1
994 bin(:,:,:) = 1
995 bin2d(:,:) = 1
996
997 end if
998
999 end subroutine da_create_bins
1000
1001 subroutine da_eigen_plot(vn, ev, main_title)
1002
1003 implicit none
1004
1005 integer, parameter :: solid_line = 65535, & ! PATTERN = 1111111111111111
1006 thick_dash = 21845, & ! PATTERN = 0101010101010101
1007 thin_dash = 3855, & ! PATTERN = 0000111100001111
1008 solid_like = 31710 ! PATTERN = 0111101111011110
1009
1010 integer, parameter :: black = 0, &
1011 white = 1, &
1012 red = 2, &
1013 green = 3, &
1014 blue = 4, &
1015 violet = 5, &
1016 cyran = 6, &
1017 magenta = 7, &
1018 freshred = 8, &
1019 tan = 9, &
1020 yellow = 10, &
1021 gray = 11
1022
1023 character(len= 5), intent(in) :: vn
1024 type (eigen_type), intent(in) :: ev
1025 character(len= 120), intent(in) :: main_title
1026 character(len=80) :: x_title, y_title
1027 real, dimension(501) :: x, y
1028
1029 integer :: mn, plot_type, mmx, mnx, mmy, mny, magnitude, &
1030 k, ix, jy, kz, line_color
1031
1032 real :: xb, xe, yb, ye, xmin, xmax, fct
1033
1034 print '("PLOT:", a,2x,a)', vn, main_title(1:40)
1035
1036 plot_type = 1
1037
1038 ix = size(ev%l_vec, dim=3) ! dimension of bin
1039 jy = size(ev%l_vec, dim=1) ! k-level
1040 kz = size(ev%l_vec, dim=2) ! mode index
1041
1042 mmy = 6
1043
1044 mny = 2
1045
1046 yb = 1.0
1047 ye = real(kz)
1048
1049 y_title = 'VERTICAL LEVEL'
1050
1051 do k=1, 501
1052 y(k) = real(k)
1053 end do
1054
1055 ! if (ev%ok_ev) then
1056
1057 !------------- plot eigen value ---------------------
1058
1059 mmx = 1
1060 mnx = 1
1061
1062 do k=1,kz
1063 x(k) = ev%g_val(k)
1064 end do
1065
1066 x_title = vn // ' E VAL'
1067
1068 xmin = minval(x(1:kz))
1069 xmax = maxval(x(1:kz))
1070
1071 call da_norm_min_max(xmin, xmax, magnitude)
1072
1073 fct = 10.0**magnitude
1074
1075 x(1:kz) = fct * x(1:kz)
1076
1077 !mslee xb = 1.0e-15
1078 xb = 1.0e-5
1079 xe = 100.0
1080
1081 mn = 0
1082
1083 call da_line_plot(x_title, y_title, main_title, &
1084 x,y,kz,mn,xb,xe,yb,ye,3, &
1085 mmx, mnx, mmy, mny, &
1086 red, 3500, thick_dash, magnitude,0)
1087
1088 call frame
1089
1090 !-------------------plot eigenvectors ----------------
1091
1092 x_title = vn // ' E VEC'
1093
1094 mmx = 8
1095 mnx = 2
1096
1097 do k=1,kz
1098 x(k) = ev%g_vec(k,1)
1099 end do
1100
1101 xmin = minval(x(1:kz))
1102 xmax = maxval(x(1:kz))
1103
1104 call da_norm_min_max(xmin, xmax, magnitude)
1105
1106 fct = 10.0**magnitude
1107
1108 xb =-80.0
1109 xe = 80.0
1110
1111 line_color = red
1112
1113 do mn=1,5
1114 fct = 10.0**magnitude *ev%g_val(mn) / ev%g_val(1)
1115
1116 do k=1,kz
1117 !test x(k) = ev%g_vec(k,mn)*fct
1118 x(k) = ev%g_vec(k,mn)*10.0**magnitude
1119 end do
1120
1121 call da_line_plot(x_title, y_title, main_title, &
1122 x,y,kz,mn,xb,xe,yb,ye,plot_type, &
1123 mmx, mnx, mmy, mny, &
1124 line_color, 6000-1000*mn, thick_dash, magnitude,0)
1125
1126 line_color = line_color + 1
1127
1128 if(line_color == cyran ) line_color = line_color + 1
1129 if(line_color > 7) line_color = red
1130 end do
1131
1132 call frame
1133
1134 end subroutine da_eigen_plot
1135
1136 subroutine da_eigen_plot2 (vn, ev, main_title, iopt)
1137 implicit none
1138
1139 integer, parameter :: solid_line = 65535, & ! PATTERN = 1111111111111111
1140 thick_dash = 21845, & ! PATTERN = 0101010101010101
1141 thin_dash = 3855, & ! PATTERN = 0000111100001111
1142 solid_like = 31710 ! PATTERN = 0111101111011110
1143
1144 integer, parameter :: black = 0, &
1145 white = 1, &
1146 red = 2, &
1147 green = 3, &
1148 blue = 4, &
1149 violet = 5, &
1150 cyran = 6, &
1151 magenta = 7, &
1152 freshred = 8, &
1153 tan = 9, &
1154 yellow = 10, &
1155 gray = 11
1156
1157 character(len= 5), intent(in) :: vn
1158 type (eigen_type), intent(in) :: ev
1159 character(len= 120), intent(in) :: main_title
1160 character(len=80) :: x_title, y_title
1161 real, dimension(501) :: x, y
1162
1163 integer :: iopt
1164
1165 integer :: mn, plot_type, mmx, mnx, mmy, mny, magnitude, &
1166 k, ix, jy, kz, line_color
1167
1168 real :: xb, xe, yb, ye, ymin, ymax, fct
1169
1170 print '("PLOT LOCAL E VAL:", a,2x,a)', vn, main_title(1:40)
1171
1172 plot_type = 1
1173
1174 if (iopt.eq.99) then
1175 ix = size(ev%l_vec, dim=3) ! dimension of bin
1176 jy = size(ev%l_vec, dim=1) ! k-level
1177 kz = size(ev%l_vec, dim=2) ! mode index
1178
1179 elseif (iopt.eq.98) then
1180 ix = size(ev%power, dim=1)
1181 jy = size(ev%power, dim=2)
1182 end if
1183
1184 mmy = 6
1185 mny = 2
1186
1187 xb = 1.0
1188 xe = real(ix)
1189
1190 if (iopt .eq. 99) y_title = 'EIGEN VALUE'
1191 if (iopt .eq. 98) y_title = 'POWER SPECTRUM'
1192
1193 do k=1, ix
1194 x(k) = real(k)
1195 end do
1196
1197 !-----------------------------------------------------
1198
1199 if (iopt .eq. 99) x_title = vn // ' LATITUDE (GRID)'
1200 if (iopt .eq. 98) x_title = vn // ' WAVE NUMBER'
1201
1202 mmx = 8
1203 mnx = 2
1204
1205 do j=1,ix
1206 if (iopt.eq.99) y(j) = ev%l_val(1,j)
1207 if (iopt.eq.98) y(j) = ev%power(j,1)
1208 end do
1209
1210 ymin = minval(y(1:ix))
1211 ymax = maxval(y(1:ix))
1212
1213 call da_norm_min_max(ymin, ymax, magnitude)
1214
1215 fct = 10.0**magnitude
1216
1217 if (iopt .eq. 99) then
1218
1219 yb =0.0
1220 ye =100.0
1221
1222 if (vn.eq.'CHI_U'.or.vn.eq.'CHI'.or.vn.eq.'CHI_B') then
1223 ye=50.0
1224 elseif (vn.eq.'RH') then
1225 ye =100.0
1226 end if
1227
1228 else
1229
1230 yb = 0.0
1231 ye = 150.0
1232 if (vn.eq.'CHI_U'.or.vn.eq.'CHI'.or.vn.eq.'CHI_B'.or. &
1233 vn.eq.'RH') ye=20.0
1234 if (vn.eq.'T_U'.or.vn.eq.'T'.or.vn.eq.'T_B') &
1235 ye=50.0
1236
1237 end if
1238
1239 ! if (iopt.eq.98.and.vn.eq.'rh') ye=150.
1240 if (iopt.eq.98.and.vn.eq.'chi_b') ye=100.
1241 line_color = red
1242
1243 do mn=1,5
1244
1245 do j=1,ix
1246 !test x(k) = ev%g_vec(k,mn)*fct
1247 if (iopt.eq.99) y(j) = ev%l_val(mn,j)*10.0**magnitude
1248 if (iopt.eq.98) y(j) = ev%power(j,mn)*10.0**magnitude
1249
1250 end do
1251 call da_line_plot(x_title, y_title, main_title, &
1252 x,y,ix,mn,xb,xe,yb,ye,plot_type, &
1253 mmx, mnx, mmy, mny, &
1254 line_color, 6000-1000*mn, thick_dash, magnitude,2)
1255
1256 line_color = line_color + 1
1257
1258 if(line_color == cyran ) line_color = line_color + 1
1259 if(line_color > 7) line_color = red
1260 end do
1261
1262 call frame
1263
1264 end subroutine da_eigen_plot2
1265
1266 subroutine da_scale_plot(vn, ev, main_title, dist)
1267
1268 implicit none
1269
1270 integer, parameter :: solid_line = 65535, & ! PATTERN = 1111111111111111
1271 thick_dash = 21845, & ! PATTERN = 0101010101010101
1272 thin_dash = 3855, & ! PATTERN = 0000111100001111
1273 solid_like = 31710 ! PATTERN = 0111101111011110
1274
1275 integer, parameter :: black = 0, &
1276 white = 1, &
1277 red = 2, &
1278 green = 3, &
1279 blue = 4, &
1280 violet = 5, &
1281 cyran = 6, &
1282 magenta = 7, &
1283 freshred = 8, &
1284 tan = 9, &
1285 yellow = 10, &
1286 gray = 11
1287
1288 character(len= 5), intent(in) :: vn
1289 type (eigen_type), intent(in) :: ev
1290 character(len= 120), intent(in) :: main_title
1291 real, intent(in) :: dist
1292 character(len=80) :: x_title, y_title
1293 real, dimension(501) :: x, y
1294
1295 integer :: mn, plot_type, mmx, mnx, mmy, mny, magnitude, &
1296 k, ix, jy, kz, line_color
1297
1298 real :: xb, xe, yb, ye, xmin, xmax, fct
1299
1300 print '("PLOT SCALE LENGTH:", a,2x,a)', vn, main_title(1:40)
1301
1302 plot_type = 1
1303
1304 ix = size(ev%l_vec, dim=3) ! dimension of bin
1305 jy = size(ev%l_vec, dim=1) ! k-level
1306 kz = size(ev%l_vec, dim=2) ! mode index
1307
1308 mmy = 6
1309
1310 mny = 2
1311
1312 yb = 1.0
1313 ye = real(kz)
1314
1315 do k=1, 501
1316 y(k) = real(k)
1317 end do
1318
1319 mmx = 5
1320 mnx = 5
1321
1322 do k=1,kz
1323 x(k) = ev%sl(k)
1324 write(unit=15, fmt='(i4,f10.1,3x,f20.1)') k,x(k), dist*x(k)
1325 end do
1326 x_title = vn // ' SCALE LENGTH(KM)'
1327
1328 x(1:kz) = x(1:kz) * dist
1329 xmin = minval(x(1:kz))
1330 xmax = maxval(x(1:kz))
1331
1332 call da_norm_min_max(xmin, xmax, magnitude)
1333
1334 fct = 10.0**magnitude
1335 x(1:kz) = fct * x(1:kz)
1336
1337 xb = 0.0
1338 xe = maxval (x(1:kz))
1339
1340 mn = 0
1341
1342 y_title = 'VERTICAL MODE'
1343
1344 call da_line_plot(x_title, y_title, main_title, &
1345 x,y,kz,mn,xb,xe,yb,ye, 1, &
1346 mmx, mnx, mmy, mny, &
1347 red, 3500, thick_dash, magnitude,0)
1348 call frame
1349
1350 !-----------------------------------------------------
1351
1352 end subroutine da_scale_plot
1353
1354 subroutine da_line_plot(x_title, y_title, main_title, &
1355 z,y,n,mn,xb,xe,yb,ye,plot_type, &
1356 mmx, mix, mmy, iy, &
1357 color, width, line_pattern, magnitude,iopt)
1358
1359 implicit none
1360
1361 real, parameter :: xfb = 0.10, &
1362 xfe = 0.90, &
1363 yfb = 0.10, &
1364 yfe = 0.90
1365
1366 integer, parameter :: black = 0, &
1367 white = 1, &
1368 red = 2, &
1369 green = 3, &
1370 blue = 4, &
1371 violet = 5, &
1372 cyran = 6, &
1373 magenta = 7, &
1374 freshred = 8, &
1375 tan = 9, &
1376 yellow = 10, &
1377 gray = 11
1378
1379 integer :: iopt
1380 character(len=80), intent(in) :: x_title, y_title
1381 character(len=120), intent(in) :: main_title
1382 integer, intent(in) :: n, & ! vertical points
1383 mn, & ! mode number
1384 plot_type, &
1385 mmx, mix, mmy, iy, &
1386 color, width, line_pattern, magnitude
1387 real, dimension(501), intent(in) :: z, y
1388 real, intent(in) :: xb,xe,yb,ye
1389
1390 integer :: lx, ly
1391 real :: xc, yc, xw, yw
1392 character(len=120) :: title,title2
1393 character(len= 2) :: chr_2
1394 character(len= 3) :: chr_3
1395
1396 real, dimension(501) :: x
1397
1398 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1399
1400 call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
1401
1402 write(chr_3(1:3), fmt='(i3)') magnitude
1403
1404 if (iopt.eq.0)then
1405 title = trim(x_title) // ' MAGNIFIED BY 10**(' // chr_3 // ')'
1406 title2 = trim(y_title)
1407 elseif (iopt.eq.2) then
1408 title = trim(x_title)
1409 title2 = trim(y_title)//' MAGNIFIED BY 10**(' // chr_3 // ')'
1410 end if
1411 print '(I2," magnitude=",i6,2x,a)', mn, magnitude, chr_3
1412 print '("title=",a)', title
1413
1414 lx = len_trim( title)
1415 ly = len_trim(title2)
1416
1417 call setusv('LW',1500)
1418
1419 call gsplci(5)
1420 call gspmci(5)
1421 call gstxci(5)
1422
1423 call pwrity(0.1-0.075,0.5,trim(title2),ly,16,90,0)
1424 call pwrity(0.5,0.1-0.075, title, lx,16, 0,0)
1425 lx = len_trim( main_title)
1426 call pwrity(0.5,1.0-0.05, main_title, lx,16, 0,0)
1427
1428 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1429
1430 ! call gsplci(cyran)
1431 ! call gspmci(cyran)
1432 ! call gstxci(cyran)
1433 call gsplci(blue)
1434 call gspmci(blue)
1435 call gstxci(blue)
1436
1437 call set(xfb,xfe,yfb,yfe,xb,xe,yb,ye,plot_type)
1438
1439 call setusv('LW',2000)
1440
1441 if(plot_type > 2) then
1442 call labmod('(E7.1)','(I4)',7,4,15,15,0,0,0)
1443 else
1444 call labmod('(I3)','(I3)',3,3,15,15,0,0,0)
1445 end if
1446
1447 call halfax(mmx, mix, mmy, iy, xb, yb, 1, 1)
1448
1449 if(plot_type < 3) call line(0.0, y(1), 0.0, y(n))
1450
1451 call gsplci(color)
1452 call gspmci(color)
1453 call gstxci(color)
1454
1455 call dashdb (line_pattern)
1456
1457 call setusv('LW',width)
1458
1459 x = z
1460 do lx=1,n
1461 ! write(unit=*, fmt='(a, i5, 2e12.4)') &
1462 ! 'n,x,y=', lx, x(lx), y(lx)
1463 if(x(lx) < xb) x(lx) = xb
1464 end do
1465
1466 call curve(x,y,n)
1467
1468 if(mn > 0) then
1469
1470 write(chr_2(1:2), fmt='(i2.2)') mn
1471
1472 call pwrity(x(1), y(1)-0.005*(y(n)-y(1)), chr_2, 2, 16, 0, 0)
1473 end if
1474
1475 if (x_title(7:9) == 'VEC') then
1476 call line(xe-25.0, ye-mn*1.5, xe-3.0,ye-mn*1.5)
1477 write(title,'("VECT ",I2)') mn
1478 call setusv('LW',2000)
1479 call pwrity(xe-15.0, ye-mn*1.5+0.5, title(1:7),7,10,0,0)
1480 end if
1481
1482 end subroutine da_line_plot
1483
1484 subroutine da_norm_min_max(xmin, xmax, magnitude)
1485
1486 implicit none
1487
1488 real, intent(inout) :: xmin, xmax
1489 integer, intent(out) :: magnitude
1490
1491
1492 real :: ymin, ymax
1493
1494 magnitude = 0
1495
1496 if(abs(xmin) < 1.0e-30 .and. abs(xmax) < 1.0e-30) return
1497
1498
1499 do
1500 if(abs(xmin) > 100.0 .or. abs(xmax) > 100.0) then
1501 magnitude = magnitude - 1
1502 xmin = 0.1*xmin
1503 xmax = 0.1*xmax
1504 else
1505 exit
1506 end if
1507 end do
1508
1509 do
1510 if(abs(xmin) < 10.0 .and. abs(xmax) < 10.0) then
1511 magnitude = magnitude + 1
1512 xmin = 10.0*xmin
1513 xmax = 10.0*xmax
1514 else
1515 exit
1516 end if
1517 end do
1518
1519 xmin = -100.0
1520 xmax = 100.0
1521
1522 end subroutine da_norm_min_max
1523
1524 subroutine da_balance_ratio_plot (ps, chi, temp, nj, nk,main_title)
1525
1526 implicit none
1527
1528 integer, parameter :: solid_line = 65535, & ! PATTERN = 1111111111111111
1529 thick_dash = 21845, & ! PATTERN = 0101010101010101
1530 thin_dash = 3855, & ! PATTERN = 0000111100001111
1531 solid_like = 31710 ! PATTERN = 0111101111011110
1532
1533 integer, parameter :: black = 0, &
1534 white = 1, &
1535 red = 2, &
1536 green = 3, &
1537 blue = 4, &
1538 violet = 5, &
1539 cyran = 6, &
1540 magenta = 7, &
1541 freshred = 8, &
1542 tan = 9, &
1543 yellow = 10, &
1544 gray = 11
1545
1546 character(len= 120), intent(in) :: main_title
1547 integer , intent(in) :: nj, nk
1548 real (kind=4), intent(in) :: ps(nj)
1549 real (kind=4), intent(in) :: chi(nk), temp(nk)
1550
1551 character(len=80) :: x_title, y_title
1552 real, dimension(501) :: x, y
1553
1554 integer :: iopt
1555
1556 integer :: plot_type, mmx, mnx, mmy, mny, k, j, mn, magnitude
1557
1558 real :: xb, xe, yb, ye, ymin, ymax, fct
1559
1560 print '("PLOT RATIO OF BALANCE PART: nj,nk", 2I5,a)', nj,nk, main_title(1:40)
1561
1562 plot_type = 1
1563 magnitude = 1
1564 mn = 0
1565
1566 ! 1.0 Surface pressure
1567
1568 mmy = 10
1569 mny = 1
1570
1571 yb = 0.0
1572 ye = 10.0
1573
1574 y_title = '<XB X> / <X X>'
1575
1576 mmx = nj/10 + 1
1577 mnx = 5
1578
1579 xb = 0.0
1580 xe = real(mmx*10)
1581
1582 x_title = ' Y-LATITUDE (GRID)'
1583
1584 do k=1, nj
1585 x(k) = real(k)
1586 end do
1587
1588 do j=1,nj
1589 y(j) = ps(j) * 10.0**magnitude
1590 end do
1591
1592 call da_line_plot(x_title, y_title, main_title, &
1593 x,y,nj, mn,xb,xe,yb,ye,plot_type, &
1594 mmx, mnx, mmy, mny, &
1595 red, 3500, thick_dash, magnitude, 2)
1596
1597 call frame
1598
1599 ! 2.0 Chi and temperature
1600
1601 mmy = 6
1602 mny = 5
1603
1604 yb = 0.0
1605 ye = float(nk)
1606
1607 y_title = 'MODEL LEVEL'
1608
1609 mmx = 10
1610 mnx = 5
1611
1612 xb = 0.0
1613 xe = real(10)
1614
1615 x_title = '<XB X> / <X X>'
1616
1617 do k=1, nk
1618 y(k) = real(k)
1619 end do
1620
1621 do k=1,nk
1622 x(k) = chi(k) * 10.0**magnitude
1623 end do
1624
1625 call da_line_plot(x_title, y_title, main_title, &
1626 x,y,nk, mn,xb,xe,yb,ye,plot_type, &
1627 mmx, mnx, mmy, mny, &
1628 blue, 3500, thick_dash, magnitude, 0)
1629
1630 do k=1,nk
1631 x(k) = temp(k) * 10.0**magnitude
1632 end do
1633
1634 call da_line_plot(x_title, y_title, main_title, &
1635 x,y,nk, mn,xb,xe,yb,ye,plot_type, &
1636 mmx, mnx, mmy, mny, &
1637 red, 3500, thick_dash, magnitude, 0)
1638
1639 call frame
1640
1641 end subroutine da_balance_ratio_plot
1642
1643 subroutine da_covariance_plot (vn, ps, nj, main_title)
1644
1645 implicit none
1646
1647 integer, parameter :: solid_line = 65535, & ! PATTERN = 1111111111111111
1648 thick_dash = 21845, & ! PATTERN = 0101010101010101
1649 thin_dash = 3855, & ! PATTERN = 0000111100001111
1650 solid_like = 31710 ! PATTERN = 0111101111011110
1651
1652 integer, parameter :: black = 0, &
1653 white = 1, &
1654 red = 2, &
1655 green = 3, &
1656 blue = 4, &
1657 violet = 5, &
1658 cyran = 6, &
1659 magenta = 7, &
1660 freshred = 8, &
1661 tan = 9, &
1662 yellow = 10, &
1663 gray = 11
1664
1665 character(len= 5), intent(in) :: vn
1666 character(len= 120), intent(in) :: main_title
1667 integer , intent(in) :: nj
1668 real (kind=4), intent(in) :: ps(nj)
1669
1670 character(len=80) :: x_title, y_title
1671 real, dimension(501) :: x, y
1672
1673 integer :: mmx, mnx, mmy, mny, k, j, magnitude
1674
1675 real :: xb, xe, yb, ye, ymin, ymax, fct
1676
1677 print '("PLOT COVARIANCE: nj",I5,a,2x,a)', nj, vn, main_title(1:40)
1678 magnitude = -3
1679
1680 y_title = trim(vn)//' COVARIANCE'
1681
1682 mmx = nj/10 + 1
1683 xe = real(mmx*10)
1684 x_title = ' Y-LATITUDE (GRID)'
1685
1686 do k=1, nj
1687 x(k) = real(k)
1688 end do
1689
1690 do j=1,nj
1691 y(j) = ps(j) * 10.0**magnitude
1692 end do
1693
1694 call da_line_plot(x_title, y_title, main_title, &
1695 x,y,nj, 0, 0.0, xe, 0.0, 4.0, 1, &
1696 mmx, 5, 4, 5, &
1697 red, 3500, thick_dash, magnitude, 2)
1698
1699 call frame
1700
1701 end subroutine da_covariance_plot
1702
1703 end program da_plot_gen_be