;*************************************************
; demonstrates the plotting of magnetic variables
;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "./gsn_csm.ncl"
load "$HAO_NCL_LIB/hao_contrib.ncl"
load "$HAO_NCL_LIB/halfzp2full.ncl"
;************************************************
begin
;************************************************
; before the data can be plotted, there is a 
; sequence of things that have to be done to it:
;************************************************
  in    = addfile("/chewy/d/murphys/files/smarch_migsdi_c_1.nc","r") 
  x     = mag2latlon(in->PHIM3D)

  t2plt = 2                     ; time index to plot
  l2plt = -5                    ; level value to plot
;************************************************
; create plot
;************************************************
  wks = gsn_open_wks("ps" ,"polar")           ; open a ps file
  gsn_define_colormap(wks,"BlWhRe") 

  res                          = True         ; plot mods desired
  res@gsnPolar                 = "NH"         ; choose hemisphere
  res@mpMinLatF                =  60          ; mpMaxLatF for NH
  res@cnFillOn                 = True         ; add color
  res@gsnSpreadColors          = True         ; use full range of color map
  res@lbLabelAutoStride        = True         ; nice lb labels

  res@mpGridAndLimbOn          = False        ; turn off lat/lon lines

  title = "Perim lat = -30"+" UT = "+in->ut(t2plt)+\
  " zplev = "+ in->lev({l2plt})+" (mag coords)"

  res@tiMainString             = title

; these resources must exist when plotting on magnetic coordinates
  res@mpOutlineOn              = False
  res@mpFillOn                 = False
  res@gsnAddCyclic             = False
  
  res@gsnPolarUT               = in->ut(t2plt); turn on local time
  res@gsnPolarTime             = True         ; turn on local time
  x&lon = x&lon-71                            ; so we can have mlt labels

  plot = gsn_csm_contour_map_polar(wks,x(t2plt,{l2plt},:,:), res)
end