;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("mag2latlon")
function mag2latlon(x:numeric)
; the gsn templates do not recognize the named dimensions
; of mlat,mlon,mlev and their cousins. This function therfore 
; changes the named dimensions of the data on magnetic coordinates. 
; This is only done to avoid NCL error messages, it does not change 
; the actual values of the coordinates. The function also assigns an 
; attribute "mag" to the variable so that the data can still be identified 
; as being on the magnetic grid for plotting purposes. 
local size
begin
  size = dimsizes(dimsizes(x))      ; four or three dimensions

; rename dimensions of x so ncl can recognize them
  if(size.eq.2)then
    x!0="lat"
    x!1="lon"
  end if
  if(size.eq.3)then
    x!1="lat"
    x!2="lon"
  end if
  if(size.eq.4)then
    x!1="lev"
    x!2="lat"
    x!3="lon"
  end if
; add units attribute 
  x&lat@units="degrees_north"
  x&lon@units="degrees_east"
  x@mag      = True
  return(x)
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
;================================================;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
; ================================================;
begin
;=================================================;
; open file and read in data
;=================================================;
  filename ="mag_timegcm_doy339_solF_qnm.nc" ; input file
  print(filename)
  itime    = 2  ; time index
  perimLat = 45.
;  
  f    = addfile(filename,"r")
  dBp_m  = f->dB_phi(itime,:,:) 
  dBl_m  = f->dB_lam(itime,:,:)    
  dBd_m  = f->dB_dwn(itime,:,:)     
  dBa_m  = f->B_abs_mag(itime,:,:)   
  kqp    = f->kqphi(itime,:,:)    
  kql    = f->kqlam(itime,:,:)  
  mtime  = f->mtime(itime,:)            
;=================================================;
; PLOT 2
;=================================================;  
  plot = new(5,graphic)                          ; create a plot array

  wks  = gsn_open_wks("ps","db_pol")    ; open a ncgm file
  gsn_define_colormap(wks,"gui_default")  ; choose a colormap

  res                     = True          ; plot mods desired  
  res@gsnDraw             = False         ; don't draw
  res@gsnFrame            = False         ; don't advance frame
  res@cnFillOn            = True          ; turn on color
  res@gsnSpreadColors     = True          ; use full range of colormap  
  res@gsnPolar            = "NH"          ; specify the hemisphere
; these resources must exist when plotting on magnetic coordinates
  res@mpOutlineOn         = False
  res@mpFillOn            = False
  res@gsnAddCyclic        = False
  res@mpMinLatF           = perimLat

  
  resv                     = True         ; plot mods desired  
  resv@gsnDraw             = False        ; don't draw
  resv@gsnFrame            = False        ; don't advance frame  
  resv@vcRefMagnitudeF	   = 0.1	  ; define vector ref mag
  resv@vcRefLengthF	   = 0.05	  ; define length of vec ref
  resv@vcMinDistanceF	   = 0.017
; these resources must exist when plotting on magnetic coordinates
  resv@mpOutlineOn         = False
  resv@mpFillOn            = False
  resv@gsnAddCyclic        = False
  resv@mpMinLatF           = perimLat


  x     = mag2latlon(dBp_m)
  plot(0) = gsn_csm_contour_map_polar(wks,x,res)  
  x     = mag2latlon(dBl_m)
  plot(1) = gsn_csm_contour_map_polar(wks,x,res)   
  x     = mag2latlon(dBd_m)  
  plot(2) = gsn_csm_contour_map_polar(wks,x,res)  
  x     = mag2latlon(kqp)      
  y     = mag2latlon(kql)      
  plot(3) = gsn_csm_vector_map_polar(wks, x, y, resv)   
  x     = mag2latlon(dBa_m)   
  plot(4) = gsn_csm_contour_map_polar(wks,x,res) 

;************************************************
  resP             = True                         ; modify the panel plot
  resP@gsnMaximize = True
  resP@txString    = filename + " UT= " + mtime(1)+":"+mtime(2)
  gsn_panel(wks,plot,(/3,2/),resP)               ; now draw as one plot
  
  end