source: palm/trunk/SCRIPTS/NCL/profiles.ncl @ 154

Last change on this file since 154 was 154, checked in by letzel, 16 years ago
  • new NCL scripts added to trunk/SCRIPTS/NCL
  • previous NCL scripts moved to trunk/SCRIPTS/NCL/Archive
File size: 16.9 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 
5begin
6
7   ; ***************************************************
8   ; read parameter_list
9   ; ***************************************************
10
11   if (isfilepresent("~/.ncl_preferences")) then
12      parameter = asciiread("~/.ncl_preferences",63,"string")
13      delete(parameter@_FillValue)
14   else
15      print(" ")
16      print("Please copy '.ncl_preferences' into your $home dircetory")
17      print(" ")
18      exit
19   end if
20 
21   ; ***************************************************
22   ; set up default parameter values and strings if not assigned in prompt or parameter list
23   ; ***************************************************
24 
25   if ( .not. isvar("file_in") ) then                   ; path+name of input file     
26      if (parameter(7) .EQ. "input file") then
27         print(" ")
28         print("Please provide input file 'file_in = ' either in prompt or parameter_list")
29         print(" ")
30         exit
31      else
32         file_in = parameter(7)
33      end if     
34   end if
35   if ( .not. isvar("format_out") ) then                ; format of output file
36      format_out = "x11"
37      if (parameter(9) .NE. "x11") then
38         format_out = parameter(9) 
39      end if
40   end if
41   if ( .not. isvar("file_out") ) then                  ; path+name of output file
42      file_out = "test"
43      if (parameter(11) .NE. "test_ts") then
44         file_out = parameter(11) 
45     end if     
46   end if
47   if ( .not. isvar("no_columns") ) then                ; number of plots in one row
48      no_columns = 1
49      if (parameter(17) .NE. "1") then
50         no_columns = stringtointeger(parameter(17)) 
51      end if
52   end if
53   if ( .not. isvar("no_lines") ) then                  ; number of plot-lines on one sheet
54      no_lines = 2
55      if (parameter(19) .NE. "2") then
56         no_lines = stringtointeger(parameter(19)) 
57      end if
58   end if
59   if ( .not. isvar("combine") ) then                   ; color of lines
60      combine = 0
61      if (parameter(23) .NE. "0") then
62         combine = stringtointeger(parameter(23))
63         if (stringtointeger(parameter(23)) .NE. 1) then
64            print(" ")
65            print("Please set 'combine' to 0 or 1")
66            print(" ")
67            exit
68         end if
69      end if
70   end if
71   if (combine .EQ. 1) then
72      if( .not. isvar("c_var") ) then
73         c_var=" "
74         if (parameter(27) .NE. "c_variables") then
75            c_var=parameter(27)
76         end if
77      end if
78   end if
79   if ( .not. isvar("black") ) then                     ; color of lines
80      black = 0
81      if (parameter(31) .NE. "0") then
82         black = stringtointeger(parameter(31))
83         if (stringtointeger(parameter(31)) .NE. 1) then
84            print(" ")
85            print("Please set 'black' to 0 or 1")
86            print(" ")
87            exit
88         end if
89      end if
90   end if
91   if ( .not. isvar("dash") ) then                      ; pattern of lines
92      dash = 0
93      if (parameter(29) .NE. "0") then
94         dash = stringtointeger(parameter(29))
95         if (stringtointeger(parameter(29)) .NE. 1) then
96            print(" ")
97            print("Please set 'dash' to 0 or 1")
98            print(" ")
99            exit
100         end if 
101      end if
102   end if
103
104   ; ***************************************************
105   ; open input file
106   ; ***************************************************
107
108   f=addfile( file_in,"r")
109   
110   vNam  = getfilevarnames(f)
111   print(" ")
112   print("Variable on netCDF file: " + vNam)
113   print(" ")
114   dim   = dimsizes(vNam)
115   z_pr  = f->zpt
116   dimz  = dimsizes(z_pr)
117   t_all = f->time
118   nt    = dimsizes(t_all)
119
120   ; ****************************************************       
121   ; start of time step and different types of mistakes that could be done
122   ; ****************************************************
123
124   if ( .not. isvar("start_time_step") ) then           
125      start_time_step = 1
126      if (parameter(13) .NE. "1") then
127         if (parameter(13) .LE. "0")
128            print(" ")
129            print("Begin with time step 1")
130            print(" ")
131            exit
132         end if
133         if (stringtointeger(parameter(13)) .GE. nt)
134            print(" ")
135            print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt-1))
136            print(" ")
137            exit
138         end if
139         start_time_step = stringtointeger(parameter(13)) 
140      end if
141   else
142      if (start_time_step .LE. 0)
143         print(" ")
144         print("Begin with time step 1")
145         print(" ")
146         exit
147      end if
148      if (start_time_step .GE. nt)
149         print(" ")
150         print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt-1))
151         print(" ")
152         exit
153      end if
154   end if
155
156   ; ****************************************************
157   ; end of time step and different types of mistakes that could be done
158   ; ****************************************************
159
160   if ( .not. isvar("end_time_step") ) then             
161      end_time_step = nt-1
162      if (parameter(15) .NE. "nt-1") then
163         if (parameter(15) .LE. "0")
164            print(" ")
165            print("'end_time_step' = "+parameter(15)+ " is too small; 'end_time_step' should be at least 1 ")
166            print(" ")
167            exit
168         end if
169         if (stringtointeger(parameter(15)) .GE. nt)
170            print(" ")
171            print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt-1))
172            print(" ")
173            exit
174         end if
175         if (stringtointeger(parameter(15)) .LT. stringtointeger(parameter(13)) )
176            print(" ")
177            print("'end_time_step' = "+ parameter(15) +" is lower than 'start_time_step' = "+parameter(13))
178            print(" ")
179            exit
180         end if
181         end_time_step = stringtointeger(parameter(15)) 
182      end if   
183   else
184      if (end_time_step .LE. 0)
185         print(" ")
186         print("'end_time_step' = "+end_time_step+ " is too small; 'end_time_step' should be at least 1 ")
187         print(" ")
188         exit
189      end if
190      if (end_time_step .GE. nt)
191         print(" ")
192         print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt-1))
193         print(" ")
194         exit
195      end if
196      if (end_time_step .LT. start_time_step)
197            print(" ")
198            print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
199            print(" ")
200            exit
201      end if
202   end if
203
204   ; ****************************************************
205   ; set up legend and colors
206   ; ****************************************************
207   
208   legend_label=new(nt-1,double)
209   do p=start_time_step,end_time_step
210      legend_label(p-start_time_step)=t_all(p)
211   end do
212   
213   np = end_time_step-start_time_step 
214   if ( black .eq. 0 ) then
215      res=True     
216      res@xyLineColors = ispan(2,237,235/np)
217   end if
218
219   ; ***************************************************
220   ; set up recourses
221   ; ***************************************************
222
223   res@gsnDraw                 = False
224   res@gsnFrame                = False
225   res@gsnPaperOrientation     = "portrait"
226   ;  res@gsnmaximize             = True
227   res@gsnPaperWidth           = 8.27
228   res@gsnPaperHeight          = 11.69
229   res@gsnPaperMargin          = 0.79
230   res@txFont                  = "helvetica"
231   res@tiMainFont              = "helvetica"
232   res@tiXAxisFont             = "helvetica"
233   res@tiYAxisFont             = "helvetica"
234   res@tmXBLabelFont           = "helvetica"
235   res@tmYLLabelFont           = "helvetica"
236   res@lgLabelFont             = "helvetica"
237   res@tmLabelAutoStride       = True
238   res@pmLegendDisplayMode     = "Always"
239   res@pmLegendSide            = "Top"
240   res@xyExplicitLegendLabels  = legend_label
241   res@pmLegendParallelPosF    = 1.4
242   res@pmLegendOrthogonalPosF  = -1.085
243   res@pmLegendWidthF          = 0.12
244   res@pmLegendHeightF         = 0.3
245   res@lgLabelFontHeightF     = .02
246   ;  res@XBLabelConstantSpacingF = 1.0
247
248   if ( dash .eq. 0 ) then
249      res@xyMonoDashPattern       = True 
250   end if
251
252   resP                        = True
253   resP@txFont                 = "helvetica"
254   resP@txString               = f@title
255   resP@txFuncCode             = "~"
256   resP@txFontHeightF          = 0.014
257
258   ; ***************************************************
259   ; set up colors and recourses for combined plot
260   ; ***************************************************
261   
262   if (combine .EQ. 1) then
263
264      if (.not. isvar("number_comb")) then
265         m=0
266         m=stringtointeger(parameter(25))
267         if(m .EQ. 0) then
268            print(" ")
269            print("Please indicate the number of variables you would like to combine ('number_comb')")
270            print(" ")
271            exit
272         end if
273      end if
274   
275      ores=True
276     
277      colors=new(m*(nt-1),integer)
278      do j=0,m-1
279         colors(j+j*(nt-2):(j+1)*(nt-1)-1) = ispan(2,237,235/np)
280      end do   
281      ores@xyLineColors = colors
282
283      dash_oplot=new(m*(nt-1),integer)
284      do j=0,m-1
285         dash_oplot(j+j*(nt-2):(j+1)*(nt-1)-1)=j
286      end do
287      ores@xyDashPatterns = dash_oplot
288
289      ores                         = True
290      ores@gsnDraw                 = False
291      ores@gsnFrame                = False
292      ores@gsnPaperOrientation     = "portrait"
293      ;  ores@gsnmaximize             = True
294      ores@gsnPaperWidth           = 8.27
295      ores@gsnPaperHeight          = 11.69
296      ores@gsnPaperMargin          = 0.79
297      ores@txFont                  = "helvetica"
298      ores@tiMainFont              = "helvetica"
299      ores@tiXAxisFont             = "helvetica"
300      ores@tiYAxisFont             = "helvetica"
301      ores@tmXBLabelFont           = "helvetica"
302      ores@tmYLLabelFont           = "helvetica"
303      ores@lgLabelFont             = "helvetica"
304      ores@tmLabelAutoStride       = True
305      ores@pmLegendDisplayMode     = "Always"
306      ores@pmLegendSide            = "Top"
307      ores@pmLegendParallelPosF    = 1.4
308      ores@pmLegendOrthogonalPosF  = -1.085
309      ores@pmLegendWidthF          = 0.15
310      ores@pmLegendHeightF         = 0.60
311      ores@lgLabelFontHeightF     = .02
312      ;  ores@XBLabelConstantSpacingF = 1.0
313 
314   end if
315
316   ; ***************************************************
317   ; set up graphics for plot
318   ; ***************************************************
319
320   plot  = new(dim,graphic)
321
322   wks=gsn_open_wks(format_out,file_out)
323   gsn_define_colormap(wks,"rainbow+white")
324 
325   ; ***************************************************
326   ; indicate plot number
327   ; ***************************************************
328
329   if (combine .EQ. 1) then
330      n = 1
331   else
332      n = 0
333   end if
334 
335   ; ***************************************************
336   ; set up minimum and maximum height
337   ; ***************************************************
338
339   if (.not. isvar("min_z"))
340      min_z=0
341      if (stringtointeger(parameter(33)) .NE. 0) then
342         if (stringtointeger(parameter(33)) .GE. max(z_pr) ) then
343            print(" ")
344            print("Minimum of height ('min_z'="+stringtointeger(parameter(33))+") is greater than available heights (="+max(z_pr)+")")
345            print(" ")
346            exit
347         end if
348         if (stringtointeger(parameter(33)) .LT. 0 ) then
349            print(" ")
350            print("Begin minimum of height 'min_z' with 0")
351            print(" ")
352            exit
353         end if
354         min_z=stringtointeger(parameter(33))
355      end if
356   else
357      if (min_z .GE. max(z_pr) ) then
358         print(" ")
359         print("Minimum of height ('min_z'="+min_z+") is greater than available heights (="+max(z_pr)+")")
360         print(" ")
361         exit
362      end if
363      if (min_z .LT. 0 ) then
364         print(" ")
365         print("Begin minimum of height 'min_z' with 0")
366         print(" ")
367         exit
368      end if
369   end if
370
371   if (.not. isvar("max_z"))
372      max_z=max(z_pr)
373      if ((parameter(35)) .NE. "max(z_pr)") then
374         if (stringtofloat(parameter(35)) .GE. max(z_pr) ) then
375            print(" ")
376            print("Maximum of height ('max_z'="+parameter(35)+") is greater than available heights (="+max(z_pr)+")")
377            print(" ")
378            exit
379         end if
380         if (stringtointeger(parameter(35)) .LE. 0 ) then
381            print(" ")
382            print("Maximum of height 'max_z' should be at least 1")
383            print(" ")
384            exit
385         end if
386         max_z=stringtointeger(parameter(35))
387      end if
388   else
389      if (max_z .GE. max(z_pr) ) then
390         print(" ")
391         print("Maximum of height ('max_z'="+max_z+") is greater than available heights (="+max(z_pr)+")")
392         print(" ")
393         exit
394      end if
395      if (max_z .LE. 0 ) then
396         print(" ")
397         print("Maximum of height 'max_z' should be at least 1")
398         print(" ")
399         exit
400      end if
401   end if
402
403   ; ***************************************************
404   ; read data and create plots
405   ; ***************************************************
406     
407   do ti = start_time_step, end_time_step
408      if( t_all(ti) .lt. 10^36) then
409         start_time_step = ti
410         break
411      end if
412   end do 
413   
414   if (combine .EQ. 1) then
415      data_o=new((/m,end_time_step,dimz/),float)
416      data_o!0 = "e"
417      legend_label_oplot=new(m*(nt-1),string)
418      mini=new(m,float)
419      maxi=new(m,float)
420      e=-1
421   end if 
422   
423   do varn = 0, dim-1   
424     
425      if ( isStrSubset( vNam(varn), "NORM") .or. isStrSubset( vNam(varn), "time") )
426         check = False
427      else
428         if (.not. isvar("var")) then
429            check = True
430            if (parameter(21) .NE. "variables") then
431               var=parameter(21)
432               check = isStrSubset( var,vNam(varn)+"," )
433            end if
434         else         
435            check = isStrSubset( var,vNam(varn)+"," )
436         end if
437      end if
438     
439      if (combine .EQ. 1) then
440         if (c_var .EQ. " ") then
441            print(" ")
442            print("Please indicate the variables you would like to combine ('c_var')")
443            print(" ")
444            exit
445         end if               
446         com=isStrSubset(c_var,vNam(varn)+",")     
447         if (com) then
448            e=e+1
449            check = False
450            temp = f->$vNam(varn)$           
451            data_o!1 = temp!0
452            data_o!2 = temp!1           
453            data_o(e,:,:) = temp(start_time_step:end_time_step,:)
454            mini(e) = min(data_o(e,:,:))
455            maxi(e) = max(data_o(e,:,:))
456           
457            do j=e,e
458               do p=start_time_step,end_time_step
459                  legend_label_oplot(e*end_time_step+p-start_time_step)=" "+flt2string(doubletofloat(t_all(p)))+" "+vNam(varn)
460               end do
461            end do
462           
463            delete(temp)
464            varn = varn + 1
465   
466         end if       
467      end if 
468 
469      if(check) then
470
471         temp = f->$vNam(varn)$     
472         data = temp(start_time_step:end_time_step,:)     
473         print("         plot of " + vNam(varn))
474         data!1 = "z"
475 
476         res@gsnLeftString      = "PROFILE plot of"
477         res@gsnRightString     = vNam(varn)
478         res@txFontHeightF      = 0.02
479         res@tiXAxisFontHeightF = 0.02
480         res@tiYAxisFontHeightF = 0.02
481         res@tiXAxisString      = " "   
482         res@tiYAxisString      = "Height [z]"
483         res@trXMinF            = min(data(:,:))
484         res@trXMaxF            = max(data(:,:))
485         res@trYMinF            = min_z
486         res@trYMaxF            = max_z   
487       
488         plot(n) = gsn_csm_xy(wks,data,z_pr,res) 
489     
490         delete(temp)
491         delete(data)
492         n = n + 1
493         varn = varn + 1
494       
495      end if     
496
497   end do
498
499   if (combine .EQ. 1) then
500      data_all=new((/m*end_time_step,dimz/),float)
501      do j=0,e
502         data_all(start_time_step-1+j*end_time_step:(j+1)*end_time_step-1,:) = data_o(j,:,:)
503      end do 
504
505      ores@gsnLeftString      = "combined PROFILE plot of"
506      ores@gsnRightString     = c_var
507      ores@txFontHeightF      = 0.02
508      ores@tiXAxisFontHeightF = 0.02
509      ores@tiYAxisFontHeightF = 0.02
510      ores@tiXAxisString      = " "   
511      ores@tiYAxisString      = "Height [z]"
512      ores@trXMinF            = min(mini(:))
513      ores@trXMaxF            = max(maxi(:))
514      ores@trYMinF            = min_z
515      ores@trYMaxF            = max_z   
516      ores@xyExplicitLegendLabels  = legend_label_oplot
517
518      print(" ")
519      print("combined plot of " + c_var)
520 
521      plot(0) = gsn_csm_xy(wks,data_all,z_pr,ores)
522     
523   end if
524
525   ; ***************************************************
526   ; merge plots onto one page
527   ; ***************************************************
528 
529   do i = 0, n-1, no_lines*no_columns
530      if( (i+no_lines*no_columns) .gt. (n-1)) then
531        gsn_panel(wks,plot(i:n-1),(/no_lines,no_columns/),resP)
532      else
533        gsn_panel(wks,plot(i:i+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
534      end if
535   end do
536
537   print(" ")
538   print("Output to: " + file_out +"."+ format_out)
539   print(" ")
540
541end
Note: See TracBrowser for help on using the repository browser.