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

Last change on this file since 1752 was 1559, checked in by heinze, 10 years ago

Changes to allow for using NCL version 6.2.1 and higher. Backward compatibility is also ensured.

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