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

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