;
pro defmapdata,map,data,ixlev,ixtime
;
; Define 2d lon x lat array map.data from field data read from file.
; Set periodic points if necessary (geographic only)
; Also set map.lats and map.lons.
;
field = map.field
ifield=map.ifield
nlon = field.nlon & nlat = field.nlat
dims = size(data,/dimensions)
ndims = n_elements(dims)
if (dims[0] ne field.nlon) then begin
  print,'>>> WARNING defmapdata: 1st dimension is not nlons: dims[0]=',$
    dims[0],' but nlons=',field.nlon
endif
if (dims[1] ne field.nlat) then begin
  print,'>>> WARNING defmapdata: 2nd dimension is not nlats: dims[1]=',$
    dims[1],' but nlats=',field.nlat
endif
;
;  If plotting on altitude then use map data from getaltzdata routine. If ZP then get map data here.
;
IF map.plotz GE 1 THEN BEGIN

  mapdata = *map.data
 
ENDIF ELSE BEGIN
;
;  Set up variables used to convert WACCM data to 'standard' log pressure level
;
IF map.ftype EQ 'WACCM' THEN BEGIN

  levs = *field.levs
  nLevs = N_ELEMENTS(levs)
  newLevs = map.levs
  newNLevs = N_ELEMENTS(newLevs)

  IF ixlev GT 0 THEN BEGIN

    newLevSelected = newLevs[ixLev]
    iLev = IXFIND_NEAREST(levs, newlevSelected)
    levSelected = levs[iLev]
    IF iLev GT 0 THEN levSelectedLower = levs[iLev-1]
    IF iLev LT nLevs-1 THEN levSelectedUpper = levs[iLev+1]

    IF newLevSelected LT levSelected AND iLev GT 0 THEN BEGIN
  
      levDiffLower = levSelected - newLevSelected
      levDiffUpper = newLevSelected - levSelectedLower
      levDiffAll =  levSelected - levSelectedLower
      levRatLower = levDiffLower / levDiffAll
      levRatUpper = levDiffUpper / levDiffAll

    ENDIF ELSE IF newLevSelected GT levSelected AND iLev LT nLevs-1 THEN BEGIN
  
      levDiffLower = newLevSelected - levSelected
      levDiffUpper = levSelectedUpper - newLevSelected
      levDiffAll = levSelectedUpper - levSelected
      levRatLower = levDiffLower / levDiffAll
      levRatUpper = levDiffUpper / levDiffAll
   
    ENDIF

  ENDIF
  
  IF ixlev EQ 0 THEN iLev = 0
  IF ixlev EQ newNLevs-1 THEN iLev = nLevs-1
  
ENDIF
;
; Get data:
;
mapdata = fltarr(nlon,nlat)
case ndims of
  4: begin
    if dims[2] ne field.nlev then $
      print,'>>> WARNING defmapdata: 3rd dimension of 4d data ',$
        'should be nlev=',field.nlev,' but is ',dims[2] 
    if dims[3] ne field.ntime and field.ntime ne 1 then $
      print,'>>> WARNING defmapdata: 4th dimension of data ',$
        'should be ntime=',field.ntime,' but is ',dims[3] 

    IF map.ftype EQ 'WACCM' THEN BEGIN
    
      IF ixlev GT 0 THEN BEGIN
        
	IF iLev GT 0 AND iLev LT nLevs-1 THEN BEGIN 
	
	  IF isdensity(map.field.name) THEN BEGIN
	  
	    dataLower = ALOG10(data[*,*,ILev-1,ixtime])
	    dataLev = ALOG10(data[*,*,ILev,ixtime])
	    dataUpper = ALOG10(data[*,*,ILev+1,ixtime])
	    
	  ENDIF ELSE BEGIN
	  
	    dataLower = data[*,*,ILev-1,ixtime]
	    dataLev = data[*,*,ILev,ixtime]
	    dataUpper = data[*,*,ILev+1,ixtime]

          ENDELSE	  
    
          IF newLevSelected LT levSelected THEN mapdata = float(dataLower)*levRatUpper + $
                                                                       float(dataLev)*levRatLower
  
          IF newLevSelected GT levSelected THEN mapdata = float(dataLev)*levRatUpper + $
	                                                               float(dataUpper)*levRatLower

	  IF isdensity(map.field.name) THEN mapdata = 10^(mapdata)

; To handle top level        
        ENDIF ELSE BEGIN
	
	  mapdata = float(data[*,*,iLev,ixtime])
	  	
        ENDELSE	    

; To handle bottom level
      ENDIF ELSE BEGIN
	
	mapdata = float(data[*,*,iLev,ixtime])
	  	
      ENDELSE

    ENDIF ELSE BEGIN
     
      mapdata = float(data[*,*,ixlev,ixtime])

    ENDELSE	  
    
  end
  3: begin
    if dims[2] ne field.ntime then $
      print,'>>> WARNING defmapdata: 3rd dimension of 3d data ',$
        'should be ntime=',field.ntime,' but is ',dims[2] 
    mapdata = float(data[*,*,ixtime])
    
  end
  2: begin
    mapdata = float(data[*,*])
    
  end
  else: begin
    print,'>>> defmapdata: var data has ',ndims,' dimensions --',$
      ' cannot make map data.'
    mapdata = fltarr(nlon,nlat)
    mapdata = 0.
  end
endcase

ENDELSE ; data for altitude or zp level

lons = *field.lons
dlon = lons[1]-lons[0] ; assume regular grid in longitude
;
; Periodic point for geographic fields (map.data will now have 
;   field.nlon+1 longitudes)
;
; 5/8/08 btf:
; Periodic points cause problems in pltce, but are necessary
;   in pltst. So do not do it here, do it only in pltst.
;
set_periodic = 0
 if ((lons[nlon-1]+dlon eq abs(lons[0])) and set_periodic gt 0) then begin
   newdata = fltarr(nlon+1,nlat)
   newdata[0:nlon-1,*] = mapdata[*,*]
   newdata[nlon,*] = mapdata[0,*]
   mapdata = newdata
   lons = fltarr(nlon+1)
   lons[0:nlon-1] = *field.lons
   lons[nlon] = 180.
   nlon = nlon+1
 endif else begin
 endelse
;
; Change latitudes if this is for a polar projection:
;
lats = *field.lats
if map.projection eq 'Polar Stereographic' then begin
  lat0 = *field.lats
  nlat0 = n_elements(lat0)
  polarlats,map,lat,hem     ; see polarlats.pro
  nlat = n_elements(lat)
  poldata = fltarr(nlon,nlat)
  if hem eq 'S' then begin
    poldata[*,*] = mapdata[*,0:nlat-1]
    lats = fltarr(nlat)
    lats = lat0[0:nlat-1]
  endif else begin ; north
    lats = fltarr(nlat)
    for j=0,nlat-1 do begin
      poldata[*,j] = mapdata[*,nlat0-1-j]
      lats[j] = lat0[nlat0-1-j]
    endfor
  endelse
  mapdata = poldata
endif
;
; Take log10 of densities:
;
if (map.log10 eq 'density fields' and isdensity(field.name) gt 0) or $
   (map.log10 eq 'this frame only') or $
   (map.log10 eq 'all fields') then begin
  missing_value = float(map.missing_value)
  log10f,mapdata,missing_value
endif
;
; Finally, update the map structure:
;
map.data = ptr_new(mapdata)
map.lons = ptr_new(lons)
map.lats = ptr_new(lats)

end