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

Last change on this file since 1374 was 1248, checked in by heinze, 11 years ago

NCL function getfilevarnames changed with NCL version 6.1.1 and higher -> adaptation required

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