source: palm/trunk/SCRIPTS/NCL/Archive/crosssections.ncl @ 509

Last change on this file since 509 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: 4.8 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.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    if ( raw_data&z .eq. -1 ) then
90      level = "-average"
91    else
92      level = "=" + raw_data&z + "m"
93    end if
94  else
95    if ( ys .eq. ye ) then
96      data = raw_data(0,:,0,:)
97      x_axis = "x"
98      y_axis = "z"
99      plane  = "y"
100      if ( raw_data&y .eq. -1 ) then
101        level = "-average"
102      else
103        level = "=" + raw_data&y + "m"
104      end if
105    else
106      if ( xs .eq. xe ) then
107        data = raw_data(0,:,:,0)
108        x_axis = "y"
109        y_axis = "z"
110        plane  = "x"
111        if ( raw_data&x .eq. -1 ) then
112          level = "-average"
113        else
114          level = "=" + raw_data&x + "m"
115        end if
116      end if
117    end if
118  end if
119  delete( raw_data )
120;
121; set up resources
122  cs_res                         = True
123  cs_res@gsnMaximize             = True
124  cs_res@gsnPaperOrientation     = "portrait"
125  cs_res@gsnPaperWidth           = 8.27
126  cs_res@gsnPaperHeight          = 11.69
127  cs_res@gsnPaperMargin          = 0.79
128  cs_res@tiMainFuncCode          = "~"
129  cs_res@tiMainFontHeightF       = 0.015
130  cs_res@tiMainString            = f@title
131  cs_res@tiXAxisString           = x_axis + " [m]"
132  cs_res@tiYAxisString           = y_axis + " [m]"
133;  cs_res@gsnLeftString           = ""   ; gsn_csm_* scripts use default data
134  cs_res@gsnCenterString         = "t=" + time + "s  " + plane + level
135;  cs_res@gsnRightString          = ""   ; gsn_csm_* scripts use default data
136  cs_res@tmXBMode                ="Automatic"
137  cs_res@tmYLMode                ="Automatic"
138  if ( mode .eq. "Fill" ) then
139    cs_res@cnFillOn                = True
140    cs_res@gsnSpreadColors         = True
141    cs_res@cnFillMode              = fill_mode
142    cs_res@lbOrientation           = "Vertical"  ; vertical label bar
143    cs_res@cnLinesOn               = False
144    cs_res@cnLineLabelsOn          = False
145  end if
146;
147; data output
148  plot_x11 = gsn_csm_contour(wks_x11,data,cs_res)
149  if ( isvar("fo") ) then
150    plot_pdf = gsn_csm_contour(wks_pdf,data,cs_res)
151    plot_eps = gsn_csm_contour(wks_eps,data,cs_res)
152    plot_ps  = gsn_csm_contour(wks_ps, data,cs_res)
153  end if
154end
Note: See TracBrowser for help on using the repository browser.