;
;-----------------------------------------------------------------------
pro proctuv,info,ncdata,field,nlev,lbcname
;
; Process TN, UN, or VN (new history format). Use TLBC, ULBC, or VLBC 
; from the file for lower boundary, then interpolate to interface levels.
;
; First get lower boundary lbc:
;
  fields = *info.fields
  for i=0,info.nflds-1 do begin
    if (fields[i].name eq lbcname) then begin
      if ptr_valid(fields[i].data) then begin ; lbc var has been read
        lbc = fields[i].data
;       print,'procfield: found ',lbcname,' (it had been read): min,max=',$
;         min(lbc),',',max(lbc)
      endif else begin
        varget,info,fields[i],lbc
;       print,'procfield: found ',lbcname,' (it had NOT been read): min,max=',$
;         min(lbc),',',max(lbc)
      endelse
    endif
  endfor
;
; Now use lbc for lower boundary, and interp remainder to interfaces:
;
  ncdata[*,*,0,*] = lbc[*,*,*] ; bottom boundary
  for k=nlev-2,1,-1 do begin   ; interp body
    ncdata[*,*,k,*] = 0.5*(ncdata[*,*,k,*]+ncdata[*,*,k-1,*])
  endfor
;
; 7/16/09: Using nlev-2 value in top nlev-1 level for now because
;   its difficult to ignore missing_values in all plot types          
;   (e.g., congrid in pltlon)
;
  ncdata[*,*,nlev-1,*] = ncdata[*,*,nlev-2,*]
; ncdata[*,*,nlev-1,*] = field.missing_value ; top gets missing data

  file = *info.pfile ; file info structure
  *field.levs = file.ilevs ; TN is now at interfaces
  field.levname = "interface levels";
end
;-----------------------------------------------------------------------
pro interp_body,data,itop
;
; Interpolate body from "half-levels" to "full-levels":
;
for k=itop,1,-1 do begin ; in f90 tgcmproc: do k=itop,2,-1
  data[*,*,k,*] = 0.5*(data[*,*,k,*]+data[*,*,k-1,*])
endfor
end
;-----------------------------------------------------------------------
pro extrap_bot,data,fname
;
; Extrapolate to bottom boundary:
;
; In the following comments, primed numbers (e.g., 2') mean the
;   half level, as received from the history (original ncdata).
; Extrapolation of bottom boundary requires use of half-levels
;   1 and 2 (f(1)=1.5*f(1')-.5f(2')), but half-level 2 (2')
;   was overwritten by interp_body (called before this routine).
;   Thus, we retrieve half-level 2 as follows:
;   f(2) = 0.5*(f(2')+f(1') ==> f(2') = 2.*f(2)-f(1')
;   This is then substituted for f(2') in the extrapolation to
;   give f(1) = 2.*f(1')-f(2)
;
data[*,*,0,*] = 2.*data[*,*,0,*]-data[*,*,1,*]
;
if isdensity(fname) gt 0 then begin
  data0 = data[*,*,0,*]
  index = where(data0 lt 0.,count) ; indices where lbc is < 0
  if count gt 0 then begin
    print,'>> Field ',fname,' # lbc < 0. = ',count,' (out of ',$
      n_elements(data[*,*,0,*]),')'
    data1 = reform(data[*,*,1,*],n_elements(data[*,*,1,*])) ; level above lbc
    data0[index] = data1[index]
    data[*,*,0,*] = data0
;   print,'Field ',fname,' min,max lbc after  check for < 0: ',$
;     min(data[*,*,0,*]),max(data[*,*,0,*])
  endif
endif ; isdensity
end
;-----------------------------------------------------------------------
pro extrap_top,data,nlev,fname
;
; Extrapolate to top (nlev) level:
;
data[*,*,nlev-1,*] = 1.5*data[*,*,nlev-2,*]-0.5*data[*,*,nlev-3,*]
;
if isdensity(fname) gt 0 then begin
  data0 = data[*,*,nlev-1,*]
  index = where(data0 lt 0.,count) ; indices where lbc is < 0
  if count gt 0 then begin
    print,'>> Field ',fname,' # ubc < 0. = ',count,' (out of ',$
      n_elements(data[*,*,nlev-1,*]),')'
    data1 = reform(data[*,*,nlev-2,*],n_elements(data[*,*,nlev-2,*]))
    data0[index] = data1[index]
    data[*,*,nlev-1,*] = data0
;   print,'Field ',fname,' min,max ubc after  check for < 0: ',$
;     min(data[*,*,0,*]),max(data[*,*,nlev-1,*])
  endif ; count > 0
endif ; isdensity
end
;-----------------------------------------------------------------------
pro procfield,info,ncdata,field,z_hist
;
; No processing necessary on f90 tgcmproc output files:
file = *info.pfile ; file info structure
if file.tgcmproc_file gt 0 then begin
; print,'Procfield: this is a tgcmproc processed file -- returning.'
; print,'file.tgcmproc_file=',file.tgcmproc_file
  field.data = ptr_new(ncdata)
  return
endif
;
; No processing necessary on magnetic fields:
if field.grid_type eq 'magnetic' then begin
; print,'Procfield: field ',field.name,' is on magnetic grid,',$
;   ' so no pre-processing is necessary.'
  field.data = ptr_new(ncdata)
  return
endif
size = size(ncdata) ; #dims, dims(#dims), type, total #elements
;print,'size(ncdata)=',size
;
; No processing on fields with fewer than 4 dimensions:
if size[0] lt 4 then begin
; print,'Procfield: field ',field.name,' has ',size[0],' dimensions',$
;   ' and will not be pre-processed.'
  field.data = ptr_new(ncdata)
  return
endif
;
; Check dimensions:
if size[1] ne field.nlon then print,'>>> WARNING procfield: ',$
  'field ',field.name,' 1st dimension ne nlon: size[1]=',size[1],$
  ' nlon=',field.nlon
if size[2] ne field.nlat then print,'>>> WARNING procfield: ',$
  'field ',field.name,' 2nd dimension ne nlat: size[2]=',size[2],$
  ' nlat=',field.nlat
if size[3] ne field.nlev then print,'>>> WARNING procfield: ',$
  'field ',field.name,' 3rd dimension ne nlev: size[3]=',size[3],$
  ' nlat=',field.nlev
nlon = size[1] & nlat = size[2] & nlev = size[3] & ntime = size[4]
;
;  Handle WACCM geopotential height which is in meters on input
;
IF info.ftype EQ 'WACCM' THEN BEGIN

  IF field.name EQ 'Z3' THEN BEGIN
  
    ncdata = ncdata * 1.e-3  ; m to km
    field.units = 'KM'
;   print,'Procfield: convert Z3 from m to km: min,max=',min(ncdata),max(ncdata)

  ENDIF  

ENDIF ELSE BEGIN
;
; 10/28/05: Check for old or new history:
;
if file.oldhist eq 0 then begin
  print,'Procfield: this is a "new" format history..' 
  case field.name of
;
; Geopotential height:
    'Z': begin
      ncdata = ncdata * 1.e-5  ; cm to km
      field.units = 'KM'
;     print,'Procfield: convert Z from cm to km: min,max=',min(ncdata),max(ncdata)
    end
;
; Geopotential height with varying gravity:
    'ZG': begin
      ncdata = ncdata * 1.e-5  ; cm to km
      field.units = 'KM'
;     print,'Procfield: convert ZG from cm to km: min,max=',min(ncdata),max(ncdata)
    end
;
; T,U,V: proctuv (see above) uses LBC for lower boundary, and interpolates
; to interface levels:
    'TN': begin
      proctuv,info,ncdata,field,nlev,'TLBC'
    end
    'UN': begin
      proctuv,info,ncdata,field,nlev,'ULBC'
    end
    'VN': begin
      proctuv,info,ncdata,field,nlev,'VLBC'
    end
;
; No processing required:
    else: begin
;     print,'Note: procfield did nothing to field ',field.name
    end
  endcase ; field.name
  field.data = ptr_new(ncdata)
  return
endif ; new hist
;
print,'Processing field ',field.name,'...'
;
; Processing is field specific:
;
case field.name of
  'TN': begin ; (lon,lat,lev,time)
; Bottom boundary is stored in top slot -- restore it to bottom level:
    ncdata[*,*,0,*] = ncdata[*,*,nlev-1,*] 
    interp_body,ncdata,nlev-2
    ncdata[*,*,nlev-1,*] = ncdata[*,*,nlev-2,*]
  end
  'UN': begin ; (lon,lat,lev,time)
; Bottom boundary is stored in top slot -- restore it to bottom level:
    ncdata[*,*,0,*] = ncdata[*,*,nlev-1,*] 
    interp_body,ncdata,nlev-2
    ncdata[*,*,nlev-1,*] = ncdata[*,*,nlev-2,*]
; Convert from cm/s to m/s:
    ncdata = ncdata*.01
    field.units = 'M/S'
  end
  'VN': begin ; (lon,lat,lev,time)
; Bottom boundary is stored in top slot -- restore it to bottom level:
    ncdata[*,*,0,*] = ncdata[*,*,nlev-1,*] 
    interp_body,ncdata,nlev-2
    ncdata[*,*,nlev-1,*] = ncdata[*,*,nlev-2,*]
; Convert from cm/s to m/s:
    ncdata = ncdata*.01
    field.units = 'M/S'
  end
  'O2': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'OX': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'N4S': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'NOZ': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'CO': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'CO2': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'H2O': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'H2': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'HOX': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'OP': begin
    ncdata[*,*,nlev-1,*] = sqrt(ncdata[*,*,nlev-2,*]^3 / $
      ncdata[*,*,nlev-3,*])
    for k=nlev-2,1,-1 do begin
      ncdata[*,*,k,*] = sqrt(ncdata[*,*,k,*]*ncdata[*,*,k-1,*])
    endfor
    ncdata[*,*,0,*] = ncdata[*,*,0,*]^2 / ncdata[*,*,1,*]
  end
  'CH4': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'AR': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'HE': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'NAT': begin
    interp_body,ncdata,nlev-1
    extrap_bot,ncdata,field.name
  end
  'O21D': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'NO2': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'NO': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'O3': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'O1': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'OH': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'HO2': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'H': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'N2D': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'TI': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'TE': begin
    extrap_top,ncdata,nlev,field.name
    interp_body,ncdata,nlev-2
    extrap_bot,ncdata,field.name
  end
  'O2P': begin
    ncdata[*,*,nlev-1,*] = sqrt(ncdata[*,*,nlev-2,*]^3 / $
      ncdata[*,*,nlev-3,*])
    for k=nlev-2,1,-1 do begin
      ncdata[*,*,k,*] = sqrt(ncdata[*,*,k,*]*ncdata[*,*,k-1,*])
    endfor
    ncdata[*,*,0,*] = ncdata[*,*,0,*]^2 / ncdata[*,*,1,*]
  end
;
; Convert omega to vertical velocity using unprocessed Z from history:
  'W': begin
    zhist = *z_hist
    pzps = fltarr(nlon,nlat,nlev,ntime)
    for k=1,nlev-2 do begin
      pzps[*,*,k,*] = zhist[*,*,k+1,*]-zhist[*,*,k-1,*]
    endfor 
    pzps[*,*,0,*] = 4.*zhist[*,*,1,*]-3.*zhist[*,*,0,*]-zhist[*,*,2,*]
    pzps[*,*,nlev-1,*] = 3.*zhist[*,*,nlev-1,*]-4.*zhist[*,*,nlev-2,*]+$
      zhist[*,*,nlev-3,*]
    ncdata = ncdata*pzps*.01 ; cm/s to m/s
    field.units = 'M/S'
  end
;
; Geopotential height:
  'Z': begin
    zmax = max(ncdata[*,*,0,*]) ; max at bottom
    if zmax*1.e-5 gt -1. and zmax*1.e-5 lt 1. then $
      ncdata = ncdata * 1.e-5 + 97. else $
      ncdata = ncdata * 1.e-5  ; cm to km
    field.units = 'KM'
  end
;
; Geopotential height (WACCM):
;  'Z3': begin
;    ncdata = ncdata * 1.e-3  ; m to km
;    field.units = 'KM'
;   print,'Procfield: convert Z3 from m to km: min,max=',min(ncdata),max(ncdata)
;  end
;
; No processing required:
  else: begin
    print,'Note: procfield did nothing to field ',field.name
  end
endcase

ENDELSE ; Not WACCM input file
;
; Set pointer: 
field.data = ptr_new(ncdata)
end