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

Last change on this file since 188 was 183, checked in by letzel, 16 years ago
  • mrun: (optionally) link local input
  • diffusion_s: bugfix for calculation of horizontal fluxes at vertical walls
  • NCL scripts in trunk/SCRIPTS/NCL updated
File size: 17.4 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 
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
[183]15      if (isfilepresent("~/palm/current_version/trunk/SCRIPTS/NCL/.ncl_preferences")) then
16         parameter = asciiread("~/palm/current_version/trunk/SCRIPTS/NCL/.ncl_preferences",129,"string")
[176]17         delete(parameter@_FillValue)
18      else
19         print(" ")
[183]20         print("'.ncl_preferences' does not exist in '~/palm/current_version/trunk/SCRIPTS/NCL/'")
[176]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
[178]356   else
357      res@trXLog = False
358      res@trYLog = False
[176]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      np=dimt
375      res@lgTitleString = "Time [h]"
376   else
377      plot  = new(dim*dimt,graphic)
378      np=dimz
379      res@lgTitleString = "Height [m]" 
380      do p=0,dimz-1
381         legend_label_zu(p)=round(zu(p),3)
382         legend_label_zw(p)=round(zw(p),3)
383      end do
384   end if
385
386   if ( black .eq. 0 ) then 
387      res@xyLineColors = ispan(2,237,235/np)
388   end if
389   if ( dash .eq. 0 ) then
390      res@xyMonoDashPattern       = True
391   end if
392
393   wks=gsn_open_wks(format_out,file_out)
394   gsn_define_colormap(wks,"rainbow+white")
395
396   n=0
397   do varn =dim-1,0,1
[178]398     
[176]399      check = True
400
401      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
402            check = False
403      end if
404
405      if (.not. isvar("var")) then
406         if (parameter(21) .NE. "variables") then
407            var=parameter(21)
408            check = isStrSubset( var,","+vNam(varn)+"," )
409         end if
410      else         
411         check = isStrSubset( var,","+vNam(varn)+"," )
412      end if
413
414      if(check) then
415
416         temp = f->$vNam(varn)$(start_time_step:end_time_step,0:dimz-1,:)
417
418         temp=temp/norm  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
419   
420         a=getvardims(temp)
421         b=dimsizes(a)           
422         do i=0,b-1           
423            if (isStrSubset( a(i),"zu_sp" ))then
424               legend_label_z=legend_label_zu
425            else
426               if (isStrSubset( a(i),"zw_sp" ))then
427                  legend_label_z=legend_label_zw     
428               end if
429            end if
430         end do 
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                     end if
450                  end do
451               end do
452               res@trXMinF = min(x_axis)
453               res@trXMaxF = max(x_axis)
454               res@gsnLeftString      = vNam(varn)
[183]455               res@gsnRightString     = "Height = "+z(p)+"m"
[176]456               if (norm .NE. 1)then
457                  res@tiYAxisString      = vNam(varn)+" / "+norm
458               else
459                  res@tiYAxisString      = vNam(varn)
460               end if
461               res@xyExplicitLegendLabels  = legend_label
462               plot(n)  = gsn_csm_xy(wks,x_axis,temp(:,p,:),res)
463               n=n+1
464            end do
465         else
466            if (sort .EQ. "layer")
467               do p=dimt-1,0,1
468                  do q=0,dimz-1
469                     do r=0,dimsizes(x_axis)-1
470                        if (temp(p,q,r) .EQ. 0)then
471                           st=p+start_time_step
472                           print(" ")
473                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis cannot be used")
474                           print(" ")
475                           res@trYLog = False
476                        end if
477                     end do
478                  end do
479                  res@trXMinF = min(x_axis)
480                  res@trXMaxF = max(x_axis)
481                  res@gsnLeftString      = vNam(varn)
[183]482                  res@gsnRightString     = "Time = "+legend_label(p)+"h"
[176]483                  if (norm .NE. 1)then
484                     res@tiYAxisString      = vNam(varn)+" / "+norm
485                  else
486                     res@tiYAxisString      = vNam(varn)
487                  end if
488                  res@xyExplicitLegendLabels  = legend_label_z
489                  plot(n)  = gsn_csm_xy(wks,x_axis,temp(p,:,:),res)
490                  n=n+1
491               end do
492            else
493               print(" ")
494               print("Please choose 'sort' either equal 'time' or 'height'")
495               print(" ")
496               exit
497            end if
498         end if
499         delete(temp)
[178]500         delete(x_axis)
[176]501      end if
502   end do
503
504   if (n .EQ. 0) then
505      print(" ")
506      print("The variables 'var=°"+var+"°' do not exist on your input file")
507      print(" ")
508      exit
509   end if
510 
511   ; ***************************************************
512   ; merge plots onto one page
513   ; ***************************************************
514
515   resP                        = True
516   resP@txFont                 = "helvetica"
517   resP@txString               = f@title
518   resP@txFuncCode             = "~"
519   resP@txFontHeightF          = 0.014
520
521   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
[183]522      gsn_panel(wks,plot,(/n,1/),resP)
[176]523   else   
524      do i = 0,n-1, no_lines*no_columns
525         if( (i+no_lines*no_columns) .gt. (n-1)) then
[183]526            gsn_panel(wks,plot(i:n-1),(/no_lines,no_columns/),resP)
[176]527         else
[183]528           gsn_panel(wks,plot(i:i+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
[176]529         end if
530      end do
531   end if
532
533   print(" ")
534   print("Output to: " + file_out +"."+ format_out)
535   print(" ")
536   
537end
Note: See TracBrowser for help on using the repository browser.