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

Last change on this file since 4358 was 2837, checked in by gronemeier, 3 years ago

palmplot: config file of ncl scripts can be chosen in shell command

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