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

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

bugfix: checking for NCL version number

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
[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   
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.