source: palm/trunk/SCRIPTS/NCL/timeseries.ncl @ 2569

Last change on this file since 2569 was 2030, checked in by gronemeier, 8 years ago

bugfix: checking for NCL version number

File size: 52.2 KB
Line 
1load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
3load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
4load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
5
6;***************************************************
7; Checking the kind of the script
8;***************************************************
9
10function which_script()
11local script
12begin
13   script = "timeseries"
14   return(script)
15end
16
17;***************************************************
18; Retrieving the NCL version used
19;***************************************************
20   
21ncl_version_ch = systemfunc("ncl -V")
22ncl_version    = stringtofloat(ncl_version_ch)
23
24;***************************************************
25; Function for checking file existence in dependence
26; on NCL version
27;***************************************************
28
29function file_exists(version:string,file_name:string)
30begin
31   if( version .GE. "6.2.1" ) then
32      existing = fileexists(file_name)
33   else
34      existing = isfilepresent(file_name)
35   end if
36   return(existing)
37end
38 
39;***************************************************
40; load .ncl.config or .ncl.config.default
41;***************************************************
42   
43file_name ="$PALM_BIN/../../.ncl.config"
44existing_file = file_exists(ncl_version_ch,file_name)
45
46if (existing_file) then
47   loadscript("$PALM_BIN/../../.ncl.config")
48else
49   file_name = "$PALM_BIN/NCL/.ncl.config.default"
50   existing_file = file_exists(ncl_version_ch,file_name)
51   if (existing_file) then
52      loadscript( "$PALM_BIN/NCL/.ncl.config.default")
53   else
54      palm_bin_path = getenv("PALM_BIN")
55      print(" ")
56      print("Neither the personal configuration file '.ncl.config' exists in")
57      print("~/palm/current_version")
58      print("nor the default configuration file '.ncl.config.default' "+\
59            "exists in")
60      print(palm_bin_path + "/NCL")
61      print(" ")
62      exit
63   end if
64end if
65
66begin
67
68   ;***************************************************
69   ; Retrieving the double quote character
70   ;***************************************************
71   
72   dq=str_get_dq()
73
74   ;***************************************************
75   ; set up default parameter values and strings
76   ;***************************************************
77   
78   if (file_1 .EQ. "File in") then
79      print(" ")
80      print("Declare input file 'file_1=' in '.ncl.config' or prompt")
81      print(" ")
82      exit
83   else
84      file_in = file_1   
85   end if
86
87   if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND.   \
88       format_out .NE. "eps" .AND. format_out .NE. "ps" .AND.    \
89       format_out .NE. "epsi" .AND. format_out .NE. "ncgm" .AND. \
90       format_out .NE. "png")then
91      print(" ")
92      print("'format_out = "+format_out+"' is invalid and set to'x11'")
93      print(" ")
94      format_out="x11"
95   end if 
96
97   if (ncl_version .LE. 5.1 .AND. format_out .EQ. "png") then
98      print(" ")
99      print("Output of png files not available")
100      print("png output is avaiable with NCL version 5.2.0 and higher ")
101      print("NCL version used: " + ncl_version_ch)
102      print(" ")
103      exit
104   end if
105
106   if (over .NE. 0 .AND. over .NE. 1) then
107      print(" ")
108      print("'over'= "+over+" is invalid and set to 0")
109      print(" ")
110      over = 0
111   end if 
112   
113   if (norm_t .EQ. 0) then
114      print(" ")
115      print("Normalising with 0 is not allowed, 'norm_t' is set to 1.0")
116      print(" ")
117      norm_t = 1.0
118   end if
119 
120
121   ;***************************************************
122   ; open input file
123   ;***************************************************
124   
125   file_in_1 = False
126   if (isStrSubset(file_in, ".nc"))then
127      start_f = -2
128      end_f = -2
129      file_in_1 = True     
130   end if
131
132   if (start_f .EQ. -1)then
133      print(" ")
134      print("'start_f' must be one of the cyclic numbers (at least 0) of "+\
135            "your input file(s)")
136      print(" ") 
137      exit
138   end if
139   if (end_f .EQ. -1)then
140      print(" ")
141      print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
142            "your input file(s)")
143      print(" ") 
144      exit
145   end if
146
147   files=new(end_f-start_f+1,string)
148
149   if (file_in_1)then
150      if (isfilepresent(file_in))then
151         files(0)=file_in
152      else
153         print(" ")
154         print("1st input file: '"+file_in+"' does not exist")
155         print(" ")
156         exit
157      end if
158   else   
159      if (start_f .EQ. 0)then
160         if (isfilepresent(file_in+".nc"))then
161            files(0)=file_in+".nc"
162            do i=1,end_f
163               if (isfilepresent(file_in+"."+i+".nc"))then   
164                  files(i)=file_in+"."+i+".nc"
165               else
166                  print(" ")
167                  print("Input file: '"+file_in+"."+i+".nc' does not exist")
168                  print(" ")
169                  exit 
170               end if       
171            end do         
172         else
173            print(" ")
174            print("Input file: '"+file_in+".nc' does not exist")
175            print(" ")
176            exit
177         end if
178      else
179         do i=start_f,end_f
180            if (isfilepresent(file_in+"."+i+".nc"))then   
181               files(i-start_f)=file_in+"."+i+".nc"
182            else
183               print(" ")
184               print("Input file: '"+file_in+"."+i+".nc' does not exist")
185               print(" ")
186               exit 
187            end if
188         end do
189      end if
190   end if
191
192   f=addfiles(files,"r")
193   f_att=addfile(files(0),"r")
194   ListSetType(f,"cat")
195
196   vNam  =getfilevarnames(f_att)
197   if( ncl_version .GE. 6.1 ) then
198       vNam = vNam(::-1)
199   end if
200   vType = getfilevartypes(f_att,vNam)
201   vNam_static = vNam
202
203   if ((all(vType .eq. "double"))) then ;distinction if data is double or float
204      check_vType = True
205   else
206      check_vType = False
207   end if
208
209   print(" ")
210   print("Variables in input file:")
211   print("- "+ vNam)
212   print(" ")
213   dim = dimsizes(vNam)
214   if (dim .EQ. 0) then
215      print(" ")
216      print("There is no data on file")
217      print(" ")
218   end if
219
220   t_all   = f[:]->time
221   nt      = dimsizes(t_all)
222
223   if (nt .EQ. 1)then
224      delta_t=t_all(nt-1)/nt
225   else
226      delta_t_array = new(nt-1,double)
227
228      do i=0,nt-2
229         delta_t_array(i) = t_all(i+1)-t_all(i)
230      end do
231
232      delta_t = min(delta_t_array)
233      delete(delta_t_array)
234   end if
235
236   if (nt .LE. 1) then
237      print(" ")
238      print("Input file contains only one time step -> " +\
239            "plot of timeseries is not possible ")
240      print(" ")
241      exit
242   end if
243
244   ;****************************************************       
245   ; start of time step and different types of mistakes that could be done
246   ;****************************************************
247
248   if (start_time_step .EQ. -1.) then           
249      start_time_step=t_all(0)/3600     
250   else
251      if (start_time_step .GE. t_all(nt-1)/3600)
252         print(" ")
253         print("'start_time_step' = "+ start_time_step +\
254               "h is equal or greater than last time step = " + \
255               t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
256         print(" ")
257         print("Select another 'start_time_step'")
258         print(" ")
259         exit
260      end if
261      if (start_time_step .LT. t_all(0)/3600)
262         print(" ")
263         print("'start_time_step' = "+ start_time_step +\
264               "h is lower than first time step = " + t_all(0)+"s = "+\
265               t_all(0)/3600+"h")
266         print(" ")
267         print("Select another 'start_time_step'")
268         print(" ")
269         exit
270      end if
271   end if
272   do i=0,nt-2     
273      if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
274          start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
275         st=i
276         break
277       else
278         st=0
279      end if
280   end do
281   if (start_time_step .GE. t_all(nt-1)-delta_t/2 .AND. \
282       start_time_step .LT. t_all(nt-1)) then
283      st=nt-2   
284    else
285      st=0
286   end if
287     
288   ;****************************************************
289   ; end of time step and different types of mistakes that could be done
290   ;****************************************************
291
292   if (end_time_step .EQ. -1.) then             
293      end_time_step = t_all(nt-1)/3600 
294   else
295      if (end_time_step .LE. t_all(0)/3600)
296         print(" ")
297         print("'end_time_step' = "+end_time_step+ \
298              "h is lower or equal than first time step = " + \
299               t_all(0)+"s = "+t_all(0)/3600+"h")
300         print(" ")
301         print("Select another 'end_time_step'")
302         print(" ")
303         exit
304      end if
305      if (end_time_step .GT. t_all(nt-1)/3600)
306         print(" ")
307         print("'end_time_step' = "+ end_time_step +\
308               "h is greater than last time step = " + t_all(nt-1)+"s = "+\
309               t_all(nt-1)/3600+"h")
310         print(" ")
311         print("Select another 'end_time_step'") 
312         print(" ")
313         exit
314      end if
315      if (end_time_step .LE. start_time_step/3600)
316         print(" ")
317         print("'end_time_step' = "+ end_time_step +\
318               "h is equal or lower than 'start_time_step' = "+\
319               start_time_step+"h")
320         print(" ")
321         print("Select another 'start_time_step' or 'end_time_step'")
322         print(" ")
323         exit
324      end if
325   end if
326   do i=0,nt-1     
327      if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
328          end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
329         et=i
330         break
331       else
332         et=0
333      end if
334   end do
335 
336   delete(start_time_step)
337   start_time_step=round(st,3)
338   delete(end_time_step)
339   end_time_step=round(et,3)
340
341   print(" ")
342   print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+\
343         t_all(start_time_step)+" s => index = "+start_time_step)
344   print("                     till "+t_all(end_time_step)/3600+" h = "+\
345         t_all(end_time_step)+" s => index = "+end_time_step)
346   print(" ")
347
348   t = t_all(start_time_step:end_time_step)/norm_t
349 
350   ; ***************************************************
351   ; set up recourses
352   ; ***************************************************
353
354   res                         = True
355   res@gsnDraw                 = False
356   res@gsnFrame                = False
357   res@tmXBMode                = True
358   res@tmYLMode                = True
359   res@txFont                  = "helvetica"
360   res@tiMainFont              = "helvetica"
361   res@tiXAxisFont             = "helvetica"
362   res@tiYAxisFont             = "helvetica"
363   res@tmXBLabelFont           = "helvetica"
364   res@tmYLLabelFont           = "helvetica"
365   res@xyLineColors            = (/237/)
366   res@trXMaxF                 = t(start_time_step)
367   res@trXMinF                 = t(end_time_step)
368   
369   res@lgLabelFontHeightF     = 0.02
370
371   resP                            = True
372   resP@gsnMaximize                = True
373   resP@gsnPanelXWhiteSpacePercent = 4.0
374   resP@gsnPanelYWhiteSpacePercent = 4.0
375   resP@txFont                     = "helvetica"
376   resP@txString                   = f_att@title+" time series "
377   resP@txFuncCode                 = "~"
378   if (format_out .NE. "x11") then
379      if (no_rows .GE. 5) then
380         resP@txFontHeightF        = 0.0105
381      else
382         resP@txFontHeightF        = 0.015
383      end if
384   else
385      resP@txFontHeightF           = 0.015
386   end if
387   res@tmXBMinorPerMajor = 4
388   res@tmYLMinorPerMajor = 4
389
390   res@vpWidthF=4
391
392   txres = True
393
394   ; ***************************************************
395   ; read data and create plots
396   ; ***************************************************
397   
398   if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
399      format_out@wkPaperSize = "A4"
400   end if
401   if ( format_out .EQ. "png" ) then
402      format_out@wkWidth  = 1000
403      format_out@wkHeight = 1000
404   end if
405
406   wks_ps = gsn_open_wks(format_out,file_out)           
407   gsn_define_colormap(wks_ps,"rainbow+white")
408   plot_ps=new(dim,graphic)
409                                 
410   n=0
411   if (check_vType) then
412      minE  =  100000.d
413      maxE  = -100000.d
414      minus =  100000.d
415      maxus = -100000.d
416      minu  =  100000.d
417      maxu  = -100000.d
418      minz  =  100000.d
419      maxz  = -100000.d
420      minw  =  100000.d
421      maxw  = -100000.d
422      minp  =  100000.d
423      maxp  = -100000.d
424      mins  =  100000.d
425      maxs  = -100000.d
426   else
427      minE  =  100000.
428      maxE  = -100000.
429      minus =  100000.
430      maxus = -100000.
431      minu  =  100000.
432      maxu  = -100000.
433      minz  =  100000.
434      maxz  = -100000.
435      minw  =  100000.
436      maxw  = -100000.
437      minp  =  100000.
438      maxp  = -100000.
439      mins  =  100000.
440      maxs  = -100000.
441   end if
442
443   if (check_vType) then
444      data   = new((/dim,(end_time_step-start_time_step)+1/),double)
445      data_0 = new((end_time_step-start_time_step)+1,double)
446      mini   = new(dim,double)
447      maxi   = new(dim,double)
448   else
449      data   = new((/dim,(end_time_step-start_time_step)+1/),float)
450      data_0 = new((end_time_step-start_time_step)+1,float)
451      mini   = new(dim,float)
452      maxi   = new(dim,float)
453   end if
454   unit   = new(dim,string)   
455   data_0 = 0.0
456 
457   
458   if (over .EQ. 1) then
459      plot_E       = gsn_csm_xy(wks_ps,t,data_0(:),res)
460      plot_Es      = gsn_csm_xy(wks_ps,t,data_0(:),res)
461      plot_us      = gsn_csm_xy(wks_ps,t,data_0(:),res)
462      plot_ws      = gsn_csm_xy(wks_ps,t,data_0(:),res)
463      plot_umax    = gsn_csm_xy(wks_ps,t,data_0(:),res)
464      plot_vmax    = gsn_csm_xy(wks_ps,t,data_0(:),res)
465      plot_wmax    = gsn_csm_xy(wks_ps,t,data_0(:),res)
466      plot_z_i_wpt = gsn_csm_xy(wks_ps,t,data_0(:),res)
467      plot_z_i_pt  = gsn_csm_xy(wks_ps,t,data_0(:),res)
468      plot_wpptp0  = gsn_csm_xy(wks_ps,t,data_0(:),res)
469      plot_wpptp   = gsn_csm_xy(wks_ps,t,data_0(:),res)
470      plot_wpt     = gsn_csm_xy(wks_ps,t,data_0(:),res)
471      plot_pt_0_   = gsn_csm_xy(wks_ps,t,data_0(:),res)
472      plot_pt_zp_  = gsn_csm_xy(wks_ps,t,data_0(:),res)
473   end if
474
475
476   if (var .NE. "all") then
477      vNam_temp = new((/dim/),string)
478      comma     = 0
479      var_char  = stringtocharacter(var)
480      no_char   = dimsizes(var_char)
481       
482      do j = 0, no_char-1
483         if(var_char(j) .eq. ",")
484            comma = comma + 1
485         end if
486      end do
487
488      if(comma .eq. 0 .or. comma .eq. 1)
489          print(" ")
490          print("The variables 'var="+var+"'" )
491          print("do not exist on your input file;")
492          print("be sure to have one comma before and after each variable")
493          print(" ")
494          exit
495      end if
496
497      indices = new((/comma/), integer)
498      comma   = 0
499       
500      do j = 0, no_char-1
501         if(var_char(j) .eq. ",")
502           indices(comma) = j
503           comma = comma + 1
504         end if
505      end do
506
507      do j = 0, comma -2
508         vNam_temp(j) = charactertostring(\
509                                   var_char(indices(j)+1:indices(j+1)-1))
510      end do
511     
512      delete(vNam)
513      vNam = new((/comma-1/),string)
514      do varn = 0, comma - 2
515         vNam(varn) = vNam_temp(comma - 2 - varn)
516      end do
517     
518      do j = 0, comma-2
519        count_check = 0
520        do varn = 0, dim-1
521           if(vNam(j) .ne. vNam_static(varn))
522             count_check=count_check+1
523           end if
524        end do   
525
526        if(count_check .eq. dim)
527            print(" ")
528            print("The variables 'var="+var+"'" )
529            print("do not exist on your input file;")
530            print("be sure to have one comma before and after each variable")
531            print(" ")
532            exit
533        end if
534      end do
535
536      dim = comma-1
537    end if
538
539
540   do varn = dim-1, 0, 1
541
542      if( isStrSubset (vNam(varn), "time") )
543         check = False
544      else
545         check = True
546      end if
547     
548      if(check) then
549
550         data_all = f[:]->$vNam(varn)$
551         data_att = f_att->$vNam(varn)$
552         unit(varn) = data_att@units
553 
554     
555         data(varn,:)=data_all(start_time_step:end_time_step)
556         
557         if (over .EQ. 1) then
558
559            mini(varn) = min(data(varn,:))
560            maxi(varn) = max(data(varn,:))
561           
562            if (vNam(varn) .EQ. "E" .OR. vNam(varn) .EQ. "Es" .OR. \
563                vNam(varn) .EQ. "E*") then
564               if (mini(varn) .EQ. maxi(varn)) then
565                  if (min(data(varn,:)) .EQ. 0)then
566                     mini(varn)= mini(varn)-0.1
567                     maxi(varn)= maxi(varn)+0.1
568                  end if
569                  if (min(data(varn,:)) .LT. 0)then
570                     mini(varn)= mini(varn)-1.+(mini(varn))/2
571                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
572                  end if
573                  if (min(data(varn,:)) .GT. 0)then
574                     mini(varn)= mini(varn)-1.-(mini(varn))/2
575                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
576                  end if
577               end if
578               minE=min((/minE,mini(varn)/))
579               maxE=max((/maxE,maxi(varn)/))
580            end if
581
582            if (vNam(varn) .EQ. "us" .OR. vNam(varn) .EQ. "ws"\
583                .OR. vNam(varn) .EQ. "u*" .OR. vNam(varn) .EQ. "w*") then
584               if (mini(varn) .EQ. maxi(varn)) then
585                  if (min(data(varn,:)) .EQ. 0)then
586                     mini(varn)= mini(varn)-0.1
587                     maxi(varn)= maxi(varn)+0.1
588                  end if
589                  if (min(data(varn,:)) .LT. 0)then
590                     mini(varn)= mini(varn)-1.+(mini(varn))/2
591                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
592                  end if
593                  if (min(data(varn,:)) .GT. 0)then
594                     mini(varn)= mini(varn)-1.-(mini(varn))/2
595                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
596                  end if
597               end if
598               minus=min((/minus,mini(varn)/))
599               maxus=max((/maxus,maxi(varn)/))
600            end if
601
602            if (vNam(varn) .EQ. "umax" .OR. vNam(varn) .EQ. "vmax" .OR. \
603                vNam(varn) .EQ. "wmax") then
604               if (mini(varn) .EQ. maxi(varn)) then
605                  if (mini(varn) .EQ. 0)then
606                     mini(varn)= mini(varn)-0.1
607                     maxi(varn)= maxi(varn)+0.1
608                  end if
609                  if (mini(varn) .LT. 0)then
610                     mini(varn)= mini(varn)-1.+(mini(varn))/2
611                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
612                  end if
613                  if (mini(varn) .GT. 0)then
614                     mini(varn)= mini(varn)-1.-(mini(varn))/2
615                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
616                  end if
617               end if
618               minu=min((/minu,mini(varn)/))
619               maxu=max((/maxu,maxi(varn)/))
620            end if
621
622            if (vNam(varn) .EQ. "z_i_wpt" .OR. vNam(varn) .EQ. "z_i_pt") then
623               if (mini(varn) .EQ. maxi(varn)) then
624                  if (min(data(varn,:)) .EQ. 0)then
625                     mini(varn)= mini(varn)-0.1
626                     maxi(varn)= maxi(varn)+0.1
627                  end if
628                  if (min(data(varn,:)) .LT. 0)then
629                     mini(varn)= mini(varn)-1.+(mini(varn))/2
630                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
631                  end if
632                  if (min(data(varn,:)) .GT. 0)then
633                     mini(varn)= mini(varn)-1.-(mini(varn))/2
634                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
635                  end if
636               end if
637               minz=min((/minz,mini(varn)/))
638               maxz=max((/maxz,maxi(varn)/))
639            end if
640
641            if (vNam(varn) .EQ. "wpptp0" .OR. vNam(varn) .EQ. "wpptp" .OR. \
642                vNam(varn) .EQ. "wpt" .OR. vNam(varn) .EQ. \
643                "w"+dq+"pt"+dq+"0" .OR. vNam(varn) .EQ. "w"+dq+"pt"+dq ) then
644               if (mini(varn) .EQ. maxi(varn)) then
645                  if (min(data(varn,:)) .EQ. 0)then
646                     mini(varn)= mini(varn)-0.1
647                     maxi(varn)= maxi(varn)+0.1
648                  end if
649                  if (min(data(varn,:)) .LT. 0)then
650                     mini(varn)= mini(varn)-1.+(mini(varn))/2
651                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
652                  end if
653                  if (min(data(varn,:)) .GT. 0)then
654                     mini(varn)= mini(varn)-1.-(mini(varn))/2
655                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
656                  end if
657               end if
658               minw=min((/minw,mini(varn)/))
659               maxw=max((/maxw,maxi(varn)/))
660            end if
661
662            if (vNam(varn) .EQ. "pt_0_" .OR. vNam(varn) .EQ. "pt_zp_"\
663               .OR. vNam(varn) .EQ. "pt(0)" .OR. vNam(varn) .EQ. "pt(zp)") then
664               if (mini(varn) .EQ. maxi(varn)) then
665                  if (min(data(varn,:)) .EQ. 0)then
666                     mini(varn)= mini(varn)-0.1
667                     maxi(varn)= maxi(varn)+0.1
668                  end if
669                  if (min(data(varn,:)) .LT. 0)then
670                     mini(varn)= mini(varn)-1.+(mini(varn))/2
671                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
672                  end if
673                  if (min(data(varn,:)) .GT. 0)then
674                     mini(varn)= mini(varn)-1.-(mini(varn))/2
675                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
676                  end if
677               end if
678               minp=min((/minp,mini(varn)/))
679               maxp=max((/maxp,maxi(varn)/))
680            end if
681           
682         end if
683      end if
684   end do
685 
686  if (isStrSubset(data@long_name," SR " ))then
687      print(" ")
688      print("If you have outputs of statistic regions you cannot overlay "+\
689            "variables;")
690      print("'over' is set to 0" )
691      print(" ")
692      over = 0
693   end if
694
695   do varn = dim-1, 0, 1
696     
697      if( isStrSubset (vNam(varn), "time") )
698         check = False
699      else
700         check = True
701      end if     
702   
703      if (var .NE. "all") then
704         check = isStrSubset( var,","+vNam(varn)+"," )
705      end if
706 
707      if(check) then
708       
709        if (over .EQ. 1) then
710 
711            res@gsnLeftString       = "overlayed plot"
712            res@gsnRightString      = unit(varn)
713            if (norm_t .NE. 1.)then
714               res@tiXAxisString        = "t ("+unit_t+")"
715            else
716               res@tiXAxisString        = "t (s)"
717            end if
718            res@tiYAxisString = " "
719            res@tiXAxisFontHeightF   = font_size
720            res@txFontHeightF        = font_size
721            res@tiYAxisFontHeightF   = font_size           
722
723            if (vNam(varn) .EQ. "E")
724               E=0
725               res@xyLineColors     = (/237/)
726               res@xyLineLabelFontHeightF = 0.05
727               res@xyLineLabelFontColor   = 237
728               res@trYMaxF           = maxE
729               res@trYMinF           = minE
730               plot_E = gsn_csm_xy(wks_ps,t,data(varn,:),res)
731            end if
732            if (vNam(varn) .EQ. "Es" .OR. vNam(varn) .EQ. "E*")
733               Es=0
734               res@xyLineColors            = (/144/)
735               res@xyLineLabelFontHeightF = 0.05
736               res@xyLineLabelFontColor   = 144
737               plot_Es = gsn_csm_xy(wks_ps,t,data(varn,:),res)
738            end if
739             
740            if (vNam(varn) .EQ. "us" .OR. vNam(varn) .EQ. "u*")
741               us=0
742               res@xyLineColors            = (/237/)
743               res@xyLineLabelFontHeightF = 0.05
744               res@xyLineLabelFontColor   = 237
745               res@trYMaxF           = maxus
746               res@trYMinF           = minus
747               plot_us = gsn_csm_xy(wks_ps,t,data(varn,:),res)
748            end if             
749            if (vNam(varn) .EQ. "ws" .OR. vNam(varn) .EQ. "w*")
750               ws=0
751               res@xyLineColors            = (/144/)
752               res@xyLineLabelFontHeightF = 0.05
753               res@xyLineLabelFontColor   = 144
754               plot_ws = gsn_csm_xy(wks_ps,t,data(varn,:),res)
755            end if
756             
757            if (vNam(varn) .EQ. "umax")
758               u=0
759               res@xyLineColors            = (/237/)
760               res@xyLineLabelFontHeightF = 0.05
761               res@xyLineLabelFontColor   = 237
762               res@trYMaxF           = maxu
763               res@trYMinF           = minu
764               plot_umax = gsn_csm_xy(wks_ps,t,data(varn,:),res)
765            end if
766            if (vNam(varn) .EQ. "vmax")
767               v=0
768               res@xyLineColors            = (/144/)
769               res@xyLineLabelFontHeightF = 0.05
770               res@xyLineLabelFontColor   = 144
771               plot_vmax = gsn_csm_xy(wks_ps,t,data(varn,:),res)
772            end if
773            if (vNam(varn) .EQ. "wmax")
774               w=0
775               res@xyLineColors            = (/80/)
776               res@xyLineLabelFontHeightF = 0.05
777               res@xyLineLabelFontColor   = 80
778               plot_wmax = gsn_csm_xy(wks_ps,t,data(varn,:),res)
779            end if
780   
781            if (vNam(varn) .EQ. "z_i_wpt")
782               zw=0
783               res@xyLineColors            = (/237/)
784               res@xyLineLabelFontHeightF = 0.05
785               res@xyLineLabelFontColor   = 237
786               res@trYMaxF           = maxz
787               res@trYMinF           = minz
788               plot_z_i_wpt = gsn_csm_xy(wks_ps,t,data(varn,:),res)
789            end if
790            if (vNam(varn) .EQ. "z_i_pt")
791               z=0
792               res@xyLineColors            = (/144/) 
793               res@xyLineLabelFontHeightF = 0.05
794               res@xyLineLabelFontColor   = 144
795               plot_z_i_pt = gsn_csm_xy(wks_ps,t,data(varn,:),res)
796            end if
797
798            if (vNam(varn) .EQ. "wpptp0" .OR. \
799                vNam(varn) .EQ. "w"+dq+"pt"+dq+"0" )
800               w0=0
801               res@xyLineColors            = (/237/)
802               res@xyLineLabelFontHeightF = 0.05
803               res@xyLineLabelFontColor   = 237
804               res@trYMaxF           = maxw
805               res@trYMinF           = minw
806               plot_wpptp0 = gsn_csm_xy(wks_ps,t,data(varn,:),res)
807            end if
808            if (vNam(varn) .EQ. "wpptp" .OR. vNam(varn) .EQ. "w"+dq+"pt"+dq)
809               wp=0
810               res@xyLineColors            = (/144/)
811               res@xyLineLabelFontHeightF = 0.05
812               res@xyLineLabelFontColor   = 144
813               plot_wpptp = gsn_csm_xy(wks_ps,t,data(varn,:),res) 
814            end if
815            if (vNam(varn) .EQ. "wpt")
816               wt=0
817               res@xyLineColors            = (/80/)
818               res@xyLineLabelFontHeightF = 0.05
819               res@xyLineLabelFontColor   = 80
820               plot_wpt = gsn_csm_xy(wks_ps,t,data(varn,:),res)
821            end if
822
823            if (vNam(varn) .EQ. "pt_0_" .OR. vNam(varn) .EQ. "pt(0)")
824               p=0
825               res@xyLineColors            = (/237/)
826               res@xyLineLabelFontHeightF = 0.05
827               res@xyLineLabelFontColor   = 237
828               res@trYMaxF           = maxp
829               res@trYMinF           = minp
830               plot_pt_0_ = gsn_csm_xy(wks_ps,t,data(varn,:),res) 
831            end if
832            if (vNam(varn) .EQ. "pt_zp_" .OR. vNam(varn) .EQ. "pt(zp)")
833               pz=0
834               res@xyLineColors            = (/144/)
835               res@xyLineLabelFontHeightF = 0.05
836               res@xyLineLabelFontColor   = 144
837               plot_pt_zp_ = gsn_csm_xy(wks_ps,t,data(varn,:),res)
838            end if
839
840         end if
841      end if
842   end do
843
844
845   do varn = dim-1, 0, 1
846   
847      if( isStrSubset (vNam(varn), "time") )
848         check = False
849      else
850         check = True
851      end if   
852 
853      if (var.NE. "all") then
854         check = isStrSubset( var,","+vNam(varn)+"," )
855      end if
856   
857      if(check) then
858       
859         if (over .EQ. 1) then       
860       
861            if (vNam(varn) .EQ. "E" .AND. Es .NE. 1) then
862               E=1   
863               overlay(plot_E,plot_Es)
864               n=n+1
865               plot_ps(n) = plot_E   
866
867               ; ***************************************************
868               ; legend for overlaid plot
869               ; ***************************************************
870     
871               lgres                    = True
872               lgMonoDashIndex          = False
873               lgres@lgLabelFont        = "helvetica"   
874               lgres@lgLabelFontHeightF = .1           
875               lgres@vpWidthF           = 0.4           
876               lgres@vpHeightF          = 0.4         
877               lgres@lgDashIndexes      = (/0,0,0/)
878               lgres@lgLineColors       = (/237,144,80/)
879               lbid = gsn_create_legend(wks_ps,2,(/"E","E*"/),lgres)       
880
881               amres = True
882               amres@amParallelPosF   = 0.6                 
883               amres@amOrthogonalPosF = -0.2           
884               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)           
885            end if
886            if ((vNam(varn) .EQ. "Es" .OR. vNam(varn) .EQ. "E*") .AND. \
887                E .NE. 1) then
888               Es=1
889               overlay(plot_E,plot_Es)
890               n=n+1
891               plot_ps(n) = plot_E 
892
893               ; ***************************************************
894               ; legend for overlaid plot
895               ; ***************************************************
896     
897               lgres                    = True
898               lgMonoDashIndex          = False
899               lgres@lgLabelFont        = "helvetica"   
900               lgres@lgLabelFontHeightF = .1           
901               lgres@vpWidthF           = 0.4           
902               lgres@vpHeightF          = 0.4         
903               lgres@lgDashIndexes      = (/0,0,0/)
904               lgres@lgLineColors       = (/237,144,80/)
905               lbid = gsn_create_legend(wks_ps,2,(/"E","E*"/),lgres)       
906
907               amres = True
908               amres@amParallelPosF   = 0.6                 
909               amres@amOrthogonalPosF = -0.2           
910               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)         
911            end if
912
913            if ((vNam(varn) .EQ. "us" .OR. vNam(varn) .EQ. "u*") .AND. \
914                 ws .NE. 1) then
915               us=1
916               overlay(plot_us,plot_ws)
917               n=n+1
918               plot_ps(n) = plot_us
919
920               ; ***************************************************
921               ; legend for overlaid plot
922               ; ***************************************************
923     
924               lgres                    = True
925               lgMonoDashIndex          = False
926               lgres@lgLabelFont        = "helvetica"   
927               lgres@lgLabelFontHeightF = .1           
928               lgres@vpWidthF           = 0.4           
929               lgres@vpHeightF          = 0.4         
930               lgres@lgDashIndexes      = (/0,0,0/)
931               lgres@lgLineColors       = (/237,144,80/)
932               lbid = gsn_create_legend(wks_ps,2,(/"u*","w*"/),lgres)       
933
934               amres = True
935               amres@amParallelPosF   = 0.6                 
936               amres@amOrthogonalPosF = -0.2           
937               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
938            end if
939            if ((vNam(varn) .EQ. "ws" .OR. vNam(varn) .EQ. "w*") .AND. \
940                 us .NE. 1) then
941               ws=1
942               overlay(plot_us,plot_ws)
943               n=n+1
944               plot_ps(n) = plot_us
945
946               ; ***************************************************
947               ; legend for overlaid  plot
948               ; ***************************************************
949     
950               lgres                    = True
951               lgMonoDashIndex          = False
952               lgres@lgLabelFont        = "helvetica"   
953               lgres@lgLabelFontHeightF = .1           
954               lgres@vpWidthF           = 0.4           
955               lgres@vpHeightF          = 0.4         
956               lgres@lgDashIndexes      = (/0,0,0/)
957               lgres@lgLineColors       = (/237,144,80/)
958               lbid = gsn_create_legend(wks_ps,2,(/"u*","w*"/),lgres)       
959
960               amres = True
961               amres@amParallelPosF   = 0.6                 
962               amres@amOrthogonalPosF = -0.2           
963               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
964            end if
965         
966            if (vNam(varn) .EQ. "umax" .AND. v .NE. 1)
967               if (w .NE. 1) then
968                  u=1         
969                  overlay(plot_umax,plot_vmax)
970                  overlay(plot_umax,plot_wmax)
971               
972                  n=n+1
973                  plot_ps(n) = plot_umax
974
975                  ; ***************************************************
976                  ; legend for overlaid plot
977                  ; ***************************************************
978     
979                  lgres                    = True
980                  lgMonoDashIndex          = False
981                  lgres@lgLabelFont        = "helvetica"   
982                  lgres@lgLabelFontHeightF = .1           
983                  lgres@vpWidthF           = 0.4           
984                  lgres@vpHeightF          = 0.4         
985                  lgres@lgDashIndexes      = (/0,0,0/)
986                  lgres@lgLineColors       = (/237,144,80/)
987                  lbid = gsn_create_legend(\
988                                  wks_ps,3,(/"umax","vmax","wmax"/),lgres)
989
990                  amres = True
991                  amres@amParallelPosF   = 0.6             
992                  amres@amOrthogonalPosF = -0.2           
993                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
994               end if
995            end if
996            if (vNam(varn) .EQ. "vmax" .AND. u .NE. 1)
997               if (w .NE. 1) then
998                  v=1 
999                  overlay(plot_umax,plot_vmax)
1000                  overlay(plot_umax,plot_wmax)
1001                  n=n+1
1002                  plot_ps(n) = plot_umax
1003
1004                  ; ***************************************************
1005                  ; legend for overlaid plot
1006                  ; ***************************************************
1007     
1008                  lgres                    = True
1009                  lgMonoDashIndex          = False
1010                  lgres@lgLabelFont        = "helvetica"   
1011                  lgres@lgLabelFontHeightF = .1         
1012                  lgres@vpWidthF           = 0.4           
1013                  lgres@vpHeightF          = 0.4         
1014                  lgres@lgDashIndexes      = (/0,0,0/)
1015                  lgres@lgLineColors       = (/237,144,80/)
1016                  lbid = gsn_create_legend(\
1017                                      wks_ps,3,(/"umax","vmax","wmax"/),lgres)
1018                  amres = True
1019                  amres@amParallelPosF   = 0.6             
1020                  amres@amOrthogonalPosF = -0.2           
1021                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
1022               end if
1023            end if
1024            if (vNam(varn) .EQ. "wmax" .AND. v .NE. 1)
1025               if(u .NE. 1) then 
1026                  w=1       
1027                  overlay(plot_umax,plot_vmax)
1028                  overlay(plot_umax,plot_wmax)
1029                  n=n+1
1030                  plot_ps(n) = plot_umax
1031
1032                  ; ***************************************************
1033                  ; legend for overlaid plot
1034                  ; ***************************************************
1035     
1036                  lgres                    = True
1037                  lgMonoDashIndex          = False
1038                  lgres@lgLabelFont        = "helvetica"   
1039                  lgres@lgLabelFontHeightF = .1           
1040                  lgres@vpWidthF           = 0.4           
1041                  lgres@vpHeightF          = 0.4         
1042                  lgres@lgDashIndexes      = (/0,0,0/)
1043                  lgres@lgLineColors       = (/237,144,80/)
1044                  lbid = gsn_create_legend(\
1045                                       wks_ps,3,(/"umax","vmax","wmax"/),lgres)
1046                  amres = True
1047                  amres@amParallelPosF   = 0.6             
1048                  amres@amOrthogonalPosF = -0.2           
1049                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
1050               end if
1051            end if
1052           
1053            if (vNam(varn) .EQ. "z_i_wpt" .AND. z .NE. 1) then
1054               zw=1       
1055               overlay(plot_z_i_wpt,plot_z_i_pt)
1056               n=n+1
1057               plot_ps(n) = plot_z_i_wpt
1058       
1059               ; ***************************************************
1060               ; legend for overlaid plot
1061               ; ***************************************************
1062     
1063               lgres                    = True
1064               lgMonoDashIndex          = False
1065               lgres@lgLabelFont        = "helvetica"   
1066               lgres@lgLabelFontHeightF = .1           
1067               lgres@vpWidthF           = 0.4           
1068               lgres@vpHeightF          = 0.4         
1069               lgres@lgDashIndexes      = (/0,0,0/)
1070               lgres@lgLineColors       = (/237,144,80/)
1071               lbid = gsn_create_legend(wks_ps,2,(/"z_i_wpt","z_i_pt"/),lgres)
1072
1073               amres = True
1074               amres@amParallelPosF   = 0.6                 
1075               amres@amOrthogonalPosF = -0.2           
1076               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)         
1077            end if
1078            if (vNam(varn) .EQ. "z_i_pt" .AND. zw .NE. 1) then
1079               z=1     
1080               overlay(plot_z_i_wpt,plot_z_i_pt)
1081               n=n+1
1082               plot_ps(n) = plot_z_i_wpt
1083
1084               ; ***************************************************
1085               ; legend for overlaid plot
1086               ; ***************************************************
1087     
1088               lgres                    = True
1089               lgMonoDashIndex          = False
1090               lgres@lgLabelFont        = "helvetica"   
1091               lgres@lgLabelFontHeightF = .1           
1092               lgres@vpWidthF           = 0.4           
1093               lgres@vpHeightF          = 0.4         
1094               lgres@lgDashIndexes      = (/0,0,0/)
1095               lgres@lgLineColors       = (/237,144,80/)
1096               lbid = gsn_create_legend(wks_ps,2,(/"z_i_wpt","z_i_pt"/),lgres)
1097
1098               amres = True
1099               amres@amParallelPosF   = 0.6                 
1100               amres@amOrthogonalPosF = -0.2           
1101               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)           
1102            end if   
1103         
1104            if ((vNam(varn) .EQ. "wpptp0" .OR. vNam(varn) .EQ. "w"+dq+"pt"+dq+"0") .AND. wp .NE. 1)
1105               if (wt .NE. 1) then
1106                  w0=1
1107                  overlay(plot_wpptp0,plot_wpptp)
1108                  overlay(plot_wpptp0,plot_wpt)
1109                  n=n+1
1110                  plot_ps(n) = plot_wpptp0
1111
1112                  ; ***************************************************
1113                  ; legend for overlaid plot
1114                  ; ***************************************************
1115     
1116                  lgres                    = True
1117                  lgMonoDashIndex          = False
1118                  lgres@lgLabelFont        = "helvetica"   
1119                  lgres@lgLabelFontHeightF = .1           
1120                  lgres@vpWidthF           = 0.4           
1121                  lgres@vpHeightF          = 0.4         
1122                  lgres@lgDashIndexes      = (/0,0,0/)
1123                  lgres@lgLineColors       = (/237,144,80/)
1124                  lbid = gsn_create_legend(\
1125                                    wks_ps,3,(/"wpptp0","wpptp","wpt"/),lgres)
1126                  amres = True
1127                  amres@amParallelPosF   = 0.6             
1128                  amres@amOrthogonalPosF = -0.2           
1129                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
1130               end if           
1131            end if
1132            if ((vNam(varn) .EQ. "wpptp" .OR. \
1133                 vNam(varn) .EQ. "w"+dq+"pt"+dq) .AND. w0 .NE. 1)
1134               if (wt .NE. 1) then
1135                 wp=1
1136                 overlay(plot_wpptp0,plot_wpptp)
1137                 overlay(plot_wpptp0,plot_wpt)
1138                 n=n+1
1139                 plot_ps(n) = plot_wpptp0
1140                 
1141                 ; ***************************************************
1142                  ; legend for overlaid plot
1143                  ; ***************************************************
1144     
1145                  lgres                    = True
1146                  lgMonoDashIndex          = False
1147                  lgres@lgLabelFont        = "helvetica"   
1148                  lgres@lgLabelFontHeightF = .1           
1149                  lgres@vpWidthF           = 0.4           
1150                  lgres@vpHeightF          = 0.4         
1151                  lgres@lgDashIndexes      = (/0,0,0/)
1152                  lgres@lgLineColors       = (/237,144,80/)
1153                  lbid = gsn_create_legend(\
1154                                     wks_ps,3,(/"wpptp0","wpptp","wpt"/),lgres)
1155
1156                  amres = True
1157                  amres@amParallelPosF   = 0.6             
1158                  amres@amOrthogonalPosF = -0.2           
1159                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres) 
1160               end if           
1161            end if
1162            if (vNam(varn) .EQ. "wpt" .AND. wp .NE. 1)
1163               if (w0 .NE. 1) then
1164                  wt=1
1165                  overlay(plot_wpptp0,plot_wpptp)
1166                  overlay(plot_wpptp0,plot_wpt)
1167                  n=n+1
1168                  plot_ps(n) = plot_wpptp0
1169
1170                  ; ***************************************************
1171                  ; legend for overlaid plot
1172                  ; ***************************************************
1173     
1174                  lgres                    = True
1175                  lgMonoDashIndex          = False
1176                  lgres@lgLabelFont        = "helvetica"   
1177                  lgres@lgLabelFontHeightF = .1           
1178                  lgres@vpWidthF           = 0.4           
1179                  lgres@vpHeightF          = 0.4         
1180                  lgres@lgDashIndexes      = (/0,0,0/)
1181                  lgres@lgLineColors       = (/237,144,80/)
1182                  lbid = gsn_create_legend(\
1183                                  wks_ps,3,(/"wpptp0","wpptp","wpt"/),lgres)
1184                  amres = True
1185                  amres@amParallelPosF   = 0.6             
1186                  amres@amOrthogonalPosF = -0.2           
1187                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
1188               end if
1189            end if
1190
1191            if ((vNam(varn) .EQ. "pt_0_" .OR. vNam(varn) .EQ. "pt(0)") .AND. \
1192                pz .NE. 1) then
1193               p=1     
1194               overlay(plot_pt_0_,plot_pt_zp_)
1195               n=n+1
1196               plot_ps(n) = plot_pt_0_ 
1197     
1198               ; ***************************************************
1199               ; legend for overlaid plot
1200               ; ***************************************************
1201     
1202               lgres                    = True
1203               lgMonoDashIndex          = False
1204               lgres@lgLabelFont        = "helvetica"   
1205               lgres@lgLabelFontHeightF = .1           
1206               lgres@vpWidthF           = 0.4           
1207               lgres@vpHeightF          = 0.4         
1208               lgres@lgDashIndexes      = (/0,0,0/)
1209               lgres@lgLineColors       = (/237,144,80/)
1210               lbid = gsn_create_legend(wks_ps,2,(/"pt_0_","pt_zp_"/),lgres)
1211
1212               amres = True
1213               amres@amParallelPosF   = 0.6                 
1214               amres@amOrthogonalPosF = -0.2           
1215               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
1216            end if
1217            if ((vNam(varn) .EQ. "pt_zp_" .OR. vNam(varn) .EQ. "pt(zp)") .AND.\
1218                 p .NE. 1) then
1219               pz=1       
1220               overlay(plot_pt_0_,plot_pt_zp_)
1221               n=n+1
1222               plot_ps(n) = plot_pt_0_   
1223
1224               ; ***************************************************
1225               ; legend for overlaid plot
1226               ; ***************************************************
1227     
1228               lgres                    = True
1229               lgMonoDashIndex          = False
1230               lgres@lgLabelFont        = "helvetica"   
1231               lgres@lgLabelFontHeightF = .1           
1232               lgres@vpWidthF           = 0.4           
1233               lgres@vpHeightF          = 0.4         
1234               lgres@lgDashIndexes      = (/0,0,0/)
1235               lgres@lgLineColors       = (/237,144,80/)
1236               lbid = gsn_create_legend(wks_ps,2,(/"pt_0_","pt_zp_"/),lgres)
1237
1238               amres = True
1239               amres@amParallelPosF   = 0.6                 
1240               amres@amOrthogonalPosF = -0.2           
1241               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
1242            end if
1243
1244            if(vNam(varn) .NE. "pt_zp_" .AND. vNam(varn) .NE. "pt(zp)" .AND. \
1245               vNam(varn) .NE. "pt_0_" .AND. vNam(varn) .NE. "pt(0)" .AND.   \
1246               vNam(varn) .NE. "wpt" .AND. vNam(varn) .NE. "wpptp" .AND.     \
1247               vNam(varn) .NE. "w"+dq+"pt"+dq .AND.                          \
1248               vNam(varn) .NE. "wpptp0" .AND.                                \
1249               vNam(varn) .NE. "w"+dq+"pt"+dq+"0" .AND.                      \
1250               vNam(varn) .NE. "z_i_pt" .AND.                                \
1251               vNam(varn) .NE. "z_i_wpt" .AND. vNam(varn) .NE. "wmax" .AND.  \
1252               vNam(varn) .NE. "vmax" .AND. vNam(varn) .NE. "umax" .AND.     \
1253               vNam(varn) .NE. "ws" .AND. vNam(varn) .NE. "w*" .AND.         \
1254               vNam(varn) .NE. "us" .AND. vNam(varn) .NE. "u*" .AND.         \
1255               vNam(varn) .NE. "Es" .AND. vNam(varn) .NE. "E*" .AND.         \
1256               vNam(varn) .NE. "E") then
1257
1258               n=n+1
1259               res@xyLineColors   = (/237/)
1260               res@xyLabelMode    = False
1261               res@gsnLeftString  = vNam(varn)
1262               res@gsnRightString = unit(varn)
1263               res@trYMaxF        = max(data(varn,:))
1264               res@trYMinF        = min(data(varn,:))
1265               if (min(data(varn,:)) .EQ. max(data(varn,:))) then
1266                  if (min(data(varn,:)) .EQ. 0)then
1267                     res@trYMaxF = max(data(varn,:))-0.1
1268                     res@trYMinF = min(data(varn,:))+0.1
1269                  end if
1270                  if (min(data(varn,:)) .LT. 0)then
1271                     res@trYMaxF = max(data(varn,:))+(min(data(varn,:)))/2
1272                     res@trYMinF = min(data(varn,:))-(max(data(varn,:)))/2
1273                  end if
1274                  if (min(data(varn,:)) .GT. 0)then
1275                     res@trYMaxF = max(data(varn,:))-(min(data(varn,:)))/2
1276                     res@trYMinF = min(data(varn,:))+(max(data(varn,:)))/2
1277                  end if
1278               end if
1279               plot_ps(n) = gsn_csm_xy(wks_ps,t,data(varn,:),res) 
1280            end if
1281
1282            ; ***************************************************
1283            ; transparent legend as placeholder for single plots
1284            ; ***************************************************
1285            lgresmono                    = True
1286            lgMonoDashIndex              = False
1287            lgresmono@lgMonoLineColor    = True
1288            lgresmono@lgLabelFont        = "helvetica"   
1289            lgresmono@lgLabelFontHeightF = .1           
1290            lgresmono@vpWidthF           = 0.4           
1291            lgresmono@vpHeightF          = 0.4         
1292            lgresmono@lgLineColor        = 0
1293            lgresmono@lgPerimOn          = False
1294            lbid = gsn_create_legend(wks_ps,1,(/""/),lgresmono)       
1295
1296            amres = True
1297            amres@amParallelPosF   = 0.6                   
1298            amres@amOrthogonalPosF = -0.2           
1299            annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
1300       
1301         else
1302           
1303            print("plot of " + vNam(varn))
1304           
1305            n=n+1
1306            res@xyLineColors            = (/237/)
1307            res@gsnLeftString       = vNam(varn)
1308            res@gsnRightString      = unit(varn)
1309       
1310            if (norm_t .NE. 1.)then
1311               res@tiXAxisString        = "t ("+unit_t+")"
1312            else
1313               res@tiXAxisString        = "t (s)"
1314            end if
1315            res@tiYAxisString = " "
1316            res@tiXAxisFontHeightF   = font_size
1317            res@txFontHeightF        = font_size
1318            res@tiYAxisFontHeightF   = font_size
1319            res@tmXBLabelFontHeightF = font_size
1320            res@tmYLLabelFontHeightF = font_size
1321            res@trYMaxF        = max(data(varn,:))
1322            res@trYMinF        = min(data(varn,:))
1323            if (min(data(varn,:)) .EQ. max(data(varn,:))) then
1324               if (min(data(varn,:)) .EQ. 0)then
1325                     res@trYMaxF = max(data(varn,:))-0.1
1326                     res@trYMinF = min(data(varn,:))+0.1
1327                  end if
1328                  if (min(data(varn,:)) .LT. 0)then
1329                     res@trYMaxF = max(data(varn,:))+(min(data(varn,:)))/2
1330                     res@trYMinF = min(data(varn,:))-(max(data(varn,:)))/2
1331                  end if
1332                  if (min(data(varn,:)) .GT. 0)then
1333                     res@trYMaxF = max(data(varn,:))-(min(data(varn,:)))/2
1334                     res@trYMinF = min(data(varn,:))+(max(data(varn,:)))/2
1335                  end if
1336            end if
1337            plot_ps(n) = gsn_csm_xy(wks_ps,t,data(varn,:),res)
1338           
1339         end if     
1340      end if
1341   end do
1342 
1343   ; ***************************************************
1344   ; merge plots onto one page
1345   ; ***************************************************
1346
1347   no_frames = 0
1348 
1349   if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
1350        n .gt. no_rows*no_columns) then
1351      gsn_panel(wks_ps,plot_ps(1:n),(/n,1/),resP)
1352      print(" ")
1353      print("Outputs to .eps or .epsi have only one frame")
1354      print(" ")
1355   else
1356      do np = 1,n,no_rows*no_columns   
1357         if ( np + no_rows*no_columns .gt. n) then   
1358            gsn_panel(wks_ps, plot_ps(np:n),(/no_rows,no_columns/),resP)
1359            no_frames = no_frames + 1           
1360         else
1361            gsn_panel(wks_ps, plot_ps(np:np+no_rows*no_columns-1),\
1362                                                   (/no_rows,no_columns/),resP)
1363            no_frames = no_frames + 1
1364         end if
1365      end do
1366   end if
1367
1368   if (format_out .EQ. "png" ) then
1369     png_output = new((/no_frames/), string)
1370     j = 0
1371
1372     if (no_frames .eq. 1) then
1373        if (ncl_version .GE. 6 ) then
1374           png_output(0) = file_out+".png"
1375        else
1376           png_output(0) = file_out+".00000"+1+".png"
1377        end if
1378        ;using imagemagick's convert for reducing the white
1379        ;space around the plot
1380        cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
1381              png_output(0) + " " + png_output(0)
1382       system(cmd)
1383     else
1384
1385       do i=0, no_frames-1
1386         j = i + 1
1387         if (j .LE. 9) then
1388           png_output(i) = file_out+".00000"+j+".png"
1389         end if
1390         if (j .GT. 9 .AND. j .LE. 99) then
1391           png_output(i) = file_out+".0000"+j+".png"
1392         end if
1393         if (j .GT. 99 .AND. j .LE. 999) then
1394           png_output(i) = file_out+".000"+j+".png"
1395         end if
1396         if (j .GT. 999) then
1397           png_output(i) = file_out+".00"+j+".png"
1398         end if
1399
1400         ;using imagemagick's convert for reducing the white
1401         ;space around the plot
1402         cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
1403                png_output(i) + " " + png_output(i)
1404         system(cmd)
1405       end do
1406     end if
1407
1408     print(" ")
1409     print("Output to: "+ png_output)
1410     print(" ")
1411   else
1412     print(" ")
1413     print("Output to: " + file_out +"."+ format_out)
1414     print(" ")
1415   end if
1416
1417end
Note: See TracBrowser for help on using the repository browser.