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

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