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

Last change on this file since 746 was 585, checked in by heinze, 14 years ago

Bugfix: enable plot of data if it is of kind double instead of kind float

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