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

Last change on this file since 190 was 190, checked in by letzel, 16 years ago
  • NCL scripts in trunk/SCRIPTS/NCL updated to version 2.0 (.ncl_preferences replaced by ncl_preferences.ncl)
File size: 15.7 KB
Line 
1load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
3load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
4
5;***************************************************
6; load ncl_preferences.ncl
7;***************************************************
8   
9if (isfilepresent("~/ncl_preferences.ncl")) then
10   loadscript("~/ncl_preferences.ncl")
11else
12  if (isfilepresent("~/palm/current_version/trunk/SCRIPTS/NCL/ncl_preferences.ncl")) then
13     loadscript( "~/palm/current_version/trunk/SCRIPTS/NCL/ncl_preferences.ncl")
14  else
15      print(" ")
16      print("'ncl_preferences.ncl' does not exist in $home or $home/palm/current_version/trunk/SCRIPTS/NCL/")
17      print(" ")
18      exit
19   end if
20end if
21   
22begin
23     
24   if (cross_sections .NE. 0 .OR. profiles .NE. 0 .OR. timeseries .NE. 0 .OR. spectra .NE. 1)then
25      print(" ")
26      print("Please specify the used script in 'ncl_preferences.ncl' (Line 7-10)")
27      print(" ")
28      print("Set 'spectra' equal 1 and the other variables equal 0")
29      print(" ")
30      exit
31   end if
32
33   ;***************************************************
34   ; set up default parameter values and strings
35   ;***************************************************
36 
37   if (file_1 .EQ. "File in") then
38      print(" ")
39      print("Please provide 1st input file 'file_1=' in 'ncl_preferences.ncl'")
40      print(" ")
41      exit
42   else
43      file_in = file_1   
44   end if
45   if (.not. isfilepresent(file_in)) then
46      print(" ")
47      print("1st input file: '"+file_in+"' does not exist")
48      print(" ")
49      exit
50   end if
51
52   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
53      print(" ")
54      print("'format_out = "+format_out+"' is invalid and set to'x11'")
55      print(" ")
56      format_out="x11"
57   end if
58   
59   if (logx .NE. 0 .AND. logx .NE. 1)then
60      print(" ")
61      print("'logx'= "+logx+" is invalid and set to 1")
62      print(" ")
63      logx = 1
64   end if 
65   
66   if (logy .NE. 0 .AND. logy .NE. 1)then
67      print(" ")
68      print("'logy'= "+logy+" is invalid and set to 1")
69      print(" ")
70      logy = 1
71   end if   
72 
73   if (normy .EQ. 0.) then
74      print(" ")
75      print("You cannot normalise the y-axis with 0, 'normy' is set to 1.0")
76      print(" ")
77      normy = 1.0
78   end if
79   if (normx .EQ. 0.) then
80      print(" ")
81      print("You cannot normalise the x-axis with 0, 'normx' is set to 1.0")
82      print(" ")
83      normx= 1.0
84   end if
85   
86   if (sort .NE. "height" .AND. sort .NE. "time") then 
87      print(" ")
88      print("'sort'= "+sort+" is invalid and set to 'height'")
89      print(" ")
90      sort = "height" 
91   end if
92   
93   if (black .NE. 0 .AND. black .NE. 1)then
94      print(" ")
95      print("'black'= "+black+" is invalid and set to 0")
96      print(" ")
97      black = 0
98   end if
99 
100   if (dash .NE. 0 .AND. dash .NE. 1)then
101      print(" ")
102      print("'dash'= "+dash+" is invalid and set to 0")
103      print(" ")
104      dash = 0
105   end if
106
107   ;***************************************************
108   ; open input file
109   ;***************************************************
110
111   f=addfile(file_in,"r")
112   
113   vNam = getfilevarnames(f)
114   print(" ")
115   print("Variables in input file:")
116   print("- "+ vNam)
117   print(" ")
118   dim = dimsizes(vNam)
119   vDim = getfiledimsizes(f)
120 
121   t_all = f->time
122   nt    = dimsizes(t_all)
123   delta_t=t_all(nt-1)/nt
124   
125   k_x=f->k_x
126   dimx=dimsizes(k_x)
127   k_y=f->k_y
128   dimy=dimsizes(k_y)
129   
130   
131   dim_level=dimsizes(height_level)
132
133   do i=0,dim-1
134      if (vNam(i) .EQ. "zu_sp")then
135         zu=f->zu_sp     
136         if (height_level(0) .EQ. -1)then
137            dimz=dimsizes(zu)
138         else
139            if (dim_level .GT. dimsizes(zu))then
140               print(" ")
141               print("'height_level' has more elements than available height levels in input file (= "+dimz+")")
142               print(" ")
143               exit
144            else
145               zuh=new(dim_level,double)
146               do le=0,dim_level-1
147                  zuh(le)=zu(height_level(le))
148               end do
149               dimz=dim_level
150            end if   
151         end if 
152      else
153         if (vNam(i) .EQ. "zw_sp")then
154            zw=f->zw_sp
155            if (height_level(0) .EQ. -1)then             
156               dimz=dimsizes(zw)
157            else
158               if (dim_level .GT. dimsizes(zw))then
159                  print(" ")
160                  print("'height_level' has more elements than available height levels in input file (= "+dimz+")")
161                  print(" ")
162                  exit
163               else
164                  zwh=new(dim_level,double)
165                  do le=0,dim_level-1
166                     zwh(le)=zw(height_level(le))
167                  end do
168                  dimz=dim_level
169               end if   
170            end if
171         end if
172      end if
173   end do
174
175   ;****************************************************       
176   ; start of time step and different types of mistakes that could be done
177   ;****************************************************
178   
179   if (start_time_step .EQ. -1.d) then         
180      start_time_step=t_all(0)/3600
181   else
182      if (start_time_step .GT. t_all(nt-1)/3600)then
183         print(" ")
184         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")
185         print(" ")
186         print("Please select another 'start_time_step'")
187         print(" ")
188         exit
189      end if
190      if (start_time_step .LT. t_all(0)/3600)then
191         print(" ")
192         print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")       
193         exit
194      end if
195   end if
196
197   do i=0,nt-1     
198      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
199         st=i
200         break
201      end if
202   end do
203   
204   if (.not. isvar("st"))then
205      print(" ")
206      print("'start_time_step' = "+ start_time_step +"h is invalid")
207      print(" ")
208      print("Please select another 'start_time_step'")
209      print(" ")
210      exit
211   end if
212   
213   ;****************************************************
214   ; end of time step and different types of mistakes that could be done
215   ;****************************************************
216
217   if (end_time_step .EQ. -1.d) then           
218      end_time_step = t_all(nt-1)/3600
219   else
220      if (end_time_step .GT. t_all(nt-1)/3600)then
221         print(" ")
222         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")
223         print(" ")
224         print("Please select another 'end_time_step'") 
225         print(" ")
226         exit
227      end if
228      if (end_time_step .LT. start_time_step/3600)then
229         print(" ")
230         print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
231         print(" ")
232         print("Please select another 'start_time_step' or 'end_time_step'")
233         print(" ")
234         exit
235      end if
236   end if
237
238   do i=0,nt-1     
239      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
240         et=i
241         break
242      end if
243   end do
244   
245   if (.not. isvar("et"))then
246      print(" ")
247      print("'end_time_step' = "+ end_time_step +"h is invalid")
248      print(" ")
249      print("Please select another 'end_time_step'")
250      print(" ")
251      exit
252   end if
253
254   delete(start_time_step)
255   start_time_step=round(st,3)
256   delete(end_time_step)
257   end_time_step=round(et,3)
258
259   print(" ")
260   print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+t_all(start_time_step)+" s => index = "+start_time_step)
261   print("                     till "+t_all(end_time_step)/3600+" h = "+t_all(end_time_step)+" s => index = "+end_time_step)
262   print(" ")
263
264   dimt = end_time_step-start_time_step+1
265 
266   ;***************************************************
267   ; set up recourses
268   ;***************************************************
269     
270   res = True
271   res@gsnDraw                 = False
272   res@gsnFrame                = False
273   res@gsnPaperOrientation     = "portrait"
274   res@gsnPaperWidth           = 8.27
275   res@gsnPaperHeight          = 11.69
276   res@gsnPaperMargin          = 0.79
277   res@txFont                  = "helvetica"
278   res@tiMainFont              = "helvetica"
279   res@tiXAxisFont             = "helvetica"
280   res@tiYAxisFont             = "helvetica"
281   res@tmXBLabelFont           = "helvetica"
282   res@tmYLLabelFont           = "helvetica"
283   res@lgLabelFont             = "helvetica"
284   res@tmLabelAutoStride       = True
285   res@pmLegendDisplayMode     = "Always"
286   res@pmLegendSide            = "Top"
287   res@pmLegendParallelPosF    = 1.2
288   res@pmLegendOrthogonalPosF  = -1.0
289   res@pmLegendWidthF          = 0.12
290   if (sort .EQ. "time")then
291      res@pmLegendHeightF         = 0.04*(end_time_step-start_time_step+1)
292   else
293      res@pmLegendHeightF         = 0.04*dimz
294   end if
295   res@lgLabelFontHeightF     = .025
296   res@lgTitleFontHeightF     = .025
297   res@txFontHeightF      = 0.025
298   res@tiXAxisFontHeightF = 0.025
299   res@tiYAxisFontHeightF = 0.025
300   
301   if (logx .EQ. 1) then
302      res@trXLog = True
303   else
304      res@trXLog = False
305   end if
306   if (logy .EQ. 1)then
307      res@trYLog = True
308   else 
309      res@trYLog = False
310   end if
311
312   legend_label=new(dimt,double)
313   legend_label_zu=new(dimz,double)
314   legend_label_zw=new(dimz,double)
315   legend_label_z=new(dimz,double)
316   do p=start_time_step,end_time_step
317      if (t_all(p)/3600 .LT. 1) then
318         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,2,True)
319      else
320         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,0,True)
321      end if
322   end do 
323   if (sort .EQ. "time")
324      plot  = new(dim*dimz,graphic)
325      np=dimt
326      res@lgTitleString = "Time [h]"
327   else
328      plot  = new(dim*dimt,graphic)
329      np=dimz
330      res@lgTitleString = "Height [m]" 
331      do p=0,dimz-1
332         if (height_level(0) .EQ. -1)then
333            legend_label_zu(p)=round(zu(p),3)
334            legend_label_zw(p)=round(zw(p),3)
335         else
336            legend_label_zu(p)=round(zuh(p),3)
337            legend_label_zw(p)=round(zwh(p),3)
338         end if
339      end do
340   end if
341   
342   step=round(235/(np-1),3)
343   if (black .eq. 0 ) then 
344      res@xyLineColors = ispan(2,237,step)
345   end if
346   if ( dash .eq. 0 ) then
347      res@xyMonoDashPattern = True
348   end if
349
350   wks=gsn_open_wks(format_out,file_out)
351   gsn_define_colormap(wks,"rainbow+white")
352
353   temp=new((/dimt,dimz,dimx/),float)
354
355   n=0
356   do varn =dim-1,0,1
357     
358      check = True
359
360      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
361            check = False
362      end if
363
364      if (var .NE. "all") then
365         check = isStrSubset( var,","+vNam(varn)+"," )
366      end if
367
368      if(check) then
369     
370         temp = f->$vNam(varn)$(start_time_step:end_time_step,0:dimz-1,:)
371           
372         a=getvardims(temp)
373         b=dimsizes(a)
374
375         if (height_level(0) .NE. -1)then
376            do te=0,dimz-1
377               temp(:,te,:) = f->$vNam(varn)$(start_time_step:end_time_step,height_level(te),:)       
378            end do
379         end if
380
381         temp=temp/(normy*normx)  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
382           
383         do i=0,b-1           
384            if (isStrSubset( a(i),"zu_sp" ))then
385               legend_label_z=legend_label_zu
386               if (height_level(0) .NE. -1)then
387                  z=zuh
388               else
389                  z=zu
390               end if
391            else
392               if (isStrSubset( a(i),"zw_sp" ))then
393                  legend_label_z=legend_label_zw
394                  if (height_level(0) .NE. -1)then
395                     z=zwh
396                  else
397                     z=zw
398                  end if   
399               end if
400            end if
401         end do 
402
403         if (isStrSubset(vNam(varn),"x"))then
404            x_axis = f->k_x
405            x_axis = x_axis/normx
406            if (normx .NE. 1.)then
407               res@tiXAxisString = "k_x / "+normx
408            else
409               res@tiXAxisString = "k_x"
410            end if
411         else
412            x_axis = f->k_y
413            x_axis = x_axis/normx
414            if (normx .NE. 1.)then
415               res@tiXAxisString = "k_y / "+normx
416            else
417               res@tiXAxisString = "k_y"
418            end if
419         end if
420       
421         if (sort .EQ. "time")
422            do p=dimz-1,0,1
423               if (logy .EQ. 1)then 
424                  do q=0,dimt-1
425                     do r=0,dimsizes(x_axis)-1
426                        if (temp(q,p,r) .EQ. 0)then
427                           st=p+start_time_step
428                           print(" ")
429                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and height "+z(p)+" cannot be used")
430                           print(" ")
431                           res@trYLog = False
432                        end if
433                     end do
434                  end do
435               end if
436               res@trXMinF = min(x_axis)
437               res@trXMaxF = max(x_axis)
438               res@gsnLeftString      = vNam(varn)
439               res@gsnRightString     = "Height = "+z(p)+"m"
440               if (normy .NE. 1.)then
441                  res@tiYAxisString      = vNam(varn)+" / "+normy
442               else
443                  res@tiYAxisString      = vNam(varn)
444               end if
445               res@xyExplicitLegendLabels  = legend_label
446               plot(n)  = gsn_csm_xy(wks,x_axis,temp(:,p,:),res)
447               n=n+1
448            end do
449         else
450            if (sort .EQ. "height")
451               do p=0,dimt-1           
452                  do q=0,dimz-1
453                     do r=0,dimsizes(x_axis)-1
454                        if (temp(p,q,r) .EQ. 0)then
455                           st=p+start_time_step
456                           print(" ")
457                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and time "+legend_label(p)+" cannot be used")
458                           print(" ")
459                           res@trYLog = False
460                        end if
461                     end do
462                  end do
463                  res@trXMinF = min(x_axis)
464                  res@trXMaxF = max(x_axis)
465                  res@gsnLeftString      = vNam(varn)
466                  res@gsnRightString     = "Time = "+legend_label(p)+"h"
467                  if (normy .NE. 1.)then
468                     res@tiYAxisString      = vNam(varn)+" / "+normy
469                  else
470                     res@tiYAxisString      = vNam(varn)
471                  end if
472                  res@xyExplicitLegendLabels  = legend_label_z
473                  plot(n)  = gsn_csm_xy(wks,x_axis,temp(p,:,:),res)
474                  n=n+1
475               end do
476            end if
477         end if
478         delete(temp)
479         delete(x_axis)
480      end if
481   end do
482
483   if (n .EQ. 0) then
484      print(" ")
485      print("The variables 'var="+var+"' do not exist on your input file;")
486      print("be sure to have one comma berfore and after each variable")
487      print(" ")
488      exit
489   end if
490 
491   ; ***************************************************
492   ; merge plots onto one page
493   ; ***************************************************
494
495   resP                        = True
496   resP@txFont                 = "helvetica"
497   resP@txString               = f@title
498   resP@txFuncCode             = "~"
499   resP@txFontHeightF          = 0.014
500
501   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
502      gsn_panel(wks,plot,(/n,1/),resP)
503   else   
504      do i = 0,n-1, no_lines*no_columns
505         if( (i+no_lines*no_columns) .gt. (n-1)) then
506            gsn_panel(wks,plot(i:n-1),(/no_lines,no_columns/),resP)
507         else
508           gsn_panel(wks,plot(i:i+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
509         end if
510      end do
511   end if
512
513   print(" ")
514   print("Output to: " + file_out +"."+ format_out)
515   print(" ")
516   
517end
Note: See TracBrowser for help on using the repository browser.