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

Last change on this file since 405 was 324, checked in by weinreis, 15 years ago

modified spectra.ncl and profiles.ncl

File size: 19.6 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   else
405      color = 2
406   end if
407   if ( dash .eq. 0 ) then
408      res@xyMonoDashPattern = True
409   end if
410
411   wks=gsn_open_wks(format_out,file_out)
412   gsn_define_colormap(wks,"rainbow+white")
413
414   n=0
415   do varn =dim-1,0,1
416     
417      check = True
418
419      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
420            check = False
421      end if
422
423      if (var .NE. "all") then
424         check = isStrSubset( var,","+vNam(varn)+"," )
425      end if
426
427      if(check) then
428
429         temp = f[:]->$vNam(varn)$
430         data = temp(start_time_step:end_time_step,0:dimz-1,:)
431
432         temp_att = f_att->$vNam(varn)$
433         a=getvardims(temp_att)
434         b=dimsizes(a)
435
436         if (height_level(0) .NE. -1)then
437            do te=0,dimz-1
438            print(height_level(te))
439               data(:,te,:) = temp(start_time_step:end_time_step,height_level(te),:)           
440            end do
441         end if
442
443         data=data/(norm_y*norm_x)  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
444           
445         do i=0,b-1           
446            if (isStrSubset( a(i),"zu_sp" ))then
447               legend_label_z=legend_label_zu
448               if (height_level(0) .NE. -1)then
449                  z=zuh
450               else
451                  z=zu
452               end if
453            else
454               if (isStrSubset( a(i),"zw_sp" ))then
455                  legend_label_z=legend_label_zw
456                  if (height_level(0) .NE. -1)then
457                     z=zwh
458                  else
459                     z=zw
460                  end if   
461               end if
462            end if
463         end do 
464         
465         min_x=new(dimz,double)
466         max_x=new(dimz,double)
467         min_y=new(dimz,float)
468         max_y=new(dimz,float)
469         plot_h  = new(dimz,graphic)
470         
471         if (isStrSubset(vNam(varn),"x"))then
472            x_axis = new((/dimz,dimx/),double)
473            do q=0,dimz-1
474               x_axis(q,:) = f_att->k_x
475               x_axis = x_axis/norm_x
476            end do
477            if (norm_x .NE. 1.)then
478               res@tiXAxisString = "k~B~x~N~ ["+unit_x+"]"
479            else
480               if (norm_height .EQ. 1)then
481                  res@tiXAxisString = "k~B~x~N~ x z [1]"
482               else
483                  res@tiXAxisString = "k~B~x~N~ [1/m]"
484               end if
485            end if
486            dim_r=dimx
487         else
488            x_axis=new((/dimz,dimy/),double)
489            do q=0,dimz-1
490               x_axis(q,:) = f_att->k_y
491               x_axis = x_axis/norm_x
492            end do
493            if (norm_x .NE. 1.)then
494               res@tiXAxisString = "k~B~x~N~ ["+unit_x+"]"
495            else
496               if (norm_height .EQ. 1)then
497                  res@tiXAxisString = "k~B~x~N~ x z [1]"
498               else
499                  res@tiXAxisString = "k~B~x~N~ [1/m]"
500               end if
501            end if
502            dim_r=dimy
503         end if
504       
505         if (sort .EQ. "time")
506            res@xyLineColors = color
507            res@pmLegendDisplayMode     = "Always"
508            res@pmLegendSide            = "Top"
509            res@pmLegendParallelPosF    = 1.2
510            res@pmLegendOrthogonalPosF  = -1.0
511            res@pmLegendWidthF          = 0.12
512            res@pmLegendHeightF         = 0.04*(end_time_step-start_time_step+1)
513            do p=dimz-1,0,1
514               if (log_y .EQ. 1)then 
515                  do q=0,dimt-1
516                     do r=0,dim_r-1
517                        if (data(q,p,r) .EQ. 0)then
518                           st=p+start_time_step
519                           print(" ")
520                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and height "+z(p)+" cannot be used")
521                           print(" ")
522                           res@trYLog = False
523                        end if
524                     end do
525                  end do
526               end if
527               res@gsnLeftString      = vNam(varn)
528               res@gsnRightString     = "Height = "+z(p)+"m"               
529               res@tiYAxisString      = "["+unit_y+"]"           
530               res@xyExplicitLegendLabels  = legend_label             
531               if (norm_height .EQ. 1)then
532                  data(:,p,:)=data(:,p,:)*doubletofloat(z(p))
533                  x_axis(p,:) = x_axis(p,:)*z(p)
534               end if
535               res@trXMinF = min(x_axis(p,:))
536               res@trXMaxF = max(x_axis(p,:))
537               plot(n)  = gsn_csm_xy(wks,x_axis(p,:),data(:,p,:),res)
538               n=n+1
539            end do
540         else
541            if (sort .EQ. "height")
542               do p=0,dimt-1           
543                  do q=0,dimz-1
544                     do r=0,dim_r-1
545                        if (data(p,q,r) .EQ. 0)then
546                           st=p+start_time_step
547                           print(" ")
548                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and time "+legend_label(p)+" cannot be used")
549                           print(" ")
550                           res@trYLog = False
551                        end if
552                     end do
553                     if (norm_height .EQ. 1 .AND. p .EQ. 0)then
554                        data(p,q,:) = data(p,q,:)*doubletofloat(legend_label_z(q))
555                        x_axis(q,:) = x_axis(q,:)*doubletofloat(legend_label_z(q))
556                     end if
557                     max_y(q)=max(data(p,q,:))
558                     min_y(q)=min(data(p,q,:))
559                     min_x(q)=min(x_axis(q,:))
560                     max_x(q)=max(x_axis(q,:))                                   
561                  end do
562                  do q=0,dimz-1
563                     res@xyLineColor = color(q)
564                     if (dash .EQ. 1)then
565                        res@xyDashPattern = q
566                     end if
567                     if (q .EQ. 0)then
568                        res@tiYAxisString      = "["+unit_y+"]"                 
569                        res@gsnLeftString      = vNam(varn)
570                        res@gsnRightString     = "Time = "+legend_label(p)+"h"
571                        res@trXMinF = min(min_x)
572                        res@trXMaxF = max(max_x)
573                        res@trYMinF = min(min_y)
574                        res@trYMaxF = max(max_y)
575                       
576                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
577                       
578                        lgres = True
579                        if (dash .EQ. 0)then
580                           lgres@lgMonoDashIndex = True
581                        else
582                           lgres@lgDashIndexes   = ispan(0,dimz-1,1)
583                        end if
584                        if (black .EQ. 1)then
585                           lgres@lgMonoLineColors = True
586                        else
587                           lgres@lgLineColors = color
588                        end if
589                        lgres@lgTitleString = "Height [m]" 
590                        lgres@lgLabelFont        = "helvetica"
591                        lgres@lgLabelFontHeightF     = font_size_legend
592                        lgres@lgTitleFontHeightF     = font_size     
593                        lgres@vpWidthF           = 0.12           
594                        lgres@vpHeightF          = font_size_legend/5*dimz           
595 
596                        lbid = gsn_create_legend(wks,dimz,legend_label_z,lgres)       
597
598                        amres = True
599                        amres@amParallelPosF   = 0.75               
600                        amres@amOrthogonalPosF = 0.15           
601                        annoid1 = gsn_add_annotation(plot_h(q),lbid,amres)
602                     else
603                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
604                        overlay(plot_h(0),plot_h(q))
605                     end if
606                  end do             
607                  plot(n)=plot_h(0)
608                  n=n+1
609               end do
610            end if
611         end if
612         delete(data)
613         delete(temp)
614         delete(x_axis)
615         delete(min_x)
616         delete(max_x)
617         delete(min_y)
618         delete(max_y)
619         delete(plot_h)
620      end if
621   end do
622
623   if (n .EQ. 0) then
624      print(" ")
625      print("The variables 'var="+var+"' do not exist on your input file;")
626      print("be sure to have one comma berfore and after each variable")
627      print(" ")
628      exit
629   end if
630 
631   ; ***************************************************
632   ; merge plots onto one page
633   ; ***************************************************
634
635   resP                        = True
636   resP@txFont                 = "helvetica"
637   resP@txString               = f_att@title
638   resP@txFuncCode             = "~"
639   resP@txFontHeightF          = 0.015
640
641   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
642      gsn_panel(wks,plot,(/n,1/),resP)
643      print(" ")
644      print("Outputs to .eps or .epsi have only one frame")
645      print(" ")
646   else   
647      do i = 0,n-1, no_rows*no_columns
648         if( (i+no_rows*no_columns) .gt. (n-1)) then
649            gsn_panel(wks,plot(i:n-1),(/no_rows,no_columns/),resP)
650         else
651           gsn_panel(wks,plot(i:i+no_rows*no_columns-1),(/no_rows,no_columns/),resP)
652         end if
653      end do
654   end if
655
656   print(" ")
657   print("Output to: " + file_out +"."+ format_out)
658   print(" ")
659   
660end
Note: See TracBrowser for help on using the repository browser.