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

Last change on this file since 1762 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
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
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
[418]46   loadscript("$PALM_BIN/../../.ncl.config")
[190]47else
[1559]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
[418]53      palm_bin_path = getenv("PALM_BIN")
[190]54      print(" ")
[418]55      print("Neither the personal configuration file '.ncl.config' exists in")
56      print("~/palm/current_version")
[566]57      print("nor the default configuration file '.ncl.config.default' "+\
58            "exists in")
[418]59      print(palm_bin_path + "/NCL")
[190]60      print(" ")
61      exit
[176]62   end if
[1559]63end if
[218]64
[190]65begin
[176]66
[190]67   ;***************************************************
[526]68   ; Retrieving the double quote character
69   ;***************************************************
70   
71   dq=str_get_dq()
72
73   ;***************************************************
[190]74   ; set up default parameter values and strings
75   ;***************************************************
76 
77   if (file_1 .EQ. "File in") then
78      print(" ")
[418]79      print("Declare input file 'file_1=' in '.ncl.config' or prompt")
[190]80      print(" ")
81      exit
[176]82   else
83      file_in = file_1   
84   end if
85
[566]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
[190]90      print(" ")
91      print("'format_out = "+format_out+"' is invalid and set to'x11'")
92      print(" ")
93      format_out="x11"
94   end if
[532]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
[190]104   
[218]105   if (log_x .NE. 0 .AND. log_x .NE. 1)then
[190]106      print(" ")
[218]107      print("'log_x'= "+log_x+" is invalid and set to 1")
[190]108      print(" ")
[218]109      log_x = 1
[190]110   end if 
111   
[218]112   if (log_y .NE. 0 .AND. log_y .NE. 1)then
[190]113      print(" ")
[218]114      print("'log_y'= "+log_y+" is invalid and set to 1")
[190]115      print(" ")
[218]116      log_y = 1
[190]117   end if   
118 
[218]119   if (norm_y .EQ. 0.) then
[190]120      print(" ")
[218]121      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
[190]122      print(" ")
[218]123      norm_y = 1.0
[176]124   end if
[218]125   if (norm_x .EQ. 0.) then
[190]126      print(" ")
[218]127      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
[190]128      print(" ")
[218]129      norm_x= 1.0
[176]130   end if
[190]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" 
[176]137   end if
[190]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(" ")
[176]143      black = 0
144   end if
[190]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(" ")
[176]150      dash = 0
151   end if
152
[190]153   ;***************************************************
154   ; open input file
155   ;***************************************************
[218]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
[176]163
[218]164   if (start_f .EQ. -1)then
165      print(" ")
[566]166      print("'start_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   if (end_f .EQ. -1)then
172      print(" ")
[566]173      print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
174            "your input file(s)")
[218]175      print(" ") 
176      exit
177   end if
178
179   files=new(end_f-start_f+1,string)
[176]180   
[218]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)
[1248]229   if( ncl_version .GE. 6.1 ) then
230      vNam = vNam(::-1)
231   end if
[585]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
[176]240   print(" ")
[190]241   print("Variables in input file:")
242   print("- "+ vNam)
[176]243   print(" ")
244   dim = dimsizes(vNam)
[218]245   vDim = getfiledimsizes(f_att)
[176]246 
[218]247   t_all = f[:]->time
[176]248   nt    = dimsizes(t_all)
[827]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
[176]262   
[218]263   k_x=f_att->k_x
[190]264   dimx=dimsizes(k_x)
[218]265   k_y=f_att->k_y
[190]266   dimy=dimsizes(k_y)
267   
268   
269   dim_level=dimsizes(height_level)
270
[176]271   do i=0,dim-1
272      if (vNam(i) .EQ. "zu_sp")then
[218]273         zu=f_att->zu_sp         
[190]274         if (height_level(0) .EQ. -1)then
275            dimz=dimsizes(zu)
276         else
277            if (dim_level .GT. dimsizes(zu))then
278               print(" ")
[566]279               print("'height_level' has more elements than available "+\
280                     "height levels in input file (= "+dimsizes(zu)+")")
[190]281               print(" ")
282               exit
283            else
284               zuh=new(dim_level,double)
285               do le=0,dim_level-1
[566]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
[190]296                  zuh(le)=zu(height_level(le))
297               end do
298               dimz=dim_level
299            end if   
300         end if 
[176]301      else
302         if (vNam(i) .EQ. "zw_sp")then
[218]303            zw=f_att->zw_sp
[190]304            if (height_level(0) .EQ. -1)then             
305               dimz=dimsizes(zw)
306            else
307               if (dim_level .GT. dimsizes(zw))then
308                  print(" ")
[566]309                  print("'height_level' has more elements than available "+\
310                        "height levels in input file (= "+dimsizes(zw)+")")
[190]311                  print(" ")
312                  exit
313               else
314                  zwh=new(dim_level,double)
315                  do le=0,dim_level-1
[566]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
[190]326                     zwh(le)=zw(height_level(le))
327                  end do
328                  dimz=dim_level
329               end if   
330            end if
[176]331         end if
332      end if
333   end do
[190]334
335   ;****************************************************       
[176]336   ; start of time step and different types of mistakes that could be done
[190]337   ;****************************************************
[176]338   
[190]339   if (start_time_step .EQ. -1.d) then         
[176]340      start_time_step=t_all(0)/3600
[190]341   else
[176]342      if (start_time_step .GT. t_all(nt-1)/3600)then
343         print(" ")
[566]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")
[176]346         print(" ")
[218]347         print("Select another 'start_time_step'")
[176]348         print(" ")
349         exit
350      end if
351      if (start_time_step .LT. t_all(0)/3600)then
352         print(" ")
[566]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")       
[176]355         exit
356      end if
357   end if
358
359   do i=0,nt-1     
[566]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
[176]362         st=i
363         break
[769]364      else
365         st=0
[176]366      end if
367   end do
368   
[190]369   ;****************************************************
[176]370   ; end of time step and different types of mistakes that could be done
[190]371   ;****************************************************
[176]372
[190]373   if (end_time_step .EQ. -1.d) then           
[176]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(" ")
[566]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")
[176]380         print(" ")
[218]381         print("Select another 'end_time_step'") 
[176]382         print(" ")
383         exit
384      end if
385      if (end_time_step .LT. start_time_step/3600)then
386         print(" ")
[566]387         print("'end_time_step' = "+ end_time_step +"h is lower than "+\
388               "'start_time_step' = "+start_time_step+"h")
[176]389         print(" ")
[218]390         print("Select another 'start_time_step' or 'end_time_step'")
[176]391         print(" ")
392         exit
393      end if
394   end if
395
396   do i=0,nt-1     
[566]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
[176]399         et=i
400         break
[769]401       else
402         et=0
[176]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(" ")
[566]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)
[176]416   print(" ")
417
418   dimt = end_time_step-start_time_step+1
419 
[190]420   ;***************************************************
[176]421   ; set up recourses
[190]422   ;***************************************************
[176]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   
[218]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
[176]448      res@trXLog = True
[190]449   else
450      res@trXLog = False
451   end if
[218]452   if (log_y .EQ. 1)then
[176]453      res@trYLog = True
[190]454   else 
[178]455      res@trYLog = False
[176]456   end if
457
[566]458   legend_label=new(dimt,string)
[176]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
[566]463      legend_label(p-start_time_step)=sprintf("%6.2f", t_all(p)/3600)
[176]464   end do 
465   if (sort .EQ. "time")
466      plot  = new(dim*dimz,graphic)
467      np=dimt
[534]468      res@lgTitleString = "Time (h)"
[176]469   else
470      plot  = new(dim*dimt,graphic)
471      np=dimz
[218]472     
[176]473      do p=0,dimz-1
[190]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
[176]481      end do
482   end if
[194]483
484   if (black .eq. 0 ) then
485      if (np .EQ. 1)then
[218]486         color = 237
[194]487      else   
488         step=round(235/(np-1),3)
[218]489         color = ispan(2,237,step)
[194]490      end if
[324]491   else
492      color = 2
[176]493   end if
494   if ( dash .eq. 0 ) then
[190]495      res@xyMonoDashPattern = True
[176]496   end if
497
[532]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
[176]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
[178]511     
[176]512      check = True
513
[566]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
[176]519            check = False
520      end if
521
[190]522      if (var .NE. "all") then
[176]523         check = isStrSubset( var,","+vNam(varn)+"," )
524      end if
525
[190]526      if(check) then
[218]527
528         temp = f[:]->$vNam(varn)$
[566]529         data = temp(start_time_step:end_time_step,0:dimz-1,:)
[218]530
531         temp_att = f_att->$vNam(varn)$
532         a=getvardims(temp_att)
[190]533         b=dimsizes(a)
[534]534         delete(temp_att)
[176]535
[190]536         if (height_level(0) .NE. -1)then
537            do te=0,dimz-1
[566]538               data(:,te,:) = temp(start_time_step:end_time_step,\
539                                                    height_level(te),:)       
[190]540            end do
541         end if
542
[566]543         data=data/(norm_y*norm_x)
[190]544           
[176]545         do i=0,b-1           
546            if (isStrSubset( a(i),"zu_sp" ))then
547               legend_label_z=legend_label_zu
[190]548               if (height_level(0) .NE. -1)then
549                  z=zuh
550               else
551                  z=zu
552               end if
[176]553            else
554               if (isStrSubset( a(i),"zw_sp" ))then
[190]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   
[176]561               end if
562            end if
563         end do 
[218]564         
[585]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
[218]572         min_x=new(dimz,double)
573         max_x=new(dimz,double)
[585]574         
[218]575         plot_h  = new(dimz,graphic)
576         
[176]577         if (isStrSubset(vNam(varn),"x"))then
[218]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
[534]584               res@tiXAxisString = "k~B~x~N~ ("+unit_x+")"
[190]585            else
[218]586               if (norm_height .EQ. 1)then
[1527]587                  res@tiXAxisString = "k~B~x~N~ x z (1)"
[218]588               else
[534]589                  res@tiXAxisString = "k~B~x~N~ (1/m)"
[218]590               end if
[190]591            end if
[218]592            dim_r=dimx
[176]593         else
[218]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
[1527]600               res@tiXAxisString = "k~B~y~N~ ("+unit_x+")"
[190]601            else
[218]602               if (norm_height .EQ. 1)then
[1527]603                  res@tiXAxisString = "k~B~y~N~ x z (1)"
[218]604               else
[1527]605                  res@tiXAxisString = "k~B~y~N~ (1/m)"
[218]606               end if
[190]607            end if
[218]608            dim_r=dimy
[176]609         end if
610       
611         if (sort .EQ. "time")
[218]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
[566]618            res@pmLegendHeightF         = 0.04*\
619                                          (end_time_step-start_time_step+1)
[190]620            do p=dimz-1,0,1
[218]621               if (log_y .EQ. 1)then 
[190]622                  do q=0,dimt-1
[218]623                     do r=0,dim_r-1
624                        if (data(q,p,r) .EQ. 0)then
[190]625                           st=p+start_time_step
626                           print(" ")
[566]627                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is "+\
628                                 "equal 0. Logarithmic scale for y-axis "+\
629                                 "and height "+z(p)+" cannot be used")
[190]630                           print(" ")
631                           res@trYLog = False
632                        end if
633                     end do
[176]634                  end do
[190]635               end if
[176]636               res@gsnLeftString      = vNam(varn)
[218]637               res@gsnRightString     = "Height = "+z(p)+"m"               
[534]638               res@tiYAxisString      = "("+unit_y+")"           
[218]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)
[176]647               n=n+1
648            end do
649         else
[190]650            if (sort .EQ. "height")
651               do p=0,dimt-1           
[176]652                  do q=0,dimz-1
[218]653                     do r=0,dim_r-1
654                        if (data(p,q,r) .EQ. 0)then
[176]655                           st=p+start_time_step
656                           print(" ")
[566]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")
[176]660                           print(" ")
661                           res@trYLog = False
662                        end if
663                     end do
[218]664                     if (norm_height .EQ. 1 .AND. p .EQ. 0)then
[566]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))
[218]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,:))
[566]673                     max_x(q)=max(x_axis(q,:))
[218]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
[566]681                        res@tiYAxisString      = "("+unit_y+")"
[218]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                       
[566]689                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),\
690                                                         data(p,q,:),res)
[218]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
[566]703                        lgres@lgTitleString      = "Height (m)" 
[218]704                        lgres@lgLabelFont        = "helvetica"
[566]705                        lgres@lgLabelFontHeightF = font_size_legend*6
706                        lgres@lgTitleFontHeightF = font_size     
[218]707                        lgres@vpWidthF           = 0.12           
[566]708                        lgres@vpHeightF          = font_size_legend*(dimz+3)
[218]709 
[566]710                        lbid = gsn_create_legend(wks,dimz,legend_label_z,lgres)
[218]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
[566]717                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),\
718                                                        data(p,q,:),res)
[218]719                        overlay(plot_h(0),plot_h(q))
720                     end if
721                  end do             
722                  plot(n)=plot_h(0)
723                  n=n+1
[176]724               end do
725            end if
[190]726         end if
[218]727         delete(data)
728         delete(temp)
[178]729         delete(x_axis)
[218]730         delete(min_x)
731         delete(max_x)
732         delete(min_y)
733         delete(max_y)
734         delete(plot_h)
[176]735      end if
736   end do
737
738   if (n .EQ. 0) then
739      print(" ")
[190]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")
[176]742      print(" ")
743      exit
744   end if
745 
746   ; ***************************************************
747   ; merge plots onto one page
748   ; ***************************************************
749
[983]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
[176]758
[534]759   no_frames = 0
760
[566]761   if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
762                                          n .gt. no_rows*no_columns) then
[183]763      gsn_panel(wks,plot,(/n,1/),resP)
[250]764      print(" ")
765      print("Outputs to .eps or .epsi have only one frame")
766      print(" ")
[176]767   else   
[218]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)
[534]771            no_frames = no_frames + 1
[176]772         else
[566]773            gsn_panel(wks,plot(i:i+no_rows*no_columns-1),\
774                                           (/no_rows,no_columns/),resP)
[534]775            no_frames = no_frames + 1   
[176]776         end if
777      end do
778   end if
779
[532]780   if (format_out .EQ. "png" ) then
781     png_output = new((/no_frames/), string)
782     j = 0
783
[957]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)
[532]794       system(cmd)
[957]795     else
[532]796
[957]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
[532]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
[176]828   
[190]829end
Note: See TracBrowser for help on using the repository browser.