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

Last change on this file since 1060 was 983, checked in by hoffmann, 12 years ago

Bugfix in netcdf.f90 and improvements in palmplot.

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