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

load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

; cross-sections
; last change: $Id: crosssections.ncl 127 2007-10-23 11:05:25Z letzel $

begin
;
; set default value(s) for shell script variables assigned on command line
  if ( .not. isvar("cm") ) then     ; colormap
    cm = "ncview_default"
  end if
  if ( .not. isvar("di") ) then     ; input directory (with final /)
    di = ""
  end if
  if ( .not. isvar("d") ) then      ; output directory (with final /)
    d = di
  end if
  if ( .not. isvar("fi") ) then     ; base name of input file (without suffix)
    fi = "example_xy"
  end if
  if ( .not. isvar("fill_mode") ) then ; "AreaFill", "RasterFill" or "CellFill"
    fill_mode = "AreaFill"
  end if
  if ( .not. isvar("fo") ) then     ; base name of output files (without suffix)
    fo = ""
  end if
  if ( .not. isvar("mode") ) then   ; output mode ("Fill" or "Line")
    mode = "Fill"
  end if
  if ( .not. isvar("t") ) then      ; output time step
    t = 0
  end if
  if ( .not. isvar("var") ) then    ; variable to be output
    var = "u_xy"
  end if
  if ( .not. isvar("xs") ) then     ; output x-coordinate range start (in m)
    xs = -1e+38
  end if
  if ( .not. isvar("xe") ) then     ; output x-coordinate range end (in m)
    xe = 1e+38
  end if
  if ( .not. isvar("ys") ) then     ; output y-coordinate range start (in m)
    ys = -1e+38
  end if
  if ( .not. isvar("ye") ) then     ; output y-coordinate range end (in m)
    ye = 1e+38
  end if
  if ( .not. isvar("zs") ) then     ; output z-coordinate range start (in m)
    zs = -1e+38
  end if
  if ( .not. isvar("ze") ) then     ; output z-coordinate range end (in m)
    ze = 1e+38
  end if
;
; open input file
  f  = addfile( di + fi  + ".nc", "r" )
;
; open workstation(s) and set colormap
  wks_x11 = gsn_open_wks("x11","cross-section") ; X11 workstation
  gsn_define_colormap(wks_x11,cm)
  if ( isvar("fo") ) then
    wks_pdf = gsn_open_wks("pdf",d+fo)  ; optional workstations
    gsn_define_colormap(wks_pdf,cm)
    wks_eps = gsn_open_wks("eps",d+fo)  ; for output on file
    gsn_define_colormap(wks_eps,cm)
    wks_ps  = gsn_open_wks("ps",d+fo)
    gsn_define_colormap(wks_ps,cm)
  end if
;
; read input data using 'coordinate subscripting'
; NCL uses the closest corresponding values (in case of two equally distant 
; values NCL chooses the smaller one)
  raw_data = f->$var$(t:t,{zs:ze},{ys:ye},{xs:xe})
  raw_data!0 = "t"
  raw_data!1 = "z"
  raw_data!2 = "y"
  raw_data!3 = "x"
  time = raw_data&t
;
; reduce variable dimensions from 4D to 2D according to output ranges
  if ( zs .eq. ze ) then
    data = raw_data(0,0,:,:)
    x_axis = "x"
    y_axis = "y"
    plane  = "z"
    if ( raw_data&z .eq. -1 ) then
      level = "-average"
    else
      level = "=" + raw_data&z + "m"
    end if
  else
    if ( ys .eq. ye ) then
      data = raw_data(0,:,0,:)
      x_axis = "x"
      y_axis = "z"
      plane  = "y"
      if ( raw_data&y .eq. -1 ) then
        level = "-average"
      else
        level = "=" + raw_data&y + "m"
      end if
    else
      if ( xs .eq. xe ) then
        data = raw_data(0,:,:,0)
        x_axis = "y"
        y_axis = "z"
        plane  = "x"
        if ( raw_data&x .eq. -1 ) then
          level = "-average"
        else
          level = "=" + raw_data&x + "m"
        end if
      end if
    end if
  end if
  delete( raw_data )
;
; set up resources
  cs_res                         = True
  cs_res@gsnMaximize             = True 
  cs_res@gsnPaperOrientation     = "portrait"
  cs_res@gsnPaperWidth           = 8.27
  cs_res@gsnPaperHeight          = 11.69
  cs_res@gsnPaperMargin          = 0.79
  cs_res@tiMainFuncCode          = "~"
  cs_res@tiMainFontHeightF       = 0.015
  cs_res@tiMainString            = f@title
  cs_res@tiXAxisString           = x_axis + " [m]"
  cs_res@tiYAxisString           = y_axis + " [m]"
;  cs_res@gsnLeftString           = ""   ; gsn_csm_* scripts use default data
  cs_res@gsnCenterString         = "t=" + time + "s  " + plane + level
;  cs_res@gsnRightString          = ""   ; gsn_csm_* scripts use default data
  cs_res@tmXBMode                ="Automatic"
  cs_res@tmYLMode                ="Automatic"
  if ( mode .eq. "Fill" ) then
    cs_res@cnFillOn                = True
    cs_res@gsnSpreadColors         = True
    cs_res@cnFillMode              = fill_mode
    cs_res@lbOrientation           = "Vertical"  ; vertical label bar
    cs_res@cnLinesOn               = False
    cs_res@cnLineLabelsOn          = False
  end if
;
; data output
  plot_x11 = gsn_csm_contour(wks_x11,data,cs_res)
  if ( isvar("fo") ) then
    plot_pdf = gsn_csm_contour(wks_pdf,data,cs_res)
    plot_eps = gsn_csm_contour(wks_eps,data,cs_res)
    plot_ps  = gsn_csm_contour(wks_ps, data,cs_res)
  end if
end