PRO interpwdata, slice
;
;  This routine interpolates WACCM history file data to the vertical pressure level grid which TIME-GCM uses.  It
;  will check for 1-D data with vertical level as the dimension or 4-D data with vertical level as the 3rd dimension
;  of the data.
;
;  Get arrays needed from slice
;
  field = slice.field 
  data = *slice.data
  levs = *field.levs
  nLevs = N_ELEMENTS(levs)
  newLevs = *field.levsPlot
  newNLevs = N_ELEMENTS(newLevs)
;
;  Get the index of the minimum and maximim vertical levels to use in original WACCM data  
;  
  levMin = slice.levmin
  iLevMin = IXFIND_NEAREST(levs,levMin)
  levMax = slice.levmax
  iLevMax = IXFIND_NEAREST(levs,levMax)
  IF iLevMin GT 0 THEN iLevMin = iLevMin
  IF iLevMax LT nLevs-1 THEN iLevMax = iLevMax
;
;  Get number of dimensions in input data and find vertical level dimension if more than one dimension
; 
  nDatDims = SIZE(data,/N_DIMENSIONS)
  
  IF nDatDims GT 1 THEN BEGIN
  
    varDims = SIZE(data,/DIMENSIONS)
    levDim = WHERE(varDims eq nLevs)
    IF levDim LT 0 THEN nLevs = nLevs + 1
    levDim = WHERE(varDims eq nLevs)
    
  ENDIF ELSE IF nDatDims EQ 1 THEN BEGIN
  
    varDims = N_ELEMENTS(data)
    levDim = 0
  
  ENDIF  
;
;  Check to make sure vertical level dimension found in input data
;  
  IF levDim LT 0 THEN BEGIN
 
    PRINT, 'INTERPWDATA: Could not find lev dimension in data '
    RETURN

  ENDIF  
;
;  Get the corresponding minimum and maximum levels in new vertical level array using original vertical min/max
;
   indexNewLevMin = IXFIND_NEAREST(newLevs,levmin)
   indexNewLevMax = IXFIND_NEAREST(newLevs,levmax)
   
   IF indexNewLevMax[0] EQ -1 THEN indexNewLevMax = newNLevs - 1
  
   iNewLevMin = 0
   iNewLevMax = newNLevs - 1
   nINLMin = N_ELEMENTS(indexNewLevMin)
;
;  New vertical min/max within original min/max
;
   IF indexNewLevMin[nINLMin-1] GT iNewLevMin AND indexNewLevMax[0] LT iNewLevMax AND indexNewLevMax[0] GE 0 THEN BEGIN
  
     iNewLevMin = indexNewLevMin[nINLMin-1]
     iNewLevMax = indexNewLevMax[0]
;
;  New vertical min above original vertical min
;   
   ENDIF ELSE IF indexNewLevMin[nINLMin-1] GT iNewLevMin AND indexNewLevMax[0] GE iNewLevMax THEN BEGIN

     iNewLevMin = indexNewLevMin[nINLMin-1]
;
;  New vertical max above original vertical max
;      
   ENDIF ELSE IF indexNewLevMin[nINLMin-1] LE iNewLevMin AND indexNewLevMax[0] LT iNewLevMax AND indexNewLevMax[0] GE 0 THEN BEGIN
  
     iNewLevMax = indexNewLevMax[0]
   
   ENDIF	 
;
;  Form subsetted 'standard' vertical level array
;
   tempLevs = newLevs[iNewLevMin:iNewLevMax]
  
   newLevs = tempLevs
   newNLevs = N_ELEMENTS(newLevs)

;  ENDELSE
;
;   Check the number of dimensions in the data field to handle as 1-D or 4-D.
;
;   Do case of one dimensional data
;
  IF nDatDims EQ 1 THEN BEGIN
  
    IF levDim EQ 0 THEN BEGIN
;
;  Form output data array
;    
      newData = FLTARR(newNLevs)
;
;  Get subset of input data and levels
;   
      iData = data[iLevMin:iLevMax]
      levsIn = levs[iLevMin:iLevMax]
;
;  Subset only non-missing data and levels
;
      indexGood = WHERE(iData NE FLOAT(slice.missing_value))
      iDataG = iData[indexGood]
      levsInG = levsIn[indexGood]
;
;  Get the output levels
;	     
      iNewLevGoodH = newNLevs - 1
      iNewLevGoodL = 0
      IF indexGood[0] GE 0 THEN BEGIN
      
        iLevGoodHigh = indexGood[N_ELEMENTS(indexGood)-1]
        iLevGoodLow = indexGood[0]
	iNewLevGoodHAll = WHERE(newLevs LT levsIn[iLevGoodHigh])
	iNewLevGoodLAll = WHERE(newLevs GT levsIn[iLevGoodLow])
	IF iNewLevGoodHAll[0] GE 0 THEN iNewLevGoodH = iNewLevGoodHAll[N_ELEMENTS(iNewLevGoodHAll)-1]
	IF iNewLevGoodLAll[0] GE 0 THEN iNewLevGoodL = iNewLevGoodLAll[0]
	  	  
      ENDIF 
;
;  Interpolate data to 'standard' output levels
;      
      loglin = 1
      IF  isdensity(slice.field.name) GT 0 THEN loglin = 0

      vfnd, levsInG, iDataG, newLevs, oData, cut=0, linear=loglin
	
      newData = oData
;
;  If data is a subset of output vertical range then extrapolate data to fill
;
      IF iNewLevGoodH GT 0 AND iNewLevGoodH LT newNLevs-1 THEN $
         newData[iNewLevGoodH+1:newNLevs-1] = newData[iNewLevGoodH] + (newData[iNewLevGoodH] - newData[iNewLevGoodH-1])
      
      IF iNewLevGoodL GT 1 THEN $      
         newData[0:iNewLevGoodL-1] = newData[iNewLevGoodL] + (newData[iNewLevGoodL] - newData[iNewLevGoodL+1])

    ENDIF ELSE BEGIN
    
      PRINT, 'INTERPWDATA: Levels is not the dimension of data as expected'
      RETURN
      
    ENDELSE 
;      
;  Do case of 2 dimensions.  Assumes levels is the second dimension
;
  ENDIF ELSE IF nDatDims EQ 2 THEN BEGIN
  
    IF levDim EQ 1 THEN BEGIN
    
      nDim1 = varDims[0]
;
;  Form output data array which will be on 'standard' vertical grid.  
;  Also, get 'nonstandard' input vertical levels
;          
      newData = FLTARR(nDim1,newNLevs)
      levsIn = levs[iLevMin:iLevMax]

      FOR iDim1=0,nDim1-1 DO BEGIN
;
;  Get input data on 'nonstandard' input vertical levels and find where data are not missing
;                
        iData = REFORM(data[iDim1,iLevMin:iLevMax])
         
	indexGoodL = -1
	indexGoodH = -1
      
        FOR iG = N_ELEMENTS(iData)-1, 0, -1 DO IF (iData(iG) NE FLOAT(slice.missing_value)) THEN indexGoodL = iG
        FOR iG = 0, N_ELEMENTS(iData)-1 DO IF (iData(iG) NE FLOAT(slice.missing_value)) THEN indexGoodH = iG
	
        IF indexGoodL NE -1 AND indexGoodH NE -1 THEN BEGIN

          iDataG = iData[indexGoodL:indexGoodH]
          levsInG = levsIn[indexGoodL:indexGoodH] 
	     
          iMinLevGoodH = newNLevs - 1
          iMinLevGoodL = 0
      
          iLevGoodHigh = indexGoodH
	  iLevGoodLow = indexGoodL
	  iNewLevGoodHAll = WHERE(newLevs LE levsIn[iLevGoodHigh])
	  iNewLevGoodLAll = WHERE(newLevs GE levsIn[iLevGoodLow])
	  IF iNewLevGoodHAll[0] GE 0 THEN iNewLevGoodH = iNewLevGoodHAll[N_ELEMENTS(iNewLevGoodHAll)-1]
	  IF iNewLevGoodLAll[0] GE 0 THEN iNewLevGoodL = iNewLevGoodLAll[0]
;	
;  Interpolate data to 'standard' vertical grid and then check for data outside original data min/max
;	
	  loglin = 1
          IF  isdensity(slice.field.name) GT 0 THEN loglin = 0
	
          vfnd, levsInG, iDataG, newLevs, oData, cut=1, linear=loglin
	
          newData[iDim1,*] = oData
;
;  Fill in last level(s) with extrapolated good value(s)
;
          IF iNewLevGoodH GT 0 AND iNewLevGoodH LT newNLevs-1 THEN BEGIN
	  
	    newData[iDim1,iNewLevGoodH+1:newNLevs-1] = slice.missing_value

          ENDIF
;
;  Fill in first level(s) with first good value
;
          IF iNewLevGoodL GT 1 THEN BEGIN	  
	    
	   newData[iDim1,0:iNewLevGoodL-1] = slice.missing_value

          ENDIF	  
	    
        ENDIF ELSE BEGIN
	
	  newData[iDim1,*] = slice.missing_value 
	
	ENDELSE 
	
      ENDFOR
 
    ENDIF ELSE BEGIN
    
      PRINT, 'INTERPWDATA: Levels is not the second dimension of data as expected'
      RETURN
      
    ENDELSE        
;
;  4-D data with vertical levels assumed as the 3rd dimension
;  
   ENDIF ELSE IF nDatDims EQ 4 THEN BEGIN
  
    IF levDim EQ 2 THEN BEGIN
    
      nDim1 = varDims[0]
      nDim2 = varDims[1]
      nDim4 = varDims[3]
      
      newData = FLTARR(nDim1,nDim2,newNLevs,nDim4)
   
      levsIn = levs[iLevMin:iLevMax]
      
      FOR iDim1=0,nDim1-1 DO BEGIN
        FOR iDim2=0,nDim2-1 DO BEGIN
          FOR iDim4=0,nDim4-1 DO BEGIN
      
            iData = TRANSPOSE(data[iDim1,iDim2,iMin:iMax,iDim4])
      
            indexGood = WHERE(iData NE FLOAT(slice.missing_value))
            iDataG = iData[indexGood]
            levsInG = levsIn[indexGood]      
      
            iNewLevGoodH = newNLevs - 1
            iNewLevGoodL = 0
            IF indexGood[0] GE 0 THEN BEGIN
      
              iLevGoodHigh = indexGood[N_ELEMENTS(indexGood)-1]
	      iLevGoodLow = indexGood[0]
	      iNewLevGoodHAll = WHERE(newLevs LE levsIn[iLevGoodHigh])
	      iNewLevGoodLAll = WHERE(newLevs GE levsIn[iLevGoodLow])
	      IF iNewLevGoodHAll[0] GE 0 THEN iNewLevGoodH = iNewLevGoodHAll[N_ELEMENTS(iNewLevGoodHAll)-1]
	      IF iNewLevGoodLAll[0] GE 0 THEN iNewLevGoodL = iNewLevGoodLAll[0]

            ENDIF 

	    loglin = 1
            IF  isdensity(slice.field.name) GT 0 THEN loglin = 0
	
            vfnd, levsInG, iDataG, newLevs, oData, cut=1, linear=loglin
	
            newData[iDim1,iDim2,*,iDim4] = oData
            IF iNewLevGoodL GE 0 AND iNewLevGoodL LT newNLevs-1 THEN $
                 newData[iDim1,iDim2,0:iNewLevGoodL-1,iDim4] = newData[iDim1,iDim2,iNewLevGoodL,iDim4]
            IF iNewLevGoodH GT 0 AND iNewLevGoodH LT newNLevs-1 THEN $
	         newData[iDim1,iDim2,iNewLevGoodH+1:newNLevs-1,iDim4] = newData[iDim1,iDim2,iNewLevGoodH,iDim4]
	
	  ENDFOR
	ENDFOR	
      ENDFOR
    
    ENDIF ELSE BEGIN
    
      PRINT, 'INTERPWDATA: Levels not the third dimension of data as expected'
      RETURN
      
    ENDELSE        
;
;  Check for common errors
;  
 ENDIF ELSE IF nDatDims EQ -1 THEN BEGIN

    PRINT, 'INTERPWDATA: Did not find levels dimension in data'
    RETURN
  
 ENDIF ELSE BEGIN
  
    PRINT, 'INTERPWDATA: Problem determining input data dimensions '
    RETURN
    
  ENDELSE
;
;  Update slice structure data and levels and field structure
;
  slice.data = PTR_NEW(newData)
  slice.levs = PTR_NEW(newLevs)
  fminmax,*slice.data,slice.fmin,slice.fmax,slice.missing_value
  
  field = slice.field   
  
END