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

Last change on this file since 176 was 176, checked in by steinfeld, 16 years ago

NCL scripts for spectra added.

File size: 17.8 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   end if
357
358   legend_label=new(dimt,double)
359   legend_label_zu=new(dimz,double)
360   legend_label_zw=new(dimz,double)
361   legend_label_z=new(dimz,double)
362   do p=start_time_step,end_time_step
363      if (t_all(p)/3600 .LT. 1) then
364         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,2,True)
365      else
366         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,0,True)
367      end if
368   end do 
369   if (sort .EQ. "time")
370      plot  = new(dim*dimz,graphic)
371      plot_ = new(dim*dimz,graphic)
372      np=dimt
373      res@lgTitleString = "Time [h]"
374   else
375      plot  = new(dim*dimt,graphic)
376      plot_ = new(dim*dimt,graphic)
377      np=dimz
378      res@lgTitleString = "Height [m]" 
379      do p=0,dimz-1
380         legend_label_zu(p)=round(zu(p),3)
381         legend_label_zw(p)=round(zw(p),3)
382      end do
383   end if
384
385   if ( black .eq. 0 ) then 
386      res@xyLineColors = ispan(2,237,235/np)
387   end if
388   if ( dash .eq. 0 ) then
389      res@xyMonoDashPattern       = True
390   end if
391
392   wks=gsn_open_wks(format_out,file_out)
393   gsn_define_colormap(wks,"rainbow+white")
394
395   n=0
396   do varn =dim-1,0,1
397
398      check = True
399
400      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
401            check = False
402      end if
403
404      if (.not. isvar("var")) then
405         if (parameter(21) .NE. "variables") then
406            var=parameter(21)
407            check = isStrSubset( var,","+vNam(varn)+"," )
408         end if
409      else         
410         check = isStrSubset( var,","+vNam(varn)+"," )
411      end if
412
413      if(check) then
414
415         temp = f->$vNam(varn)$(start_time_step:end_time_step,0:dimz-1,:)
416
417         temp=temp/norm  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
418   
419         a=getvardims(temp)
420         b=dimsizes(a)           
421         do i=0,b-1           
422            if (isStrSubset( a(i),"zu_sp" ))then
423               legend_label_z=legend_label_zu
424            else
425               if (isStrSubset( a(i),"zw_sp" ))then
426                  legend_label_z=legend_label_zw     
427               end if
428            end if
429         end do 
430
431         if (isStrSubset(vNam(varn),"x"))then
432            x_axis = f->k_x
433            res@tiXAxisString = "k_x"
434         else
435            x_axis = f->k_y
436            res@tiXAxisString = "k_y"
437         end if
438       
439         if (sort .EQ. "time")
440            do p=dimz-1,0,1 
441               do q=0,dimt-1
442                  do r=0,dimsizes(x_axis)-1
443                     if (temp(q,p,r) .EQ. 0)then
444                        st=p+start_time_step
445                        print(" ")
446                        print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis cannot be used")
447                        print(" ")
448                        res@trYLog = False
449                     else
450                        res@trYLog = True
451                     end if
452                  end do
453               end do
454               res@trYMinF = min(temp(:,p,:))
455               res@trYMaxF = max(temp(:,p,:))
456               res@trXMinF = min(x_axis)
457               res@trXMaxF = max(x_axis)
458               res@gsnLeftString      = vNam(varn)
459               res@gsnRightString     = "Height: "+z(p)+"m"
460               if (norm .NE. 1)then
461                  res@tiYAxisString      = vNam(varn)+" / "+norm
462               else
463                  res@tiYAxisString      = vNam(varn)
464               end if
465               res@xyExplicitLegendLabels  = legend_label
466               plot(n)  = gsn_csm_xy(wks,x_axis,temp(:,p,:),res)
467               n=n+1
468            end do
469         else
470            if (sort .EQ. "layer")
471               do p=dimt-1,0,1
472                  do q=0,dimz-1
473                     do r=0,dimsizes(x_axis)-1
474                        if (temp(p,q,r) .EQ. 0)then
475                           st=p+start_time_step
476                           print(" ")
477                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis cannot be used")
478                           print(" ")
479                           res@trYLog = False
480                        else
481                           res@trYLog = True
482                        end if
483                     end do
484                  end do
485                  res@trYMinF = min(temp(p,:,:))
486                  res@trYMaxF = max(temp(p,:,:))
487                  res@trXMinF = min(x_axis)
488                  res@trXMaxF = max(x_axis)
489                  res@gsnLeftString      = vNam(varn)
490                  res@gsnRightString     = "Time: "+legend_label(p)+"h"
491                  if (norm .NE. 1)then
492                     res@tiYAxisString      = vNam(varn)+" / "+norm
493                  else
494                     res@tiYAxisString      = vNam(varn)
495                  end if
496                  res@xyExplicitLegendLabels  = legend_label_z
497                  plot(n)  = gsn_csm_xy(wks,x_axis,temp(p,:,:),res)
498                  n=n+1
499               end do
500            else
501               print(" ")
502               print("Please choose 'sort' either equal 'time' or 'height'")
503               print(" ")
504               exit
505            end if
506         end if
507         delete(temp)
508      end if
509   end do
510
511   if (n .EQ. 0) then
512      print(" ")
513      print("The variables 'var=°"+var+"°' do not exist on your input file")
514      print(" ")
515      exit
516   end if
517 
518   ; ***************************************************
519   ; merge plots onto one page
520   ; ***************************************************
521
522   resP                        = True
523   resP@txFont                 = "helvetica"
524   resP@txString               = f@title
525   resP@txFuncCode             = "~"
526   resP@txFontHeightF          = 0.014
527
528   do m=0,n-1
529      plot_(m)=plot(n-1-m)
530   end do
531
532   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
533      gsn_panel(wks,plot_,(/n,1/),resP)
534   else   
535      do i = 0,n-1, no_lines*no_columns
536         if( (i+no_lines*no_columns) .gt. (n-1)) then
537            gsn_panel(wks,plot_(i:n-1),(/no_lines,no_columns/),resP)
538         else
539           gsn_panel(wks,plot_(i:i+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
540         end if
541      end do
542   end if
543
544   print(" ")
545   print("Output to: " + file_out +"."+ format_out)
546   print(" ")
547   
548end
Note: See TracBrowser for help on using the repository browser.