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

Last change on this file since 2644 was 2563, checked in by Giersch, 7 years ago

Restart runs with the usage of the wind turbine model are possible now. Further small at reading/writing restart data

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 .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   
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*10)
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.