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

Last change on this file since 3126 was 2837, checked in by gronemeier, 7 years ago

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

File size: 24.4 KB
RevLine 
[176]1load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
3load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
4
[190]5;***************************************************
[1559]6; Checking the kind of the script
[190]7;***************************************************
[194]8
9function which_script()
10local script
11begin
[1559]12   script = "spectra"
[194]13   return(script)
14end
[1559]15
16;***************************************************
17; Retrieving the NCL version used
18;***************************************************
[190]19   
[1559]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
[2030]30   if( version .GE. "6.2.1" ) then
[1559]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
[2837]42existing_file = file_exists(ncl_version_ch,file_config)
[1559]43if (existing_file) then
[2837]44   loadscript(file_config)
[190]45else
[2837]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
[218]56
[2837]57
[190]58begin
[176]59
[190]60   ;***************************************************
[526]61   ; Retrieving the double quote character
62   ;***************************************************
63   
64   dq=str_get_dq()
65
66   ;***************************************************
[190]67   ; set up default parameter values and strings
68   ;***************************************************
69 
70   if (file_1 .EQ. "File in") then
71      print(" ")
[418]72      print("Declare input file 'file_1=' in '.ncl.config' or prompt")
[190]73      print(" ")
74      exit
[176]75   else
76      file_in = file_1   
77   end if
78
[566]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
[190]83      print(" ")
84      print("'format_out = "+format_out+"' is invalid and set to'x11'")
85      print(" ")
86      format_out="x11"
87   end if
[532]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
[190]97   
[218]98   if (log_x .NE. 0 .AND. log_x .NE. 1)then
[190]99      print(" ")
[218]100      print("'log_x'= "+log_x+" is invalid and set to 1")
[190]101      print(" ")
[218]102      log_x = 1
[190]103   end if 
104   
[218]105   if (log_y .NE. 0 .AND. log_y .NE. 1)then
[190]106      print(" ")
[218]107      print("'log_y'= "+log_y+" is invalid and set to 1")
[190]108      print(" ")
[218]109      log_y = 1
[190]110   end if   
111 
[218]112   if (norm_y .EQ. 0.) then
[190]113      print(" ")
[218]114      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
[190]115      print(" ")
[218]116      norm_y = 1.0
[176]117   end if
[218]118   if (norm_x .EQ. 0.) then
[190]119      print(" ")
[218]120      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
[190]121      print(" ")
[218]122      norm_x= 1.0
[176]123   end if
[190]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" 
[176]130   end if
[190]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(" ")
[176]136      black = 0
137   end if
[190]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(" ")
[176]143      dash = 0
144   end if
145
[190]146   ;***************************************************
147   ; open input file
148   ;***************************************************
[218]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
[176]156
[218]157   if (start_f .EQ. -1)then
158      print(" ")
[566]159      print("'start_f' must be one of the cyclic numbers (at least 0) of "+\
160            "your input file(s)")
[218]161      print(" ") 
162      exit
163   end if
164   if (end_f .EQ. -1)then
165      print(" ")
[566]166      print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
167            "your input file(s)")
[218]168      print(" ") 
169      exit
170   end if
171
172   files=new(end_f-start_f+1,string)
[176]173   
[218]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)
[1248]222   if( ncl_version .GE. 6.1 ) then
223      vNam = vNam(::-1)
224   end if
[585]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
[176]233   print(" ")
[190]234   print("Variables in input file:")
235   print("- "+ vNam)
[176]236   print(" ")
237   dim = dimsizes(vNam)
[218]238   vDim = getfiledimsizes(f_att)
[176]239 
[218]240   t_all = f[:]->time
[176]241   nt    = dimsizes(t_all)
[827]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
[176]255   
[218]256   k_x=f_att->k_x
[190]257   dimx=dimsizes(k_x)
[218]258   k_y=f_att->k_y
[190]259   dimy=dimsizes(k_y)
260   
261   
262   dim_level=dimsizes(height_level)
263
[176]264   do i=0,dim-1
265      if (vNam(i) .EQ. "zu_sp")then
[218]266         zu=f_att->zu_sp         
[190]267         if (height_level(0) .EQ. -1)then
268            dimz=dimsizes(zu)
269         else
270            if (dim_level .GT. dimsizes(zu))then
271               print(" ")
[566]272               print("'height_level' has more elements than available "+\
273                     "height levels in input file (= "+dimsizes(zu)+")")
[190]274               print(" ")
275               exit
276            else
277               zuh=new(dim_level,double)
278               do le=0,dim_level-1
[566]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
[190]289                  zuh(le)=zu(height_level(le))
290               end do
291               dimz=dim_level
292            end if   
293         end if 
[176]294      else
295         if (vNam(i) .EQ. "zw_sp")then
[218]296            zw=f_att->zw_sp
[190]297            if (height_level(0) .EQ. -1)then             
298               dimz=dimsizes(zw)
299            else
300               if (dim_level .GT. dimsizes(zw))then
301                  print(" ")
[566]302                  print("'height_level' has more elements than available "+\
303                        "height levels in input file (= "+dimsizes(zw)+")")
[190]304                  print(" ")
305                  exit
306               else
307                  zwh=new(dim_level,double)
308                  do le=0,dim_level-1
[566]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
[190]319                     zwh(le)=zw(height_level(le))
320                  end do
321                  dimz=dim_level
322               end if   
323            end if
[176]324         end if
325      end if
326   end do
[190]327
328   ;****************************************************       
[176]329   ; start of time step and different types of mistakes that could be done
[190]330   ;****************************************************
[176]331   
[190]332   if (start_time_step .EQ. -1.d) then         
[176]333      start_time_step=t_all(0)/3600
[190]334   else
[176]335      if (start_time_step .GT. t_all(nt-1)/3600)then
336         print(" ")
[566]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")
[176]339         print(" ")
[218]340         print("Select another 'start_time_step'")
[176]341         print(" ")
342         exit
343      end if
344      if (start_time_step .LT. t_all(0)/3600)then
345         print(" ")
[566]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")       
[176]348         exit
349      end if
350   end if
351
352   do i=0,nt-1     
[566]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
[176]355         st=i
356         break
[769]357      else
358         st=0
[176]359      end if
360   end do
361   
[190]362   ;****************************************************
[176]363   ; end of time step and different types of mistakes that could be done
[190]364   ;****************************************************
[176]365
[190]366   if (end_time_step .EQ. -1.d) then           
[176]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(" ")
[566]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")
[176]373         print(" ")
[218]374         print("Select another 'end_time_step'") 
[176]375         print(" ")
376         exit
377      end if
378      if (end_time_step .LT. start_time_step/3600)then
379         print(" ")
[566]380         print("'end_time_step' = "+ end_time_step +"h is lower than "+\
381               "'start_time_step' = "+start_time_step+"h")
[176]382         print(" ")
[218]383         print("Select another 'start_time_step' or 'end_time_step'")
[176]384         print(" ")
385         exit
386      end if
387   end if
388
389   do i=0,nt-1     
[566]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
[176]392         et=i
393         break
[769]394       else
395         et=0
[176]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(" ")
[566]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)
[176]409   print(" ")
410
411   dimt = end_time_step-start_time_step+1
412 
[190]413   ;***************************************************
[176]414   ; set up recourses
[190]415   ;***************************************************
[176]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   
[218]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
[176]441      res@trXLog = True
[190]442   else
443      res@trXLog = False
444   end if
[218]445   if (log_y .EQ. 1)then
[176]446      res@trYLog = True
[190]447   else 
[178]448      res@trYLog = False
[176]449   end if
450
[566]451   legend_label=new(dimt,string)
[176]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
[566]456      legend_label(p-start_time_step)=sprintf("%6.2f", t_all(p)/3600)
[176]457   end do 
458   if (sort .EQ. "time")
459      plot  = new(dim*dimz,graphic)
460      np=dimt
[534]461      res@lgTitleString = "Time (h)"
[176]462   else
463      plot  = new(dim*dimt,graphic)
464      np=dimz
[218]465     
[176]466      do p=0,dimz-1
[190]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
[176]474      end do
475   end if
[194]476
477   if (black .eq. 0 ) then
478      if (np .EQ. 1)then
[218]479         color = 237
[194]480      else   
481         step=round(235/(np-1),3)
[218]482         color = ispan(2,237,step)
[194]483      end if
[324]484   else
485      color = 2
[176]486   end if
487   if ( dash .eq. 0 ) then
[190]488      res@xyMonoDashPattern = True
[176]489   end if
490
[532]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
[176]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
[178]504     
[176]505      check = True
506
[566]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
[176]512            check = False
513      end if
514
[190]515      if (var .NE. "all") then
[176]516         check = isStrSubset( var,","+vNam(varn)+"," )
517      end if
518
[190]519      if(check) then
[218]520
521         temp = f[:]->$vNam(varn)$
[566]522         data = temp(start_time_step:end_time_step,0:dimz-1,:)
[218]523
524         temp_att = f_att->$vNam(varn)$
525         a=getvardims(temp_att)
[190]526         b=dimsizes(a)
[534]527         delete(temp_att)
[176]528
[190]529         if (height_level(0) .NE. -1)then
530            do te=0,dimz-1
[566]531               data(:,te,:) = temp(start_time_step:end_time_step,\
532                                                    height_level(te),:)       
[190]533            end do
534         end if
535
[566]536         data=data/(norm_y*norm_x)
[190]537           
[176]538         do i=0,b-1           
539            if (isStrSubset( a(i),"zu_sp" ))then
540               legend_label_z=legend_label_zu
[190]541               if (height_level(0) .NE. -1)then
542                  z=zuh
543               else
544                  z=zu
545               end if
[176]546            else
547               if (isStrSubset( a(i),"zw_sp" ))then
[190]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   
[176]554               end if
555            end if
556         end do 
[218]557         
[585]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
[218]565         min_x=new(dimz,double)
566         max_x=new(dimz,double)
[585]567         
[218]568         plot_h  = new(dimz,graphic)
569         
[176]570         if (isStrSubset(vNam(varn),"x"))then
[218]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
[534]577               res@tiXAxisString = "k~B~x~N~ ("+unit_x+")"
[190]578            else
[218]579               if (norm_height .EQ. 1)then
[1527]580                  res@tiXAxisString = "k~B~x~N~ x z (1)"
[218]581               else
[534]582                  res@tiXAxisString = "k~B~x~N~ (1/m)"
[218]583               end if
[190]584            end if
[218]585            dim_r=dimx
[176]586         else
[218]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
[1527]593               res@tiXAxisString = "k~B~y~N~ ("+unit_x+")"
[190]594            else
[218]595               if (norm_height .EQ. 1)then
[1527]596                  res@tiXAxisString = "k~B~y~N~ x z (1)"
[218]597               else
[1527]598                  res@tiXAxisString = "k~B~y~N~ (1/m)"
[218]599               end if
[190]600            end if
[218]601            dim_r=dimy
[176]602         end if
603       
604         if (sort .EQ. "time")
[218]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
[566]611            res@pmLegendHeightF         = 0.04*\
612                                          (end_time_step-start_time_step+1)
[190]613            do p=dimz-1,0,1
[218]614               if (log_y .EQ. 1)then 
[190]615                  do q=0,dimt-1
[218]616                     do r=0,dim_r-1
617                        if (data(q,p,r) .EQ. 0)then
[190]618                           st=p+start_time_step
619                           print(" ")
[566]620                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is "+\
621                                 "equal 0. Logarithmic scale for y-axis "+\
622                                 "and height "+z(p)+" cannot be used")
[190]623                           print(" ")
624                           res@trYLog = False
625                        end if
626                     end do
[176]627                  end do
[190]628               end if
[176]629               res@gsnLeftString      = vNam(varn)
[218]630               res@gsnRightString     = "Height = "+z(p)+"m"               
[534]631               res@tiYAxisString      = "("+unit_y+")"           
[218]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)
[176]640               n=n+1
641            end do
642         else
[190]643            if (sort .EQ. "height")
644               do p=0,dimt-1           
[176]645                  do q=0,dimz-1
[218]646                     do r=0,dim_r-1
647                        if (data(p,q,r) .EQ. 0)then
[176]648                           st=p+start_time_step
649                           print(" ")
[566]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")
[176]653                           print(" ")
654                           res@trYLog = False
655                        end if
656                     end do
[218]657                     if (norm_height .EQ. 1 .AND. p .EQ. 0)then
[566]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))
[218]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,:))
[566]666                     max_x(q)=max(x_axis(q,:))
[218]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
[566]674                        res@tiYAxisString      = "("+unit_y+")"
[218]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)
[2563]680                        res@trYMaxF = max(max_y*10)
[218]681                       
[566]682                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),\
683                                                         data(p,q,:),res)
[218]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
[566]696                        lgres@lgTitleString      = "Height (m)" 
[218]697                        lgres@lgLabelFont        = "helvetica"
[566]698                        lgres@lgLabelFontHeightF = font_size_legend*6
699                        lgres@lgTitleFontHeightF = font_size     
[218]700                        lgres@vpWidthF           = 0.12           
[566]701                        lgres@vpHeightF          = font_size_legend*(dimz+3)
[218]702 
[566]703                        lbid = gsn_create_legend(wks,dimz,legend_label_z,lgres)
[218]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
[566]710                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),\
711                                                        data(p,q,:),res)
[218]712                        overlay(plot_h(0),plot_h(q))
713                     end if
714                  end do             
715                  plot(n)=plot_h(0)
716                  n=n+1
[176]717               end do
718            end if
[190]719         end if
[218]720         delete(data)
721         delete(temp)
[178]722         delete(x_axis)
[218]723         delete(min_x)
724         delete(max_x)
725         delete(min_y)
726         delete(max_y)
727         delete(plot_h)
[176]728      end if
729   end do
730
731   if (n .EQ. 0) then
732      print(" ")
[190]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")
[176]735      print(" ")
736      exit
737   end if
738 
739   ; ***************************************************
740   ; merge plots onto one page
741   ; ***************************************************
742
[983]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
[176]751
[534]752   no_frames = 0
753
[566]754   if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
755                                          n .gt. no_rows*no_columns) then
[183]756      gsn_panel(wks,plot,(/n,1/),resP)
[250]757      print(" ")
758      print("Outputs to .eps or .epsi have only one frame")
759      print(" ")
[176]760   else   
[218]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)
[534]764            no_frames = no_frames + 1
[176]765         else
[566]766            gsn_panel(wks,plot(i:i+no_rows*no_columns-1),\
767                                           (/no_rows,no_columns/),resP)
[534]768            no_frames = no_frames + 1   
[176]769         end if
770      end do
771   end if
772
[532]773   if (format_out .EQ. "png" ) then
774     png_output = new((/no_frames/), string)
775     j = 0
776
[957]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)
[532]787       system(cmd)
[957]788     else
[532]789
[957]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
[532]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
[176]821   
[190]822end
Note: See TracBrowser for help on using the repository browser.