;
pro defwinddata,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
vector=map.vector
lats=*map.lats
lons=*map.lons
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

;
; 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] 
    mapdata = float(data[*,*,ixlev,ixtime])
  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

;
; Periodic point for geographic fields (map.data will now have 
;   field.nlon+1 longitudes)
;
lons = *field.lons
dlon = lons[1]-lons[0] ; assume regular grid in longitude
if lons[nlon-1]+dlon eq abs(lons[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[*]   <-idl thought the 2nd * was referring to a pointer
  lons[0:nlon-1] = *field.lons
  lons[nlon] = 180.
  nlon = nlon+1
  ;print,'defmapdata: added periodic point.'
endif else begin
  ;print,'defmapdata: did NOT need to set periodic point.'
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

dims=size(mapdata)
if map.projection eq 'Cylindrical Equidistant' then begin
  lonce=3
  latce=3
  for i=0,dims[1]-1 do begin
    if i mod lonce ne 0 then $
      mapdata[i,*]=map.missing_value
  endfor
  for j=0,dims[2]-1 do begin
    if j mod latce ne 0 then $
      mapdata[*,j]=map.missing_value
  endfor
endif else if map.projection eq 'Polar Stereographic' then begin
  lat1=lats(dims[2]-1)
  for i=0,dims[1]-1 do begin
    latst=1
    for j=0,dims[2]-1 do begin
      if abs(map.stperim) lt 10. then begin 
	lonst=1
        if abs(lats[j]) lt 20 then begin
	  lonst=1
        endif else if abs(lats[j]) lt 30 then begin
	  lonst=2
        endif else if abs(lats[j]) lt 40 then begin
	  lonst=2
        endif else if abs(lats[j]) lt 50 then begin
	  lonst=3
        endif else if abs(lats[j]) lt 60 then begin
	  lonst=3
        endif else if abs(lats[j]) lt 70 then begin
	  lonst=4
        endif else if abs(lats[j]) lt 80 then begin
	  lonst=4
        endif else if abs(lats[j]) le 90 then begin
	  lonst=15
        endif
      endif else if abs(map.stperim) lt 20. then begin
	lonst=1
        if abs(lats[j]) lt 20 then begin
	  lonst=1
        endif else if abs(lats[j]) lt 30 then begin
	  lonst=2
        endif else if abs(lats[j]) lt 40 then begin
	  lonst=2
        endif else if abs(lats[j]) lt 50 then begin
	  lonst=3
        endif else if abs(lats[j]) lt 60 then begin
	  lonst=3
        endif else if abs(lats[j]) lt 70 then begin
	  lonst=4
        endif else if abs(lats[j]) lt 80 then begin
	  lonst=4
        endif else if abs(lats[j]) le 90 then begin
	  lonst=15
        endif
      endif else if abs(map.stperim) lt 30. then begin
	lonst=1
        if abs(lats[j]) lt 30 then begin
	  lonst=2
        endif else if abs(lats[j]) lt 40 then begin
	  lonst=2
        endif else if abs(lats[j]) lt 50 then begin
	  lonst=3
        endif else if abs(lats[j]) lt 60 then begin
	  lonst=3
        endif else if abs(lats[j]) lt 70 then begin
	  lonst=4
        endif else if abs(lats[j]) lt 80 then begin
	  lonst=4
        endif else if abs(lats[j]) le 90 then begin
	  lonst=15
        endif
      endif else if abs(map.stperim) lt 40. then begin
	lonst=1
        if abs(lats[j]) lt 40 then begin
	  lonst=2
        endif else if abs(lats[j]) lt 50 then begin
	  lonst=3
        endif else if abs(lats[j]) lt 60 then begin
	  lonst=3
        endif else if abs(lats[j]) lt 70 then begin
	  lonst=4
        endif else if abs(lats[j]) lt 80 then begin
	  lonst=4
        endif else if abs(lats[j]) le 90 then begin
	  lonst=15
        endif
      endif else if abs(map.stperim) lt 50. then begin
	lonst=1
        if abs(lats[j]) lt 50 then begin
	  lonst=1
        endif else if abs(lats[j]) lt 60 then begin
	  lonst=2
        endif else if abs(lats[j]) lt 70 then begin
	  lonst=4
        endif else if abs(lats[j]) lt 80 then begin
	  lonst=4
        endif else if abs(lats[j]) le 90 then begin
	  lonst=15
        endif
      endif else if abs(map.stperim) lt 60. then begin
	lonst=1
        if abs(lats[j]) lt 60 then begin
	  lonst=2
        endif else if abs(lats[j]) lt 70 then begin
	  lonst=4
        endif else if abs(lats[j]) lt 80 then begin
	  lonst=4
        endif else if abs(lats[j]) le 90 then begin
	  lonst=15
        endif
      endif else if abs(map.stperim) lt 70. then begin
	lonst=1
        if abs(lats[j]) lt 70 then begin
	  lonst=4
        endif else if abs(lats[j]) lt 80 then begin
	  lonst=4
        endif else if abs(lats[j]) le 90 then begin
	  lonst=14
        endif
      endif else if abs(map.stperim) lt 80. then begin
	lonst=1
        if abs(lats[j]) lt 80 then begin
	  lonst=3
        endif else if abs(lats[j]) le 90 then begin
	  lonst=7
        endif
      endif else if abs(map.stperim) lt 90. then begin
	lonst=1
        if abs(lats[j]) le 90 then begin
	  lonst=3
        endif
      endif
      if i mod lonst ne 0 or j mod latst ne 0 then begin
        mapdata[i,j]=map.missing_value
      endif 
    endfor
  endfor
  if map.fixpolar eq 'SLT' then begin
    uthrs = float(map.mtime[1])+float(map.mtime[2])/60. ; ut (hrs)
    dlon=abs(lons[0]-lons[1])
    ishift = fix((uthrs-12.)*15./dlon)
    mapdata = shift(mapdata,ishift,0)
  endif

endif else if map.projection eq 'Satellite Projection' then begin
  lat1=lats(dims[2]-1)
  for i=0,dims[1]-1 do begin
    latsv=2
    for j=0,dims[2]-1 do begin
      lonsv=2
      if i mod lonsv ne 0 or j mod latsv ne 0 then begin
        mapdata[i,j]=map.missing_value
      endif 
    endfor
  endfor
endif else print,'>>Unknown projection:',map.projection

nlats=n_elements(*field.lats)
if nlats ne dims[2] then begin
  veclats=fltarr(dims[2])
  for j=0, dims[2]-1 do begin
    veclats[j]=lats[j]
  endfor
  *vector.latvecs=veclats
endif else *vector.latvecs=*field.lats
*vector.lonvecs=*field.lons

;
; Finally, update the map structure:
;
map.data = ptr_new(mapdata)
map.lons = ptr_new(lons)
map.lats = ptr_new(lats)
end