1 | load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
---|
2 | load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
---|
3 | load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
---|
4 | |
---|
5 | begin |
---|
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 | |
---|
541 | end |
---|