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

Last change on this file since 510 was 418, checked in by heinze, 15 years ago

ncl_preferences is renamed .ncl.config.default

File size: 19.8 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;***************************************************
[418]6; load .ncl.config or .ncl.config.default
[190]7;***************************************************
[194]8
9function which_script()
10local script
11begin
12   script="spectra"
13   return(script)
14end
[190]15   
[418]16if (isfilepresent("$PALM_BIN/../../.ncl.config")) then
17   loadscript("$PALM_BIN/../../.ncl.config")
[190]18else
[418]19  if (isfilepresent("$PALM_BIN/NCL/.ncl.config.default")) then
20     loadscript( "$PALM_BIN/NCL/.ncl.config.default")
[190]21  else
[418]22      palm_bin_path = getenv("PALM_BIN")
[190]23      print(" ")
[418]24      print("Neither the personal configuration file '.ncl.config' exists in")
25      print("~/palm/current_version")
26      print("nor the default configuration file '.ncl.config.default' exists in")
27      print(palm_bin_path + "/NCL")
[190]28      print(" ")
29      exit
[176]30   end if
[418]31end if   
[218]32
[190]33begin
[176]34
[190]35   ;***************************************************
36   ; set up default parameter values and strings
37   ;***************************************************
38 
39   if (file_1 .EQ. "File in") then
40      print(" ")
[418]41      print("Declare input file 'file_1=' in '.ncl.config' or prompt")
[190]42      print(" ")
43      exit
[176]44   else
45      file_in = file_1   
46   end if
47
[190]48   if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND. format_out .NE. "eps" .AND. format_out .NE. "ps" .AND. format_out .NE. "epsi" .AND. format_out .NE. "ncgm")then
49      print(" ")
50      print("'format_out = "+format_out+"' is invalid and set to'x11'")
51      print(" ")
52      format_out="x11"
53   end if
54   
[218]55   if (log_x .NE. 0 .AND. log_x .NE. 1)then
[190]56      print(" ")
[218]57      print("'log_x'= "+log_x+" is invalid and set to 1")
[190]58      print(" ")
[218]59      log_x = 1
[190]60   end if 
61   
[218]62   if (log_y .NE. 0 .AND. log_y .NE. 1)then
[190]63      print(" ")
[218]64      print("'log_y'= "+log_y+" is invalid and set to 1")
[190]65      print(" ")
[218]66      log_y = 1
[190]67   end if   
68 
[218]69   if (norm_y .EQ. 0.) then
[190]70      print(" ")
[218]71      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
[190]72      print(" ")
[218]73      norm_y = 1.0
[176]74   end if
[218]75   if (norm_x .EQ. 0.) then
[190]76      print(" ")
[218]77      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
[190]78      print(" ")
[218]79      norm_x= 1.0
[176]80   end if
[190]81   
82   if (sort .NE. "height" .AND. sort .NE. "time") then 
83      print(" ")
84      print("'sort'= "+sort+" is invalid and set to 'height'")
85      print(" ")
86      sort = "height" 
[176]87   end if
[190]88   
89   if (black .NE. 0 .AND. black .NE. 1)then
90      print(" ")
91      print("'black'= "+black+" is invalid and set to 0")
92      print(" ")
[176]93      black = 0
94   end if
[190]95 
96   if (dash .NE. 0 .AND. dash .NE. 1)then
97      print(" ")
98      print("'dash'= "+dash+" is invalid and set to 0")
99      print(" ")
[176]100      dash = 0
101   end if
102
[190]103   ;***************************************************
104   ; open input file
105   ;***************************************************
[218]106   
107   file_in_1 = False
108   if (isStrSubset(file_in, ".nc"))then
109      start_f = -2
110      end_f = -2
111      file_in_1 = True     
112   end if
[176]113
[218]114   if (start_f .EQ. -1)then
115      print(" ")
116      print("'start_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
117      print(" ") 
118      exit
119   end if
120   if (end_f .EQ. -1)then
121      print(" ")
122      print("'end_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
123      print(" ") 
124      exit
125   end if
126
127   files=new(end_f-start_f+1,string)
[176]128   
[218]129   if (file_in_1)then
130      if (isfilepresent(file_in))then
131         files(0)=file_in
132      else
133         print(" ")
134         print("1st input file: '"+file_in+"' does not exist")
135         print(" ")
136         exit
137      end if
138   else   
139      if (start_f .EQ. 0)then
140         if (isfilepresent(file_in+".nc"))then
141            files(0)=file_in+".nc"
142            do i=1,end_f
143               if (isfilepresent(file_in+"."+i+".nc"))then   
144                  files(i)=file_in+"."+i+".nc"
145               else
146                  print(" ")
147                  print("Input file: '"+file_in+"."+i+".nc' does not exist")
148                  print(" ")
149                  exit 
150               end if       
151            end do         
152         else
153            print(" ")
154            print("Input file: '"+file_in+".nc' does not exist")
155            print(" ")
156            exit
157         end if
158      else
159         do i=start_f,end_f
160            if (isfilepresent(file_in+"."+i+".nc"))then   
161               files(i-start_f)=file_in+"."+i+".nc"
162            else
163               print(" ")
164               print("Input file: '"+file_in+"."+i+".nc' does not exist")
165               print(" ")
166               exit 
167            end if
168         end do
169      end if
170   end if
171
172   f=addfiles(files,"r")
173   f_att=addfile(files(0),"r")
174   ListSetType(f,"cat")
175   
176   vNam = getfilevarnames(f_att)
[176]177   print(" ")
[190]178   print("Variables in input file:")
179   print("- "+ vNam)
[176]180   print(" ")
181   dim = dimsizes(vNam)
[218]182   vDim = getfiledimsizes(f_att)
[176]183 
[218]184   t_all = f[:]->time
[176]185   nt    = dimsizes(t_all)
186   delta_t=t_all(nt-1)/nt
187   
[218]188   k_x=f_att->k_x
[190]189   dimx=dimsizes(k_x)
[218]190   k_y=f_att->k_y
[190]191   dimy=dimsizes(k_y)
192   
193   
194   dim_level=dimsizes(height_level)
195
[176]196   do i=0,dim-1
197      if (vNam(i) .EQ. "zu_sp")then
[218]198         zu=f_att->zu_sp         
[190]199         if (height_level(0) .EQ. -1)then
200            dimz=dimsizes(zu)
201         else
202            if (dim_level .GT. dimsizes(zu))then
203               print(" ")
204               print("'height_level' has more elements than available height levels in input file (= "+dimz+")")
205               print(" ")
206               exit
207            else
208               zuh=new(dim_level,double)
209               do le=0,dim_level-1
210                  zuh(le)=zu(height_level(le))
211               end do
212               dimz=dim_level
213            end if   
214         end if 
[176]215      else
216         if (vNam(i) .EQ. "zw_sp")then
[218]217            zw=f_att->zw_sp
[190]218            if (height_level(0) .EQ. -1)then             
219               dimz=dimsizes(zw)
220            else
221               if (dim_level .GT. dimsizes(zw))then
222                  print(" ")
223                  print("'height_level' has more elements than available height levels in input file (= "+dimz+")")
224                  print(" ")
225                  exit
226               else
227                  zwh=new(dim_level,double)
228                  do le=0,dim_level-1
229                     zwh(le)=zw(height_level(le))
230                  end do
231                  dimz=dim_level
232               end if   
233            end if
[176]234         end if
235      end if
236   end do
[190]237
238   ;****************************************************       
[176]239   ; start of time step and different types of mistakes that could be done
[190]240   ;****************************************************
[176]241   
[190]242   if (start_time_step .EQ. -1.d) then         
[176]243      start_time_step=t_all(0)/3600
[190]244   else
[176]245      if (start_time_step .GT. t_all(nt-1)/3600)then
246         print(" ")
247         print("'start_time_step' = "+ start_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
248         print(" ")
[218]249         print("Select another 'start_time_step'")
[176]250         print(" ")
251         exit
252      end if
253      if (start_time_step .LT. t_all(0)/3600)then
254         print(" ")
255         print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")       
256         exit
257      end if
258   end if
259
260   do i=0,nt-1     
[190]261      if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[176]262         st=i
263         break
264      end if
265   end do
266   
[190]267   if (.not. isvar("st"))then
268      print(" ")
269      print("'start_time_step' = "+ start_time_step +"h is invalid")
270      print(" ")
[218]271      print("Select another 'start_time_step'")
[190]272      print(" ")
273      exit
274   end if
275   
276   ;****************************************************
[176]277   ; end of time step and different types of mistakes that could be done
[190]278   ;****************************************************
[176]279
[190]280   if (end_time_step .EQ. -1.d) then           
[176]281      end_time_step = t_all(nt-1)/3600
282   else
283      if (end_time_step .GT. t_all(nt-1)/3600)then
284         print(" ")
285         print("'end_time_step' = "+ end_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
286         print(" ")
[218]287         print("Select another 'end_time_step'") 
[176]288         print(" ")
289         exit
290      end if
291      if (end_time_step .LT. start_time_step/3600)then
292         print(" ")
293         print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
294         print(" ")
[218]295         print("Select another 'start_time_step' or 'end_time_step'")
[176]296         print(" ")
297         exit
298      end if
299   end if
300
301   do i=0,nt-1     
[190]302      if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[176]303         et=i
304         break
305      end if
306   end do
[190]307   
308   if (.not. isvar("et"))then
309      print(" ")
310      print("'end_time_step' = "+ end_time_step +"h is invalid")
311      print(" ")
[218]312      print("Select another 'end_time_step'")
[190]313      print(" ")
314      exit
315   end if
[176]316
317   delete(start_time_step)
318   start_time_step=round(st,3)
319   delete(end_time_step)
320   end_time_step=round(et,3)
321
322   print(" ")
323   print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+t_all(start_time_step)+" s => index = "+start_time_step)
324   print("                     till "+t_all(end_time_step)/3600+" h = "+t_all(end_time_step)+" s => index = "+end_time_step)
325   print(" ")
326
327   dimt = end_time_step-start_time_step+1
328 
[190]329   ;***************************************************
[176]330   ; set up recourses
[190]331   ;***************************************************
[176]332     
333   res = True
334   res@gsnDraw                 = False
335   res@gsnFrame                = False
336   res@gsnPaperOrientation     = "portrait"
337   res@gsnPaperWidth           = 8.27
338   res@gsnPaperHeight          = 11.69
339   res@gsnPaperMargin          = 0.79
340   res@txFont                  = "helvetica"
341   res@tiMainFont              = "helvetica"
342   res@tiXAxisFont             = "helvetica"
343   res@tiYAxisFont             = "helvetica"
344   res@tmXBLabelFont           = "helvetica"
345   res@tmYLLabelFont           = "helvetica"
346   res@lgLabelFont             = "helvetica"
347   res@tmLabelAutoStride       = True
348   
[218]349   res@lgLabelFontHeightF     = font_size_legend
350   res@lgTitleFontHeightF     = font_size
351   res@txFontHeightF      = font_size
352   res@tiXAxisFontHeightF = font_size
353   res@tiYAxisFontHeightF = font_size
354   res@tmXBLabelFontHeightF = font_size
355   res@tmYLLabelFontHeightF = font_size
356   
357   res@tmXBMinorPerMajor = 4
358   res@tmYLMinorPerMajor = 4
359   
360   if (log_x .EQ. 1) then
[176]361      res@trXLog = True
[190]362   else
363      res@trXLog = False
364   end if
[218]365   if (log_y .EQ. 1)then
[176]366      res@trYLog = True
[190]367   else 
[178]368      res@trYLog = False
[176]369   end if
370
371   legend_label=new(dimt,double)
372   legend_label_zu=new(dimz,double)
373   legend_label_zw=new(dimz,double)
374   legend_label_z=new(dimz,double)
375   do p=start_time_step,end_time_step
376      if (t_all(p)/3600 .LT. 1) then
377         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,2,True)
378      else
379         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,0,True)
380      end if
381   end do 
382   if (sort .EQ. "time")
383      plot  = new(dim*dimz,graphic)
384      np=dimt
385      res@lgTitleString = "Time [h]"
386   else
387      plot  = new(dim*dimt,graphic)
388      np=dimz
[218]389     
[176]390      do p=0,dimz-1
[190]391         if (height_level(0) .EQ. -1)then
392            legend_label_zu(p)=round(zu(p),3)
393            legend_label_zw(p)=round(zw(p),3)
394         else
395            legend_label_zu(p)=round(zuh(p),3)
396            legend_label_zw(p)=round(zwh(p),3)
397         end if
[176]398      end do
399   end if
[194]400
401   if (black .eq. 0 ) then
402      if (np .EQ. 1)then
[218]403         color = 237
[194]404      else   
405         step=round(235/(np-1),3)
[218]406         color = ispan(2,237,step)
[194]407      end if
[324]408   else
409      color = 2
[176]410   end if
411   if ( dash .eq. 0 ) then
[190]412      res@xyMonoDashPattern = True
[176]413   end if
414
415   wks=gsn_open_wks(format_out,file_out)
416   gsn_define_colormap(wks,"rainbow+white")
417
418   n=0
419   do varn =dim-1,0,1
[178]420     
[176]421      check = True
422
423      if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "zu_sp") .OR. isStrSubset( vNam(varn), "zw_sp") .OR. isStrSubset( vNam(varn), "k_x") .OR. isStrSubset( vNam(varn), "k_y")) then
424            check = False
425      end if
426
[190]427      if (var .NE. "all") then
[176]428         check = isStrSubset( var,","+vNam(varn)+"," )
429      end if
430
[190]431      if(check) then
[218]432
433         temp = f[:]->$vNam(varn)$
434         data = temp(start_time_step:end_time_step,0:dimz-1,:)
435
436         temp_att = f_att->$vNam(varn)$
437         a=getvardims(temp_att)
[190]438         b=dimsizes(a)
[176]439
[190]440         if (height_level(0) .NE. -1)then
441            do te=0,dimz-1
[218]442            print(height_level(te))
443               data(:,te,:) = temp(start_time_step:end_time_step,height_level(te),:)           
[190]444            end do
445         end if
446
[218]447         data=data/(norm_y*norm_x)  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
[190]448           
[176]449         do i=0,b-1           
450            if (isStrSubset( a(i),"zu_sp" ))then
451               legend_label_z=legend_label_zu
[190]452               if (height_level(0) .NE. -1)then
453                  z=zuh
454               else
455                  z=zu
456               end if
[176]457            else
458               if (isStrSubset( a(i),"zw_sp" ))then
[190]459                  legend_label_z=legend_label_zw
460                  if (height_level(0) .NE. -1)then
461                     z=zwh
462                  else
463                     z=zw
464                  end if   
[176]465               end if
466            end if
467         end do 
[218]468         
469         min_x=new(dimz,double)
470         max_x=new(dimz,double)
471         min_y=new(dimz,float)
472         max_y=new(dimz,float)
473         plot_h  = new(dimz,graphic)
474         
[176]475         if (isStrSubset(vNam(varn),"x"))then
[218]476            x_axis = new((/dimz,dimx/),double)
477            do q=0,dimz-1
478               x_axis(q,:) = f_att->k_x
479               x_axis = x_axis/norm_x
480            end do
481            if (norm_x .NE. 1.)then
482               res@tiXAxisString = "k~B~x~N~ ["+unit_x+"]"
[190]483            else
[218]484               if (norm_height .EQ. 1)then
485                  res@tiXAxisString = "k~B~x~N~ x z [1]"
486               else
487                  res@tiXAxisString = "k~B~x~N~ [1/m]"
488               end if
[190]489            end if
[218]490            dim_r=dimx
[176]491         else
[218]492            x_axis=new((/dimz,dimy/),double)
493            do q=0,dimz-1
494               x_axis(q,:) = f_att->k_y
495               x_axis = x_axis/norm_x
496            end do
497            if (norm_x .NE. 1.)then
498               res@tiXAxisString = "k~B~x~N~ ["+unit_x+"]"
[190]499            else
[218]500               if (norm_height .EQ. 1)then
501                  res@tiXAxisString = "k~B~x~N~ x z [1]"
502               else
503                  res@tiXAxisString = "k~B~x~N~ [1/m]"
504               end if
[190]505            end if
[218]506            dim_r=dimy
[176]507         end if
508       
509         if (sort .EQ. "time")
[218]510            res@xyLineColors = color
511            res@pmLegendDisplayMode     = "Always"
512            res@pmLegendSide            = "Top"
513            res@pmLegendParallelPosF    = 1.2
514            res@pmLegendOrthogonalPosF  = -1.0
515            res@pmLegendWidthF          = 0.12
516            res@pmLegendHeightF         = 0.04*(end_time_step-start_time_step+1)
[190]517            do p=dimz-1,0,1
[218]518               if (log_y .EQ. 1)then 
[190]519                  do q=0,dimt-1
[218]520                     do r=0,dim_r-1
521                        if (data(q,p,r) .EQ. 0)then
[190]522                           st=p+start_time_step
523                           print(" ")
524                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and height "+z(p)+" cannot be used")
525                           print(" ")
526                           res@trYLog = False
527                        end if
528                     end do
[176]529                  end do
[190]530               end if
[176]531               res@gsnLeftString      = vNam(varn)
[218]532               res@gsnRightString     = "Height = "+z(p)+"m"               
533               res@tiYAxisString      = "["+unit_y+"]"           
534               res@xyExplicitLegendLabels  = legend_label             
535               if (norm_height .EQ. 1)then
536                  data(:,p,:)=data(:,p,:)*doubletofloat(z(p))
537                  x_axis(p,:) = x_axis(p,:)*z(p)
538               end if
539               res@trXMinF = min(x_axis(p,:))
540               res@trXMaxF = max(x_axis(p,:))
541               plot(n)  = gsn_csm_xy(wks,x_axis(p,:),data(:,p,:),res)
[176]542               n=n+1
543            end do
544         else
[190]545            if (sort .EQ. "height")
546               do p=0,dimt-1           
[176]547                  do q=0,dimz-1
[218]548                     do r=0,dim_r-1
549                        if (data(p,q,r) .EQ. 0)then
[176]550                           st=p+start_time_step
551                           print(" ")
[190]552                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and time "+legend_label(p)+" cannot be used")
[176]553                           print(" ")
554                           res@trYLog = False
555                        end if
556                     end do
[218]557                     if (norm_height .EQ. 1 .AND. p .EQ. 0)then
558                        data(p,q,:) = data(p,q,:)*doubletofloat(legend_label_z(q))
559                        x_axis(q,:) = x_axis(q,:)*doubletofloat(legend_label_z(q))
560                     end if
561                     max_y(q)=max(data(p,q,:))
562                     min_y(q)=min(data(p,q,:))
563                     min_x(q)=min(x_axis(q,:))
564                     max_x(q)=max(x_axis(q,:))                                   
565                  end do
566                  do q=0,dimz-1
567                     res@xyLineColor = color(q)
568                     if (dash .EQ. 1)then
569                        res@xyDashPattern = q
570                     end if
571                     if (q .EQ. 0)then
572                        res@tiYAxisString      = "["+unit_y+"]"                 
573                        res@gsnLeftString      = vNam(varn)
574                        res@gsnRightString     = "Time = "+legend_label(p)+"h"
575                        res@trXMinF = min(min_x)
576                        res@trXMaxF = max(max_x)
577                        res@trYMinF = min(min_y)
578                        res@trYMaxF = max(max_y)
579                       
580                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
581                       
582                        lgres = True
583                        if (dash .EQ. 0)then
584                           lgres@lgMonoDashIndex = True
585                        else
586                           lgres@lgDashIndexes   = ispan(0,dimz-1,1)
587                        end if
588                        if (black .EQ. 1)then
589                           lgres@lgMonoLineColors = True
590                        else
591                           lgres@lgLineColors = color
592                        end if
593                        lgres@lgTitleString = "Height [m]" 
594                        lgres@lgLabelFont        = "helvetica"
595                        lgres@lgLabelFontHeightF     = font_size_legend
596                        lgres@lgTitleFontHeightF     = font_size     
597                        lgres@vpWidthF           = 0.12           
598                        lgres@vpHeightF          = font_size_legend/5*dimz           
599 
600                        lbid = gsn_create_legend(wks,dimz,legend_label_z,lgres)       
601
602                        amres = True
603                        amres@amParallelPosF   = 0.75               
604                        amres@amOrthogonalPosF = 0.15           
605                        annoid1 = gsn_add_annotation(plot_h(q),lbid,amres)
606                     else
607                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
608                        overlay(plot_h(0),plot_h(q))
609                     end if
610                  end do             
611                  plot(n)=plot_h(0)
612                  n=n+1
[176]613               end do
614            end if
[190]615         end if
[218]616         delete(data)
617         delete(temp)
[178]618         delete(x_axis)
[218]619         delete(min_x)
620         delete(max_x)
621         delete(min_y)
622         delete(max_y)
623         delete(plot_h)
[176]624      end if
625   end do
626
627   if (n .EQ. 0) then
628      print(" ")
[190]629      print("The variables 'var="+var+"' do not exist on your input file;")
630      print("be sure to have one comma berfore and after each variable")
[176]631      print(" ")
632      exit
633   end if
634 
635   ; ***************************************************
636   ; merge plots onto one page
637   ; ***************************************************
638
639   resP                        = True
640   resP@txFont                 = "helvetica"
[218]641   resP@txString               = f_att@title
[176]642   resP@txFuncCode             = "~"
[218]643   resP@txFontHeightF          = 0.015
[176]644
[418]645   if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. n .gt. no_rows*no_columns) then
[183]646      gsn_panel(wks,plot,(/n,1/),resP)
[250]647      print(" ")
648      print("Outputs to .eps or .epsi have only one frame")
649      print(" ")
[176]650   else   
[218]651      do i = 0,n-1, no_rows*no_columns
652         if( (i+no_rows*no_columns) .gt. (n-1)) then
653            gsn_panel(wks,plot(i:n-1),(/no_rows,no_columns/),resP)
[176]654         else
[218]655           gsn_panel(wks,plot(i:i+no_rows*no_columns-1),(/no_rows,no_columns/),resP)
[176]656         end if
657      end do
658   end if
659
660   print(" ")
661   print("Output to: " + file_out +"."+ format_out)
662   print(" ")
663   
[190]664end
Note: See TracBrowser for help on using the repository browser.