source: palm/trunk/SCRIPTS/NCL/timeseries.ncl @ 157

Last change on this file since 157 was 157, checked in by letzel, 14 years ago
  • NCL scripts in trunk/SCRIPTS/NCL updated for vector plots
File size: 8.2 KB
Line 
1load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
3load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
4load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
5
6begin
7 
8   ; ***************************************************
9   ; read parameter_list
10   ; ***************************************************
11   
12   if (isfilepresent("~/.ncl_preferences")) then
13      parameter = asciiread("~/.ncl_preferences",73,"string")
14      delete(parameter@_FillValue)
15   else
16      print(" ")
17      print("Please copy '.ncl_preferences' into your $home dircetory")
18      print(" ")
19      exit
20   end if
21
22   ; ***************************************************
23   ; set up default parameter values and strings if not assigned in prompt or parameter list
24   ; ***************************************************
25   
26   if ( .not. isvar("file_in") ) then                   ; path+name of input file     
27      if (parameter(7) .EQ. "input file") then
28         print(" ")
29         print("Please provide input file 'file_in = ' either in prompt or parameter_list")
30         print(" ")
31         exit
32      else
33         file_in = parameter(7)
34      end if     
35   end if
36   if ( .not. isvar("format_out") ) then                ; format of output file
37      format_out = "x11"
38      if (parameter(9) .NE. "x11") then
39         format_out = parameter(9) 
40      end if
41   end if
42   if ( .not. isvar("file_out") ) then                  ; path+name of output file
43      file_out = "test"
44      if (parameter(11) .NE. "test") then
45         file_out = parameter(11) 
46      end if     
47   end if
48   if ( .not. isvar("no_columns") ) then                ; number of plots in one row
49      no_columns = 1
50      if (parameter(17) .NE. "1") then
51         no_columns = stringtointeger(parameter(17)) 
52      end if
53   end if
54   if ( .not. isvar("no_lines") ) then                  ; number of plot-lines on one sheet
55      no_lines = 2
56      if (parameter(19) .NE. "2") then
57         no_lines = stringtointeger(parameter(19)) 
58      end if
59   end if 
60   if ( .not. isvar("var") ) then                       ; variable name
61      check = True
62   end if
63 
64   ; ***************************************************
65   ; open input file
66   ; ***************************************************
67
68   f  = addfile(file_in , "r" )
69
70   vNam  = getfilevarnames(f)
71   print(" ")
72   print("Variable on netCDF file: " + vNam)
73   print(" ")
74   dim   = dimsizes(vNam)
75   t_all = f->time
76   nt    = dimsizes(t_all)
77
78   ; ****************************************************       
79   ; start of time step and different types of mistakes that could be done
80   ; ****************************************************
81
82   if ( .not. isvar("start_time_step") ) then           
83      start_time_step = 1
84      if (parameter(13) .LE. "1") then
85         if (parameter(13) .EQ. "0")
86            print(" ")
87            print("Begin with time step 1")
88            print(" ")
89            exit
90         end if
91         if (stringtointeger(parameter(13)) .GE. nt)
92            print(" ")
93            print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt-1))
94            print(" ")
95            exit
96         end if
97         start_time_step = stringtointeger(parameter(13)) 
98      end if
99   else
100      if (start_time_step .LE. 0)
101         print(" ")
102         print("Begin with time step 1")
103         print(" ")
104         exit
105      end if
106      if (start_time_step .GE. nt)
107         print(" ")
108         print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt-1))
109         print(" ")
110         exit
111      end if
112   end if
113
114   ; ****************************************************
115   ; end of time step and different types of mistakes that could be done
116   ; ****************************************************
117
118   if ( .not. isvar("end_time_step") ) then             
119      end_time_step = nt-1
120      if (parameter(15) .NE. "nt") then
121         if (parameter(15) .LE. "0")
122            print(" ")
123            print("'end_time_step' = "+parameter(15)+ " is too small; 'end_time_step' should be at least 1 ")
124            print(" ")
125            exit
126         end if
127         if (stringtointeger(parameter(15)) .GE. nt)
128            print(" ")
129            print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt-1)) 
130            print(" ")
131            exit
132         end if
133         if (stringtointeger(parameter(15)) .LT. stringtointeger(parameter(13)) )
134            print(" ")
135            print("'end_time_step' = "+ parameter(15) +" is lower than 'start_time_step' = "+parameter(13))
136            print(" ")
137            exit
138         end if
139         end_time_step = stringtointeger(parameter(15)) 
140      end if   
141   else
142      if (end_time_step .LE. 0)
143         print(" ")
144         print("'end_time_step' = "+end_time_step+ " is too small; 'end_time_step' should be at least 1 ")
145         print(" ")
146         exit
147      end if
148      if (end_time_step .GE. nt)
149         print(" ")
150         print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt-1))
151         print(" ")
152         exit
153      end if
154      if (end_time_step .LT. start_time_step)
155            print(" ")
156            print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
157            print(" ")
158            exit
159      end if
160   end if 
161
162   t = f->time(start_time_step:end_time_step)
163
164 
165   ; ***************************************************
166   ; set up recourses
167   ; ***************************************************
168
169   res                         = True
170   res@gsnDraw                 = False
171   res@gsnFrame                = False
172   res@gsnPaperOrientation     = "portrait"
173   res@gsnPaperWidth           = 8.27
174   res@gsnPaperHeight          = 11.69
175   res@gsnPaperMargin          = 0.79
176   res@tmXBMode                = True
177   res@tmYLMode                = True
178   res@txFont                  = "helvetica"
179   res@tiMainFont              = "helvetica"
180   res@tiXAxisFont             = "helvetica"
181   res@tiYAxisFont             = "helvetica"
182   res@tmXBLabelFont           = "helvetica"
183   res@tmYLLabelFont           = "helvetica" 
184   res@xyLineColors            = (/237/)
185   resP                        = True
186   resP@txFont                 = "helvetica"
187   resP@txString               = f@title+" time series "
188   resP@txFuncCode             = "~"
189   resP@txFontHeightF          = 0.015
190
191   res@vpWidthF=4
192
193
194   ; ***************************************************
195   ; read data and create plots
196   ; ***************************************************
197   
198   wks_ps = gsn_open_wks(format_out,file_out)           
199   gsn_define_colormap(wks_ps,"rainbow+white")
200   plot_ps=new(dim,graphic)                             
201   n=0
202 
203   do varn = 0, dim-1 
204      if( isvar("var") ) then                   
205         check = isStrSubset( var,vNam(varn)+"," )
206      end if
207      if (parameter(21) .NE. "variables") then
208         var = parameter(21)
209         check = isStrSubset( var,vNam(varn)+"," )
210      end if
211   
212      if( isStrSubset (vNam(varn), "time") )
213         check = False
214      end if
215 
216      if(check) then
217 
218         n=n+1
219         data = f ->$vNam(varn)$         
220         print("plot of " + vNam(varn))
221
222         res@gsnRightString       = vNam(varn)
223         res@tiXAxisString        = " time [s] "
224         res@tiXAxisFontHeightF   = 0.07
225         res@txFontHeightF        = 0.07
226         res@tiYAxisFontHeightF   = 0.07
227     
228         plot_ps(n) = gsn_csm_xy(wks_ps,t,data(start_time_step:end_time_step),res)
229   
230      end if
231 
232   end do
233
234   ; ***************************************************
235   ; merge plots onto one page
236   ; ***************************************************
237
238   do np = 1,n,no_lines*no_columns
239   
240      if ( np + no_lines*no_columns .gt. n) then   
241         gsn_panel(wks_ps, plot_ps(np:n),(/no_lines,no_columns/),resP)
242      else
243         gsn_panel(wks_ps, plot_ps(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
244      end if
245     
246   end do
247
248   print(" ")
249   print("Output to: " + file_out +"."+ format_out)
250   print(" ")
251 
252end
Note: See TracBrowser for help on using the repository browser.