;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("slt2lon")
function slt2lon(ut:numeric,lt:numeric)
local lon,ind_lt_neg180,ind_ge_pos180,flon
; takes a slt (lt) and ut and returns the 
; corresponding longitude. 
begin
  
  lon           = (lt-ut)*15
  ind_lt_neg180 = ind(lon.lt.-180)
  ind_ge_pos180 = ind(lon.ge.180)
  if(.not.any(ismissing(ind_lt_neg180))) then
    lon(ind_lt_neg180) = lon(ind_lt_neg180) + 360.
  end if
  if(.not.any(ismissing(ind_ge_pos180))) then
    lon(ind_ge_pos180) = lon(ind_ge_pos180) - 360.
  end if

  if(typeof(lon).eq."double")then
    flon = doubletofloat(lon)
    return(flon)
  else
    return(lon)
  end if
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("lon2slt")
function lon2slt(ut:numeric,lon:numeric)
local slt,ind_lt_00,ind_ge_24,fslt
; takes a lon and ut and returns the 
; corresponding solar local time. 
begin
  slt       = ut + lon/15
  ind_lt_00 = ind(slt.lt.0)
  ind_ge_24 = ind(slt.ge.24)
  if(.not.any(ismissing(ind_lt_00))) then
    slt(ind_lt_00) = slt(ind_lt_00) + 24
  end if
  if(.not.any(ismissing(ind_ge_24))) then
    slt(ind_ge_24) = slt(ind_ge_24) - 24
  end if

  if(typeof(slt).eq."double")then
    fslt = doubletofloat(slt)
    return(fslt)
  else
    return(slt)
  end if
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("mlon2mlt")
function mlon2mlt(ut:numeric,mlon:numeric)
local mlt,ind_lt_00,ind_ge_24,fmlt
; given a magnetic local time and ut, 
; calculates the corresponding magnetic
; longitude.
begin
; note on ce plots that the mlt at 71 should = ut in top text
  mlt       = ut + (mlon-71)/15
  ind_lt_00 = ind(mlt.lt.0)
  ind_ge_24 = ind(mlt.ge.24)
  if(.not.any(ismissing(ind_lt_00))) then
    mlt(ind_lt_00) = mlt(ind_lt_00) + 24
  end if
  if(.not.any(ismissing(ind_ge_24))) then
    mlt(ind_ge_24) = mlt(ind_ge_24) - 24
  end if

  if(typeof(mlt).eq."double")then
    fmlt = doubletofloat(mlt)
    return(fmlt)
  else
    return(mlt)
  end if
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("mlt2mlon")
function mlt2mlon(ut:numeric,lt:numeric)
local lon,ind_lt_neg180,ind_ge_pos180,flon
; takes a slt (lt) and ut and returns the 
; corresponding magnetic longitude. 
begin
  
  lon           = (lt-ut)*15 + 71
  ind_lt_neg180 = ind(lon.lt.-180)
  ind_ge_pos180 = ind(lon.ge.180)
  if(.not.any(ismissing(ind_lt_neg180))) then
    lon(ind_lt_neg180) = lon(ind_lt_neg180) + 360.
  end if
  if(.not.any(ismissing(ind_ge_pos180))) then
    lon(ind_ge_pos180) = lon(ind_ge_pos180) - 360.
  end if

  if(typeof(lon).eq."double")then
    flon = doubletofloat(lon)
    return(flon)
  else
    return(lon)
  end if
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("mlt_axis_lnoon")
function mlt_axis_lnoon(ut:numeric)
local lnoon,lon_vals,ind_lt_neg180,ind_ge_180
; given ut, calculates the longitude of magnetic
; local noon and then creates a set of longitude
; values around noon for use in the magnetic
; local time axis. 
begin
  lnoon                =  15*(12-ut)+71.
  lon_vals             =  (/lnoon-180,lnoon-120,lnoon-60,lnoon,lnoon+60,\
  lnoon+120,lnoon+180/)
  ind_lt_neg180             = ind(lon_vals.lt.-180)
  ind_ge_180                = ind(lon_vals.ge.180)
  if(.not.any(ismissing(ind_lt_neg180))) then
    lon_vals(ind_lt_neg180) = lon_vals(ind_lt_neg180) + 360
  end if
  if(.not.any(ismissing(ind_ge_180))) then
    lon_vals(ind_ge_180)   = lon_vals(ind_ge_180) - 360
  end if
  return(lon_vals)
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("set_lt_axis_utplot")
function set_lt_axis_utplot(ut[*]:numeric,values[*]:numeric,x:numeric,\
type:string,isMag:logical,res:logical)
local labs,ind_lt_00,ind_ge_24
; creates the local time axis on ut plots. 
begin
  res@tmXUseBottom     = False       ; keep top axis indp of bot
  res@tmXTOn           = True        ; tickmarks on top axis
  res@tmXTLabelsOn     = True        ; tickmark labs on top axis
  res@tmXTMode         = "Explicit"  ; tickmarks independent
  res@tmXTValues       =  values
  res@tiXAxisSide      = "top"        ; move x title to top
  res@tiXAxisFontHeightF = 0.017      ; xaxis font height
  res@tmXTLabelFontHeightF= 0.018
  res@gsnLeftString    = " "
  res@gsnRightString   = " "


  if(isMag.eq.True)then
    if(type.eq."lon")then
      res@tiXAxisString    = "Mag Local Time (hrs)" ; LT labels
      labs =  sprintf("%5.1f",mlon2mlt(ut,x))
    end if
    if(type.eq."slt")then
      res@tiXAxisString    = "Mag Longitude (deg)" ; LT labels
      labs = mlt2mlon(ut,x) 
    end if
  else
    if(type.eq."lon")then
      res@tiXAxisString    = "Local Time (hrs)" ; LT labels
      labs =  sprintf("%5.1f",lon2slt(ut,x))
    end if
    if(type.eq."slt")then
      res@tiXAxisString    = "Longitude (deg)" ; LT labels
      labs = slt2lon(ut,x) 
    end if
  end if

  res@tmXTLabels  = labs
  return(res)
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("set_lt_axis_cnplot")
function set_lt_axis_cnplot(ut:numeric,ismag:logical,res:logical)
; used for zpHeight vs. lat and zpHeight vs. lon templates
local lon_vals
begin
  res@tiXAxisFontHeightF = 0.017     ; xaxis font height
  res@gsnLeftString    = " "         ; turn off left string
  res@gsnRightString   = " "         ; trun off right string
  res@tiXAxisSide      = "top"       ; (lt axis on top for ce/polar/cont)
  res@tmXUseBottom     = False       ; keep top axis indpendent of bottom
  res@tmXTOn           = True        ; tickmarks on top axis ON
  res@tmXTLabelsOn     = True        ; tickmark labels on top axis ON
  res@tmXTMode         = "Explicit"  ; explicit labels on top axis
  if(isatt(res,"tmXTValues"))then
    delete(res@tmXTValues)
  end if
  if(isatt(res,"tmXTLabels"))then
    delete(res@tmXTLabels)
  end if
; MAGNETIC VARIBLES (local time axis on top)
  if(ismag)then
    lon_vals             =  mlt_axis_lnoon(ut)
    res@tiXAxisString    =  "Magnetic Local Time (hrs)" 
    res@tmXTValues       =  lon_vals
    res@tmXTLabels       =  sprintf("%5.1f",mlon2mlt(ut,lon_vals))
; not magnetic
  else
   
    lon_vals             =  (/-180.,-120.,-60.,0.,60.,120.,180./)
    res@tiXAxisString    =  "Local Time (hrs)" ; LT labels
    res@tmXTValues       =  lon_vals
    res@tmXTLabels       =  lon2slt(ut,lon_vals) 
  end if
  return(res)
end
;*********************************
undef("set_lt_axis_ceplot")
function set_lt_axis_ceplot(ut:numeric,ismag:logical,res:logical)
local lnoon,lon_vals,labs,ind_lt_180,ind_ge_180
; local time axis on ce plots
begin
  res@tiXAxisFontHeightF = 0.017     ; xaxis font height
  res@gsnLeftString    = " "         ; turn off left string
  res@gsnRightString   = " "         ; trun off right string
  
; MAGNETIC VARIBLES (lt axis on bottom)
  if(ismag)then
    lnoon                =  15*(12-ut)+71.
    lon_vals             =  (/lnoon-180,lnoon-120,lnoon-60,lnoon,lnoon+60,\
    lnoon+120,lnoon+180/)
    ind_lt_neg180             = ind(lon_vals.lt.-180)
    ind_ge_180                = ind(lon_vals.ge.180)
    if(.not.any(ismissing(ind_lt_neg180))) then
      lon_vals(ind_lt_neg180) = lon_vals(ind_lt_neg180) + 360
    end if
    if(.not.any(ismissing(ind_ge_180))) then
      lon_vals(ind_ge_180)   = lon_vals(ind_ge_180) - 360
    end if

    res@pmLabelBarOrthogonalPosF  = 0.13
    res@tiXAxisString        = "Magnetic Local Time (hrs)" 
    res@tmXBMode             = "Explicit"  ; explicit labels on top axis
    res@tmXBValues           =  lon_vals   ; location of slt/mlt labels
    res@tmXBLabels           =  sprintf("%5.2f",mlon2mlt(ut,lon_vals))
    if(isatt(res,"cnFillOn"))then
      res@tiXAxisOffsetYF    = .14     ; move xaxis title up
    end if
; not magnetic
  else
    lon_vals               =  (/-180.,-120.,-60.,0.,60.,120.,180./)
    res@tiXAxisString      =  "Local Time (hrs)" ; LT labels
    res@tiXAxisSide        = "top"       ; (lt axis on top for ce/polar/cont)
    res@tmXUseBottom       = False       ; keep top axis indpendent of bottom
    res@tmXTOn             = True        ; tickmarks on top axis ON
    res@tmXTLabelsOn       = True        ; tickmark labels on top axis ON
    res@tmXTMode           = "Explicit"  ; explicit labels on top axis
    res@tmXTValues         =  lon_vals   ; location of slt/mlt labels
    res@tmXTLabels         =  lon2slt(ut,lon_vals)
    res@tiXAxisOffsetYF    = .025        ; move xaxis title up
    res@tiMainOffsetYF     = .03         ; move main  title up
  end if
  return(res)
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("set_lt_axis_polarplot")
function set_lt_axis_polarplot(ut:numeric,ismag:logical,res:logical)
; local time axis on polar plots
begin
  res@tiXAxisFontHeightF = 0.017        ; xaxis font height
  res@tiXAxisSide        = "top"        ; move x title to top (can be only 1)
  res@tiXAxisOffsetYF    = .025         ; move xaxis title up
  res@tiMainOffsetYF     = .03          ; move main  title up
  res@gsnPolarUT         = ut           ; turn on local time
  res@gsnPolarTime       = True         ; turn on local time
  
; MAGNETIC VARIBLES 
  if(ismag)then
    res@tiXAxisString      =  "Magnetic Local Time (hrs)" 
  else
    res@tiXAxisString      =  "Local Time (hrs)" 
  end if
  return(res)
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("h_axis")
function h_axis(z:numeric,res:logical)
local tmp,zavg,zfloat,nlev,labels,lev_int,kk,yback,yfwd,kback,kfwd,minor_labels,min_lev,wb,wf,value,vint
begin 
; this function creates a YR height axis for the zpHeight vs. xxx plots

; we take either a meridonal,zonal, or gloabal average of z before the
; calculation. z is assumed to have already been subsampled in the
; call

  if(z@units .eq. "KM")then
    sc = 1.0                ; tgcmproc files are already interpolated
  else
    sc = 100000
  end if


  if(dimsizes(dimsizes(z)).eq.1)then ; ut or profile plot don't take dim_avg
    zavg         = z/sc
  else
;    zavg         = z(:,0)
    zavg         = dim_avg(z)/sc
  end if

  if(isdouble(zavg))then 
    zfloat     = doubletofloat(zavg)
    zfloat!0   ="lev"
    zfloat&lev = doubletofloat(z&lev)
  else
    zfloat     = zavg
  end if
  nlev  = dimsizes(zfloat)
;*********************************************
; MAJOR TICKMARK LABELS: interpolate to find zplev values 
;*********************************************
  labels  = (/20,30,40,50,60,70,80,90,100,150,200,300,400,500,600,700,800/)
  lev_int = new(dimsizes(labels),"float")
  do i = 0, dimsizes(labels)-1
    if(labels(i).lt.zfloat(0))then  ; label < lowest z value
      lev_int(i) = -999.
      continue
    end if
    if(labels(i).gt.zfloat(nlev-1))then  ; label > highest z value
      lev_int(i) = -999.
      continue
    end if

    kk = 0
    do k=0,dimsizes(zfloat)-2
      if(labels(i).ge.zfloat(k).and.labels(i).lt.zfloat(k+1))then 
        kk=k
        break
      end if
    end do
    if ( (labels(i) - zfloat(kk)) .eq. 0.0) then ; exact match
      lev_int(i) = zfloat&lev(kk)
    else
      yback = zfloat(kk)
      yfwd  = zfloat(kk+1)
      kfwd  = z&lev(kk+1)
      kback = z&lev(kk)
      if(ismissing(yfwd))then ; missing value in column
        lev_int=kback
      else
        wb    = fabs((yfwd-labels(i))/(yfwd-yback))
        wf    = fabs((yback-labels(i))/(yfwd-yback))
        value = (kback*wb+kfwd*wf)*1.
        if(typeof(value).eq."double")then
          vint = doubletofloat(value)
          lev_int(i) = vint
        else
          lev_int(i) = value
        end if
      end if
    end if
  end do
;************************************************
; MINOR TICKMARK LABELS
;************************************************
  minor_labels = (/25,35,45,55,65,75,85,95,110,120,130,140,160,170,180,\
  190,225,250,275,325,350,375,425,450,475,525,550,575,625,650,675,725,\
  750,775/)
  min_lev = new(dimsizes(minor_labels),"float")
  do i = 0, dimsizes(minor_labels)-1
    if(minor_labels(i).lt.zfloat(0))then  ; label < lowest z value
      min_lev(i) = -999.
      continue
    end if
    if(minor_labels(i).gt.zfloat(nlev-1))then  ; label > highest z value
      min_lev(i) = -999.
      continue
    end if

    kk = 0
    do k=0,dimsizes(zfloat)-2
      if(minor_labels(i).ge.zfloat(k).and.minor_labels(i).lt.zfloat(k+1))then 
        kk=k
        break
      end if
    end do
    if ( (minor_labels(i) - zfloat(kk)) .eq. 0.0) then ; exact match
      min_lev(i) = zfloat&lev(kk)
    else
      yback = zfloat(kk)
      yfwd  = zfloat(kk+1)
      kfwd  = z&lev(kk+1)
      kback = z&lev(kk)
      if(ismissing(yfwd))then ; missing value in column
        min_lev=kback
      else
        wb    = fabs((yfwd-minor_labels(i))/(yfwd-yback))
        wf    = fabs((yback-minor_labels(i))/(yfwd-yback))
        value = (kback*wb+kfwd*wf)*1.
        if(typeof(value).eq."double")then
          vint = doubletofloat(value)
          min_lev(i) = vint
        else
          min_lev(i) = value
        end if
      end if
    end if
  end do
  
  res@tmYUseLeft         = False           ; Keep left axis indep of right.
  res@tmYRLabelsOn       = True            ; have tick mark labels
  res@tmYROn             = True            ; have tick marks
  res@tmYRMode           = "Explicit"      ; label independently
  res@tmYRValues         = lev_int         ; location of labels
  res@tmYRLabels         = labels          ; labels for height axis
  res@tmLabelAutoStride  = True            ; nice tick mark labels
  res@tmYRMinorOn        = True            ; turn on minor tickmarks
  res@tmYRMinorValues    = min_lev         ; location of labels
  return(res)
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("remove_h_axis")
function remove_h_axis(res:logical)
begin 
; since the addition of a h_axis is done inside the main variable loop,
; it must be removed in the event of a magnetic variable. The only
; way to do this is to delete the actual resources
  
  if(isatt(res,"tmYRMode"))then  ; h_axis set from previous variable
    delete(res@tmYUseLeft)
    delete(res@tmYRLabelsOn)
    delete(res@tmYROn)
    delete(res@tmYRMode)
    delete(res@tmYRValues)
    delete(res@tmYRLabels)
    delete(res@tmLabelAutoStride)
    delete(res@tmYRMinorOn)
    delete(res@tmYRMinorValues)
  end if
  return(res)
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM