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

Last change on this file since 178 was 178, checked in by steinfeld, 14 years ago

Updated version of the NCL-script for spectras

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