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

Last change on this file since 530 was 526, checked in by heinze, 15 years ago

Adjustment of the NCL scripts and palmplot to allow for the use of special characters in NetCDF variable names

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