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

  sigmas=new((/18/),float)
  sigmas=fspan(1.,18.,18)
  lev_p           = (/200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950,1000/)
  lev_p@units     = "hPa"                        ; required for vinth2p
  lev_p!0         = "lev_p"                  ; variable/dim name
  lev_p&lev_p     =  lev_p                   ; create coordinate variable
  lev_p@long_name = "pressure"               ; attach some attributes

 data = addfile("/ptmp/notaro/streamfn/ha.foam1-5lpj.0k0kfixveg.mone6351034.nc","r")
 lat = data->lat
 lon = data->lon
 lev = data->lev
 hyam = data->hyam
 hybm = data->hybm
 P0mb = data->P0*0.01
 ps =   data->PS(6,:,:)
 W =    data->OMEGA(6,:,:,:)
 U =    data->U(6,:,:,:)
 V =    data->V(6,:,:,:)

  vp   = vinth2p (V, hyam,hybm, lev_p ,ps, 2, P0mb, 1, False )
  psi = zonal_mpsi(vp,lat,lev_p*100.,ps)   ; calculate the zonal mean msf
  psi = psi * 1.e-11
  psi!0   = "lev"                             ; name the coordinates since
  psi!1   = "lat"                             ; zmmsf does not copy them.
  psi&lev = lev_p
  psi&lat = lat                               ; cp lat to "lat"

  wp   = vinth2p (W, hyam,hybm, lev_p ,ps, 2, P0mb, 1, False )
  vpzon=dim_avg(vp)
  wpzon=dim_avg(wp)
  vpzon!0="lev"
  vpzon!1="lat"
  vpzon&lev=lev_p
  vpzon&lat=lat
  copy_VarCoords(vpzon,wpzon)
  wpzon=wpzon*50.  ; just to amplify the magnitude of omega compared to v

;;;;;;;;;;;;;;;;;;;;;;;;;
  wks  = gsn_open_wks("ps","streamfunction")    ; open a ncgm file
  gsn_define_colormap(wks,"gui_default")  ; choose a colormap
  plot=new(1,graphic)

  res                     = True          ; plot mods desired
  res@gsnDraw=False
  res@gsnFrame=False
  res@cnFillOn            = False          ; turn on color
  res@cnLinesOn = True
  res@cnLineLabelsOn = True
  res@cnInfoLabelOn=False
  res@tiXAxisString="Latitude"
  res@tiYAxisString="Pressure"
  res@tmXTOn=False
  res@tiMainFontHeightF = 0.014
  res@tiMainString = "July Zonal Mean Merid Stream Fn & Merid/Vert Wind"
  res@lbOrientation        = "vertical"
  res@pmLabelBarOrthogonalPosF = 0.15
  res@vpWidthF            = 0.5           ; change aspect ratio of plot
  res@vpHeightF           = 0.45
  res@cnLineLabelPlacementMode = "computed"
  res@cnLineThicknessF=1.1
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels=(/-4.5,-4.2,-3.9,-3.6,-3.3,-3.0,-2.7,-2.4,-2.1,-1.8,-1.5,-1.2,-0.9,-0.6,-0.3,-0.1, \
                    .1,.3,.6,.9,1.2,1.5,1.8,2.1,2.4,2.7,3.0,3.3,3.6,3.9,4.2,4.5/)
  res@cnMonoLineDashPattern = False
  res@cnLineDashPatterns = (/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
  res@vcRefMagnitudeF = 0.5                ; define vector ref mag
  res@vcRefLengthF    = 0.01              ; define length of vec ref
  res@vcGlyphStyle    = "CurlyVector"      ; turn on curley vectors
  res@vcMinMagnitudeF = 0.2
  res@vcMinDistanceF  = 0.01               ; thin out vectors
  res@vcFillArrowsOn           = True  ; Fill the vector arrows
  res@vcMonoFillArrowFillColor = True ; in different colors
  res@vcFillArrowFillColor = "gray"
  res@vcFillArrowEdgeColor     = 1     ; Draw the edges in black.
  res@vcFillArrowWidthF        = 0.02 ; Make vectors thinner.
  res@vcGlyphStyle="FillArrow"
  res@vcRefAnnoOn=False
  res@vcFillArrowHeadInteriorXF=0.05
  plot(0)=gsn_csm_pres_hgt_vector(wks,psi,vpzon(:,0:38:2),wpzon(:,0:38:2),res)

  pres             = True
  pres@gsnDraw=True
  pres@gsnFrame=True
  pres@gsnMaximize = True
  pres@gsnPaperWidth=8.5
  pres@gsnPaperHeight=11.0
  pres@gsnPaperOrientation="portrait"
  pres@gsnPaperMargin=0.02
  pres@gsnPanelLeft = 0.02
  pres@gsnPanelRight = 0.98
  pres@gsnPanelBottom=0.03
  pres@gsnPanelTop=0.98
  pres@gsnPanelYWhiteSpacePercent = 2
  pres@gsnPanelRowSpec=True
  gsn_panel(wks,plot,(/1/),pres)

end