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

Last change on this file since 1746 was 1559, checked in by heinze, 10 years ago

Changes to allow for using NCL version 6.2.1 and higher. Backward compatibility is also ensured.

File size: 52.2 KB
RevLine 
[194]1load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
[154]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"
[1559]5
[190]6;***************************************************
[1559]7; Checking the kind of the script
[190]8;***************************************************
[194]9
10function which_script()
11local script
12begin
[1559]13   script = "timeseries"
[194]14   return(script)
15end
[418]16
[1559]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 .EQ. "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
[418]47   loadscript("$PALM_BIN/../../.ncl.config")
[190]48else
[1559]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
[418]54      palm_bin_path = getenv("PALM_BIN")
[190]55      print(" ")
[418]56      print("Neither the personal configuration file '.ncl.config' exists in")
57      print("~/palm/current_version")
[566]58      print("nor the default configuration file '.ncl.config.default' "+\
59            "exists in")
[418]60      print(palm_bin_path + "/NCL")
[190]61      print(" ")
62      exit
[154]63   end if
[1559]64end if
65
[190]66begin
[154]67
[190]68   ;***************************************************
[526]69   ; Retrieving the double quote character
70   ;***************************************************
71   
72   dq=str_get_dq()
73
74   ;***************************************************
[190]75   ; set up default parameter values and strings
76   ;***************************************************
[218]77   
[190]78   if (file_1 .EQ. "File in") then
79      print(" ")
[418]80      print("Declare input file 'file_1=' in '.ncl.config' or prompt")
[190]81      print(" ")
82      exit
[175]83   else
84      file_in = file_1   
[154]85   end if
[175]86
[566]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
[190]91      print(" ")
92      print("'format_out = "+format_out+"' is invalid and set to'x11'")
93      print(" ")
94      format_out="x11"
[154]95   end if 
[175]96
[532]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
[190]106   if (over .NE. 0 .AND. over .NE. 1) then
107      print(" ")
108      print("'over'= "+over+" is invalid and set to 0")
109      print(" ")
[161]110      over = 0
[218]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
[190]119 
[175]120
[190]121   ;***************************************************
[154]122   ; open input file
[190]123   ;***************************************************
[218]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
[154]131
[218]132   if (start_f .EQ. -1)then
133      print(" ")
[566]134      print("'start_f' must be one of the cyclic numbers (at least 0) of "+\
135            "your input file(s)")
[218]136      print(" ") 
137      exit
138   end if
139   if (end_f .EQ. -1)then
140      print(" ")
[566]141      print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
142            "your input file(s)")
[218]143      print(" ") 
144      exit
145   end if
[154]146
[218]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
[585]196   vNam  =getfilevarnames(f_att)
[1248]197   if( ncl_version .GE. 6.1 ) then
198       vNam = vNam(::-1)
199   end if
[585]200   vType = getfilevartypes(f_att,vNam)
201   vNam_static = vNam
[218]202
[585]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
[154]209   print(" ")
[190]210   print("Variables in input file:")
211   print("- "+ vNam)
[154]212   print(" ")
[161]213   dim = dimsizes(vNam)
214   if (dim .EQ. 0) then
215      print(" ")
[310]216      print("There is no data on file")
[161]217      print(" ")
218   end if
219
[310]220   t_all   = f[:]->time
221   nt      = dimsizes(t_all)
[218]222
[827]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
[769]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
[190]244   ;****************************************************       
[154]245   ; start of time step and different types of mistakes that could be done
[190]246   ;****************************************************
[154]247
[190]248   if (start_time_step .EQ. -1.) then           
249      start_time_step=t_all(0)/3600     
[154]250   else
[162]251      if (start_time_step .GE. t_all(nt-1)/3600)
[154]252         print(" ")
[566]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")
[154]256         print(" ")
[218]257         print("Select another 'start_time_step'")
[162]258         print(" ")
[154]259         exit
260      end if
[162]261      if (start_time_step .LT. t_all(0)/3600)
[154]262         print(" ")
[566]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")
[154]266         print(" ")
[218]267         print("Select another 'start_time_step'")
[162]268         print(" ")
[154]269         exit
270      end if
271   end if
[162]272   do i=0,nt-2     
[566]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
[162]275         st=i
276         break
[769]277       else
278         st=0
[162]279      end if
280   end do
[566]281   if (start_time_step .GE. t_all(nt-1)-delta_t/2 .AND. \
282       start_time_step .LT. t_all(nt-1)) then
[162]283      st=nt-2   
[769]284    else
285      st=0
[162]286   end if
287     
[190]288   ;****************************************************
[154]289   ; end of time step and different types of mistakes that could be done
[190]290   ;****************************************************
[154]291
[190]292   if (end_time_step .EQ. -1.) then             
293      end_time_step = t_all(nt-1)/3600 
[154]294   else
[162]295      if (end_time_step .LE. t_all(0)/3600)
[154]296         print(" ")
[566]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")
[154]300         print(" ")
[218]301         print("Select another 'end_time_step'")
[162]302         print(" ")
[154]303         exit
304      end if
[162]305      if (end_time_step .GT. t_all(nt-1)/3600)
[154]306         print(" ")
[566]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")
[154]310         print(" ")
[218]311         print("Select another 'end_time_step'") 
[162]312         print(" ")
[154]313         exit
314      end if
[162]315      if (end_time_step .LE. start_time_step/3600)
[161]316         print(" ")
[566]317         print("'end_time_step' = "+ end_time_step +\
318               "h is equal or lower than 'start_time_step' = "+\
319               start_time_step+"h")
[161]320         print(" ")
[218]321         print("Select another 'start_time_step' or 'end_time_step'")
[162]322         print(" ")
[161]323         exit
[154]324      end if
[162]325   end if
326   do i=0,nt-1     
[566]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
[162]329         et=i
330         break
[769]331       else
332         et=0
[162]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)
[154]340
[175]341   print(" ")
[566]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)
[175]346   print(" ")
347
[218]348   t = t_all(start_time_step:end_time_step)/norm_t
[154]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"
[161]364   res@tmYLLabelFont           = "helvetica"
365   res@xyLineColors            = (/237/)
[758]366   res@trXMaxF                 = t(start_time_step)
367   res@trXMinF                 = t(end_time_step)
[161]368   
[218]369   res@lgLabelFontHeightF     = 0.02
[161]370
[983]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
[218]387   res@tmXBMinorPerMajor = 4
388   res@tmYLMinorPerMajor = 4
389
[154]390   res@vpWidthF=4
391
[161]392   txres = True
[154]393
394   ; ***************************************************
395   ; read data and create plots
396   ; ***************************************************
397   
[532]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
[154]406   wks_ps = gsn_open_wks(format_out,file_out)           
407   gsn_define_colormap(wks_ps,"rainbow+white")
[161]408   plot_ps=new(dim,graphic)
409                                 
[154]410   n=0
[585]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
[161]442
[585]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)   
[161]455   data_0 = 0.0
[585]456 
[161]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
[162]474
[161]475
[310]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
[317]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
[310]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
[566]508         vNam_temp(j) = charactertostring(\
509                                   var_char(indices(j)+1:indices(j+1)-1))
[310]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
[317]524        end do   
525
[310]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
[161]542      if( isStrSubset (vNam(varn), "time") )
543         check = False
544      else
545         check = True
[310]546      end if
[161]547     
548      if(check) then
[162]549
[218]550         data_all = f[:]->$vNam(varn)$
551         data_att = f_att->$vNam(varn)$
552         unit(varn) = data_att@units
[310]553 
554     
[161]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           
[566]562            if (vNam(varn) .EQ. "E" .OR. vNam(varn) .EQ. "Es" .OR. \
563                vNam(varn) .EQ. "E*") then
[161]564               if (mini(varn) .EQ. maxi(varn)) then
565                  if (min(data(varn,:)) .EQ. 0)then
[218]566                     mini(varn)= mini(varn)-0.1
567                     maxi(varn)= maxi(varn)+0.1
[161]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
[526]582            if (vNam(varn) .EQ. "us" .OR. vNam(varn) .EQ. "ws"\
583                .OR. vNam(varn) .EQ. "u*" .OR. vNam(varn) .EQ. "w*") then
[161]584               if (mini(varn) .EQ. maxi(varn)) then
585                  if (min(data(varn,:)) .EQ. 0)then
[218]586                     mini(varn)= mini(varn)-0.1
587                     maxi(varn)= maxi(varn)+0.1
[161]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
[566]602            if (vNam(varn) .EQ. "umax" .OR. vNam(varn) .EQ. "vmax" .OR. \
603                vNam(varn) .EQ. "wmax") then
[161]604               if (mini(varn) .EQ. maxi(varn)) then
605                  if (mini(varn) .EQ. 0)then
[218]606                     mini(varn)= mini(varn)-0.1
607                     maxi(varn)= maxi(varn)+0.1
[161]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
[218]625                     mini(varn)= mini(varn)-0.1
626                     maxi(varn)= maxi(varn)+0.1
[161]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
[566]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
[161]644               if (mini(varn) .EQ. maxi(varn)) then
645                  if (min(data(varn,:)) .EQ. 0)then
[218]646                     mini(varn)= mini(varn)-0.1
647                     maxi(varn)= maxi(varn)+0.1
[161]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
[526]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
[161]664               if (mini(varn) .EQ. maxi(varn)) then
665                  if (min(data(varn,:)) .EQ. 0)then
[218]666                     mini(varn)= mini(varn)-0.1
667                     maxi(varn)= maxi(varn)+0.1
[161]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
[310]685 
686  if (isStrSubset(data@long_name," SR " ))then
[162]687      print(" ")
[566]688      print("If you have outputs of statistic regions you cannot overlay "+\
689            "variables;")
[310]690      print("'over' is set to 0" )
[162]691      print(" ")
[310]692      over = 0
[175]693   end if
[162]694
[310]695   do varn = dim-1, 0, 1
[161]696     
[154]697      if( isStrSubset (vNam(varn), "time") )
698         check = False
[161]699      else
700         check = True
701      end if     
[190]702   
703      if (var .NE. "all") then
[161]704         check = isStrSubset( var,","+vNam(varn)+"," )
705      end if
[154]706 
707      if(check) then
[310]708       
[161]709        if (over .EQ. 1) then
[154]710 
[161]711            res@gsnLeftString       = "overlayed plot"
712            res@gsnRightString      = unit(varn)
[218]713            if (norm_t .NE. 1.)then
[534]714               res@tiXAxisString        = "t ("+unit_t+")"
[218]715            else
[534]716               res@tiXAxisString        = "t (s)"
[218]717            end if
718            res@tiYAxisString = " "
719            res@tiXAxisFontHeightF   = font_size
720            res@txFontHeightF        = font_size
721            res@tiYAxisFontHeightF   = font_size           
[154]722
[161]723            if (vNam(varn) .EQ. "E")
724               E=0
725               res@xyLineColors     = (/237/)
726               res@xyLineLabelFontHeightF = 0.05
727               res@xyLineLabelFontColor   = 237
[310]728               res@trYMaxF           = maxE
729               res@trYMinF           = minE
[161]730               plot_E = gsn_csm_xy(wks_ps,t,data(varn,:),res)
731            end if
[526]732            if (vNam(varn) .EQ. "Es" .OR. vNam(varn) .EQ. "E*")
[161]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             
[526]740            if (vNam(varn) .EQ. "us" .OR. vNam(varn) .EQ. "u*")
[161]741               us=0
742               res@xyLineColors            = (/237/)
743               res@xyLineLabelFontHeightF = 0.05
744               res@xyLineLabelFontColor   = 237
[310]745               res@trYMaxF           = maxus
746               res@trYMinF           = minus
[161]747               plot_us = gsn_csm_xy(wks_ps,t,data(varn,:),res)
748            end if             
[526]749            if (vNam(varn) .EQ. "ws" .OR. vNam(varn) .EQ. "w*")
[161]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
[310]762               res@trYMaxF           = maxu
763               res@trYMinF           = minu
[161]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
[154]780   
[161]781            if (vNam(varn) .EQ. "z_i_wpt")
782               zw=0
783               res@xyLineColors            = (/237/)
784               res@xyLineLabelFontHeightF = 0.05
785               res@xyLineLabelFontColor   = 237
[310]786               res@trYMaxF           = maxz
787               res@trYMinF           = minz
[161]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
[566]798            if (vNam(varn) .EQ. "wpptp0" .OR. \
799                vNam(varn) .EQ. "w"+dq+"pt"+dq+"0" )
[161]800               w0=0
801               res@xyLineColors            = (/237/)
802               res@xyLineLabelFontHeightF = 0.05
803               res@xyLineLabelFontColor   = 237
[310]804               res@trYMaxF           = maxw
805               res@trYMinF           = minw
[161]806               plot_wpptp0 = gsn_csm_xy(wks_ps,t,data(varn,:),res)
807            end if
[526]808            if (vNam(varn) .EQ. "wpptp" .OR. vNam(varn) .EQ. "w"+dq+"pt"+dq)
[161]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
[526]823            if (vNam(varn) .EQ. "pt_0_" .OR. vNam(varn) .EQ. "pt(0)")
[161]824               p=0
825               res@xyLineColors            = (/237/)
826               res@xyLineLabelFontHeightF = 0.05
827               res@xyLineLabelFontColor   = 237
[310]828               res@trYMaxF           = maxp
829               res@trYMinF           = minp
[161]830               plot_pt_0_ = gsn_csm_xy(wks_ps,t,data(varn,:),res) 
831            end if
[526]832            if (vNam(varn) .EQ. "pt_zp_" .OR. vNam(varn) .EQ. "pt(zp)")
[161]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
[154]841      end if
[161]842   end do
[310]843
844
845   do varn = dim-1, 0, 1
[161]846   
847      if( isStrSubset (vNam(varn), "time") )
848         check = False
849      else
850         check = True
851      end if   
[190]852 
853      if (var.NE. "all") then
[161]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
[162]865               plot_ps(n) = plot_E   
866
867               ; ***************************************************
[566]868               ; legend for overlaid plot
[162]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/)
[526]879               lbid = gsn_create_legend(wks_ps,2,(/"E","E*"/),lgres)       
[162]880
881               amres = True
882               amres@amParallelPosF   = 0.6                 
883               amres@amOrthogonalPosF = -0.2           
884               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)           
[161]885            end if
[566]886            if ((vNam(varn) .EQ. "Es" .OR. vNam(varn) .EQ. "E*") .AND. \
887                E .NE. 1) then
[161]888               Es=1
889               overlay(plot_E,plot_Es)
890               n=n+1
[162]891               plot_ps(n) = plot_E 
892
893               ; ***************************************************
[566]894               ; legend for overlaid plot
[162]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/)
[526]905               lbid = gsn_create_legend(wks_ps,2,(/"E","E*"/),lgres)       
[162]906
907               amres = True
908               amres@amParallelPosF   = 0.6                 
909               amres@amOrthogonalPosF = -0.2           
910               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)         
[161]911            end if
912
[566]913            if ((vNam(varn) .EQ. "us" .OR. vNam(varn) .EQ. "u*") .AND. \
914                 ws .NE. 1) then
[161]915               us=1
916               overlay(plot_us,plot_ws)
917               n=n+1
918               plot_ps(n) = plot_us
[162]919
920               ; ***************************************************
[566]921               ; legend for overlaid plot
[162]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/)
[526]932               lbid = gsn_create_legend(wks_ps,2,(/"u*","w*"/),lgres)       
[162]933
934               amres = True
935               amres@amParallelPosF   = 0.6                 
936               amres@amOrthogonalPosF = -0.2           
937               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
[161]938            end if
[566]939            if ((vNam(varn) .EQ. "ws" .OR. vNam(varn) .EQ. "w*") .AND. \
940                 us .NE. 1) then
[161]941               ws=1
942               overlay(plot_us,plot_ws)
943               n=n+1
944               plot_ps(n) = plot_us
[162]945
946               ; ***************************************************
[566]947               ; legend for overlaid  plot
[162]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/)
[526]958               lbid = gsn_create_legend(wks_ps,2,(/"u*","w*"/),lgres)       
[162]959
960               amres = True
961               amres@amParallelPosF   = 0.6                 
962               amres@amOrthogonalPosF = -0.2           
963               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
[161]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)
[310]971               
[161]972                  n=n+1
973                  plot_ps(n) = plot_umax
[162]974
975                  ; ***************************************************
[566]976                  ; legend for overlaid plot
[162]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/)
[566]987                  lbid = gsn_create_legend(\
988                                  wks_ps,3,(/"umax","vmax","wmax"/),lgres)
[162]989
990                  amres = True
991                  amres@amParallelPosF   = 0.6             
992                  amres@amOrthogonalPosF = -0.2           
993                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
[161]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
[162]1003
1004                  ; ***************************************************
[566]1005                  ; legend for overlaid plot
[162]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/)
[566]1016                  lbid = gsn_create_legend(\
1017                                      wks_ps,3,(/"umax","vmax","wmax"/),lgres)
[162]1018                  amres = True
1019                  amres@amParallelPosF   = 0.6             
1020                  amres@amOrthogonalPosF = -0.2           
1021                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
[161]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
[162]1031
1032                  ; ***************************************************
[566]1033                  ; legend for overlaid plot
[162]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/)
[566]1044                  lbid = gsn_create_legend(\
1045                                       wks_ps,3,(/"umax","vmax","wmax"/),lgres)
[162]1046                  amres = True
1047                  amres@amParallelPosF   = 0.6             
1048                  amres@amOrthogonalPosF = -0.2           
1049                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
[161]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
[162]1057               plot_ps(n) = plot_z_i_wpt
1058       
1059               ; ***************************************************
[566]1060               ; legend for overlaid plot
[162]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/)
[566]1071               lbid = gsn_create_legend(wks_ps,2,(/"z_i_wpt","z_i_pt"/),lgres)
[162]1072
1073               amres = True
1074               amres@amParallelPosF   = 0.6                 
1075               amres@amOrthogonalPosF = -0.2           
1076               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)         
[161]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
[162]1082               plot_ps(n) = plot_z_i_wpt
1083
1084               ; ***************************************************
[566]1085               ; legend for overlaid plot
[162]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/)
[566]1096               lbid = gsn_create_legend(wks_ps,2,(/"z_i_wpt","z_i_pt"/),lgres)
[162]1097
1098               amres = True
1099               amres@amParallelPosF   = 0.6                 
1100               amres@amOrthogonalPosF = -0.2           
1101               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)           
[161]1102            end if   
1103         
[526]1104            if ((vNam(varn) .EQ. "wpptp0" .OR. vNam(varn) .EQ. "w"+dq+"pt"+dq+"0") .AND. wp .NE. 1)
[161]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
[162]1111
1112                  ; ***************************************************
[566]1113                  ; legend for overlaid plot
[162]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/)
[566]1124                  lbid = gsn_create_legend(\
1125                                    wks_ps,3,(/"wpptp0","wpptp","wpt"/),lgres)
[162]1126                  amres = True
1127                  amres@amParallelPosF   = 0.6             
1128                  amres@amOrthogonalPosF = -0.2           
1129                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
[161]1130               end if           
1131            end if
[566]1132            if ((vNam(varn) .EQ. "wpptp" .OR. \
1133                 vNam(varn) .EQ. "w"+dq+"pt"+dq) .AND. w0 .NE. 1)
[161]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
[162]1139                 plot_ps(n) = plot_wpptp0
1140                 
1141                 ; ***************************************************
[566]1142                  ; legend for overlaid plot
[162]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/)
[566]1153                  lbid = gsn_create_legend(\
1154                                     wks_ps,3,(/"wpptp0","wpptp","wpt"/),lgres)
[162]1155
1156                  amres = True
1157                  amres@amParallelPosF   = 0.6             
1158                  amres@amOrthogonalPosF = -0.2           
1159                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres) 
[161]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
[162]1168                  plot_ps(n) = plot_wpptp0
1169
1170                  ; ***************************************************
[566]1171                  ; legend for overlaid plot
[162]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/)
[566]1182                  lbid = gsn_create_legend(\
1183                                  wks_ps,3,(/"wpptp0","wpptp","wpt"/),lgres)
[162]1184                  amres = True
1185                  amres@amParallelPosF   = 0.6             
1186                  amres@amOrthogonalPosF = -0.2           
1187                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
[161]1188               end if
1189            end if
1190
[566]1191            if ((vNam(varn) .EQ. "pt_0_" .OR. vNam(varn) .EQ. "pt(0)") .AND. \
1192                pz .NE. 1) then
[161]1193               p=1     
1194               overlay(plot_pt_0_,plot_pt_zp_)
1195               n=n+1
[162]1196               plot_ps(n) = plot_pt_0_ 
1197     
1198               ; ***************************************************
[566]1199               ; legend for overlaid plot
[162]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/)
[566]1210               lbid = gsn_create_legend(wks_ps,2,(/"pt_0_","pt_zp_"/),lgres)
[162]1211
1212               amres = True
1213               amres@amParallelPosF   = 0.6                 
1214               amres@amOrthogonalPosF = -0.2           
1215               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
[161]1216            end if
[566]1217            if ((vNam(varn) .EQ. "pt_zp_" .OR. vNam(varn) .EQ. "pt(zp)") .AND.\
1218                 p .NE. 1) then
[161]1219               pz=1       
1220               overlay(plot_pt_0_,plot_pt_zp_)
1221               n=n+1
[162]1222               plot_ps(n) = plot_pt_0_   
1223
1224               ; ***************************************************
[566]1225               ; legend for overlaid plot
[162]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/)
[566]1236               lbid = gsn_create_legend(wks_ps,2,(/"pt_0_","pt_zp_"/),lgres)
[162]1237
1238               amres = True
1239               amres@amParallelPosF   = 0.6                 
1240               amres@amOrthogonalPosF = -0.2           
1241               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
[161]1242            end if
[162]1243
[526]1244            if(vNam(varn) .NE. "pt_zp_" .AND. vNam(varn) .NE. "pt(zp)" .AND. \
[566]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.         \
[526]1256               vNam(varn) .NE. "E") then
[162]1257
[161]1258               n=n+1
1259               res@xyLineColors   = (/237/)
1260               res@xyLabelMode    = False
1261               res@gsnLeftString  = vNam(varn)
1262               res@gsnRightString = unit(varn)
[310]1263               res@trYMaxF        = max(data(varn,:))
1264               res@trYMinF        = min(data(varn,:))
[161]1265               if (min(data(varn,:)) .EQ. max(data(varn,:))) then
1266                  if (min(data(varn,:)) .EQ. 0)then
[310]1267                     res@trYMaxF = max(data(varn,:))-0.1
1268                     res@trYMinF = min(data(varn,:))+0.1
[161]1269                  end if
1270                  if (min(data(varn,:)) .LT. 0)then
[310]1271                     res@trYMaxF = max(data(varn,:))+(min(data(varn,:)))/2
1272                     res@trYMinF = min(data(varn,:))-(max(data(varn,:)))/2
[161]1273                  end if
1274                  if (min(data(varn,:)) .GT. 0)then
[310]1275                     res@trYMaxF = max(data(varn,:))-(min(data(varn,:)))/2
1276                     res@trYMinF = min(data(varn,:))+(max(data(varn,:)))/2
[161]1277                  end if
1278               end if
1279               plot_ps(n) = gsn_csm_xy(wks_ps,t,data(varn,:),res) 
1280            end if
[310]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)
[161]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)
[310]1309       
[218]1310            if (norm_t .NE. 1.)then
[534]1311               res@tiXAxisString        = "t ("+unit_t+")"
[218]1312            else
[534]1313               res@tiXAxisString        = "t (s)"
[218]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
[310]1321            res@trYMaxF        = max(data(varn,:))
1322            res@trYMinF        = min(data(varn,:))
[161]1323            if (min(data(varn,:)) .EQ. max(data(varn,:))) then
1324               if (min(data(varn,:)) .EQ. 0)then
[310]1325                     res@trYMaxF = max(data(varn,:))-0.1
1326                     res@trYMinF = min(data(varn,:))+0.1
[161]1327                  end if
1328                  if (min(data(varn,:)) .LT. 0)then
[310]1329                     res@trYMaxF = max(data(varn,:))+(min(data(varn,:)))/2
1330                     res@trYMinF = min(data(varn,:))-(max(data(varn,:)))/2
[161]1331                  end if
1332                  if (min(data(varn,:)) .GT. 0)then
[310]1333                     res@trYMaxF = max(data(varn,:))-(min(data(varn,:)))/2
1334                     res@trYMinF = min(data(varn,:))+(max(data(varn,:)))/2
[161]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
[154]1341   end do
[161]1342 
[154]1343   ; ***************************************************
1344   ; merge plots onto one page
1345   ; ***************************************************
[532]1346
1347   no_frames = 0
[162]1348 
[566]1349   if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
1350        n .gt. no_rows*no_columns) then
[161]1351      gsn_panel(wks_ps,plot_ps(1:n),(/n,1/),resP)
[250]1352      print(" ")
1353      print("Outputs to .eps or .epsi have only one frame")
1354      print(" ")
[161]1355   else
[218]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)
[532]1359            no_frames = no_frames + 1           
[161]1360         else
[566]1361            gsn_panel(wks_ps, plot_ps(np:np+no_rows*no_columns-1),\
1362                                                   (/no_rows,no_columns/),resP)
[532]1363            no_frames = no_frames + 1
[161]1364         end if
1365      end do
1366   end if
[154]1367
[532]1368   if (format_out .EQ. "png" ) then
1369     png_output = new((/no_frames/), string)
1370     j = 0
1371
[969]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)
[532]1382       system(cmd)
[969]1383     else
[532]1384
[969]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
[532]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
[969]1416
[154]1417end
Note: See TracBrowser for help on using the repository browser.