source: palm/trunk/SCRIPTS/NCL/spectra.ncl @ 1262

Last change on this file since 1262 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: 24.0 KB
RevLine 
[176]1load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
3load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
4
[190]5;***************************************************
[418]6; load .ncl.config or .ncl.config.default
[190]7;***************************************************
[194]8
9function which_script()
10local script
11begin
12   script="spectra"
13   return(script)
14end
[190]15   
[418]16if (isfilepresent("$PALM_BIN/../../.ncl.config")) then
17   loadscript("$PALM_BIN/../../.ncl.config")
[190]18else
[418]19  if (isfilepresent("$PALM_BIN/NCL/.ncl.config.default")) then
20     loadscript( "$PALM_BIN/NCL/.ncl.config.default")
[190]21  else
[418]22      palm_bin_path = getenv("PALM_BIN")
[190]23      print(" ")
[418]24      print("Neither the personal configuration file '.ncl.config' exists in")
25      print("~/palm/current_version")
[566]26      print("nor the default configuration file '.ncl.config.default' "+\
27            "exists in")
[418]28      print(palm_bin_path + "/NCL")
[190]29      print(" ")
30      exit
[176]31   end if
[418]32end if   
[218]33
[190]34begin
[176]35
[190]36   ;***************************************************
[532]37   ; Retrieving the NCL version used
38   ;***************************************************
39   
40   ncl_version_ch = systemfunc("ncl -V")
41   ncl_version    = stringtofloat(ncl_version_ch)
42
43   ;***************************************************
[526]44   ; Retrieving the double quote character
45   ;***************************************************
46   
47   dq=str_get_dq()
48
49   ;***************************************************
[190]50   ; set up default parameter values and strings
51   ;***************************************************
52 
53   if (file_1 .EQ. "File in") then
54      print(" ")
[418]55      print("Declare input file 'file_1=' in '.ncl.config' or prompt")
[190]56      print(" ")
57      exit
[176]58   else
59      file_in = file_1   
60   end if
61
[566]62   if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND.   \
63       format_out .NE. "eps" .AND. format_out .NE. "ps" .AND.    \
64       format_out .NE. "epsi" .AND. format_out .NE. "ncgm" .AND. \
65       format_out .NE. "png")then
[190]66      print(" ")
67      print("'format_out = "+format_out+"' is invalid and set to'x11'")
68      print(" ")
69      format_out="x11"
70   end if
[532]71
72   if (ncl_version .LE. 5.1 .AND. format_out .EQ. "png") then
73      print(" ")
74      print("Output of png files not available")
75      print("png output is avaiable with NCL version 5.2.0 and higher ")
76      print("NCL version used: " + ncl_version_ch)
77      print(" ")
78      exit
79   end if
[190]80   
[218]81   if (log_x .NE. 0 .AND. log_x .NE. 1)then
[190]82      print(" ")
[218]83      print("'log_x'= "+log_x+" is invalid and set to 1")
[190]84      print(" ")
[218]85      log_x = 1
[190]86   end if 
87   
[218]88   if (log_y .NE. 0 .AND. log_y .NE. 1)then
[190]89      print(" ")
[218]90      print("'log_y'= "+log_y+" is invalid and set to 1")
[190]91      print(" ")
[218]92      log_y = 1
[190]93   end if   
94 
[218]95   if (norm_y .EQ. 0.) then
[190]96      print(" ")
[218]97      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
[190]98      print(" ")
[218]99      norm_y = 1.0
[176]100   end if
[218]101   if (norm_x .EQ. 0.) then
[190]102      print(" ")
[218]103      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
[190]104      print(" ")
[218]105      norm_x= 1.0
[176]106   end if
[190]107   
108   if (sort .NE. "height" .AND. sort .NE. "time") then 
109      print(" ")
110      print("'sort'= "+sort+" is invalid and set to 'height'")
111      print(" ")
112      sort = "height" 
[176]113   end if
[190]114   
115   if (black .NE. 0 .AND. black .NE. 1)then
116      print(" ")
117      print("'black'= "+black+" is invalid and set to 0")
118      print(" ")
[176]119      black = 0
120   end if
[190]121 
122   if (dash .NE. 0 .AND. dash .NE. 1)then
123      print(" ")
124      print("'dash'= "+dash+" is invalid and set to 0")
125      print(" ")
[176]126      dash = 0
127   end if
128
[190]129   ;***************************************************
130   ; open input file
131   ;***************************************************
[218]132   
133   file_in_1 = False
134   if (isStrSubset(file_in, ".nc"))then
135      start_f = -2
136      end_f = -2
137      file_in_1 = True     
138   end if
[176]139
[218]140   if (start_f .EQ. -1)then
141      print(" ")
[566]142      print("'start_f' must be one of the cyclic numbers (at least 0) of "+\
143            "your input file(s)")
[218]144      print(" ") 
145      exit
146   end if
147   if (end_f .EQ. -1)then
148      print(" ")
[566]149      print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
150            "your input file(s)")
[218]151      print(" ") 
152      exit
153   end if
154
155   files=new(end_f-start_f+1,string)
[176]156   
[218]157   if (file_in_1)then
158      if (isfilepresent(file_in))then
159         files(0)=file_in
160      else
161         print(" ")
162         print("1st input file: '"+file_in+"' does not exist")
163         print(" ")
164         exit
165      end if
166   else   
167      if (start_f .EQ. 0)then
168         if (isfilepresent(file_in+".nc"))then
169            files(0)=file_in+".nc"
170            do i=1,end_f
171               if (isfilepresent(file_in+"."+i+".nc"))then   
172                  files(i)=file_in+"."+i+".nc"
173               else
174                  print(" ")
175                  print("Input file: '"+file_in+"."+i+".nc' does not exist")
176                  print(" ")
177                  exit 
178               end if       
179            end do         
180         else
181            print(" ")
182            print("Input file: '"+file_in+".nc' does not exist")
183            print(" ")
184            exit
185         end if
186      else
187         do i=start_f,end_f
188            if (isfilepresent(file_in+"."+i+".nc"))then   
189               files(i-start_f)=file_in+"."+i+".nc"
190            else
191               print(" ")
192               print("Input file: '"+file_in+"."+i+".nc' does not exist")
193               print(" ")
194               exit 
195            end if
196         end do
197      end if
198   end if
199
200   f=addfiles(files,"r")
201   f_att=addfile(files(0),"r")
202   ListSetType(f,"cat")
203   
204   vNam = getfilevarnames(f_att)
[1248]205   if( ncl_version .GE. 6.1 ) then
206      vNam = vNam(::-1)
207   end if
[585]208   vType = getfilevartypes(f_att,vNam)
209
210   if ((all(vType .eq. "double"))) then ;distinction if data is double or float
211      check_vType = True
212   else
213      check_vType = False
214   end if
215
[176]216   print(" ")
[190]217   print("Variables in input file:")
218   print("- "+ vNam)
[176]219   print(" ")
220   dim = dimsizes(vNam)
[218]221   vDim = getfiledimsizes(f_att)
[176]222 
[218]223   t_all = f[:]->time
[176]224   nt    = dimsizes(t_all)
[827]225
226   if (nt .EQ. 1)then
227      delta_t=t_all(nt-1)/nt
228   else
229      delta_t_array = new(nt-1,double)
230
231      do i=0,nt-2
232         delta_t_array(i) = t_all(i+1)-t_all(i)
233      end do
234
235      delta_t = min(delta_t_array)
236      delete(delta_t_array)
237   end if
[176]238   
[218]239   k_x=f_att->k_x
[190]240   dimx=dimsizes(k_x)
[218]241   k_y=f_att->k_y
[190]242   dimy=dimsizes(k_y)
243   
244   
245   dim_level=dimsizes(height_level)
246
[176]247   do i=0,dim-1
248      if (vNam(i) .EQ. "zu_sp")then
[218]249         zu=f_att->zu_sp         
[190]250         if (height_level(0) .EQ. -1)then
251            dimz=dimsizes(zu)
252         else
253            if (dim_level .GT. dimsizes(zu))then
254               print(" ")
[566]255               print("'height_level' has more elements than available "+\
256                     "height levels in input file (= "+dimsizes(zu)+")")
[190]257               print(" ")
258               exit
259            else
260               zuh=new(dim_level,double)
261               do le=0,dim_level-1
[566]262                  if (height_level(le) .GE. dimsizes(zu)) then
263                     no_levels=dimsizes(zu)-1
264                     print(" ")
265                     print("Element "+le+" of 'height_level' is greater "  +\
266                          "than the maximum available index in input file "+\
267                          "which is "+no_levels+". Note that the first "   +\
268                          "element has the index 0.")
269                     print(" ")
270                     exit
271                  end if
[190]272                  zuh(le)=zu(height_level(le))
273               end do
274               dimz=dim_level
275            end if   
276         end if 
[176]277      else
278         if (vNam(i) .EQ. "zw_sp")then
[218]279            zw=f_att->zw_sp
[190]280            if (height_level(0) .EQ. -1)then             
281               dimz=dimsizes(zw)
282            else
283               if (dim_level .GT. dimsizes(zw))then
284                  print(" ")
[566]285                  print("'height_level' has more elements than available "+\
286                        "height levels in input file (= "+dimsizes(zw)+")")
[190]287                  print(" ")
288                  exit
289               else
290                  zwh=new(dim_level,double)
291                  do le=0,dim_level-1
[566]292                     if (height_level(le) .GE. dimsizes(zw)) then
293                        no_levels=dimsizes(zw)-1
294                        print(" ")
295                        print("Element "+le+" of 'height_level' is greater "+\
296                              "than the maximum available index in input "  +\
297                              "file which is "+no_levels+". Note that the " +\
298                              "first element has the index 0.")
299                        print(" ")
300                        exit
301                     end if
[190]302                     zwh(le)=zw(height_level(le))
303                  end do
304                  dimz=dim_level
305               end if   
306            end if
[176]307         end if
308      end if
309   end do
[190]310
311   ;****************************************************       
[176]312   ; start of time step and different types of mistakes that could be done
[190]313   ;****************************************************
[176]314   
[190]315   if (start_time_step .EQ. -1.d) then         
[176]316      start_time_step=t_all(0)/3600
[190]317   else
[176]318      if (start_time_step .GT. t_all(nt-1)/3600)then
319         print(" ")
[566]320         print("'start_time_step' = "+ start_time_step +"h is greater than "+\
321               "last time step = "+ t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
[176]322         print(" ")
[218]323         print("Select another 'start_time_step'")
[176]324         print(" ")
325         exit
326      end if
327      if (start_time_step .LT. t_all(0)/3600)then
328         print(" ")
[566]329         print("'start_time_step' = "+ start_time_step +"h is lower than "+\
330               "first time step = "+ t_all(0)+"s = "+t_all(0)/3600+"h")       
[176]331         exit
332      end if
333   end if
334
335   do i=0,nt-1     
[566]336      if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
337          start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[176]338         st=i
339         break
[769]340      else
341         st=0
[176]342      end if
343   end do
344   
[190]345   ;****************************************************
[176]346   ; end of time step and different types of mistakes that could be done
[190]347   ;****************************************************
[176]348
[190]349   if (end_time_step .EQ. -1.d) then           
[176]350      end_time_step = t_all(nt-1)/3600
351   else
352      if (end_time_step .GT. t_all(nt-1)/3600)then
353         print(" ")
[566]354         print("'end_time_step' = "+ end_time_step +"h is greater than "+\
355               "last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
[176]356         print(" ")
[218]357         print("Select another 'end_time_step'") 
[176]358         print(" ")
359         exit
360      end if
361      if (end_time_step .LT. start_time_step/3600)then
362         print(" ")
[566]363         print("'end_time_step' = "+ end_time_step +"h is lower than "+\
364               "'start_time_step' = "+start_time_step+"h")
[176]365         print(" ")
[218]366         print("Select another 'start_time_step' or 'end_time_step'")
[176]367         print(" ")
368         exit
369      end if
370   end if
371
372   do i=0,nt-1     
[566]373      if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
374          end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[176]375         et=i
376         break
[769]377       else
378         et=0
[176]379      end if
380   end do
381
382   delete(start_time_step)
383   start_time_step=round(st,3)
384   delete(end_time_step)
385   end_time_step=round(et,3)
386
387   print(" ")
[566]388   print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+\
389         t_all(start_time_step)+" s => index = "+start_time_step)
390   print("                     till "+t_all(end_time_step)/3600+" h = "+\
391         t_all(end_time_step)+" s => index = "+end_time_step)
[176]392   print(" ")
393
394   dimt = end_time_step-start_time_step+1
395 
[190]396   ;***************************************************
[176]397   ; set up recourses
[190]398   ;***************************************************
[176]399     
400   res = True
401   res@gsnDraw                 = False
402   res@gsnFrame                = False
403   res@txFont                  = "helvetica"
404   res@tiMainFont              = "helvetica"
405   res@tiXAxisFont             = "helvetica"
406   res@tiYAxisFont             = "helvetica"
407   res@tmXBLabelFont           = "helvetica"
408   res@tmYLLabelFont           = "helvetica"
409   res@lgLabelFont             = "helvetica"
410   res@tmLabelAutoStride       = True
411   
[218]412   res@lgLabelFontHeightF     = font_size_legend
413   res@lgTitleFontHeightF     = font_size
414   res@txFontHeightF      = font_size
415   res@tiXAxisFontHeightF = font_size
416   res@tiYAxisFontHeightF = font_size
417   res@tmXBLabelFontHeightF = font_size
418   res@tmYLLabelFontHeightF = font_size
419   
420   res@tmXBMinorPerMajor = 4
421   res@tmYLMinorPerMajor = 4
422   
423   if (log_x .EQ. 1) then
[176]424      res@trXLog = True
[190]425   else
426      res@trXLog = False
427   end if
[218]428   if (log_y .EQ. 1)then
[176]429      res@trYLog = True
[190]430   else 
[178]431      res@trYLog = False
[176]432   end if
433
[566]434   legend_label=new(dimt,string)
[176]435   legend_label_zu=new(dimz,double)
436   legend_label_zw=new(dimz,double)
437   legend_label_z=new(dimz,double)
438   do p=start_time_step,end_time_step
[566]439      legend_label(p-start_time_step)=sprintf("%6.2f", t_all(p)/3600)
[176]440   end do 
441   if (sort .EQ. "time")
442      plot  = new(dim*dimz,graphic)
443      np=dimt
[534]444      res@lgTitleString = "Time (h)"
[176]445   else
446      plot  = new(dim*dimt,graphic)
447      np=dimz
[218]448     
[176]449      do p=0,dimz-1
[190]450         if (height_level(0) .EQ. -1)then
451            legend_label_zu(p)=round(zu(p),3)
452            legend_label_zw(p)=round(zw(p),3)
453         else
454            legend_label_zu(p)=round(zuh(p),3)
455            legend_label_zw(p)=round(zwh(p),3)
456         end if
[176]457      end do
458   end if
[194]459
460   if (black .eq. 0 ) then
461      if (np .EQ. 1)then
[218]462         color = 237
[194]463      else   
464         step=round(235/(np-1),3)
[218]465         color = ispan(2,237,step)
[194]466      end if
[324]467   else
468      color = 2
[176]469   end if
470   if ( dash .eq. 0 ) then
[190]471      res@xyMonoDashPattern = True
[176]472   end if
473
[532]474   if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
475      format_out@wkPaperSize = "A4"
476   end if
477   if ( format_out .EQ. "png" ) then
478      format_out@wkWidth  = 1000
479      format_out@wkHeight = 1000
480   end if
481
[176]482   wks=gsn_open_wks(format_out,file_out)
483   gsn_define_colormap(wks,"rainbow+white")
484
485   n=0
486   do varn =dim-1,0,1
[178]487     
[176]488      check = True
489
[566]490      if ( isStrSubset( vNam(varn), "time") .OR.  \
491           isStrSubset( vNam(varn), "zu_sp") .OR. \
492           isStrSubset( vNam(varn), "zw_sp") .OR. \
493           isStrSubset( vNam(varn), "k_x") .OR.   \
494           isStrSubset( vNam(varn), "k_y")) then
[176]495            check = False
496      end if
497
[190]498      if (var .NE. "all") then
[176]499         check = isStrSubset( var,","+vNam(varn)+"," )
500      end if
501
[190]502      if(check) then
[218]503
504         temp = f[:]->$vNam(varn)$
[566]505         data = temp(start_time_step:end_time_step,0:dimz-1,:)
[218]506
507         temp_att = f_att->$vNam(varn)$
508         a=getvardims(temp_att)
[190]509         b=dimsizes(a)
[534]510         delete(temp_att)
[176]511
[190]512         if (height_level(0) .NE. -1)then
513            do te=0,dimz-1
[566]514               data(:,te,:) = temp(start_time_step:end_time_step,\
515                                                    height_level(te),:)       
[190]516            end do
517         end if
518
[566]519         data=data/(norm_y*norm_x)
[190]520           
[176]521         do i=0,b-1           
522            if (isStrSubset( a(i),"zu_sp" ))then
523               legend_label_z=legend_label_zu
[190]524               if (height_level(0) .NE. -1)then
525                  z=zuh
526               else
527                  z=zu
528               end if
[176]529            else
530               if (isStrSubset( a(i),"zw_sp" ))then
[190]531                  legend_label_z=legend_label_zw
532                  if (height_level(0) .NE. -1)then
533                     z=zwh
534                  else
535                     z=zw
536                  end if   
[176]537               end if
538            end if
539         end do 
[218]540         
[585]541         if (check_vType) then
542            min_y=new(dimz,double)
543            max_y=new(dimz,double)
544         else
545            min_y=new(dimz,float)
546            max_y=new(dimz,float)
547         end if
[218]548         min_x=new(dimz,double)
549         max_x=new(dimz,double)
[585]550         
[218]551         plot_h  = new(dimz,graphic)
552         
[176]553         if (isStrSubset(vNam(varn),"x"))then
[218]554            x_axis = new((/dimz,dimx/),double)
555            do q=0,dimz-1
556               x_axis(q,:) = f_att->k_x
557               x_axis = x_axis/norm_x
558            end do
559            if (norm_x .NE. 1.)then
[534]560               res@tiXAxisString = "k~B~x~N~ ("+unit_x+")"
[190]561            else
[218]562               if (norm_height .EQ. 1)then
[534]563                  res@tiXAxisString = "k~B~x~N~ x z (1"
[218]564               else
[534]565                  res@tiXAxisString = "k~B~x~N~ (1/m)"
[218]566               end if
[190]567            end if
[218]568            dim_r=dimx
[176]569         else
[218]570            x_axis=new((/dimz,dimy/),double)
571            do q=0,dimz-1
572               x_axis(q,:) = f_att->k_y
573               x_axis = x_axis/norm_x
574            end do
575            if (norm_x .NE. 1.)then
[534]576               res@tiXAxisString = "k~B~x~N~ ("+unit_x+")"
[190]577            else
[218]578               if (norm_height .EQ. 1)then
[534]579                  res@tiXAxisString = "k~B~x~N~ x z (1)"
[218]580               else
[534]581                  res@tiXAxisString = "k~B~x~N~ (1/m)"
[218]582               end if
[190]583            end if
[218]584            dim_r=dimy
[176]585         end if
586       
587         if (sort .EQ. "time")
[218]588            res@xyLineColors = color
589            res@pmLegendDisplayMode     = "Always"
590            res@pmLegendSide            = "Top"
591            res@pmLegendParallelPosF    = 1.2
592            res@pmLegendOrthogonalPosF  = -1.0
593            res@pmLegendWidthF          = 0.12
[566]594            res@pmLegendHeightF         = 0.04*\
595                                          (end_time_step-start_time_step+1)
[190]596            do p=dimz-1,0,1
[218]597               if (log_y .EQ. 1)then 
[190]598                  do q=0,dimt-1
[218]599                     do r=0,dim_r-1
600                        if (data(q,p,r) .EQ. 0)then
[190]601                           st=p+start_time_step
602                           print(" ")
[566]603                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is "+\
604                                 "equal 0. Logarithmic scale for y-axis "+\
605                                 "and height "+z(p)+" cannot be used")
[190]606                           print(" ")
607                           res@trYLog = False
608                        end if
609                     end do
[176]610                  end do
[190]611               end if
[176]612               res@gsnLeftString      = vNam(varn)
[218]613               res@gsnRightString     = "Height = "+z(p)+"m"               
[534]614               res@tiYAxisString      = "("+unit_y+")"           
[218]615               res@xyExplicitLegendLabels  = legend_label             
616               if (norm_height .EQ. 1)then
617                  data(:,p,:)=data(:,p,:)*doubletofloat(z(p))
618                  x_axis(p,:) = x_axis(p,:)*z(p)
619               end if
620               res@trXMinF = min(x_axis(p,:))
621               res@trXMaxF = max(x_axis(p,:))
622               plot(n)  = gsn_csm_xy(wks,x_axis(p,:),data(:,p,:),res)
[176]623               n=n+1
624            end do
625         else
[190]626            if (sort .EQ. "height")
627               do p=0,dimt-1           
[176]628                  do q=0,dimz-1
[218]629                     do r=0,dim_r-1
630                        if (data(p,q,r) .EQ. 0)then
[176]631                           st=p+start_time_step
632                           print(" ")
[566]633                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' "+\
634                                 "is equal 0. Logarithmic scale for y-axis "+\
635                                 "and time "+legend_label(p)+" h cannot be used")
[176]636                           print(" ")
637                           res@trYLog = False
638                        end if
639                     end do
[218]640                     if (norm_height .EQ. 1 .AND. p .EQ. 0)then
[566]641                        data(p,q,:) = data(p,q,:)*\
642                                           doubletofloat(legend_label_z(q))
643                        x_axis(q,:) = x_axis(q,:)*\
644                                           doubletofloat(legend_label_z(q))
[218]645                     end if
646                     max_y(q)=max(data(p,q,:))
647                     min_y(q)=min(data(p,q,:))
648                     min_x(q)=min(x_axis(q,:))
[566]649                     max_x(q)=max(x_axis(q,:))
[218]650                  end do
651                  do q=0,dimz-1
652                     res@xyLineColor = color(q)
653                     if (dash .EQ. 1)then
654                        res@xyDashPattern = q
655                     end if
656                     if (q .EQ. 0)then
[566]657                        res@tiYAxisString      = "("+unit_y+")"
[218]658                        res@gsnLeftString      = vNam(varn)
659                        res@gsnRightString     = "Time = "+legend_label(p)+"h"
660                        res@trXMinF = min(min_x)
661                        res@trXMaxF = max(max_x)
662                        res@trYMinF = min(min_y)
663                        res@trYMaxF = max(max_y)
664                       
[566]665                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),\
666                                                         data(p,q,:),res)
[218]667                       
668                        lgres = True
669                        if (dash .EQ. 0)then
670                           lgres@lgMonoDashIndex = True
671                        else
672                           lgres@lgDashIndexes   = ispan(0,dimz-1,1)
673                        end if
674                        if (black .EQ. 1)then
675                           lgres@lgMonoLineColors = True
676                        else
677                           lgres@lgLineColors = color
678                        end if
[566]679                        lgres@lgTitleString      = "Height (m)" 
[218]680                        lgres@lgLabelFont        = "helvetica"
[566]681                        lgres@lgLabelFontHeightF = font_size_legend*6
682                        lgres@lgTitleFontHeightF = font_size     
[218]683                        lgres@vpWidthF           = 0.12           
[566]684                        lgres@vpHeightF          = font_size_legend*(dimz+3)
[218]685 
[566]686                        lbid = gsn_create_legend(wks,dimz,legend_label_z,lgres)
[218]687
688                        amres = True
689                        amres@amParallelPosF   = 0.75               
690                        amres@amOrthogonalPosF = 0.15           
691                        annoid1 = gsn_add_annotation(plot_h(q),lbid,amres)
692                     else
[566]693                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),\
694                                                        data(p,q,:),res)
[218]695                        overlay(plot_h(0),plot_h(q))
696                     end if
697                  end do             
698                  plot(n)=plot_h(0)
699                  n=n+1
[176]700               end do
701            end if
[190]702         end if
[218]703         delete(data)
704         delete(temp)
[178]705         delete(x_axis)
[218]706         delete(min_x)
707         delete(max_x)
708         delete(min_y)
709         delete(max_y)
710         delete(plot_h)
[176]711      end if
712   end do
713
714   if (n .EQ. 0) then
715      print(" ")
[190]716      print("The variables 'var="+var+"' do not exist on your input file;")
717      print("be sure to have one comma berfore and after each variable")
[176]718      print(" ")
719      exit
720   end if
721 
722   ; ***************************************************
723   ; merge plots onto one page
724   ; ***************************************************
725
[983]726   resP                            = True
727   resP@gsnMaximize                = True
728   resP@gsnPanelXWhiteSpacePercent = 4.0
729   resP@gsnPanelYWhiteSpacePercent = 4.0
730   resP@txFont                     = "helvetica"
731   resP@txString                   = f_att@title
732   resP@txFuncCode                 = "~"
733   resP@txFontHeightF              = 0.0105
[176]734
[534]735   no_frames = 0
736
[566]737   if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
738                                          n .gt. no_rows*no_columns) then
[183]739      gsn_panel(wks,plot,(/n,1/),resP)
[250]740      print(" ")
741      print("Outputs to .eps or .epsi have only one frame")
742      print(" ")
[176]743   else   
[218]744      do i = 0,n-1, no_rows*no_columns
745         if( (i+no_rows*no_columns) .gt. (n-1)) then
746            gsn_panel(wks,plot(i:n-1),(/no_rows,no_columns/),resP)
[534]747            no_frames = no_frames + 1
[176]748         else
[566]749            gsn_panel(wks,plot(i:i+no_rows*no_columns-1),\
750                                           (/no_rows,no_columns/),resP)
[534]751            no_frames = no_frames + 1   
[176]752         end if
753      end do
754   end if
755
[532]756   if (format_out .EQ. "png" ) then
757     png_output = new((/no_frames/), string)
758     j = 0
759
[957]760     if (no_frames .eq. 1) then
761        if (ncl_version .GE. 6 ) then
762           png_output(0) = file_out+".png"
763        else
764           png_output(0) = file_out+".00000"+1+".png"
765        end if
766        ;using imagemagick's convert for reducing the white
767        ;space around the plot
768        cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
769              png_output(0) + " " + png_output(0)
[532]770       system(cmd)
[957]771     else
[532]772
[957]773       do i=0, no_frames-1
774         j = i + 1
775         if (j .LE. 9) then
776           png_output(i) = file_out+".00000"+j+".png"
777         end if
778         if (j .GT. 9 .AND. j .LE. 99) then
779           png_output(i) = file_out+".0000"+j+".png"
780         end if
781         if (j .GT. 99 .AND. j .LE. 999) then
782           png_output(i) = file_out+".000"+j+".png"
783         end if
784         if (j .GT. 999) then
785           png_output(i) = file_out+".00"+j+".png"
786         end if
787
788         ;using imagemagick's convert for reducing the white
789         ;space around the plot
790         cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
791                png_output(i) + " " + png_output(i)
792         system(cmd)
793       end do
794     end if
795
[532]796     print(" ")
797     print("Output to: "+ png_output)
798     print(" ")
799   else
800     print(" ")
801     print("Output to: " + file_out +"."+ format_out)
802     print(" ")
803   end if
[176]804   
[190]805end
Note: See TracBrowser for help on using the repository browser.