;
pro defmeriddata,line,data,linedata

  field = line.field
  fieldlevs = *field.levs
  nfieldlevs = n_elements(fieldlevs) 

  lats = *line.lats
  nlat = N_ELEMENTS(lats)
  lons = *line.lons
  nlon = N_ELEMENTS(lons)

  iLat = line.ilat
  ilon = line.ilon
  iLev = line.ilev

  IF line.ftype EQ 'WACCM' THEN BEGIN

    iFieldLev = 0
;      
;  Need to handle fields on interfaces levels(instead of midpoints)      
;
    IF iLev GT 0 THEN BEGIN

      nLevs = N_ELEMENTS(fieldlevs)
      newLevs = *line.levs
      newNLevs = N_ELEMENTS(newLevs)
      newLevSelected = newLevs[iLev]
      iFieldLev = IXFIND_NEAREST(fieldlevs, newlevSelected)
      levSelected = fieldlevs[iFieldLev]
      IF iFieldLev GT 0 THEN levSelectedLower = fieldlevs[iFieldLev-1]
      IF iFieldLev LT nLevs-1 THEN levSelectedUpper = fieldlevs[iFieldLev+1]

      IF newLevSelected LT levSelected AND iFieldLev 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 iFieldLev LT nLevs-1 THEN BEGIN
   
        levDiffLower = newLevSelected - levSelected
        levDiffUpper = levSelectedUpper - newLevSelected
        levDiffAll = levSelectedUpper - levSelected
        levRatLower = levDiffLower / levDiffAll
        levRatUpper = levDiffUpper / levDiffAll
   
      ENDIF

    ENDIF

  ENDIF

  linedata = fltarr(nlat)
  linedata[*] = 0.
;
;  Handle WACCM special case
;   
  IF line.ftype EQ 'WACCM' THEN BEGIN

    IF line.amplitude NE 0 AND field.name NE 'YAXZ' THEN BEGIN  ;Get wavenumber amplitude
 
  print, 'calculating amplitude in defmeriddata '
  
      iwvn = line.wavenumber
      ampdata = fltarr(nlon)
      ampmin = 1.0E+24
      for j=0,nlat-1 do begin
      
    	ampdata[*] = float(data[*,j,iFieldLev,line.imtime])
    	ampdata[*] = ABS(FFT(ampdata[*], -1)) * 2.
    	if ampdata[iwvn] LT ampmin AND ampdata[iwvn] GT 0.0 then ampmin = ampdata[iwvn]
    	linedata[j] = ampdata[iwvn]

      endfor

      zeroind = WHERE(linedata[*] EQ 0.0)
      if zeroind[0] NE -1 then linedata[zeroind] = ampmin
 
      IF line.amplitude EQ 2 THEN linedata[*] = ALOG10(linedata[*])
;
;  Update min/max for data
; 
      line.fmin = MIN(linedata)
      line.fmax = MAX(linedata) 

    ENDIF ELSE BEGIN ;Not wavenumber amplitude

      if line.zonalmean eq 0 then begin  ; No zonal mean
    
        linedata[*] = float(data[ilon,*,iFieldLev,line.imtime])
;      
;  Need to handle fields on interfaces levels(instead of midpoints)      
;
        IF iLev GT 0 THEN BEGIN
        
	  IF iFieldLev GT 0 AND iFieldLev LT nLevs-1 THEN BEGIN 
	
	    IF isdensity(line.field.name) THEN BEGIN
	  
	      dataLower = ALOG10(data[ilon,*,IFieldLev-1,line.imtime])
	      dataLev = ALOG10(data[ilon,*,IFieldLev,line.imtime])
	      dataUpper = ALOG10(data[ilon,*,IFieldLev+1,line.imtime])
	    
	    ENDIF ELSE BEGIN
	  
	      dataLower = data[ilon,*,IFieldLev-1,line.imtime]
	      dataLev = data[ilon,*,IFieldLev,line.imtime]
	      dataUpper = data[ilon,*,IFieldLev+1,line.imtime]

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

	    IF isdensity(line.field.name) THEN linedata = 10^(linedata)
	    
          ENDIF
  
        ENDIF
      
      endif else begin ; zonal mean
      
        for k=0,nlat-1 do begin
          linedata[k] = float(mean(data[*,k,iFieldLev,line.imtime]))
        endfor
	
        dataLower = linedata
        dataLev = linedata
        dataUpper = linedata
;      
;  Need to handle fields on interfaces levels(instead of midpoints)      
;
        IF iLev GT 0 THEN BEGIN
        
	  IF iFieldLev GT 0 AND iFieldLev LT nLevs-1 THEN BEGIN 
	
	    IF isdensity(line.field.name) THEN BEGIN
	  
              FOR k=0,nlat-1 do begin
	      
	        dataLower[k] = ALOG10(mean(data[*,k,IFieldLev-1,line.imtime]))
	        dataLev[k] = ALOG10(mean(data[*,k,IFieldLev,line.imtime]))
	        dataUpper[k] = ALOG10(mean(data[*,k,IFieldLev+1,line.imtime]))
		
	      ENDFOR
	    
	    ENDIF ELSE BEGIN
	  
              FOR k=0,nlat-1 do begin
	      
	        dataLower[k] = mean(data[*,k,IFieldLev-1,line.imtime])
	        dataLev[k] = mean(data[*,k,IFieldLev,line.imtime])
	        dataUpper[k] = mean(data[*,k,IFieldLev+1,line.imtime])
		
	      ENDFOR

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

	    IF isdensity(line.field.name) THEN linedata = 10^(linedata)
	    
          ENDIF
  
        ENDIF
           
      endelse ; zonal mean
      
    ENDELSE ;Not wavenumber amplitude

  ENDIF ELSE BEGIN ;Not WACCM 
;
;  Handle wave amplitude calculation here
;      
    IF line.amplitude NE 0 AND field.name NE 'YAXZ' THEN BEGIN
  
      iwvn = line.wavenumber
      ampdata = fltarr(nlon)
      ampmin = 1.0E+24
      for j=0,nlat-1 do begin
      
    	ampdata[*] = float(data[*,j,iLev,line.imtime])
    	ampdata[*] = ABS(FFT(ampdata[*], -1)) * 2.
    	if ampdata[iwvn] LT ampmin AND ampdata[iwvn] GT 0.0 then ampmin = ampdata[iwvn]
    	linedata[j] = ampdata[iwvn]
    
      endfor

      zeroind = WHERE(linedata[*] EQ 0.0)
      if zeroind[0] NE -1 then linedata[zeroind] = ampmin
 
      IF line.amplitude EQ 2 THEN linedata[*] = ALOG10(linedata[*])

      fmin = min(linedata) & fmax = max(linedata)
      line.fmin = fmin
      line.fmax = fmax

    ENDIF ELSE BEGIN
      
      if line.zonalmean eq 0 then begin
 
          linedata[*] = float(data[ilon,*,iLev,line.imtime])

      endif else begin ; zonal mean
        for k=0,nlat-1 do begin
          linedata[k] = float(mean(data[*,k,iLev,line.imtime]))
        endfor
      endelse ; zonal mean
    
    ENDELSE;Not wavenumber amplitude    
 
  ENDELSE; WACCM or not
;
; Take log10 of densities:
;
    if (line.log10 eq 'density fields' and isdensity(field.name) gt 0) or $
       (line.log10 eq 'this frame only') or $
       (line.log10 eq 'all fields') then begin
      missing_value = float(line.missing_value)
      log10f,linedata,missing_value
    endif    
;
; Update line structure:
;
  line.data = ptr_new(linedata)
;
;  Update min/max for data
;  
  fminmax,linedata,line.fmin,line.fmax,line.missing_value

end