source: palm/tags/release-3.6/SCRIPTS/NCL/spectra.ncl

Last change on this file was 218, checked in by letzel, 14 years ago
  • User can check user parameters and deduce further quantities in user_check_parameters
  • Topography definition according to new user parameter topography_grid_convention (init_grid, modules, user_header, user_init_grid, user_parin)
  • New subdirectory trunk/EXAMPLES with two examples highlighting the topography_grid_convention
File size: 19.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("Declare input file 'file_1=' in 'ncl_preferences.ncl' or prompt")
38      print(" ")
39      exit
40   else
41      file_in = file_1   
42   end if
43
44   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
45      print(" ")
46      print("'format_out = "+format_out+"' is invalid and set to'x11'")
47      print(" ")
48      format_out="x11"
49   end if
50   
51   if (log_x .NE. 0 .AND. log_x .NE. 1)then
52      print(" ")
53      print("'log_x'= "+log_x+" is invalid and set to 1")
54      print(" ")
55      log_x = 1
56   end if 
57   
58   if (log_y .NE. 0 .AND. log_y .NE. 1)then
59      print(" ")
60      print("'log_y'= "+log_y+" is invalid and set to 1")
61      print(" ")
62      log_y = 1
63   end if   
64 
65   if (norm_y .EQ. 0.) then
66      print(" ")
67      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
68      print(" ")
69      norm_y = 1.0
70   end if
71   if (norm_x .EQ. 0.) then
72      print(" ")
73      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
74      print(" ")
75      norm_x= 1.0
76   end if
77   
78   if (sort .NE. "height" .AND. sort .NE. "time") then 
79      print(" ")
80      print("'sort'= "+sort+" is invalid and set to 'height'")
81      print(" ")
82      sort = "height" 
83   end if
84   
85   if (black .NE. 0 .AND. black .NE. 1)then
86      print(" ")
87      print("'black'= "+black+" is invalid and set to 0")
88      print(" ")
89      black = 0
90   end if
91 
92   if (dash .NE. 0 .AND. dash .NE. 1)then
93      print(" ")
94      print("'dash'= "+dash+" is invalid and set to 0")
95      print(" ")
96      dash = 0
97   end if
98
99   ;***************************************************
100   ; open input file
101   ;***************************************************
102   
103   file_in_1 = False
104   if (isStrSubset(file_in, ".nc"))then
105      start_f = -2
106      end_f = -2
107      file_in_1 = True     
108   end if
109
110   if (start_f .EQ. -1)then
111      print(" ")
112      print("'start_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
113      print(" ") 
114      exit
115   end if
116   if (end_f .EQ. -1)then
117      print(" ")
118      print("'end_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
119      print(" ") 
120      exit
121   end if
122
123   files=new(end_f-start_f+1,string)
124   
125   if (file_in_1)then
126      if (isfilepresent(file_in))then
127         files(0)=file_in
128      else
129         print(" ")
130         print("1st input file: '"+file_in+"' does not exist")
131         print(" ")
132         exit
133      end if
134   else   
135      if (start_f .EQ. 0)then
136         if (isfilepresent(file_in+".nc"))then
137            files(0)=file_in+".nc"
138            do i=1,end_f
139               if (isfilepresent(file_in+"."+i+".nc"))then   
140                  files(i)=file_in+"."+i+".nc"
141               else
142                  print(" ")
143                  print("Input file: '"+file_in+"."+i+".nc' does not exist")
144                  print(" ")
145                  exit 
146               end if       
147            end do         
148         else
149            print(" ")
150            print("Input file: '"+file_in+".nc' does not exist")
151            print(" ")
152            exit
153         end if
154      else
155         do i=start_f,end_f
156            if (isfilepresent(file_in+"."+i+".nc"))then   
157               files(i-start_f)=file_in+"."+i+".nc"
158            else
159               print(" ")
160               print("Input file: '"+file_in+"."+i+".nc' does not exist")
161               print(" ")
162               exit 
163            end if
164         end do
165      end if
166   end if
167
168   f=addfiles(files,"r")
169   f_att=addfile(files(0),"r")
170   ListSetType(f,"cat")
171   
172   vNam = getfilevarnames(f_att)
173   print(" ")
174   print("Variables in input file:")
175   print("- "+ vNam)
176   print(" ")
177   dim = dimsizes(vNam)
178   vDim = getfiledimsizes(f_att)
179 
180   t_all = f[:]->time
181   nt    = dimsizes(t_all)
182   delta_t=t_all(nt-1)/nt
183   
184   k_x=f_att->k_x
185   dimx=dimsizes(k_x)
186   k_y=f_att->k_y
187   dimy=dimsizes(k_y)
188   
189   
190   dim_level=dimsizes(height_level)
191
192   do i=0,dim-1
193      if (vNam(i) .EQ. "zu_sp")then
194         zu=f_att->zu_sp         
195         if (height_level(0) .EQ. -1)then
196            dimz=dimsizes(zu)
197         else
198            if (dim_level .GT. dimsizes(zu))then
199               print(" ")
200               print("'height_level' has more elements than available height levels in input file (= "+dimz+")")
201               print(" ")
202               exit
203            else
204               zuh=new(dim_level,double)
205               do le=0,dim_level-1
206                  zuh(le)=zu(height_level(le))
207               end do
208               dimz=dim_level
209            end if   
210         end if 
211      else
212         if (vNam(i) .EQ. "zw_sp")then
213            zw=f_att->zw_sp
214            if (height_level(0) .EQ. -1)then             
215               dimz=dimsizes(zw)
216            else
217               if (dim_level .GT. dimsizes(zw))then
218                  print(" ")
219                  print("'height_level' has more elements than available height levels in input file (= "+dimz+")")
220                  print(" ")
221                  exit
222               else
223                  zwh=new(dim_level,double)
224                  do le=0,dim_level-1
225                     zwh(le)=zw(height_level(le))
226                  end do
227                  dimz=dim_level
228               end if   
229            end if
230         end if
231      end if
232   end do
233
234   ;****************************************************       
235   ; start of time step and different types of mistakes that could be done
236   ;****************************************************
237   
238   if (start_time_step .EQ. -1.d) then         
239      start_time_step=t_all(0)/3600
240   else
241      if (start_time_step .GT. t_all(nt-1)/3600)then
242         print(" ")
243         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")
244         print(" ")
245         print("Select another 'start_time_step'")
246         print(" ")
247         exit
248      end if
249      if (start_time_step .LT. t_all(0)/3600)then
250         print(" ")
251         print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")       
252         exit
253      end if
254   end if
255
256   do i=0,nt-1     
257      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
258         st=i
259         break
260      end if
261   end do
262   
263   if (.not. isvar("st"))then
264      print(" ")
265      print("'start_time_step' = "+ start_time_step +"h is invalid")
266      print(" ")
267      print("Select another 'start_time_step'")
268      print(" ")
269      exit
270   end if
271   
272   ;****************************************************
273   ; end of time step and different types of mistakes that could be done
274   ;****************************************************
275
276   if (end_time_step .EQ. -1.d) then           
277      end_time_step = t_all(nt-1)/3600
278   else
279      if (end_time_step .GT. t_all(nt-1)/3600)then
280         print(" ")
281         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")
282         print(" ")
283         print("Select another 'end_time_step'") 
284         print(" ")
285         exit
286      end if
287      if (end_time_step .LT. start_time_step/3600)then
288         print(" ")
289         print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
290         print(" ")
291         print("Select another 'start_time_step' or 'end_time_step'")
292         print(" ")
293         exit
294      end if
295   end if
296
297   do i=0,nt-1     
298      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
299         et=i
300         break
301      end if
302   end do
303   
304   if (.not. isvar("et"))then
305      print(" ")
306      print("'end_time_step' = "+ end_time_step +"h is invalid")
307      print(" ")
308      print("Select another 'end_time_step'")
309      print(" ")
310      exit
311   end if
312
313   delete(start_time_step)
314   start_time_step=round(st,3)
315   delete(end_time_step)
316   end_time_step=round(et,3)
317
318   print(" ")
319   print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+t_all(start_time_step)+" s => index = "+start_time_step)
320   print("                     till "+t_all(end_time_step)/3600+" h = "+t_all(end_time_step)+" s => index = "+end_time_step)
321   print(" ")
322
323   dimt = end_time_step-start_time_step+1
324 
325   ;***************************************************
326   ; set up recourses
327   ;***************************************************
328     
329   res = True
330   res@gsnDraw                 = False
331   res@gsnFrame                = False
332   res@gsnPaperOrientation     = "portrait"
333   res@gsnPaperWidth           = 8.27
334   res@gsnPaperHeight          = 11.69
335   res@gsnPaperMargin          = 0.79
336   res@txFont                  = "helvetica"
337   res@tiMainFont              = "helvetica"
338   res@tiXAxisFont             = "helvetica"
339   res@tiYAxisFont             = "helvetica"
340   res@tmXBLabelFont           = "helvetica"
341   res@tmYLLabelFont           = "helvetica"
342   res@lgLabelFont             = "helvetica"
343   res@tmLabelAutoStride       = True
344   
345   res@lgLabelFontHeightF     = font_size_legend
346   res@lgTitleFontHeightF     = font_size
347   res@txFontHeightF      = font_size
348   res@tiXAxisFontHeightF = font_size
349   res@tiYAxisFontHeightF = font_size
350   res@tmXBLabelFontHeightF = font_size
351   res@tmYLLabelFontHeightF = font_size
352   
353   res@tmXBMinorPerMajor = 4
354   res@tmYLMinorPerMajor = 4
355   
356   if (log_x .EQ. 1) then
357      res@trXLog = True
358   else
359      res@trXLog = False
360   end if
361   if (log_y .EQ. 1)then
362      res@trYLog = True
363   else 
364      res@trYLog = False
365   end if
366
367   legend_label=new(dimt,double)
368   legend_label_zu=new(dimz,double)
369   legend_label_zw=new(dimz,double)
370   legend_label_z=new(dimz,double)
371   do p=start_time_step,end_time_step
372      if (t_all(p)/3600 .LT. 1) then
373         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,2,True)
374      else
375         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,0,True)
376      end if
377   end do 
378   if (sort .EQ. "time")
379      plot  = new(dim*dimz,graphic)
380      np=dimt
381      res@lgTitleString = "Time [h]"
382   else
383      plot  = new(dim*dimt,graphic)
384      np=dimz
385     
386      do p=0,dimz-1
387         if (height_level(0) .EQ. -1)then
388            legend_label_zu(p)=round(zu(p),3)
389            legend_label_zw(p)=round(zw(p),3)
390         else
391            legend_label_zu(p)=round(zuh(p),3)
392            legend_label_zw(p)=round(zwh(p),3)
393         end if
394      end do
395   end if
396
397   if (black .eq. 0 ) then
398      if (np .EQ. 1)then
399         color = 237
400      else   
401         step=round(235/(np-1),3)
402         color = ispan(2,237,step)
403      end if
404   end if
405   if ( dash .eq. 0 ) then
406      res@xyMonoDashPattern = True
407   end if
408
409   wks=gsn_open_wks(format_out,file_out)
410   gsn_define_colormap(wks,"rainbow+white")
411
412   n=0
413   do varn =dim-1,0,1
414     
415      check = True
416
417      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
418            check = False
419      end if
420
421      if (var .NE. "all") then
422         check = isStrSubset( var,","+vNam(varn)+"," )
423      end if
424
425      if(check) then
426
427         temp = f[:]->$vNam(varn)$
428         data = temp(start_time_step:end_time_step,0:dimz-1,:)
429
430         temp_att = f_att->$vNam(varn)$
431         a=getvardims(temp_att)
432         b=dimsizes(a)
433
434         if (height_level(0) .NE. -1)then
435            do te=0,dimz-1
436            print(height_level(te))
437               data(:,te,:) = temp(start_time_step:end_time_step,height_level(te),:)           
438            end do
439         end if
440
441         data=data/(norm_y*norm_x)  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
442           
443         do i=0,b-1           
444            if (isStrSubset( a(i),"zu_sp" ))then
445               legend_label_z=legend_label_zu
446               if (height_level(0) .NE. -1)then
447                  z=zuh
448               else
449                  z=zu
450               end if
451            else
452               if (isStrSubset( a(i),"zw_sp" ))then
453                  legend_label_z=legend_label_zw
454                  if (height_level(0) .NE. -1)then
455                     z=zwh
456                  else
457                     z=zw
458                  end if   
459               end if
460            end if
461         end do 
462         
463         min_x=new(dimz,double)
464         max_x=new(dimz,double)
465         min_y=new(dimz,float)
466         max_y=new(dimz,float)
467         plot_h  = new(dimz,graphic)
468         
469         if (isStrSubset(vNam(varn),"x"))then
470            x_axis = new((/dimz,dimx/),double)
471            do q=0,dimz-1
472               x_axis(q,:) = f_att->k_x
473               x_axis = x_axis/norm_x
474            end do
475            if (norm_x .NE. 1.)then
476               res@tiXAxisString = "k~B~x~N~ ["+unit_x+"]"
477            else
478               if (norm_height .EQ. 1)then
479                  res@tiXAxisString = "k~B~x~N~ x z [1]"
480               else
481                  res@tiXAxisString = "k~B~x~N~ [1/m]"
482               end if
483            end if
484            dim_r=dimx
485         else
486            x_axis=new((/dimz,dimy/),double)
487            do q=0,dimz-1
488               x_axis(q,:) = f_att->k_y
489               x_axis = x_axis/norm_x
490            end do
491            if (norm_x .NE. 1.)then
492               res@tiXAxisString = "k~B~x~N~ ["+unit_x+"]"
493            else
494               if (norm_height .EQ. 1)then
495                  res@tiXAxisString = "k~B~x~N~ x z [1]"
496               else
497                  res@tiXAxisString = "k~B~x~N~ [1/m]"
498               end if
499            end if
500            dim_r=dimy
501         end if
502       
503         if (sort .EQ. "time")
504            res@xyLineColors = color
505            res@pmLegendDisplayMode     = "Always"
506            res@pmLegendSide            = "Top"
507            res@pmLegendParallelPosF    = 1.2
508            res@pmLegendOrthogonalPosF  = -1.0
509            res@pmLegendWidthF          = 0.12
510            res@pmLegendHeightF         = 0.04*(end_time_step-start_time_step+1)
511            do p=dimz-1,0,1
512               if (log_y .EQ. 1)then 
513                  do q=0,dimt-1
514                     do r=0,dim_r-1
515                        if (data(q,p,r) .EQ. 0)then
516                           st=p+start_time_step
517                           print(" ")
518                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and height "+z(p)+" cannot be used")
519                           print(" ")
520                           res@trYLog = False
521                        end if
522                     end do
523                  end do
524               end if
525               res@gsnLeftString      = vNam(varn)
526               res@gsnRightString     = "Height = "+z(p)+"m"               
527               res@tiYAxisString      = "["+unit_y+"]"           
528               res@xyExplicitLegendLabels  = legend_label             
529               if (norm_height .EQ. 1)then
530                  data(:,p,:)=data(:,p,:)*doubletofloat(z(p))
531                  x_axis(p,:) = x_axis(p,:)*z(p)
532               end if
533               res@trXMinF = min(x_axis(p,:))
534               res@trXMaxF = max(x_axis(p,:))
535               plot(n)  = gsn_csm_xy(wks,x_axis(p,:),data(:,p,:),res)
536               n=n+1
537            end do
538         else
539            if (sort .EQ. "height")
540               do p=0,dimt-1           
541                  do q=0,dimz-1
542                     do r=0,dim_r-1
543                        if (data(p,q,r) .EQ. 0)then
544                           st=p+start_time_step
545                           print(" ")
546                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and time "+legend_label(p)+" cannot be used")
547                           print(" ")
548                           res@trYLog = False
549                        end if
550                     end do
551                     if (norm_height .EQ. 1 .AND. p .EQ. 0)then
552                        data(p,q,:) = data(p,q,:)*doubletofloat(legend_label_z(q))
553                        x_axis(q,:) = x_axis(q,:)*doubletofloat(legend_label_z(q))
554                     end if
555                     max_y(q)=max(data(p,q,:))
556                     min_y(q)=min(data(p,q,:))
557                     min_x(q)=min(x_axis(q,:))
558                     max_x(q)=max(x_axis(q,:))                                   
559                  end do
560                  do q=0,dimz-1
561                     res@xyLineColor = color(q)
562                     if (dash .EQ. 1)then
563                        res@xyDashPattern = q
564                     end if
565                     if (q .EQ. 0)then
566                        res@tiYAxisString      = "["+unit_y+"]"                 
567                        res@gsnLeftString      = vNam(varn)
568                        res@gsnRightString     = "Time = "+legend_label(p)+"h"
569                        res@trXMinF = min(min_x)
570                        res@trXMaxF = max(max_x)
571                        res@trYMinF = min(min_y)
572                        res@trYMaxF = max(max_y)
573                       
574                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
575                       
576                        lgres = True
577                        if (dash .EQ. 0)then
578                           lgres@lgMonoDashIndex = True
579                        else
580                           lgres@lgDashIndexes   = ispan(0,dimz-1,1)
581                        end if
582                        if (black .EQ. 1)then
583                           lgres@lgMonoLineColors = True
584                        else
585                           lgres@lgLineColors = color
586                        end if
587                        lgres@lgTitleString = "Height [m]" 
588                        lgres@lgLabelFont        = "helvetica"
589                        lgres@lgLabelFontHeightF     = font_size_legend
590                        lgres@lgTitleFontHeightF     = font_size     
591                        lgres@vpWidthF           = 0.12           
592                        lgres@vpHeightF          = font_size_legend/5*dimz           
593 
594                        lbid = gsn_create_legend(wks,dimz,legend_label_z,lgres)       
595
596                        amres = True
597                        amres@amParallelPosF   = 0.75               
598                        amres@amOrthogonalPosF = 0.15           
599                        annoid1 = gsn_add_annotation(plot_h(q),lbid,amres)
600                     else
601                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
602                        overlay(plot_h(0),plot_h(q))
603                     end if
604                  end do             
605                  plot(n)=plot_h(0)
606                  n=n+1
607               end do
608            end if
609         end if
610         delete(data)
611         delete(temp)
612         delete(x_axis)
613         delete(min_x)
614         delete(max_x)
615         delete(min_y)
616         delete(max_y)
617         delete(plot_h)
618      end if
619   end do
620
621   if (n .EQ. 0) then
622      print(" ")
623      print("The variables 'var="+var+"' do not exist on your input file;")
624      print("be sure to have one comma berfore and after each variable")
625      print(" ")
626      exit
627   end if
628 
629   ; ***************************************************
630   ; merge plots onto one page
631   ; ***************************************************
632
633   resP                        = True
634   resP@txFont                 = "helvetica"
635   resP@txString               = f_att@title
636   resP@txFuncCode             = "~"
637   resP@txFontHeightF          = 0.015
638
639   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
640      gsn_panel(wks,plot,(/n,1/),resP)
641   else   
642      do i = 0,n-1, no_rows*no_columns
643         if( (i+no_rows*no_columns) .gt. (n-1)) then
644            gsn_panel(wks,plot(i:n-1),(/no_rows,no_columns/),resP)
645         else
646           gsn_panel(wks,plot(i:i+no_rows*no_columns-1),(/no_rows,no_columns/),resP)
647         end if
648      end do
649   end if
650
651   print(" ")
652   print("Output to: " + file_out +"."+ format_out)
653   print(" ")
654   
655end
Note: See TracBrowser for help on using the repository browser.