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

Last change on this file since 221 was 218, checked in by letzel, 16 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
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
[190]5;***************************************************
6; load ncl_preferences.ncl
7;***************************************************
[194]8
9function which_script()
10local script
11begin
12   script="spectra"
13   return(script)
14end
[190]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
[176]26   end if
[190]27end if
[218]28
[190]29begin
[176]30
[190]31   ;***************************************************
32   ; set up default parameter values and strings
33   ;***************************************************
34 
35   if (file_1 .EQ. "File in") then
36      print(" ")
[218]37      print("Declare input file 'file_1=' in 'ncl_preferences.ncl' or prompt")
[190]38      print(" ")
39      exit
[176]40   else
41      file_in = file_1   
42   end if
43
[190]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   
[218]51   if (log_x .NE. 0 .AND. log_x .NE. 1)then
[190]52      print(" ")
[218]53      print("'log_x'= "+log_x+" is invalid and set to 1")
[190]54      print(" ")
[218]55      log_x = 1
[190]56   end if 
57   
[218]58   if (log_y .NE. 0 .AND. log_y .NE. 1)then
[190]59      print(" ")
[218]60      print("'log_y'= "+log_y+" is invalid and set to 1")
[190]61      print(" ")
[218]62      log_y = 1
[190]63   end if   
64 
[218]65   if (norm_y .EQ. 0.) then
[190]66      print(" ")
[218]67      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
[190]68      print(" ")
[218]69      norm_y = 1.0
[176]70   end if
[218]71   if (norm_x .EQ. 0.) then
[190]72      print(" ")
[218]73      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
[190]74      print(" ")
[218]75      norm_x= 1.0
[176]76   end if
[190]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" 
[176]83   end if
[190]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(" ")
[176]89      black = 0
90   end if
[190]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(" ")
[176]96      dash = 0
97   end if
98
[190]99   ;***************************************************
100   ; open input file
101   ;***************************************************
[218]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
[176]109
[218]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)
[176]124   
[218]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)
[176]173   print(" ")
[190]174   print("Variables in input file:")
175   print("- "+ vNam)
[176]176   print(" ")
177   dim = dimsizes(vNam)
[218]178   vDim = getfiledimsizes(f_att)
[176]179 
[218]180   t_all = f[:]->time
[176]181   nt    = dimsizes(t_all)
182   delta_t=t_all(nt-1)/nt
183   
[218]184   k_x=f_att->k_x
[190]185   dimx=dimsizes(k_x)
[218]186   k_y=f_att->k_y
[190]187   dimy=dimsizes(k_y)
188   
189   
190   dim_level=dimsizes(height_level)
191
[176]192   do i=0,dim-1
193      if (vNam(i) .EQ. "zu_sp")then
[218]194         zu=f_att->zu_sp         
[190]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 
[176]211      else
212         if (vNam(i) .EQ. "zw_sp")then
[218]213            zw=f_att->zw_sp
[190]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
[176]230         end if
231      end if
232   end do
[190]233
234   ;****************************************************       
[176]235   ; start of time step and different types of mistakes that could be done
[190]236   ;****************************************************
[176]237   
[190]238   if (start_time_step .EQ. -1.d) then         
[176]239      start_time_step=t_all(0)/3600
[190]240   else
[176]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(" ")
[218]245         print("Select another 'start_time_step'")
[176]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     
[190]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
[176]258         st=i
259         break
260      end if
261   end do
262   
[190]263   if (.not. isvar("st"))then
264      print(" ")
265      print("'start_time_step' = "+ start_time_step +"h is invalid")
266      print(" ")
[218]267      print("Select another 'start_time_step'")
[190]268      print(" ")
269      exit
270   end if
271   
272   ;****************************************************
[176]273   ; end of time step and different types of mistakes that could be done
[190]274   ;****************************************************
[176]275
[190]276   if (end_time_step .EQ. -1.d) then           
[176]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(" ")
[218]283         print("Select another 'end_time_step'") 
[176]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(" ")
[218]291         print("Select another 'start_time_step' or 'end_time_step'")
[176]292         print(" ")
293         exit
294      end if
295   end if
296
297   do i=0,nt-1     
[190]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
[176]299         et=i
300         break
301      end if
302   end do
[190]303   
304   if (.not. isvar("et"))then
305      print(" ")
306      print("'end_time_step' = "+ end_time_step +"h is invalid")
307      print(" ")
[218]308      print("Select another 'end_time_step'")
[190]309      print(" ")
310      exit
311   end if
[176]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 
[190]325   ;***************************************************
[176]326   ; set up recourses
[190]327   ;***************************************************
[176]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   
[218]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
[176]357      res@trXLog = True
[190]358   else
359      res@trXLog = False
360   end if
[218]361   if (log_y .EQ. 1)then
[176]362      res@trYLog = True
[190]363   else 
[178]364      res@trYLog = False
[176]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
[218]385     
[176]386      do p=0,dimz-1
[190]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
[176]394      end do
395   end if
[194]396
397   if (black .eq. 0 ) then
398      if (np .EQ. 1)then
[218]399         color = 237
[194]400      else   
401         step=round(235/(np-1),3)
[218]402         color = ispan(2,237,step)
[194]403      end if
[176]404   end if
405   if ( dash .eq. 0 ) then
[190]406      res@xyMonoDashPattern = True
[176]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
[178]414     
[176]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
[190]421      if (var .NE. "all") then
[176]422         check = isStrSubset( var,","+vNam(varn)+"," )
423      end if
424
[190]425      if(check) then
[218]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)
[190]432         b=dimsizes(a)
[176]433
[190]434         if (height_level(0) .NE. -1)then
435            do te=0,dimz-1
[218]436            print(height_level(te))
437               data(:,te,:) = temp(start_time_step:end_time_step,height_level(te),:)           
[190]438            end do
439         end if
440
[218]441         data=data/(norm_y*norm_x)  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
[190]442           
[176]443         do i=0,b-1           
444            if (isStrSubset( a(i),"zu_sp" ))then
445               legend_label_z=legend_label_zu
[190]446               if (height_level(0) .NE. -1)then
447                  z=zuh
448               else
449                  z=zu
450               end if
[176]451            else
452               if (isStrSubset( a(i),"zw_sp" ))then
[190]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   
[176]459               end if
460            end if
461         end do 
[218]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         
[176]469         if (isStrSubset(vNam(varn),"x"))then
[218]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+"]"
[190]477            else
[218]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
[190]483            end if
[218]484            dim_r=dimx
[176]485         else
[218]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+"]"
[190]493            else
[218]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
[190]499            end if
[218]500            dim_r=dimy
[176]501         end if
502       
503         if (sort .EQ. "time")
[218]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)
[190]511            do p=dimz-1,0,1
[218]512               if (log_y .EQ. 1)then 
[190]513                  do q=0,dimt-1
[218]514                     do r=0,dim_r-1
515                        if (data(q,p,r) .EQ. 0)then
[190]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
[176]523                  end do
[190]524               end if
[176]525               res@gsnLeftString      = vNam(varn)
[218]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)
[176]536               n=n+1
537            end do
538         else
[190]539            if (sort .EQ. "height")
540               do p=0,dimt-1           
[176]541                  do q=0,dimz-1
[218]542                     do r=0,dim_r-1
543                        if (data(p,q,r) .EQ. 0)then
[176]544                           st=p+start_time_step
545                           print(" ")
[190]546                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and time "+legend_label(p)+" cannot be used")
[176]547                           print(" ")
548                           res@trYLog = False
549                        end if
550                     end do
[218]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
[176]607               end do
608            end if
[190]609         end if
[218]610         delete(data)
611         delete(temp)
[178]612         delete(x_axis)
[218]613         delete(min_x)
614         delete(max_x)
615         delete(min_y)
616         delete(max_y)
617         delete(plot_h)
[176]618      end if
619   end do
620
621   if (n .EQ. 0) then
622      print(" ")
[190]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")
[176]625      print(" ")
626      exit
627   end if
628 
629   ; ***************************************************
630   ; merge plots onto one page
631   ; ***************************************************
632
633   resP                        = True
634   resP@txFont                 = "helvetica"
[218]635   resP@txString               = f_att@title
[176]636   resP@txFuncCode             = "~"
[218]637   resP@txFontHeightF          = 0.015
[176]638
639   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
[183]640      gsn_panel(wks,plot,(/n,1/),resP)
[176]641   else   
[218]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)
[176]645         else
[218]646           gsn_panel(wks,plot(i:i+no_rows*no_columns-1),(/no_rows,no_columns/),resP)
[176]647         end if
648      end do
649   end if
650
651   print(" ")
652   print("Output to: " + file_out +"."+ format_out)
653   print(" ")
654   
[190]655end
Note: See TracBrowser for help on using the repository browser.