;==============================================================================
; Function geov_extractor::number_density
;
; Returns and array of number densities (molecules/cm^3) of the current volume
; mixing ratio field (VMR).  The array return has dimensions appropriate for
; the count vector.
;
; Parameters:
;
;    offset: A vector containing the starting position of the field in the
;            netCDF file.
;
;     count: A vector containing the counts of the data values of the field
;            in the netCDF file. COUNT is a 1-based vector with an element
;            for each dimension of the data field.
;
;==============================================================================
Function geov_extractor::number_density, offset, count

  ; check that the units of the current field is 'VMR'
  vunits = self.ncfiles->get_vunits()
  units = (*vunits)[self.vid]
  if ( STRCMP( units, 'kg/kg', 3,/fold_case ) or $
       STRCMP( units, 'mmr', 3,/fold_case )   or $
       STRCMP( units, 'g/g', 3,/fold_case )  ) then begin

       result = DIALOG_MESSAGE([ 'The units of the mixing ratio are MMR.',$
                                 'This computation assumes units of VMR.',$  
                                 'Do you wish to continue?'] , /question )
       if (STRCMP( result,'no',2,/fold_case)) then return, ptr_new() ; null ptr
  endif else if ( (STRCMP( units, 'VMR', 3,/fold_case ) eq 0)       and $
       (STRCMP( units, 'mol/mol', 7,/fold_case ) eq 0) and $
       (STRCMP( units, 'mole/mole', 8,/fold_case ) eq 0) and $
       (STRCMP( units, 'cm3/cm3', 5,/fold_case ) eq 0) and $
       (STRCMP( units, 'm3/m3', 5,/fold_case ) eq 0)   ) then begin
       result = DIALOG_MESSAGE([ 'The units of the selected variable is not a ',$
                                 'mixing ratio.  The result of this extractor ',$
                                 'will be the variable times background density.',$
                                 'Do you wish to continue?'] , /question )
       if (STRCMP( result,'no',2,/fold_case)) then return, ptr_new() ; null ptr
;    message,'The variable has incorrect units for the number density calculation'
;    return, -1
  endif
 
  k = 1.380044E-16 ; boltzman's constant (erg/K)

  ; retrieve the needed pressure and temperature fields 
  pressure = self->pressure_field('hyam','hybm', offset, count )
  temperature = self->temperature_field( offset, count )

  P = *pressure
  u = *self.data
  T = *temperature

  ; assuming P,u,T are all the same size we can imploy IDL's
  ; array opperators
  density = 10.*(P*u)/( k*T )
  ptr_free,pressure
  ptr_free,temperature

  return, ptr_new(density, /NO_COPY)


End



;=================================================================================
; Function geov_extractor::int_col_dens
;
; Returns an array of number total integrated column densities (molecules/cm^2)
; of the current volume mixing ratio field (VMR). The array return has dimensions
; appropriate for the count vector with the exception that the pressure level
; dimension has been colapsed, i.e., the result is independent of pressure level
; since pressure level dimension has been integrated over.
;
; Parameters:
;
;    offset: A vector containing the starting position of the field in the
;            netCDF file.
;
;     count: A vector containing the counts of the data values of the field
;            in the netCDF file. COUNT is a 1-based vector with an element
;            for each dimension of the data field.
;
;    dobson: Set dobson keyword if Dobson Units (DU) is desired
;
;=================================================================================
Function geov_extractor::int_col_dens, offset, count, dobson=dobson, above_lev=above_lev, below_lev=below_lev, intlevels=intlevels
  @nccoord_pars

  ; check that the units of the current field is 'VMR'
  vunits = self.ncfiles->get_vunits()
  units = (*vunits)[self.vid]
  if ( STRCMP( units, 'kg/kg', 3,/fold_case ) or $
       STRCMP( units, 'mmr', 3,/fold_case )   or $
       STRCMP( units, 'g/g', 3,/fold_case )  ) then begin

       result = DIALOG_MESSAGE([ 'The units of the mixing ratio are MMR.',$
                                 'This computation assumes units of VMR.',$  
                                 'Do you wish to continue?'] , /question )
       if (STRCMP( result,'no',2,/fold_case)) then return, ptr_new()
  endif else if ( (STRCMP( units, 'VMR', 3,/fold_case ) eq 0)       and $
       (STRCMP( units, 'mol/mol', 7,/fold_case ) eq 0) and $
       (STRCMP( units, 'mole/mole', 8,/fold_case ) eq 0) and $
       (STRCMP( units, 'cm3/cm3', 5,/fold_case ) eq 0) and $
       (STRCMP( units, 'm3/m3', 5,/fold_case ) eq 0)   ) then begin
    message,'The variable has incorrect units for the number density calculation'
;    return, -1 
    return, ptr_new() 
  endif 

  result = self->variable_dimensions( self.vid )
  itim = result.timdim

  ntimes = 1
  if (itim ge 0) then ntimes = count[itim]

  ; get the ncfile coord num for each of the dimensions
  dsizes = self.ncfiles->get_dsizes()

  ;get the vertical coord index and total number of pressure levels

  ilon = result.londim
  ilat = result.latdim
  ilev = result.levdim
  lev_did = self.ncfiles->get_dimid_of_var(ilev,self.vid)
  nlevs = (*dsizes)[lev_did]

  col_offset = offset
  col_count  = count
  int_count  = count
  int_count[ilev] = 1

  nmax = 10 ; max number of times process at a time

  nloops = (ntimes/nmax)

  total_times = 0

  ; set the col_count according to the levels we want to integrate over
  if ( n_elements(above_lev) ) then begin
    ; we want to integrate only above a specified level index
    if (self.ncfiles->is_cflag_set(self.ncfiles->get_ilev(),NCCOORD_FLAG_DOWN)) then begin
      ;if above_lev lt 1 then above_lev=1
      col_count[ilev]  = above_lev + 1
      col_offset[ilev] = 0
    endif else begin
      col_count[ilev]  = nlevs - above_lev + 1
      col_offset[ilev] = above_lev 
    endelse
  endif else begin
    ; want to get entire lev range - all pressure levels
    ; we need all pressure levels since we need to integrate over the all press levels
    col_count[ilev]  = nlevs
    col_offset[ilev] = 0 
  endelse

  if ( n_elements(below_lev) ne 0 ) then begin
    if (self.ncfiles->is_cflag_set(self.ncfiles->get_ilev(),NCCOORD_FLAG_DOWN)) then begin
      col_count[ilev]  = col_count[ilev] - below_lev
      col_offset[ilev] = below_lev
    endif else begin
      col_count[ilev]  = col_count[ilev] - ( nlevs - below_lev)
    endelse
  endif 

  while ( total_times lt ntimes ) do begin

    number = (ntimes-total_times) < nmax

    col_count[itim] = number
    int_count[itim] = number
    col_offset[itim] = total_times + offset[itim]

    total_times = total_times + number

    ;need dp for all levels - use col_offset and col_count here
    dp = self->dPress_field( col_offset, col_count )

    ;also need to extract u for all levels
    u = *( self.ncfiles->get_data( self.vid, col_offset, col_count ) )

    ; assuming delP and u have the same dimensions and sizes we can employ
    ; IDL's array multiplication opperator
    delP = *dp

    ptr_free, dp

    u = reform(u,col_count)
    delp = reform(delp,col_count)
    X = u*delP

    X = reform(X,col_count)

    ; define some constants need for the column integration
    k = 1.380044E-16 ; boltzman's constant (erg/K)
    R=2.8704e6       ; gas constant (erg/g/K)
    g=980.616        ; grav acceleration (cm/sec2)

    if ( n_elements(intlevels) eq 0 ) then begin
      ; use IDL's TOTAL function the perform the summation calculation
      ; needed for the integration
      if (col_count[ilev] gt 1) then int_col = (TOTAL(X,ilev+1))*R*10./(k*g) $
      else  int_col = X*R*10./(k*g)
    endif else begin

      ; sum each column separately
      crds = lonarr(4)
      int_col = dblarr( int_count )
      c = lonarr(4,2)
      c[ilev,0] = 0
      c[ilev,1] = nlevs-1
      crds[ilev] = 0
      tmax = count[itim]-1
      jmax = count[ilon]-1
      imax = count[ilat]-1
      for t=0,tmax do begin
        for j=0,jmax do begin
          for i=0,imax do begin
            crds[itim] = t
            crds[ilat] = i
            crds[ilon] = j 
            indices = *(intlevels[ crds[0],crds[1],crds[2],crds[3] ])
            c[itim,*]=t
            c[ilon,*]=j
            c[ilat,*]=i
            X_col = X[ c[0,0]:c[0,1], c[1,0]:c[1,1], c[2,0]:c[2,1], c[3,0]:c[3,1] ]
            int_col[ crds[0],crds[1],crds[2],crds[3] ] = (TOTAL(X_col[min(indices):max(indices)]))*R*10./(k*g)

          endfor
        endfor
      endfor
    endelse

    if (keyword_set(dobson) ) then int_col = int_col/(2.687e16)

    int_col = reform( int_col, int_count )

    if (n_elements(trans_int_cols) eq 0) then trans_int_cols = transpose( int_col ) $
    else trans_int_cols = [trans_int_cols, transpose( int_col )]

  endwhile

  int_cols = transpose( trans_int_cols )

  return, ptr_new(int_cols, /NO_COPY)
End


Function geov_extractor::variable_dimensions, vid
  result = {latdim:-1, londim:-1, levdim:-1, timdim:-1 }

  vndims = (*(self.ncfiles->get_vndims()))[vid]
  d2coord = self.ncfiles->get_d2coord()
  ilat = self.ncfiles->get_ilat()
  ilon = self.ncfiles->get_ilon()
  ilev = self.ncfiles->get_ilev()
  itim = self.ncfiles->get_itim()

  if ( self.ncfiles->var_depends_on_ilevi(self.vid) ) then $
    ilev = self.ncfiles->get_ilevi()

  for idim = 0, vndims-1 do begin
    did = self.ncfiles->get_dimid_of_var(idim, vid)
    ic = (*d2coord)[did]
    
    case ic of 
      ilat: result.latdim = idim
      ilon: result.londim = idim
      ilev: result.levdim = idim
      itim: result.timdim = idim
      else:

    endcase
  endfor

  return, result

end

;==============================================================================
; Function geov_extractor::dPress_field
;
; Returns and array of delta-pressures. The array return has dimensions
; appropriate for the count vector.
;
; Parameters:
;
;    offset: A vector containing the starting position of the field in the
;            netCDF file.
;
;     count: A vector containing the counts of the data values of the field
;            in the netCDF file. COUNT is a 1-based vector with an element
;            for each dimension of the data field.
;
;===============================================================================
Function geov_extractor::dPress_field, offset, count

  ; get the array which maps the gui coord number to the coord number in the
  ; netCDF file
  d2coord = self.ncfiles->get_d2coord()

  ; determine which coordinate is the level dimension
;  ilev_id = self.ncfiles->get_dimid_of_var(self.ncfiles->get_ilev(),self.vid)
;  ilev = (*d2coord)[ilev_id]

  result = self->variable_dimensions( self.vid )

  ilon = result.londim
  ilat = result.latdim
  ilev = result.levdim
  itim = result.timdim

  ; need to add one to the level deminsion since
  ; delta-p = press(lev+1)-press(lev)
  dpcount = count
  dpcount[ilev] = count[ilev]+1

  ; use hyai and hybi to get pressure field at intermediate model grid points
  press = self->pressure_field('hyai','hybi',offset, dpcount)

  ncoords = intarr(4)
  ncoords[*] = 1
  ; assume a 4-D array - remove the redundent dimensions later
;  dp = make_array (count[0], count[1],count[2],count[3] )
  if ilat ge 0 then ncoords[ilat]=count[ilat]
  if ilon ge 0 then ncoords[ilon]=count[ilon]
  if ilev ge 0 then ncoords[ilev]=count[ilev]
  if itim ge 0 then ncoords[itim]=count[itim]
  dp = make_array (ncoords[0], ncoords[1],ncoords[2],ncoords[3] )

  pressure = *press
      
;  crds0 = [ [0,count[0]-1],[0,count[1]-1],[0,count[2]-1],[0,count[3]-1] ]
  crds0 = [ [0,ncoords[0]-1],[0,ncoords[1]-1],[0,ncoords[2]-1],[0,ncoords[3]-1] ]
  crds1 = crds0
  crds1[ *,ilev ] = crds0[ *,ilev] + 1

  press_below = pressure[ crds0[0,0]:crds0[1,0], crds0[0,1]:crds0[1,1], crds0[0,2]:crds0[1,2], crds0[0,3]:crds0[1,3]  ]
  press_above = pressure[ crds1[0,0]:crds1[1,0], crds1[0,1]:crds1[1,1], crds1[0,2]:crds1[1,2], crds1[0,3]:crds1[1,3]  ]

  dp = press_above - press_below

  ptr_free, press

  return, ptr_new(dp, /NO_COPY)
End



;==============================================================================
; Function geov_extractor::pressure_field
;
; Returns and array containing a pressure field. The array return has dimensions
; appropriate for the count vector.
;
; Parameters:
;
;       hya: A string containing the name of a vertical level coefficient
;            need to compute the pressure field
;
;       hyb: A string containing the name of a vertical level coefficient
;            need to compute the pressure field
;
;    offset: A vector containing the starting position of the field in the
;            netCDF file.
;
;     count: A vector containing the counts of the data values of the field
;            in the netCDF file. COUNT is a 1-based vector with an element
;            for each dimension of the data field.
;
;===============================================================================
Function geov_extractor::pressure_field, hya, hyb, offset, count

  P0 = 1.e5
  PS = self->surface_press_field( offset, count, psdims=psdims )

  hya_units=''

  result = self->variable_dimensions(self.vid)

  ilat = result.latdim
  ilon = result.londim
  ilev = result.levdim
  itim = result.timdim

  ; retrieve the vertical level coefficients
  hyam = self->hy_field( hya, offset, count, units=hya_units )
  hybm = self->hy_field( hyb, offset, count )

  ; assume a 4-D pressure field
  ncoords = intarr(4)
  ncoords[*] = 1
  if (ilon ge 0) then ncoords[ilon] = count[ilon]
  if (ilat ge 0) then ncoords[ilat] = count[ilat]
  if (ilev ge 0) then ncoords[ilev] = count[ilev]
  if (itim ge 0) then ncoords[itim] = count[itim]
  pressure = make_array (ncoords[0], ncoords[1],ncoords[2],ncoords[3], /float )
 
  a = *hyam
  b = *hybm
  surPress = *PS

  ptr_free, hyam
  ptr_free, hybm
  ptr_free, PS
  
  c = [ [0,ncoords[0]-1], [0,ncoords[1]-1], [0,ncoords[2]-1], [0,ncoords[3]-1] ] 
 
  nlevs = count[ilev]
  ntims = count[itim]

  if ( not strcmp(hya_units,'pa',2,/fold_case ) ) then a = a*P0

  pscoord0 = lonarr(3)
  pscoord1 = lonarr(3)

  pscoord0[*] = 0
  pscoord1[*] = 0

  if (ilat ge 0) then pscoord1[psdims.latdim] = count[ilat]-1

  if (ilon ge 0) then pscoord1[psdims.londim] = count[ilon]-1 

  press_size = intarr(4)
  spress_size = intarr(3)

  press_size[*] = 1
  spress_size[*] = 1

  if (ilev ge 0 ) then press_size[ilev] = count[ilev]
  if (ilat ge 0 ) then press_size[ilat] = count[ilat]
  if (ilon ge 0 ) then press_size[ilon] = count[ilon]
  if (itim ge 0 ) then press_size[itim] = count[itim]

  if (ilat ge 0 ) then spress_size[psdims.latdim] = count[ilat]
  if (ilon ge 0 ) then spress_size[psdims.londim] = count[ilon]
  if (itim ge 0 ) then spress_size[psdims.timdim] = count[itim]

  pressure = reform(pressure,press_size, /overwrite)
  surPress = reform(surPress,spress_size,/overwrite)

  for t = 0, ntims-1 do begin
    for k = 0, nlevs-1 do begin

      c[*,ilev] = k
      pscoord0[psdims.timdim]=t
      pscoord1[psdims.timdim]=t
      c[*,itim] = t

      pressure[ c[0,0]:c[1,0], c[0,1]:c[1,1], c[0,2]:c[1,2], c[0,3]:c[1,3] ] $
              =  a[k] + surPress[pscoord0[0]:pscoord1[0], pscoord0[1]:pscoord1[1], pscoord0[2]:pscoord1[2] ]*b[k]

    endfor
  endfor

 return, ptr_new(pressure, /NO_COPY)
;  return, pressure
End




;==============================================================================
; Function geov_extractor::temperature_field
;
; Returns and array containing a temperature field. The array return has dimensions
; appropriate for the count vector.
;
; Parameters:
;
;    offset: A vector containing the starting position of the field in the
;            netCDF file.
;
;     count: A vector containing the counts of the data values of the field
;            in the netCDF file. COUNT is a 1-based vector with an element
;            for each dimension of the data field.
;
;===============================================================================
Function geov_extractor::temperature_field, offset, count

  ; get the file id
  fid = self.ncfiles->get_fid()

  ; get the variable id and check that T exists in the netCDF file
  varid = NCDF_VARID( fid, 'T' )
  if ( varid eq -1 ) then varid = NCDF_VARID( fid, 'ST' )   ;; ECHAM naming convention
  if ( varid eq -1 ) then varid = NCDF_VARID( fid, 'st' )   ;; ECHAM naming convention
  if ( varid eq -1 ) then varid = NCDF_VARID( fid, 't' )   ;; ECHAM naming convention
  if ( varid eq -1 ) then begin
     message,'T (or st) does NOT exist in the netCDF file !!!'
     return, -1
  endif

  ; extract the temperature field from the netCDF file
  temperature = self.ncfiles->get_data( varid, offset, count )

  return, temperature
End


;==============================================================================
; Function geov_extractor::altitude_field
;
; Returns and array containing an altitude field. The array return has dimensions
; appropriate for the count vector.
;
; Parameters:
;
;    offset: A vector containing the starting position of the field in the
;            netCDF file.
;
;     count: A vector containing the counts of the data values of the field
;            in the netCDF file. COUNT is a 1-based vector with an element
;            for each dimension of the data field.
;
;===============================================================================
Function geov_extractor::altitude_field, offset, count

  ; get the file id
  fid = self.ncfiles->get_fid()

  ; get the variable id and check that Z3 exists in the netCDF file
  varid = NCDF_VARID( fid, 'Z3' )
  if ( varid eq -1 ) then begin
     message,'Z3 does NOT exist in the netCDF file !!!'
     return, -1
  endif

  ; extract the altitude field from the netCDF file
  geopalt = self.ncfiles->get_data( varid, offset, count )
  
  return, ptr_new(geopalt, /NO_COPY)
;  return, geopalt
End

;==============================================================================
; Function geov_extractor::surface_press_field
;
; Returns and array containing a surface pressure field. The array return has
; dimensions appropriate for the count vector.
;
; Parameters:
;
;    offset: A vector containing the starting position of the field in the
;            netCDF file.
;
;     count: A vector containing the counts of the data values of the field
;            in the netCDF file. COUNT is a 1-based vector with an element
;            for each dimension of the data field.
;
;===============================================================================
Function geov_extractor::surface_press_field, offset, count, psdims=psdims

  ; get the file id
  fid = self.ncfiles->get_fid()

  ; get the variable id for the surf press field and check that
  ; the field exists in the netCDF file
  psvarid = NCDF_VARID( fid, 'PS' ) ; try 'PS'
  if ( psvarid eq -1 ) then  psvarid = NCDF_VARID(fid,'APS') ;try 'APS' -ECHAM
  if ( psvarid eq -1 ) then  psvarid = NCDF_VARID(fid,'ps') ;try lowercase
  if ( psvarid eq -1 ) then  psvarid = NCDF_VARID(fid,'aps') ;try lowercase
  if ( psvarid eq -1 ) then begin
    message,'PS (or APS) does NOT exist in the netCDF file '
    return, -1
  endif

  psdims = self->variable_dimensions(psvarid)

  ; set up the offset and count arrays...

  vndim = (*(self.ncfiles->get_vndims()))[psvarid]

  psoffset = ptr_new( lonarr(vndim))
  pscount  = ptr_new( lonarr(vndim))


;  ; get the array which maps the gui coord number to the coord number in the
;  ; netCDF file
;  d2coord = self.ncfiles->get_d2coord()

;  ; get the id number of each dimension
;  ilon_did = self.ncfiles->get_dimid_of_var(self.ncfiles->get_ilon(),self.vid)
;  ilat_did = self.ncfiles->get_dimid_of_var(self.ncfiles->get_ilat(),self.vid)
;  itim_did = self.ncfiles->get_dimid_of_var(self.ncfiles->get_itim(),self.vid)

;  ; get the ncfile coord num for each of the dimensions
;  ilat = (*d2coord)[ilat_did]
;  ilon = (*d2coord)[ilon_did]
;  itim = (*d2coord)[itim_did]

  psres = self->variable_dimensions( psvarid )
  res = self->variable_dimensions( self.vid )

  ; setup the count and offset vectors so that the appropriate subset is
;;  ; extrated from the file  -- assume that the surface press field is
;;  ; ordered as (lon, lat, time)


  if (res.londim ge 0) then (*pscount)[psres.londim] = count[res.londim]
  if (res.latdim ge 0) then (*pscount)[psres.latdim] = count[res.latdim]
  if (res.timdim ge 0) then (*pscount)[psres.timdim] = count[res.timdim]

  if (res.londim ge 0) then (*psoffset)[psres.londim] = offset[res.londim]
  if (res.latdim ge 0) then (*psoffset)[psres.latdim] = offset[res.latdim]
  if (res.timdim ge 0) then (*psoffset)[psres.timdim] = offset[res.timdim]

;  (*psoffset)[0] = offset[ilon]
;  (*psoffset)[1] = offset[ilat]
;  (*psoffset)[2] = offset[itim]


  

  ; extract the surface press field from the netCDF file
  psdata = self.ncfiles->get_data( psvarid, *psoffset,*pscount)

  ptr_free, psoffset
  ptr_free, pscount

  return, psdata

  ;return, ptr_new(psdata, /NO_COPY)
End



;==============================================================================
; Function geov_extractor::hy_field
;
; Returns and array containing a vertical level coefficient field. The array
; return has dimensions appropriate for the count vector.
;
; Parameters:
;
;    offset: A vector containing the starting position of the field in the
;            netCDF file.
;
;     count: A vector containing the counts of the data values of the field
;            in the netCDF file. COUNT is a 1-based vector with an element
;            for each dimension of the data field.
;
;===============================================================================
Function geov_extractor::hy_field, hyname, offset, count, units=units
  ; get the file id
  fid = self.ncfiles->get_fid()

  ; get the variable id and check that hyname field exists in the file
  varid = NCDF_VARID( fid, hyname )
  if ( varid eq -1 ) then begin
    message, string(hyname+ ' does NOT exist in the netCDF file !!!')
    return, -1
  endif

;  if ( n_elements(units) gt 0 ) then begin

   result = NCDF_VARINQ(fid, varid)
   natts = result.natts
   attnames = strarr(natts)
   FOR i=0L,natts-1 DO attnames[i] = NCDF_ATTNAME(fid, varid, i)
   ;; get units
   wa = ( where(strlowcase(attnames) EQ 'units', cnt) )[0]

   if  (cnt ne 0) then begin 
     NCDF_ATTGET, fid, varid, 'units', value
     units = strtrim(string(value),2)
   endif 
;  endif

  ; set up the offset and count arrays...

  vndim = (*(self.ncfiles->get_vndims()))[varid]


  hyoffset = ptr_new( lonarr(vndim))
  hycount  = ptr_new( lonarr(vndim))

  dsizes = self.ncfiles->get_dsizes()

;  ; get the array which maps the gui coord number to the coord number in the
;  ; netCDF file
;  d2coord = self.ncfiles->get_d2coord()

;  ; get the id number of each dimension
;  ilev_did = self.ncfiles->get_dimid_of_var(self.ncfiles->get_ilev(),self.vid)

;  ; get the ncfile coord num for each of the dimensions
;  ilev = (*d2coord)[ilev_did]

  result = self->variable_dimensions(self.vid)
  ilev = result.levdim

  (*hycount)[0] = count[ilev]
  (*hyoffset)[0] = offset[ilev]

  ; extract the vertical coefficient data from the netCDF file
  hydata = self.ncfiles->get_data( varid, *hyoffset, *hycount)

  ptr_free, hyoffset
  ptr_free, hycount

  return, hydata
End

Function geov_extractor::check_tot_col
    @geov_pars

    ; retrieve coordinate mask for current variable selection
    pvbits= self.ncfiles->get_pvbits()   ; coordinate bits for all plotting variables
    mask = ((*pvbits)[self.gui->get_ipvar()])[0] ; use first variable selected

    if (self.colflag ne 0) then mask = mask and (not LEVBIT) ; turn off lev dependency

   ptype = self.gui->get_ptype()

   ; check that plot type is compatible with variable selected
   pstruct = geov_get_bit_struct(ptype)

   ; X axis requested for lat, lon, lev or time but plotting variable does not have that coordinate
   if (   ((pstruct.xbit ne VARBIT) and ((pstruct.xbit and mask) eq 0))  $

   ; Y axis requested for lat, lon, lev or time but plotting variable does not have that coordinate
       or ((pstruct.ybit ne VARBIT) and ((pstruct.ybit and mask) eq 0))  $

   ; user requested averaging on come coordinate but variable does not depend on that coordinate
       or ((pstruct.avebit gt 0) and ((pstruct.avebit and mask) eq 0)) ) then begin

         if (keyword_set(changeit)) then begin

             self->choose_default_plot
             return, 1

         endif else begin

             result = dialog_message(['EXTRACTOR SELECTION IS INCOMPATIBLE WITH PLOT TYPE.', $
                                      'PLEASE SELECT ANOTHER PLOT TYPE OR ANOTHER EXTRACTOR'])
             return, 0
         endelse

    endif

  return, 1 ; plot is valid for this operator
End

Function geov_extractor::troposphere_levels, offset, count, dims=dims
  @geov_pars
  @nccoord_pars

  ; get the file id
  fid = self.ncfiles->get_fid()

  ; get the variable id for the surf press field and check that
  ; the field exists in the netCDF file
  zvarid = NCDF_VARID( fid, 'TROPLEV' ) ; try 'PS'
  if ( zvarid eq -1 ) then begin
    message,'TROPLEV does NOT exist in the netCDF file '
    return, ptr_new()
  endif

  dims = self->variable_dimensions(zvarid)

  ; set up the offset and count arrays...

  vndim = (*(self.ncfiles->get_vndims()))[zvarid]

  zoffset = lonarr(vndim)
  zcount  = lonarr(vndim)

  zres = self->variable_dimensions( zvarid )
  res = self->variable_dimensions( self.vid )

  zcount[zres.londim] = count[res.londim]
  zcount[zres.latdim] = count[res.latdim]
  zcount[zres.timdim] = count[res.timdim]

  zoffset[zres.londim] = offset[res.londim]
  zoffset[zres.latdim] = offset[res.latdim]
  zoffset[zres.timdim] = offset[res.timdim]

  zdata = self.ncfiles->get_data( zvarid, zoffset, zcount )
  *zdata = reform(*zdata,zcount )

  zcords = intarr(3)
  pcords = intarr(4)
  lcords = intarr(4)

  poffset = offset
  pcount  = count

  dsizes = self.ncfiles->get_dsizes()
  ;get the vertical coord index and total number of pressure levels

  ilev = res.levdim
  lev_did = self.ncfiles->get_dimid_of_var(ilev,self.vid)
  nlevs = (*dsizes)[lev_did]
  poffset[res.levdim ] = 0
  pcount[ res.levdim ] = nlevs

  p = self->pressure_field('hyai','hybi',poffset, pcount )

  if (self.ncfiles->is_cflag_set(self.ncfiles->get_ilev(),NCCOORD_FLAG_DOWN)) then begin
    srfind = nlevs-1
  endif else begin
    srfind = 0
  endelse

  p_col = fltarr( nlevs )

  tlvscount = count
  tlvscount[res.levdim] = 1
  
  troplevs = ptrarr(tlvscount)

  for t=0,count[res.timdim]-1 do begin
    for i=0,count[res.latdim]-1 do begin
      for j=0,count[res.londim]-1 do begin

        zcords[zres.timdim]=t
        zcords[zres.latdim]=i
        zcords[zres.londim]=j

        pcords[res.timdim]=t
        pcords[res.latdim]=i
        pcords[res.londim]=j
        pcords[res.levdim]=srfind

        lcords[res.timdim]=t
        lcords[res.latdim]=i
        lcords[res.londim]=j
        lcords[res.levdim]=0

        z = (*zdata)[ zcords[0],zcords[1],zcords[2] ]
        p0 = (*p)[ pcords[0],pcords[1],pcords[2],pcords[3] ]
        p_trop = p0*exp( -z/H0 )

        for k=0,nlevs-1 do begin
          pcords[res.levdim]=k
          p_col[k] = (*p)[ pcords[0],pcords[1],pcords[2],pcords[3] ]
        endfor          
        if (p_trop gt max(p_col) or p_trop lt min(p_col)) then begin
           print,'tropopause is outside data range... p_trop,max(p_col),min(p_col) ',p_trop,max(p_col),min(p_col)
           return, ptr_new()
        endif
        ind = where( p_col ge p_trop )
        if (ind[0] eq -1) then begin
          print,'not able to find troposphere levels... '
          print,'p_col = ',p_col
          print,'p_trop = ',p_trop
          return, ptr_new()
        endif

        troplevs[ lcords[0],lcords[1],lcords[2],lcords[3] ] = ptr_new(ind)

      endfor
    endfor
  endfor

  ptr_free, zdata
  ptr_free, p

  return, ptr_new( troplevs )

End

Function geov_extractor::trop_col_dens, offset, count, dobson=dobson

  troplevs = self->troposphere_levels( offset, count, dims=trldims )

  if ( not ptr_valid(troplevs) ) then return, ptr_new()

  icd = self->int_col_dens( offset, count, intlevels=*troplevs, dobson=dobson )

  ptr_free, troplevs

  return, icd

end

Pro geov_extractor::extract
  return
End