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

Last change on this file since 4472 was 2837, checked in by gronemeier, 7 years ago

palmplot: config file of ncl scripts can be chosen in shell command

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