source: palm/trunk/SCRIPTS/NCL/Archive/crosssections_new.ncl @ 2429

Last change on this file since 2429 was 154, checked in by letzel, 17 years ago
  • new NCL scripts added to trunk/SCRIPTS/NCL
  • previous NCL scripts moved to trunk/SCRIPTS/NCL/Archive
File size: 12.6 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
6; cross-sections
7; last change: $Id: crosssections_new.ncl 127 2007-10-23 11:05:25Z letzel $
8
9begin
10;
11; set default value(s) for shell script variables assigned on command line
12  if ( .not. isvar("cm") ) then     ; colormap
13    cm = "ncview_default"
14  end if
15  if ( .not. isvar("di") ) then     ; input directory (with final /)
16    di = ""
17  end if
18  if ( .not. isvar("d") ) then      ; output directory (with final /)
19    d = di
20  end if
21  if ( .not. isvar("fi") ) then     ; base name of input file (without suffix)
22    fi = "example_xy"
23  end if
24  if ( .not. isvar("fill_mode") ) then ; "AreaFill", "RasterFill" or "CellFill"
25    fill_mode = "AreaFill"
26  end if
27  if ( .not. isvar("fo") ) then     ; base name of output files (without suffix)
28    fo = ""
29  end if
30  if ( .not. isvar("mode") ) then   ; output mode ("Fill" or "Line")
31    mode = "Fill"
32  end if
33  if ( .not. isvar("t") ) then      ; output time step
34    t = 0
35  end if
36  if ( .not. isvar("var") ) then    ; variable to be output
37    var = "u_xy"
38  end if
39  if ( .not. isvar("xs") ) then     ; output x-coordinate range start (in m)
40    xs = -1e+38
41  end if
42  if ( .not. isvar("xe") ) then     ; output x-coordinate range end (in m)
43    xe = 1e+38
44  end if
45  if ( .not. isvar("ys") ) then     ; output y-coordinate range start (in m)
46    ys = -1e+38
47  end if
48  if ( .not. isvar("ye") ) then     ; output y-coordinate range end (in m)
49    ye = 1e+38
50  end if
51  if ( .not. isvar("zs") ) then     ; output z-coordinate range start (in m)
52    zs = -1e+38
53  end if
54  if ( .not. isvar("ze") ) then     ; output z-coordinate range end (in m)
55    ze = 1e+38
56  end if
57;
58; open input file
59  f  = addfile( di + fi  + ".nc", "r" )
60;
61; open workstation(s) and set colormap
62  wks_x11 = gsn_open_wks("x11","cross-section") ; X11 workstation
63  gsn_define_colormap(wks_x11,cm)
64  if ( isvar("fo") ) then
65    wks_pdf = gsn_open_wks("pdf",d+fo)  ; optional workstations
66    gsn_define_colormap(wks_pdf,cm)
67    wks_eps = gsn_open_wks("eps",d+fo)  ; for output on file
68    gsn_define_colormap(wks_eps,cm)
69    wks_ps  = gsn_open_wks("ps",d+fo)
70    gsn_define_colormap(wks_ps,cm)
71  end if
72;
73; read input data using 'coordinate subscripting'
74; NCL uses the closest corresponding values (in case of two equally distant
75; values NCL chooses the smaller one)
76  raw_data = f->$var$(t:t,{zs:ze},{ys:ye},{xs:xe})
77  raw_data!0 = "t"
78  raw_data!1 = "z"
79  raw_data!2 = "y"
80  raw_data!3 = "x"
81  time = raw_data&t
82;
83; reduce variable dimensions from 4D to 2D according to output ranges
84  if ( zs .eq. ze ) then
85    data = raw_data(0,0,:,:)
86    x_axis = "x"
87    y_axis = "y"
88    plane  = "z"
89
90    dx =  f->x(1)
91    dy =  f->y(1)
92
93    if ( raw_data&z .eq. -1 ) then
94      level = "-average"
95    else
96      level = "=" + raw_data&z + "m"
97    end if
98  else
99    if ( ys .eq. ye ) then
100      data = raw_data(0,:,0,:)
101      x_axis = "x"
102      y_axis = "z"
103      plane  = "y"
104
105      dx =  f->x(1)
106      if (isfilevar(f, "zw_3d")) then
107        dz =  f->zw_3d(1)
108      else
109        dz =  f->zw(1)
110      end if
111
112      if ( raw_data&y .eq. -1 ) then
113        level = "-average"
114      else
115        level = "=" + raw_data&y + "m"
116      end if
117    else
118      if ( xs .eq. xe ) then
119        data = raw_data(0,:,:,0)
120        x_axis = "y"
121        y_axis = "z"
122        plane  = "x"
123
124        dy =  f->y(1)
125        if (isfilevar(f, "zw_3d")) then
126          dz =  f->zw_3d(1)
127        else
128          dz =  f->zw(1)
129        end if
130
131        if ( raw_data&x .eq. -1 ) then
132          level = "-average"
133        else
134          level = "=" + raw_data&x + "m"
135        end if
136      end if
137    end if
138  end if
139  delete( raw_data )
140;
141; set up resources
142  cs_res                         = True
143  cs_res@gsnDraw                 = False 
144  cs_res@gsnFrame                = False
145  cs_res@gsnMaximize             = True
146  cs_res@gsnPaperOrientation     = "portrait"
147  cs_res@gsnPaperWidth           = 8.27
148  cs_res@gsnPaperHeight          = 11.69
149  cs_res@gsnPaperMargin          = 0.79
150  cs_res@tiMainFuncCode          = "~"
151  cs_res@tiMainFontHeightF       = 0.015
152  cs_res@tiMainString            = f@title
153  cs_res@tiXAxisString           = x_axis + " [m]"
154  cs_res@tiYAxisString           = y_axis + " [m]"
155;  cs_res@gsnLeftString           = ""   ; gsn_csm_* scripts use default data
156  cs_res@gsnCenterString         = "t=" + time + "s  " + plane + level
157;  cs_res@gsnRightString          = ""   ; gsn_csm_* scripts use default data
158  cs_res@tmXBMode                ="Automatic"
159  cs_res@tmYLMode                ="Automatic"
160  if ( mode .eq. "Fill" ) then
161    cs_res@cnFillOn                = True
162    cs_res@gsnSpreadColors         = True
163    cs_res@cnFillMode              = fill_mode
164    cs_res@lbOrientation           = "Vertical"  ; vertical label bar
165    cs_res@cnLinesOn               = False
166    cs_res@cnLineLabelsOn          = False
167  end if
168
169  vector     = True
170  co_overlay = True
171
172  if ( .not. isvar("wv1") ) then
173    vector = False
174  end if
175  if ( .not. isvar("wv2") ) then
176    vector = False
177  end if
178  if ( .not. isvar("co") ) then
179    co_overlay = False
180  end if
181
182  if ( vector .eq. True ) then
183
184    if ( wv1 .eq. "u_xy" .or.  wv1 .eq. "u_xz" .or. wv1 .eq. "u_yz" ) then
185
186      raw_data = f->$wv1$(t:t,{zs:ze},{ys:ye},{xs-(0.5*dx):xe-(0.5*dx)})
187      raw_data!0 = "t"
188      raw_data!1 = "z"
189      raw_data!2 = "y"
190      raw_data!3 = "x"
191      time = raw_data&t
192
193    else
194
195      if ( wv1 .eq. "v_xy" .or.  wv1 .eq. "v_xz" .or. wv1 .eq. "v_yz" ) then
196
197        raw_data = f->$wv1$(t:t,{zs:ze},{ys-(0.5*dy):ye-(0.5*dy)},{xs:xe})
198        raw_data!0 = "t"
199        raw_data!1 = "z"
200        raw_data!2 = "y"
201        raw_data!3 = "x"
202        time = raw_data&t
203
204      else
205
206        if ( wv1 .eq. "w_xy" .or.  wv1 .eq. "w_xz" .or. wv1 .eq. "w_yz" ) then
207
208          raw_data = f->$wv1$(t:t,{zs+(0.5*dz):ze+(0.5*dz)},{ys:ye},{xs:xe})
209          raw_data!0 = "t"
210          raw_data!1 = "z"
211          raw_data!2 = "y"
212          raw_data!3 = "x"
213          time = raw_data&t
214
215        end if
216      end if
217    end if
218
219    if ( zs .eq. ze ) then
220      data_u = raw_data(0,0,:,:)
221    else
222      if ( ys .eq. ye ) then
223        data_u = raw_data(0,:,0,:)
224      else
225        if ( xs .eq. xe ) then
226          data_u = raw_data(0,:,:,0)
227        end if
228      end if
229    end if
230    delete( raw_data )
231
232    if ( wv2 .eq. "u_xy" .or.  wv2 .eq. "u_xz" .or. wv2 .eq. "u_yz" ) then
233
234      raw_data = f->$wv2$(t:t,{zs:ze},{ys:ye},{xs-(0.5*dx):xe-(0.5*dx)})
235      raw_data!0 = "t"
236      raw_data!1 = "z"
237      raw_data!2 = "y"
238      raw_data!3 = "x"
239      time = raw_data&t
240
241    else
242
243      if ( wv2 .eq. "v_xy" .or.  wv2 .eq. "v_xz" .or. wv2 .eq. "v_yz" ) then
244
245        raw_data = f->$wv2$(t:t,{zs:ze},{ys-(0.5*dy):ye-(0.5*dy)},{xs:xe})
246        raw_data!0 = "t"
247        raw_data!1 = "z"
248        raw_data!2 = "y"
249        raw_data!3 = "x"
250        time = raw_data&t
251
252      else
253
254        if ( wv2 .eq. "w_xy" .or.  wv2 .eq. "w_xz" .or. wv2 .eq. "w_yz" ) then
255
256          raw_data = f->$wv2$(t:t,{zs+(0.5*dz):ze+(0.5*dz)},{ys:ye},{xs:xe})
257          raw_data!0 = "t"
258          raw_data!1 = "z"
259          raw_data!2 = "y"
260          raw_data!3 = "x"
261          time = raw_data&t
262
263        end if
264      end if
265    end if
266
267
268    if ( zs .eq. ze ) then
269      data_v = raw_data(0,0,:,:)
270    else
271      if ( ys .eq. ye ) then
272        data_v = raw_data(0,:,0,:)
273      else
274        if ( xs .eq. xe ) then
275          data_v = raw_data(0,:,:,0)
276        end if
277      end if
278    end if
279    delete( raw_data )
280
281    wv_res                     = True
282    wv_res@gsnDraw             = False 
283    wv_res@gsnFrame            = False
284    wv_res@gsnRightString      = ""
285    wv_res@gsnLeftString       = ""
286    wv_res@gsnCenterString     = ""
287    wv_res@vpClipOn            = True
288
289  end if
290
291  if ( co_overlay .eq. True ) then
292
293    raw_data = f->$co$(t:t,{zs:ze},{ys:ye},{xs:xe})
294    raw_data!0 = "t"
295    raw_data!1 = "z"
296    raw_data!2 = "y"
297    raw_data!3 = "x"
298    time = raw_data&t
299
300    if ( zs .eq. ze ) then
301      data2 = raw_data(0,0,:,:)
302    else
303      if ( ys .eq. ye ) then
304        data2 = raw_data(0,:,0,:)
305      else
306        if ( xs .eq. xe ) then
307          data2 = raw_data(0,:,:,0)
308        end if
309      end if
310    end if
311    delete( raw_data )
312
313    co_res                              = True
314    co_res@cnLineThicknessF             = 1.0
315    co_res@gsnDraw                      = False
316    co_res@gsnFrame                     = False
317    co_res@gsnContourZeroLineThicknessF = 2.0
318    co_res@gsnContourNegLineDashPattern = 4
319    co_res@cnLineColor                  = "Black"
320    co_res@gsnRightString   = ""
321    co_res@gsnLeftString    = ""
322    co_res@gsnCenterString  = ""
323
324;    co_res1@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
325;    co_res1@cnMinLevelValF       = -5.                ; set min contour level
326;    co_res1@cnMaxLevelValF       =  5.                ; set max contour level
327;    co_res1@cnLevelSpacingF      =  0.5               ; set contour spacing
328     co_res@cnLineLabelsOn        = True
329     co_res@cnLineLabelPlacementMode   = "constant"
330;    co_res@cnLineLabelBackgroundColor = "white"
331
332  end if
333
334  if ( co_overlay .eq. True .and. vector .eq. True ) then
335
336    if ( isvar("fo") ) then
337
338      plot_pdf2 = gsn_csm_contour(wks_pdf,data2,co_res)
339      plot_pdf1 = gsn_csm_contour(wks_pdf,data,cs_res)
340      plot_pdf3 = gsn_csm_vector(wks_pdf,data_u,data_v,wv_res)
341      overlay(plot_pdf1,plot_pdf2)
342      overlay(plot_pdf1,plot_pdf3)
343      draw(plot_pdf1)
344      frame(wks_pdf)
345
346      plot_eps2 = gsn_csm_contour(wks_eps,data2,co_res)
347      plot_eps1 = gsn_csm_contour(wks_eps,data,cs_res)
348      plot_eps3 = gsn_csm_vector(wks_eps,data_u,data_v,wv_res)
349      overlay(plot_eps1,plot_eps2)
350      overlay(plot_eps1,plot_eps3)
351      draw(plot_eps1)
352      frame(wks_eps)
353
354      plot_ps2 = gsn_csm_contour(wks_ps,data2,co_res)
355      plot_ps1 = gsn_csm_contour(wks_ps,data,cs_res)
356      plot_ps3 = gsn_csm_vector(wks_ps,data_u,data_v,wv_res)
357      overlay(plot_ps1,plot_ps2)
358      overlay(plot_ps1,plot_ps3)
359      draw(plot_ps1)
360      frame(wks_ps)
361
362    end if
363
364    plot1 = gsn_csm_contour(wks_x11,data,cs_res)
365    plot2 = gsn_csm_contour(wks_x11,data2,co_res)
366    plot3 = gsn_csm_vector(wks_x11,data_u,data_v,wv_res)
367
368    overlay(plot1,plot2)
369    overlay(plot1,plot3)
370
371    draw(plot1)
372    frame(wks_x11)
373
374  else
375
376    if ( co_overlay .eq. True .and. vector .ne. True ) then
377
378      if ( isvar("fo") ) then
379
380        plot_pdf2 = gsn_csm_contour(wks_pdf,data2,co_res)
381        plot_pdf1 = gsn_csm_contour(wks_pdf,data,cs_res)
382        overlay(plot_pdf1,plot_pdf2)
383        draw(plot_pdf1)       
384        frame(wks_pdf)
385
386        plot_eps2 = gsn_csm_contour(wks_eps,data2,co_res)
387        plot_eps1 = gsn_csm_contour(wks_eps,data,cs_res)
388        overlay(plot_eps1,plot_eps2)
389        draw(plot_eps1)
390        frame(wks_eps)
391
392        plot_ps2 = gsn_csm_contour(wks_ps,data2,co_res)
393        plot_ps1 = gsn_csm_contour(wks_ps,data,cs_res)
394        overlay(plot_ps1,plot_ps2)
395        draw(plot_ps1)
396        frame(wks_ps)
397
398      end if
399
400      plot1 = gsn_csm_contour(wks_x11,data,cs_res)
401      plot2 = gsn_csm_contour(wks_x11,data2,co_res)
402
403      overlay(plot1,plot2)
404
405      draw(plot1)
406      frame(wks_x11)
407
408    else
409
410    if ( co_overlay .ne. True .and. vector .eq. True ) then
411
412      if ( isvar("fo") ) then
413
414        plot_pdf1 = gsn_csm_contour(wks_pdf,data,cs_res)
415        plot_pdf3 = gsn_csm_vector(wks_pdf,data_u,data_v,wv_res)
416        overlay(plot_pdf1,plot_pdf3)
417        draw(plot_pdf1)
418        frame(wks_pdf)
419
420        plot_eps1 = gsn_csm_contour(wks_eps,data,cs_res)
421        plot_eps3 = gsn_csm_vector(wks_eps,data_u,data_v,wv_res)
422        overlay(plot_eps1,plot_eps3)
423        draw(plot_eps1)
424        frame(wks_eps)
425
426        plot_ps1 = gsn_csm_contour(wks_ps,data,cs_res)
427        plot_ps3 = gsn_csm_vector(wks_ps,data_u,data_v,wv_res)
428        overlay(plot_ps1,plot_ps3)
429        draw(plot_ps1)
430        frame(wks_ps)
431
432      end if
433
434      plot1 = gsn_csm_contour(wks_x11,data,cs_res)
435      plot3 = gsn_csm_vector(wks_x11,data_u,data_v,wv_res)
436
437      overlay(plot1,plot3)
438
439      draw(plot1)
440      frame(wks_x11)
441
442    else
443
444    if ( co_overlay .ne. True .and. vector .ne. True) then
445
446      if ( isvar("fo") ) then
447
448        plot_pdf1 = gsn_csm_contour(wks_pdf,data,cs_res)
449        draw(plot_pdf1)
450        frame(wks_pdf)
451
452        plot_eps1 = gsn_csm_contour(wks_eps,data,cs_res)
453        draw(plot_eps1)
454        frame(wks_eps)
455
456        plot_ps1 = gsn_csm_contour(wks_ps,data,cs_res)
457        draw(plot_ps1)
458        frame(wks_ps)
459
460      end if
461
462      plot1 = gsn_csm_contour(wks_x11,data,cs_res)
463
464      draw(plot1)
465      frame(wks_x11)
466
467    end if
468    end if
469    end if
470
471  end if
472
473end
Note: See TracBrowser for help on using the repository browser.