;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.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
undef("isMag")
function isMag(x:numeric)
local dims
; this function checks whether the data is in magnetic 
; coordinates
begin
  dims = (/"mlev","mlat","mlong","rmlat","rmlon","rmlev"/)
  if(any(isdim(x,dims)))then
    return(True)
  else
; variables that have been fixed by mag2latlon, no longer have
; the magnetic names, so we assigned the variable an attribute to
; check for.
    if(isatt(x,"mag"))then 
      return(True)
    end if
; a non-magnetic variable
    return(False)
  end if
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("time_interval")
function time_interval(in:file)
local mtime,time,ntime,days2min,hrs2min,mins,totmin_start,totmin_end,interval
; all something vs. ut plots can have the ut axis cover hours, days, or 
; months. We change the plot labels based on a subjective switch 
; threshold. The value returned from this function is used in the label
; of the plot as well as to determine the type of labels.
begin
  mtime = (/in->mtime/)
  time  = in->time
  ntime = dimsizes(time)
  
; beginning of time record
  days2min     = mtime(0,0)*24*60
  hrs2min      = mtime(0,1)*60
  mins         = mtime(0,2)
  totmin_start = days2min+hrs2min+mins
; end of time record
  days2min     = mtime(ntime-1,0)*24*60
  hrs2min      = mtime(ntime-1,1)*60
  mins         = mtime(ntime-1,2)
  
  totmin_end   = days2min+hrs2min+mins
  interval     = totmin_end - totmin_start   ; in minutes 
; 
  if(interval/60 .gt. 72.)then
    if(interval/(60*24) .gt. 90.)then
      return("months")
    else
      return("days")
    end if
  end if
  return("hours")
end 
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("zplev_to_height")
function zplev_to_height(z[*][*][*][*]:numeric,x[*][*][*][*]:numeric,\
arg:logical)
local lev,tmp,x_tmp,z_tmp,dims,x_int,zsc
; this function interpolates from zplev to height levels
; arg@height can be an array.
begin
 
; file processed by tgcmproc may already be interpolated to height
  if(x&lev@units.eq."km")then
    return(x)
  end if

; magnetic variables can not be interpolated to height

  if(isMag(x))then
    print("fields on mag grid can not be interpolated to height, select a zplev and try again")
    exit
  end if

; on difference files, Z is a difference, and can therefore 
; not be used in interpolating to height.
  if(isatt(x,"differences"))then
    print("diff files can not be interpolated to height")
    return(x)
  end if

  print("vertically interpolating  " + x@long_name + " to km") 

; user provides levels or calculated from top/bott/dz
  if(isatt(arg,"height"))then
    lev = arg@height                      ; maps
  else
    np = (arg@top-arg@bott)/arg@dz
    if(typeof(np).eq."float")then
      npts = floattointeger(np)
    else
      npts = np
    end if 
    tmp = fspan(arg@bott,arg@top,npts)  ; slice plots
    if(typeof(x).eq."double")then
      lev = new(dimsizes(tmp),"double")
      lev = fspan(arg@bott,arg@top,npts)
    else
      lev = fspan(arg@bott,arg@top,npts)
    end if
    delete(tmp)
  end if

; scale z if not already in km. If z has gone through halfzp, then it
; is scaled there.
  if(z@units .ne."KM")then
;    zsc = z/100000
;    zsc@units = "KM"
    sc = 100000
  else
    sc = 1.
  end if

  x_tmp = x(time|:,lat|:,lon|:,lev|:)
  z_tmp = z(time|:,lat|:,lon|:,lev|:)
  dims  = dimsizes(x_tmp)

  if(isatt(arg,"height"))then           ; map plots
    x_int = new( (/dims(0),dims(1),dims(2)/),typeof(x))
    tmp   = linint1(z_tmp/sc,x_tmp,False,lev,0)
    x_int = tmp(:,:,:,0)                ; dimension reduction
  else
    x_int = new( (/dims(0),dims(1),dims(2),dimsizes(lev)/),typeof(x))
    x_int = linint1(z(time|:,lat|:,lon|:,lev|:)/sc,x_tmp,False,lev,0)
    x_int!3             = "lev"
    x_int&lev           = lev
    x_int&lev@units     = "km"
    x_int&lev@long_name = "height (km)"
  end if
  
; copy over named dimensions and coordinate variables
  if(isatt(arg,"height"))then
    do n = 0, dimsizes(dimsizes(x_int))-1
      x_int!n=x_tmp!n
      if(iscoord(x_tmp,x_tmp!n))then
        x_int&$x_int!n$ = x_tmp&$x_tmp!n$
      end if 
    end do
  else
    do n = 0, dimsizes(dimsizes(x_int))-2
      x_int!n=x_tmp!n
      if(iscoord(x_tmp,x_tmp!n))then
        x_int&$x_int!n$ = x_tmp&$x_tmp!n$
      end if 
    end do
  end if

; copy over attributes
  att_names = getvaratts(x);
  if(.not.all(ismissing(att_names)))
    do i = 0, dimsizes(att_names)-1
      x_int@$att_names(i)$ = x@$att_names(i)$
    end do
  end if
  if(isatt(arg,"height"))then
    return(x_int(time|:,lat|:,lon|:))
  else
    return(x_int(time|:,lev|:,lat|:,lon|:))
  end if
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
undef("model_times2idx")
function model_times2idx(request[3]:integer,in:file)
local test2,dims,ntimes,j
; this function takes a model time given by the user in the
; form of 106,45,3 and determines which array index of the
; time dimension it equates to. This is index is then
; used to subscript the variable to be plotted.
begin
  mtime  = (/in->mtime/)   ; avoid bad coord varb
  dims   = dimsizes(mtime)
  ntimes = dims(0)

  do j=0,ntimes-1
    test2 = mtime(j,:)
    if all((request.eq.test2))then
      return(j)
      exit
    end if
  end do
  print("That model time does not exist on the file, defaulting to index 0")
  return(0)
end
;MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM