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

Last change on this file since 2172 was 2030, checked in by gronemeier, 7 years ago

bugfix: checking for NCL version number

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
[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   
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.